En af de grundlæggende antagelser for lineær modellering er konstant, eller homogen, varians. Hvad betyder det helt præcist? Lad os simulere nogle data, der opfylder denne betingelse, for at illustrere begrebet.
Nedenfor opretter vi en sorteret vektor af tal fra 1 til 10 kaldet x, og opretter derefter en vektor af tal kaldet y, der er en funktion af x. Når vi plotter x mod y, får vi en ret linje med et skæringspunkt på 1,2 og en hældning på 2,1.
x <- seq(1,10, length.out = 100)y <- 1.2 + 2.1 * xplot(x, y)
Nu skal vi tilføje noget “støj” til vores data, så y ikke er fuldstændig bestemt af x. Det kan vi gøre ved at trække tilfældige værdier fra en teoretisk normalfordeling med middelværdi 0 og en vis varians og derefter tilføje dem til den formel, der genererer y. Funktionen rnorm i R giver os mulighed for nemt at gøre dette. Nedenfor trækker vi 100 tilfældige værdier fra en normalfordeling med middelværdi 0 og standardafvigelse 2 og gemmer dem som en vektor kaldet noise. (Husk, at standardafvigelse simpelthen er kvadratroden af variansen.) Derefter genererer vi y med støjen tilføjet. Funktionen set.seed(222) giver dig mulighed for at få de samme “tilfældige” data, hvis du vil følge med. Til sidst plotter vi dataene igen.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2)y <- 1.2 + 2.1 * x + noiseplot(x, y)
Nu har vi data, der er en kombination af en lineær, deterministisk komponent (\(y = 1,2 + 2,1x\))) og tilfældig støj trukket fra en \(N(0, 2)\)-fordeling. Dette er de grundlæggende antagelser, vi gør om vores data, når vi tilpasser en traditionel lineær model. Nedenfor bruger vi funktionen lm til at “genfinde” de “sande” værdier, som vi har brugt til at generere dataene, hvoraf der er tre:
- Skæringspunktet: 1,2
- Hældningen: 2.1
- Standardafvigelsen: 2
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-16Den (Intercept) og x (dvs. hældningen) estimerede værdi i afsnittet Koefficienter, 1,27 og 2,09, ligger ret tæt på henholdsvis 1,2 og 2,1. Reststandardfejlen på 1,926 er også ret tæt på den konstante værdi på 2. Vi producerede en “god” model, fordi vi vidste, hvordan y blev genereret, og gav lm-funktionen den “korrekte” model til at passe. Selv om vi ikke kan gøre dette i det virkelige liv, er det en nyttig øvelse til at hjælpe os med at forstå forudsætningerne for lineær modellering.
Hvad nu, hvis variansen ikke var konstant? Hvad hvis vi multiplicerede standardafvigelsen på 2 med kvadratroden af x? Som vi ser i nedenstående plot, øges den lodrette spredning af punkterne, når x øges.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2 * sqrt(x))y <- 1.2 + 2.1 * x + noiseplot(x, y)
Vi har ganget 2 med sqrt(x), fordi vi angiver standardafvigelsen. Hvis vi kunne angive varians, ville vi have ganget 4 med bare x.
Hvis vi tilpasser den samme model ved hjælp af lm, får vi en residualstandardfejl på 4,488.
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-16Vi ved, at det ikke er rigtigt, fordi vi simulerede dataene. Der var ingen konstant standardafvigelse, da vi skabte “støjen”. Hver tilfældig værdi blev trukket fra en anden normalfordeling, hver med middelværdi 0, men med en standardafvigelse, der varierede i henhold til x. Det betyder, at vores antagelse om konstant varians er overtrådt. Hvordan ville vi opdage dette i det virkelige liv?
Den mest almindelige måde er at plotte residualer mod tilpassede værdier. Dette er let at gøre i R. Du skal blot kalde plot på modelobjektet. Dette genererer fire forskellige plot til vurdering af de traditionelle modelleringsforudsætninger. Se dette blogindlæg for flere oplysninger. De plots, vi er interesserede i, er det 1. og 3. plot, som vi kan angive med which-argumentet.
plot(m2, which = 1)
plot(m2, which = 3)
I det første plot ser vi, at variabiliteten omkring 0 stiger, når vi bevæger os længere mod højre med større tilpassede værdier. I det tredje plot ser vi også stigende variabilitet, efterhånden som vi bevæger os mod højre, selv om residualerne denne gang er blevet standardiseret og transformeret til kvadratrodsskalaen. Den positive bane for den glatte røde linje indikerer en stigning i variansen.
Så nu, hvor vi har fået bekræftet, at vores antagelse om ikke-konstant varians er ugyldig, hvad kan vi så gøre? En tilgang er at logtransformere dataene. Dette virker nogle gange, hvis din responsvariabel er positiv og stærkt skævvredet. Men det er ikke rigtig tilfældet her. y er kun svagt skævvredet. (Kald hist() på y for at verificere det.) Desuden ved vi, at vores data ikke er på en logaritmisk skala.
En anden fremgangsmåde er at modellere den ikke-konstante varians. Det er emnet for dette blogindlæg.
For at gøre dette vil vi bruge funktioner fra nlme-pakken, der er inkluderet i grundinstallationen af R. Arbejdsfunktionen er gls, som står for “generalized least squares” (generaliserede mindste kvadrater). Vi bruger den på samme måde som vi ville bruge lm-funktionen, bortset fra at vi også bruger weights-argumentet sammen med en håndfuld variansfunktioner. Lad os demonstrere det med et eksempel.
Nedenfor tilpasser vi den “korrekte” model til vores data, der udviste en ikke-konstant varians. Vi indlæser nlme-pakken, så vi kan bruge gls-funktionen1. Vi angiver modelsyntaksen som før, y ~ x. Derefter bruger vi weights-argumentet til at angive variansfunktionen, i dette tilfælde varFixed, som også er en del af nlme-pakken. Dette siger, at vores variansfunktion ikke har nogen parametre og en enkelt covariat, x, hvilket er præcis, hvordan vi genererede dataene. Syntaksen, ~x, er en ensidig formel, der kan læses som “modelvarians som en funktion af x.”
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 residualResultatet giver en residualstandardfejl på 1,925, der er meget tæt på 2, den “sande” værdi, som vi brugte til at generere dataene. Værdierne for (Intercept) og hældning, 1,37 og 2,09, er også meget tæt på 1,2 og 2,1.
Gennu engang kan vi kalde plot på modelobjektet. I dette tilfælde genereres kun ét plot: standardiserede residualer versus tilpassede værdier:
plot(vm1)
Dette plot ser godt ud. Så længe vi modellerer vores varians som en funktion af x, passer den tilpassede model hverken overfit eller underfit på nogen systematisk måde (i modsætning til, da vi brugte lm til at tilpasse modellen og antog konstant varians.)
Funktionen varFixed opretter vægte for hver observation, symboliseret som \(w_i\). Ideen er, at jo højere vægt en given observation har, jo mindre er variansen af den normalfordeling, hvorfra dens støjkomponent blev udtrukket. I vores eksempel øges variansen, efterhånden som x øges. Derfor bør højere vægte være forbundet med mindre x-værdier. Dette kan opnås ved at tage den reciprokke værdi af x, dvs. \(w_i = 1/x\). Så når \(x = 2\), er \(w = 1/2\). Når \(x = 10\), \(w = 1/10\). Større x får mindre vægte.
Finalt, for at sikre at større vægte er forbundet med mindre varianser, dividerer vi den konstante varians med vægten. Matematisk udtrykt,
\
Og ved at tage kvadratroden for at udtrykke standardafvigelsen,
Så jo større nævneren er (dvs. jo større vægt og dermed mindre x), jo mindre varians og mere præcis er den observerede y.
I øvrigt kan vi også bruge lm til at vægte observationer. Den har også et weights-argument ligesom gls-funktionen. Den eneste forskel er, at vi er nødt til at være mere eksplicitte med hensyn til, hvordan vi udtrykker vægtene. I dette tilfælde er vi nødt til at angive det reciprokke af x. Bemærk, at resultatet er næsten identisk med det, vi får ved at bruge gls og varFixed.
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-16Den stærke side af nlme-pakken er, at den tillader en række forskellige variansfunktioner. varFixed-funktionen, som vi lige har illustreret, er den enkleste og noget, der kan udføres med lm. De andre variansfunktioner omfatter:
varIdentvarPowervarExpvarConstPowervarComb
Lad os undersøge hver af disse ved hjælp af simulerede data.
varIdent
Den varIdent funktion giver os mulighed for at modellere forskellige varianser pr. stratum. For at se, hvordan den fungerer, skal vi først simulere data med denne egenskab. Nedenfor bruger vi set.seed(11) i tilfælde af, at nogen ønsker at simulere de samme tilfældige data. Derefter sætter vi n lig med 400, det antal observationer, vi vil simulere. x genereres på samme måde som før. I dette eksempel medtager vi en ekstra prædiktor kaldet g, som kan opfattes som køn. Vi udtager tilfældigt 400 værdier af “m” og “f”. Dernæst simulerer vi en vektor med to standardafvigelser og gemmer den som msd. Hvis g == "m", så er standardafvigelsen \(2 \ gange 2,5\). Ellers er det \(2\) for g == "f". Vi bruger denne vektor i den næste linje til at generere y. Bemærk, hvor msd er sat ind i rnorm-funktionen. Læg også mærke til, hvordan vi genererer y. Vi medtager et samspil mellem x og y. Når g == "f", er interceptet og hældningen 1,2 og 2,1. Når g == "m", er interceptet og hældningen henholdsvis (1,2 + 2,8) og (2,1 + 2,8). Endelig placerer vi vores simulerede data i en dataramme og plotter med 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()
Bemærk de forskellige varianser for hvert niveau af g. Variansen for “m” er meget større end variansen for “f”. Den har en meget større spredning end “f”. Vi simulerede dataene på den måde. Vi indstillede variansen for “m” til at være 2,5 gange så stor som variansen for “f”.
Nu skal vi bruge gls-funktionen med varIdent i et forsøg på at genfinde disse sande værdier. Vi bruger den samme måde som før: vi definerer vores variansfunktion i weights-argumentet. Nedenfor angiver vi i form-argumentet, at formlen til modellering af variansen er betinget af g. Udtrykket ~ 1|g er en ensidig formel, der siger, at variansen varierer mellem niveauerne af g. 1 betyder blot, at vi ikke har yderligere kovariater i vores model for varians. Det er muligt at medtage en formel som ~ x|g, men det ville være forkert i dette tilfælde, da vi ikke har brugt x til at generere variansen. Når vi ser på plottet, viser det også, at mens variabiliteten i y er forskellig mellem grupperne, øges variabiliteten ikke, når x stiger.
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 residualFørst bemærker vi nederst i outputtet, at den estimerede residualstandardfejl er 2.0053974, hvilket er meget tæt på den “sande” værdi på 2. Bemærk også i afsnittet “Variansfunktion”, at vi får en estimeret værdi på 2,58127 for gruppe “m”, hvilket er meget tæt på den “sande” værdi på 2,5, som vi brugte til at generere den forskellige varians for gruppe “m”. Generelt gælder det, at når du bruger varIdent-funktionen til at estimere forskellige varianser mellem niveauer af strata, vil et af niveauerne blive sat til grundlinjen, og de andre vil blive estimeret som multipla af den residuelle standardfejl. I dette tilfælde, da “f” kommer før “m” alfabetisk, blev “f” sat til baseline, dvs. 1. Den estimerede residualstandardfejl for gruppe “f” er \(2,005397 \ gange 1\). Den anslåede residualstandardfejl for gruppe “m” er \(2,005397 \ gange 2,58127\).
Det er vigtigt at bemærke, at der kun er tale om skøn. For at få en fornemmelse af usikkerheden ved disse skøn kan vi bruge intervals-funktionen fra nlme-pakken til at få 95 % konfidensintervaller. For at reducere outputtet gemmer vi resultatet og viser udvalgte elementer i objektet.
int <- intervals(vm2)Elementet varStruct indeholder det 95% konfidensinterval for den parameter, der er estimeret i variansfunktionen. Parameteren er i dette tilfælde den resterende standardfejlmultiplikator for mandegruppen.
int$varStruct## lower est. upper## m 2.245615 2.58127 2.967095## attr(,"label")## "Variance function:"Elementet sigma indeholder det 95 % konfidensinterval for den resterende standardfejl.
int$sigma## lower est. upper ## 1.818639 2.005397 2.211335 ## attr(,"label")## "Residual standard error:"Både intervallerne er ret små og indeholder den “sande” værdi, som vi har brugt til at generere dataene.
varPower
Funktionen varPower giver os mulighed for at modellere variansen som en potens af den absolutte værdi af en kovariabel. For at hjælpe med at forklare dette vil vi endnu en gang simulere data med denne egenskab. Nedenfor er det vigtigste at bemærke sd-argumentet i rnorm-funktionen. Det er der, hvor vi tager standardafvigelsen på 2 og derefter multiplicerer den med den absolutte værdi af x hævet til potensen af 1,5. Dette svarer til den måde, hvorpå vi simulerede data for at demonstrere varFixed-funktionen ovenfor. I det tilfælde antog vi simpelthen, at eksponenten var 0,5. (Husk, at det at tage kvadratroden af et tal svarer til at hæve det til potensen 0,5). Her valgte vi vilkårligt en potens på 1,5. Når vi bruger gls med varPower, vil vi forsøge at genfinde den “sande” værdi på 1,5 ud over 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()
Vi ser, at dataene udviser den klassiske form for varians, der stiger, når prædiktoren øges. For at modellere disse data korrekt ved hjælp af gls, forsyner vi den med formlen y ~x og bruger weights-argumentet med varPower. Bemærk, at vi angiver den ensidige formel, ligesom vi gjorde med varFixed-funktionen. I denne model får vi dog et skøn for potensen.
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 residualint <- 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:"Potensen skønnes at være 1,52, hvilket er meget tæt på den “sande” værdi på 1,5. Den estimerede residualstandardfejl på 1,8729439 er også ret tæt på 2. Begge intervaller indfanger de værdier, vi har brugt til at simulere dataene. Koefficienterne i modellen er på den anden side noget dårligt estimeret. Dette er ikke overraskende i betragtning af, hvor stor variabilitet der er i y, især for \(x > 2\).
varExp
Den varExp funktion giver os mulighed for at modellere variansen som en eksponentiel funktion af en kovariabel. Endnu en gang vil vi forklare denne variansfunktion ved hjælp af simulerede data. Den eneste ændring er i sd-argumentet i rnorm-funktionen. Vi har en fast værdi på 2, som vi multiplicerer med x, der selv multipliceres med 0,5 og eksponeres. Bemærk, hvor hurtigt variansen stiger, når vi øger x. Dette skyldes den eksponentielle vækst i variansen.
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()
For at arbejde baglæns og genvinde disse værdier bruger vi varExp-funktionen i weights-argumentet i gls. Den ensidige formel ændres ikke i dette tilfælde. Den siger, at man skal modellere variansen som en funktion af x. varExp-funktionen siger, at x er blevet ganget med en vis værdi og eksponeret, så gls vil forsøge at estimere denne værdi.
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 residualint <- 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:"Det “ekspon”-estimat på 0,4845623 i afsnittet “variansfunktion” er meget tæt på vores angivne værdi på 0,5. Ligeledes er den estimerede residualstandardfejl på 2,1191918 tæt på den “sande” værdi på 2. Skønnene af modelkoefficienterne ligger også tæt på de værdier, vi brugte til at generere dataene, men læg mærke til usikkerheden i (Intercept). Hypotesetesten i resuméet kan ikke udelukke et negativt intercept. Dette er igen ikke overraskende, da y har så meget ikke-konstant varians, og da vi ikke har nogen observationer af x lig med 0. Da interceptet er den værdi, som y antager, når x er lig med 0, er vores estimerede intercept en ekstrapolation til en begivenhed, som vi ikke har observeret.
varConstPower
Funktionen varConstPower giver os mulighed for at modellere varians som en positiv konstant plus en potens af den absolutte værdi af en kovariabel. Det er en mundfuld, men det er grundlæggende det samme som varPower-funktionen, bortset fra at vi nu tilføjer en positiv konstant til den. Den følgende kode simulerer data, som varConstPower-funktionen ville være velegnet til at bruge. Bemærk, at den er identisk med den kode, vi brugte til at simulere data til varPower-delen, bortset fra at vi tilføjer 0,7 til x i abs-funktionen. Hvorfor 0,7? Ingen særlig grund. Det er bare en positiv konstant, som vi har valgt.
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()
Den korrekte måde at modellere disse data på er at bruge gls med varConstPower-funktionen. Den ensidige formel er den samme som før. Resuméet returnerer tre estimater for variansmodellen: konstanten, effekten og den resterende standardfejl. Husk, at de “sande” værdier er henholdsvis 0,7, 1,5 og 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 residualint <- 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:"Intervallerne indfanger de “sande” værdier, selv om intervallet for konstanten er ret bredt. Hvis vi ikke kendte den sande værdi, kunne konstanten tilsyneladende plausibelt være 2 eller endog 4.
varComb
Sluttelig giver varComb-funktionen os mulighed for at modellere kombinationer af to eller flere variansmodeller ved at multiplicere de tilsvarende variansfunktioner sammen. Det er klart, at dette kan rumme nogle meget komplekse variansmodeller. Vi vil simulere nogle grundlæggende data, der ville være passende at bruge med varComb-funktionen.
Det, vi gør nedenfor, er at kombinere to forskellige variansprocesser:
- En, der tillader standardafvigelse at være forskellig mellem køn (“f” = \(2\), “m” = \(2 \ gange 2.5\))
- En, der tillader standardafvigelsen at stige, når
xstiger, hvorxganges med 1,5 og eksponeres
For at hjælpe med at visualisere dataene har vi begrænset x til et interval fra 1 til 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()
Det viste plot viser stigende varians, når x stiger, men også forskelle i varians mellem kønnene. Variansen af y for gruppe “m” er meget større end variansen af y i gruppe “f”, især når x er større end 1,5.
For at modellere den datagenereringsproces, vi specificerede ovenfor, korrekt og forsøge at genfinde de sande værdier, bruger vi varComb-funktionen som en indpakning omkring to andre variansfunktioner: varIdent og varExp. Hvorfor disse to? Fordi vi har forskellige varianser mellem kønnene, og vi har varians, der stiger eksponentielt som en funktion af 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 residualSammenfatningsafsnittet indeholder to afsnit til modellering af varians. I det første vurderes multiplikatoren for gruppen “m” til at være 2,437, hvilket er meget tæt på den sande værdi på 2,5. Den eksponentielle parameter estimeres til at være 1,54, hvilket ligger ekstremt tæt på den sande værdi på 1,5, som vi brugte ved genereringen af dataene. Endelig anslås reststandardfejlen til at være ca. 1,79, hvilket ligger tæt på den sande værdi på 2. Kald intervals(vm6) viser meget snævre konfidensintervaller. Præfikserne A og B i varStruct-elementet er blot etiketter for de to forskellige variansmodeller.
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:"Unheldigvis er estimaterne af modelkoefficienterne på grund af den store eksponentielle variabilitet sørgeligt dårlige.
Hvad er pointen?
Så hvorfor simulerede vi alle disse falske data og forsøgte derefter at genfinde de “sande” værdier? Hvad er det godt for? Rimelige spørgsmål. Svaret er, at det hjælper os til at forstå, hvilke antagelser vi gør, når vi specificerer og tilpasser en model. Når vi specificerer og tilpasser følgende model …
m <- lm(y ~ x1 + x2) … siger vi, at vi mener, at y er omtrent lig med en vægtet sum af x1 og x2 (plus et intercept), hvor fejlene er tilfældige træk fra en normalfordeling med middelværdi 0 og en fast standardafvigelse. På samme måde, hvis vi specificerer og tilpasser følgende model …
m <- gls(y ~ x1 + x2, data = d, weights = varPower(form = ~x1)) … siger vi, at vi mener, at y er omtrent lig med en vægtet sum af x1 og x2 (plus et intercept), hvor fejlene er tilfældige trækninger fra en normalfordeling med middelværdi 0 og en standardafvigelse, der vokser som et multiplum af x1 forhøjet med en vis potens.
Hvis vi kan simulere data, der passer til disse modeller, har vi en bedre forståelse af de antagelser, vi gør, når vi bruger disse modeller. Forhåbentlig har du nu en bedre forståelse af, hvad du kan gøre for at modellere varians ved hjælp af nlme-pakken.
For spørgsmål eller afklaringer vedrørende denne artikel kan du kontakte UVA Library StatLab: [email protected]
Se hele samlingen af UVA Library StatLab-artikler.
Clay Ford
Statistisk forskningskonsulent
University of Virginia Library
7. april 2020
-
1.
nlme-pakken er måske bedre kendt for sinlme-funktion, som bruges til at tilpasse modeller med blandede effekter (dvs. modeller med både faste og tilfældige effekter). Dette blogindlæg demonstrerer variansfunktioner ved hjælp afgls, som ikke tilpasser tilfældige effekter. Men alt, hvad vi præsenterer i dette blogindlæg, kan bruges medlme.