x <- seq(1,10, length.out = 100)y <- 1.2 + 2.1 * xplot(x, y)

set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2)y <- 1.2 + 2.1 * x + noiseplot(x, y)

m <- lm(y ~ x)summary(m)

## ## Call:## lm(formula = y ~ x)## ## Residuals:## Min 1Q Median 3Q Max ## -5.5831 -1.2165 0.3288 1.3022 4.3714 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.27426 0.44720 2.849 0.00534 ** ## x 2.09449 0.07338 28.541 < 2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 1.926 on 98 degrees of freedom## Multiple R-squared: 0.8926, Adjusted R-squared: 0.8915 ## F-statistic: 814.6 on 1 and 98 DF, p-value: < 2.2e-16

set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2 * sqrt(x))y <- 1.2 + 2.1 * x + noiseplot(x, y)

m2 <- lm(y ~ x)summary(m2)

## ## Call:## lm(formula = y ~ x)## ## Residuals:## Min 1Q Median 3Q Max ## -15.0460 -2.4013 0.5638 2.8734 10.2996 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.313 1.042 1.26 0.211 ## x 2.096 0.171 12.26 <2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 4.488 on 98 degrees of freedom## Multiple R-squared: 0.6051, Adjusted R-squared: 0.6011 ## F-statistic: 150.2 on 1 and 98 DF, p-value: < 2.2e-16

plot(m2, which = 1)

plot(m2, which = 3)

library(nlme)vm1 <- gls(y ~ x, weights = varFixed(~x))summary(vm1)

## Generalized least squares fit by REML## Model: y ~ x ## Data: NULL ## AIC BIC logLik## 576.2928 584.0477 -285.1464## ## Variance function:## Structure: fixed weights## Formula: ~x ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 1.369583 0.6936599 1.974431 0.0511## x 2.085425 0.1504863 13.857900 0.0000## ## Correlation: ## (Intr)## x -0.838## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.8942967 -0.6293867 0.1551594 0.6758773 2.2722755 ## ## Residual standard error: 1.925342 ## Degrees of freedom: 100 total; 98 residual

plot(vm1)

m3 <- lm(y ~ x, weights = 1/x)summary(m3)

## ## Call:## lm(formula = y ~ x, weights = 1/x)## ## Weighted Residuals:## Min 1Q Median 3Q Max ## -5.5725 -1.2118 0.2987 1.3013 4.3749 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.3696 0.6937 1.974 0.0511 . ## x 2.0854 0.1505 13.858 <2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 1.925 on 98 degrees of freedom## Multiple R-squared: 0.6621, Adjusted R-squared: 0.6587 ## F-statistic: 192 on 1 and 98 DF, p-value: < 2.2e-16

varIdent

Funcția varIdent ne permite să modelăm diferite varianțe pe strat. Pentru a vedea cum funcționează, vom simula mai întâi date cu această proprietate. Mai jos folosim set.seed(11) în cazul în care cineva dorește să simuleze aceleași date aleatoare. Apoi setăm n egal cu 400, numărul de observații pe care le vom simula. x este generat la fel ca înainte. În acest exemplu, includem un predictor suplimentar numit g, care poate fi considerat ca fiind sexul. Eșantionăm aleatoriu 400 de valori ale lui „m” și „f”. Apoi simulăm un vector de două deviații standard și îl salvăm ca msd. Dacă g == "m", atunci abaterea standard este \(2 \ ori 2,5\). În caz contrar, este \(2\\) pentru g == "f". Folosim acest vector în linia următoare pentru a genera y. Observați unde este introdus msd în funcția rnorm. De asemenea, observați cum generăm y. Am inclus o interacțiune între x și y. Atunci când g == "f", intercepția și panta sunt de 1,2 și 2,1. Când g == "m", intercepția și panta sunt (1,2 + 2,8) și, respectiv, (2,1 + 2,8). În cele din urmă, plasăm datele noastre simulate într-un cadru de date și le reprezentăm grafic cu ggplot2.

set.seed(11)n <- 400x <- seq(1,10, length.out = n)g <- sample(c("m","f"), size = n, replace = TRUE)msd <- ifelse(g=="m", 2*2.5, 2)y <- 1.2 + 2.1*x - 1.5*(g == "m") + 2.8*x*(g == "m") + rnorm(n, sd = msd)d <- data.frame(y, x, g)library(ggplot2)ggplot(d, aes(x, y, color = g)) + geom_point()

Observați varianțele diferite pentru fiecare nivel de g. Varianța pentru „m” este mult mai mare decât varianța pentru „f”. Are o dispersie mult mai mare decât „f”. Am simulat datele în acest fel. Am setat varianța pentru „m” să fie de 2,5 ori mai mare decât cea a lui „f”.

Acum să folosim funcția gls cu varIdent în încercarea de a recupera aceste valori reale. Folosim același mod ca înainte: definim funcția noastră de varianță în argumentul weights. Mai jos specificăm în argumentul form că formula de modelare a varianței este condiționată de g. Expresia ~ 1|g este o formulă unilaterală care spune că varianța diferă între nivelurile lui g. 1 înseamnă doar că nu avem covariate suplimentare în modelul nostru pentru varianță. Este posibil să includem o formulă precum ~ x|g, dar aceasta ar fi incorectă în acest caz, deoarece nu am utilizat x în generarea varianței. Examinarea graficului arată, de asemenea, că în timp ce variabilitatea în y este diferită între grupuri, variabilitatea nu crește pe măsură ce x crește.

vm2 <- gls(y ~ x*g, data = d, weights = varIdent(form = ~ 1|g))summary(vm2)

## Generalized least squares fit by REML## Model: y ~ x * g ## Data: d ## AIC BIC logLik## 2081.752 2105.64 -1034.876## ## Variance function:## Structure: Different standard deviations per stratum## Formula: ~1 | g ## Parameter estimates:## f m ## 1.00000 2.58127 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 0.9686349 0.3237724 2.99172 0.0029## x 2.1222707 0.0525024 40.42239 0.0000## gm -1.9765090 0.9352152 -2.11343 0.0352## x:gm 2.9957974 0.1553551 19.28355 0.0000## ## Correlation: ## (Intr) x gm ## x -0.901 ## gm -0.346 0.312 ## x:gm 0.304 -0.338 -0.907## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.74115039 -0.67013954 0.01619031 0.69793776 2.72985748 ## ## Residual standard error: 2.005397 ## Degrees of freedom: 400 total; 396 residual

Observați mai întâi, în partea de jos a rezultatului, că eroarea standard reziduală estimată este 2.0053974, foarte aproape de valoarea „adevărată” de 2. De asemenea, observați în secțiunea „Funcția de varianță” că obținem o valoare estimată de 2,58127 pentru grupul „m”, care este foarte aproape de valoarea „adevărată” de 2,5 pe care am folosit-o pentru a genera varianța diferită pentru grupul „m”. În general, atunci când utilizați funcția varIdent pentru a estima diferite varianțe între niveluri de straturi, unul dintre niveluri va fi setat la nivelul de bază, iar celelalte vor fi estimate ca multipli ai erorii standard reziduale. În acest caz, deoarece „f” vine înaintea lui „m” din punct de vedere alfabetic, „f” a fost setat la linia de bază, sau 1. Eroarea standard reziduală estimată pentru grupul „f” este \(2,005397 \ ori 1\). Eroarea standard reziduală reziduală estimată pentru grupul „m” este \(2,005397 \ ori 2,58127\).

Este important să rețineți că acestea sunt doar estimări. Pentru a ne face o idee despre incertitudinea acestor estimări, putem folosi funcția intervals din pachetul nlme pentru a obține intervale de încredere de 95%. Pentru a reduce ieșirea, salvăm rezultatul și vizualizăm elementele selectate ale obiectului.

int <- intervals(vm2)

Elementul varStruct conține intervalul de încredere de 95% pentru parametrul estimat în funcția de varianță. În acest caz, parametrul este multiplicatorul de eroare standard reziduală pentru grupul de bărbați.

int$varStruct

## lower est. upper## m 2.245615 2.58127 2.967095## attr(,"label")## "Variance function:"

Elementul sigma conține intervalul de încredere de 95% pentru eroarea standard reziduală.

int$sigma

## lower est. upper ## 1.818639 2.005397 2.211335 ## attr(,"label")## "Residual standard error:"

Ambele intervale sunt destul de mici și conțin valoarea „adevărată” pe care am folosit-o pentru a genera datele.

varPower

Funcția varPower ne permite să modelăm varianța ca o putere a valorii absolute a unei covariate. Încă o dată, pentru a ne ajuta să explicăm acest lucru, vom simula date cu această proprietate. Mai jos, principalul lucru care trebuie observat este argumentul sd al funcției rnorm. Acesta este locul unde luăm abaterea standard de 2 și apoi o înmulțim cu valoarea absolută a x ridicată la puterea 1,5. Acest lucru este similar cu modul în care am simulat datele pentru a demonstra funcția varFixed de mai sus. În acel caz, am presupus pur și simplu că exponentul a fost 0,5. (Reamintim că luarea rădăcinii pătrate a unui număr este echivalentă cu ridicarea acestuia la puterea 0,5.) Aici am ales în mod arbitrar o putere de 1,5. Atunci când folosim gls cu varPower vom încerca să recuperăm „adevărata” valoare de 1,5 în plus față de 2.

set.seed(4)n <- 1000x <- seq(1,10,length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*abs(x)^1.5)d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()

Veziem că datele prezintă forma clasică de creștere a varianței pe măsură ce crește predictorul. Pentru a modela corect aceste date folosind gls, îi furnizăm formula y ~x și folosim argumentul weights cu varPower. Observați că specificăm formula unilaterală la fel cum am făcut-o cu funcția varFixed. Cu toate acestea, în acest model, vom obține o estimare a puterii.

vm3 <- gls(y ~ x, data = d, weights = varPower(form = ~x))summary(vm3)

## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 8840.188 8859.811 -4416.094## ## Variance function:## Structure: Power of variance covariate## Formula: ~x ## Parameter estimates:## power ## 1.520915 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 2.321502 0.4750149 4.887218 0## x 1.582272 0.2241367 7.059404 0## ## Correlation: ## (Intr)## x -0.845## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.892319514 -0.655237190 0.001653178 0.691241690 3.346273816 ## ## Residual standard error: 1.872944 ## Degrees of freedom: 1000 total; 998 residual

int <- intervals(vm3)int$varStruct

## lower est. upper## power 1.445064 1.520915 1.596765## attr(,"label")## "Variance function:"

int$sigma

## lower est. upper ## 1.650987 1.872944 2.124740 ## attr(,"label")## "Residual standard error:"

Puterea este estimată ca fiind 1,52, ceea ce este foarte aproape de valoarea „adevărată” de 1,5. Eroarea standard reziduală reziduală estimată de 1,8729439 este, de asemenea, destul de apropiată de 2. Ambele intervale surprind valorile pe care le-am folosit pentru a simula datele. Coeficienții din model, pe de altă parte, sunt oarecum slab estimați. Acest lucru nu este surprinzător, având în vedere cât de multă variabilitate există în y, în special pentru \(x > 2\).

varExp

Funcția varExp ne permite să modelăm varianța ca o funcție exponențială a unei covariate. Încă o dată vom explica această funcție de varianță folosind date simulate. Singura schimbare este în argumentul sd al funcției rnorm. Avem o valoare fixă de 2 pe care o înmulțim cu x, care este la rândul său înmulțită cu 0,5 și exponențiată. Observați cât de rapid crește varianța pe măsură ce creștem x. Acest lucru se datorează creșterii exponențiale a varianței.

set.seed(55)n <- 400x <- seq(1, 10, length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*exp(0.5*x))d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()

Pentru a lucra în sens invers și a recupera aceste valori, folosim funcția varExp în argumentul weights al lui gls. Formula unilaterală nu se modifică în acest caz. Ea spune că modelează varianța ca funcție de x. Funcția varExp spune că x a fost înmulțită cu o anumită valoare și exponențiată, astfel încât gls va încerca să estimeze acea valoare.

vm4 <- gls(y ~ x, data = d, weights = varExp(form = ~x))summary(vm4)

## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 3873.153 3889.098 -1932.576## ## Variance function:## Structure: Exponential of variance covariate## Formula: ~x ## Parameter estimates:## expon ## 0.4845623 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 1.198525 1.1172841 1.072713 0.284## x 2.073297 0.4933502 4.202486 0.000## ## Correlation: ## (Intr)## x -0.892## ## Standardized residuals:## Min Q1 Med Q3 Max ## -3.16976627 -0.69525696 0.00614222 0.63718022 2.90423217 ## ## Residual standard error: 2.119192 ## Degrees of freedom: 400 total; 398 residual

int <- intervals(vm4)int$varStruct

## lower est. upper## expon 0.4562871 0.4845623 0.5128374## attr(,"label")## "Variance function:"

int$sigma

## lower est. upper ## 1.786737 2.119192 2.513505 ## attr(,"label")## "Residual standard error:"

Stimarea „expon” de 0,4845623 din secțiunea „funcție de varianță” este foarte apropiată de valoarea noastră specificată de 0,5. De asemenea, eroarea standard reziduală reziduală estimată de 2,1191918 este aproape de valoarea „adevărată” de 2. Estimările coeficienților modelului sunt, de asemenea, apropiate de valorile pe care le-am folosit pentru a genera datele, dar observați incertitudinea din (Intercept). Testul de ipoteză din rezumat nu poate exclude o intercepție negativă. Din nou, acest lucru nu este surprinzător, deoarece y are o varianță neconstantă atât de mare, precum și faptul că nu avem observații ale x egal cu 0. Deoarece intercepția este valoarea pe care o ia y atunci când x este egal cu 0, intercepția noastră estimată este o extrapolare la un eveniment pe care nu l-am observat.

varConstPower

Funcția varConstPower ne permite să modelăm varianța ca o constantă pozitivă plus o putere a valorii absolute a unei covariate. Este un cuvânt greu de spus, dar aceasta este practic aceeași funcție ca și funcția varPower, cu excepția faptului că acum îi adăugăm o constantă pozitivă. Următorul cod simulează date pentru care funcția varConstPower ar fi potrivită pentru a fi utilizată. Observați că este identic cu codul pe care l-am folosit pentru a simula datele pentru secțiunea varPower, cu excepția faptului că adăugăm 0,7 la x în funcția abs. De ce 0,7? Nici un motiv special. Este doar o constantă pozitivă pe care am ales-o.

set.seed(4)n <- 1000x <- seq(1,10,length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*abs(0.7 + x)^1.5)d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()

Modul corect de a modela aceste date este de a folosi gls cu funcția varConstPower. Formula unilaterală este aceeași ca înainte. Rezumatul returnează trei estimări pentru modelul de varianță: constanta, puterea și eroarea standard reziduală. Reamintim că valorile „adevărate” sunt 0,7, 1,5 și, respectiv, 2.

vm5 <- gls(y ~ x, data = d, weights = varConstPower(form = ~x))summary(vm5)

## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 9319.794 9344.323 -4654.897## ## Variance function:## Structure: Constant plus power of variance covariate## Formula: ~x ## Parameter estimates:## const power ## 0.3554921 1.3593016 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 2.763788 0.7607449 3.633003 3e-04## x 1.456677 0.2915621 4.996112 0e+00## ## Correlation: ## (Intr)## x -0.824## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.846999539 -0.664551194 0.006709051 0.679895874 3.287366683 ## ## Residual standard error: 2.887761 ## Degrees of freedom: 1000 total; 998 residual

int <- intervals(vm5)int$varStruct

## lower est. upper## const 0.02658804 0.3554921 4.753063## power 1.14258677 1.3593016 1.576016## attr(,"label")## "Variance function:"

int$sigma

## lower est. upper ## 1.814757 2.887761 4.595195 ## attr(,"label")## "Residual standard error:"

Intervalele surprind valorile „adevărate”, deși intervalul pentru constantă este destul de larg. Dacă nu am cunoaște valoarea adevărată, s-ar părea că constanta ar putea fi plauzibil 2 sau chiar 4.

varComb

În cele din urmă, funcția varComb ne permite să modelăm combinații de două sau mai multe modele de varianță prin înmulțirea împreună a funcțiilor de varianță corespunzătoare. Evident, acest lucru poate acomoda unele modele de varianță foarte complexe. Vom simula câteva date de bază care ar fi adecvate pentru a fi utilizate cu funcția varComb.

Ce facem mai jos este să combinăm două procese de varianță diferite:

• Unul care permite ca abaterea standard să difere între sexe („f” = \(2\), „m” = \(2 \ ori 2.5\))
• Unul care permite ca abaterea standard să crească pe măsură ce x crește, unde x este înmulțit cu 1,5 și exponențiat

Pentru a ajuta la vizualizarea datelor, am limitat x la un interval de la 1 la 3.

set.seed(77)n <- 400x <- seq(1, 3, length.out = n)g <- sample(c("m","f"), size = n, replace = TRUE)msd <- ifelse(g=="m", 2*2.5, 2) * exp(1.5*x)y <- 1.2 + 2.1*x - 1.5*(g == "m") + 2.8*x*(g == "m") + rnorm(n, sd = msd)d <- data.frame(y, x, g)ggplot(d, aes(x, y, color = g)) + geom_point()

Pratgraficul arată creșterea varianței pe măsură ce x crește, dar și diferențe de varianță între sexe. Varianța y pentru grupul „m” este mult mai mare decât varianța y în grupul „f”, în special atunci când x este mai mare de 1,5.

Pentru a modela corect procesul de generare a datelor pe care l-am specificat mai sus și pentru a încerca să recuperăm adevăratele valori, folosim funcția varComb ca un înveliș în jurul a încă două funcții de varianță: varIdent și varExp. De ce aceste două? Pentru că avem varianțe diferite între sexe și pentru că varianța crește exponențial în funcție de x.

vm6 <- gls(y ~ x*g, data = d, weights = varComb(varIdent(form = ~ 1|g), varExp(form = ~x)))summary(vm6)

## Generalized least squares fit by REML## Model: y ~ x * g ## Data: d ## AIC BIC logLik## 4431.45 4459.32 -2208.725## ## Combination of variance functions: ## Structure: Different standard deviations per stratum## Formula: ~1 | g ## Parameter estimates:## f m ## 1.000000 2.437046 ## Structure: Exponential of variance covariate## Formula: ~x ## Parameter estimates:## expon ## 1.540608 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) -5.996337 6.539552 -0.9169339 0.3597## x 8.829971 4.849272 1.8208860 0.0694## gm -17.572238 16.548540 -1.0618603 0.2889## x:gm 16.932938 12.192793 1.3887661 0.1657## ## Correlation: ## (Intr) x gm ## x -0.972 ## gm -0.395 0.384 ## x:gm 0.387 -0.398 -0.974## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.74441479 -0.67478610 -0.00262221 0.66079254 3.14975288 ## ## Residual standard error: 1.794185 ## Degrees of freedom: 400 total; 396 residual

Secțiunea rezumat conține două secțiuni pentru modelarea varianței. Prima estimează că multiplicatorul pentru grupul „m” este de 2,437, ceea ce este foarte aproape de valoarea reală de 2,5. Parametrul exponențial este estimat la 1,54, extrem de aproape de valoarea reală de 1,5 pe care am folosit-o la generarea datelor. În cele din urmă, eroarea standard reziduală este estimată la aproximativ 1,79, aproape de valoarea reală de 2. Apelul intervals(vm6) arată intervale de încredere foarte strânse. Prefixele A și B din elementul varStruct sunt doar etichete pentru cele două modele de varianță diferite.

int <- intervals(vm6)int$varStruct

## lower est. upper## A.m 2.119874 2.437046 2.801673## B.expon 1.419366 1.540608 1.661850## attr(,"label")## "Variance function:"

int$sigma

## lower est. upper ## 1.380008 1.794185 2.332669 ## attr(,"label")## "Residual standard error:"

Din păcate, din cauza variabilității exponențiale mari, estimările coeficienților modelului sunt lamentabil de proaste.

Ce rost are?

Atunci de ce am simulat toate aceste date false și apoi am încercat să recuperăm valorile „adevărate”? La ce bun așa ceva? Întrebări corecte. Răspunsul este că ne ajută să înțelegem ce ipoteze facem atunci când specificăm și potrivim un model. Atunci când specificăm și potrivim următorul model…

m <- lm(y ~ x1 + x2)

…spunem că noi credem că y este aproximativ egal cu o sumă ponderată de x1 și x2 (plus o intercepție), erorile fiind extrageri aleatorii dintr-o distribuție normală cu media 0 și o anumită abatere standard fixă. De asemenea, dacă specificăm și adaptăm următorul model…

m <- gls(y ~ x1 + x2, data = d, weights = varPower(form = ~x1))

… spunem că noi credem că y este aproximativ egal cu o sumă ponderată de x1 și x2 (plus o interceptare), erorile fiind extrageri aleatorii dintr-o distribuție normală cu media 0 și o abatere standard care crește ca multiplu de x1 mărit cu o anumită putere.

Dacă putem simula date potrivite pentru aceste modele, atunci avem o mai bună înțelegere a ipotezelor pe care le facem atunci când folosim aceste modele. Sperăm că acum aveți o mai bună înțelegere a ceea ce puteți face pentru a modela varianța folosind pachetul nlme.

Pentru întrebări sau clarificări cu privire la acest articol, contactați UVA Library StatLab: [email protected]

Vezi întreaga colecție de articole ale UVA Library StatLab.

Clay Ford
Consilier de cercetare statistică
University of Virginia Library
Aprilie 7, 2020