Ett av de grundläggande antagandena för linjär modellering är konstant, eller homogen, varians. Vad innebär det exakt? Låt oss simulera några data som uppfyller detta villkor för att illustrera begreppet.
Nedan skapar vi en sorterad vektor av tal från 1 till 10 som kallas x och skapar sedan en vektor av tal som kallas y som är en funktion av x. När vi plottar x mot y får vi en rät linje med ett intercept på 1,2 och en lutning på 2,1.
x <- seq(1,10, length.out = 100)y <- 1.2 + 2.1 * xplot(x, y)
Nu ska vi lägga till lite ”brus” till våra data så att y inte helt bestäms av x. Det kan vi göra genom att slumpmässigt dra värden från en teoretisk normalfördelning med medelvärde 0 och en viss fastställd varians och sedan lägga till dem i formeln som genererar y. Med funktionen rnorm i R kan vi enkelt göra detta. Nedan drar vi 100 slumpmässiga värden från en normalfördelning med medelvärde 0 och standardavvikelse 2 och sparar dem som en vektor kallad noise. (Kom ihåg att standardavvikelsen helt enkelt är kvadratroten av variansen.) Sedan genererar vi y med bruset tillsatt. Med funktionen set.seed(222) kan du få fram samma ”slumpmässiga” data om du vill följa med. Slutligen plottar vi data på nytt.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2)y <- 1.2 + 2.1 * x + noiseplot(x, y)
Nu har vi data som är en kombination av en linjär, deterministisk komponent (\(y = 1,2 + 2,1x\)) och slumpmässigt brus som dras från en \(N(0, 2)\)-fördelning. Detta är de grundläggande antaganden vi gör om våra data när vi anpassar en traditionell linjär modell. Nedan använder vi funktionen lm för att ”återskapa” de ”sanna” värden som vi använde för att generera data, varav det finns tre:
- Avsnittet: 1,2
- Höjningen: 2.1
- Standardavvikelsen: 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-16Skattningarna (Intercept) och x (dvs. lutningen) i avsnittet Koefficienter, 1,27 och 2,09, ligger ganska nära 1,2 respektive 2,1. Reststandardfelet på 1,926 ligger också ganska nära det konstanta värdet 2. Vi producerade en ”bra” modell eftersom vi visste hur y genererades och gav lm-funktionen den ”korrekta” modellen att anpassa. Även om vi inte kan göra detta i verkligheten är det en användbar övning för att hjälpa oss att förstå antagandena för linjär modellering.
Nu vad händer om variansen inte var konstant? Vad händer om vi multiplicerar standardavvikelsen 2 med kvadratroten av x? Som vi ser i diagrammet nedan ökar den vertikala spridningen av punkterna när x ökar.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2 * sqrt(x))y <- 1.2 + 2.1 * x + noiseplot(x, y)
Vi multiplicerar 2 med sqrt(x) eftersom vi anger standardavvikelsen. Om vi kunde ange varians skulle vi ha multiplicerat 4 med bara x.
Om vi anpassar samma modell med lm får vi ett reststandardfel 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 vet att det inte stämmer, eftersom vi simulerade data. Det fanns ingen konstant standardavvikelse när vi skapade ”bruset”. Varje slumpmässigt värde drogs från en annan normalfördelning, var och en med medelvärde 0 men en standardavvikelse som varierade enligt x. Detta innebär att vårt antagande om konstant varians bryts. Hur skulle vi upptäcka detta i verkligheten?
Det vanligaste sättet är att plotta residualer mot anpassade värden. Detta är lätt att göra i R. Kalla bara plot på modellobjektet. Detta genererar fyra olika plottar för att bedöma de traditionella modelleringsantagandena. Se det här blogginlägget för mer information. De plottar vi är intresserade av är den första och tredje plotten, som vi kan specificera med argumentet which.
plot(m2, which = 1)
plot(m2, which = 3)
I den första plotten ser vi att variabiliteten runt 0 ökar när vi rör oss längre till höger med större anpassade värden. I den tredje plotten ser vi också att variabiliteten ökar när vi rör oss åt höger, även om residualerna denna gång har standardiserats och transformerats till kvadratrotsskalan. Den positiva banan för den jämna röda linjen indikerar en ökad varians.
Så nu när vi har bekräftat att vårt antagande om icke-konstant varians är ogiltigt, vad kan vi göra? Ett tillvägagångssätt är att logtransformera data. Detta fungerar ibland om din svarsvariabel är positiv och starkt skev. Men det är inte riktigt fallet här. y är endast svagt skevt. (Ring hist() på y för att verifiera.) Dessutom vet vi att våra data inte är på en logaritmisk skala.
En annan metod är att modellera den icke-konstanta variansen. Det är ämnet för det här blogginlägget.
För att göra detta kommer vi att använda funktioner från paketet nlme som ingår i grundinstallationen av R. Arbetshästfunktionen är gls, som står för ”generalized least squares”. Vi använder den precis som vi skulle använda lm-funktionen, förutom att vi också använder weights-argumentet tillsammans med en handfull variansfunktioner. Låt oss demonstrera med ett exempel.
Nedan anpassar vi den ”korrekta” modellen till våra data som uppvisade icke-konstant varians. Vi laddar nlme-paketet så att vi kan använda gls-funktionen1. Vi anger modellens syntax som tidigare, y ~ x. Sedan använder vi argumentet weights för att ange variansfunktionen, i det här fallet varFixed, som också ingår i paketet nlme. Detta säger att vår variansfunktion inte har några parametrar och en enda kovariat, x, vilket är exakt hur vi genererade data. Syntaxen, ~x, är en ensidig formel som kan läsas som ”modellvarians som en funktion av 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 ger ett residualstandardfel på 1,925 som ligger mycket nära 2, det ”sanna” värde vi använde för att generera data. Värdena för (Intercept) och lutning, 1,37 och 2,09, ligger också mycket nära 1,2 och 2,1.
En gång till kan vi anropa plot på modellobjektet. I det här fallet genereras endast en plott: standardiserade residualer mot anpassade värden:
plot(vm1)
Denna plott ser bra ut. Så länge vi modellerar vår varians som en funktion av x, passar den anpassade modellen varken över eller under på något systematiskt sätt (till skillnad från när vi använde lm för att anpassa modellen och antog konstant varians.)
Funktionen varFixed skapar vikter för varje observation, symboliserat som \(w_i\). Tanken är att ju högre vikt en viss observation har, desto mindre är variansen för den normalfördelning från vilken dess bruskomponent drogs. I vårt exempel ökar variansen när x ökar. Därför bör högre vikter förknippas med mindre x-värden. Detta kan åstadkommas genom att ta reciproken av x, det vill säga \(w_i = 1/x\). Så när \(x = 2\), \(w = 1/2\). När \(x = 10\), \(w = 1/10\). Större x får mindre vikter.
För att säkerställa att större vikter är förknippade med mindre varianser delar vi slutligen den konstanta variansen med vikten. Matematiskt uttryckt,
\
Och genom att ta kvadratroten för att uttrycka standardavvikelsen,
\
Så ju större nämnaren är (dvs. ju större vikt och därmed mindre x), desto mindre är variansen och desto mer exakt är den observerade y.
Incidentellt kan vi använda lm för att vikta observationer också. Även den har ett weights argument som gls funktionen. Den enda skillnaden är att vi måste vara mer explicita om hur vi uttrycker vikterna. I det här fallet måste vi ange det reciproka av x. Lägg märke till att resultatet är nästan identiskt med vad vi får med gls och 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-16Kraften hos nlme-paketet är att det tillåter en mängd olika variansfunktioner. Funktionen varFixed som vi just illustrerade är den enklaste och något som kan göras med lm. De andra variansfunktionerna inkluderar:
varIdentvarPowervarExpvarConstPowervarComb
Låt oss utforska var och en av dessa med hjälp av simulerade data.
varIdent
Funktionen varIdent gör att vi kan modellera olika varianser per stratum. För att se hur den fungerar ska vi först simulera data med denna egenskap. Nedan använder vi set.seed(11) ifall någon vill simulera samma slumpmässiga data. Sedan ställer vi in n lika med 400, antalet observationer som vi kommer att simulera. x genereras på samma sätt som tidigare. I det här exemplet inkluderar vi ytterligare en prediktor som heter g och som kan betraktas som kön. Vi tar ett slumpmässigt urval av 400 värden för ”m” och ”f”. Därefter simulerar vi en vektor med två standardavvikelser och sparar som msd. Om g == "m" är standardavvikelsen \(2 \ gånger 2,5\). Annars är det \(2\) för g == "f". Vi använder denna vektor i nästa rad för att generera y. Lägg märke till var msd sätts in i rnorm-funktionen. Lägg också märke till hur vi genererar y. Vi inkluderar en interaktion mellan x och y. När g == "f" är interceptet och lutningen 1,2 och 2,1. När g == "m" används är interceptet och lutningen (1,2 + 2,8) respektive (2,1 + 2,8). Slutligen placerar vi våra simulerade data i en dataram och plottar 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()
Märk de olika varianserna för varje nivå av g. Variansen för ”m” är mycket större än variansen för ”f”. Den har mycket större spridning än ”f”. Vi simulerade data på det sättet. Vi ställde in variansen för ”m” till att vara 2,5 gånger större än variansen för ”f”.
Nu ska vi använda gls-funktionen med varIdent för att försöka återskapa dessa sanna värden. Vi använder samma sätt som tidigare: definiera vår variansfunktion i weights-argumentet. Nedan specificerar vi i argumentet form att formeln för modellering av variansen är villkorad av g. Uttrycket ~ 1|g är en ensidig formel som säger att variansen skiljer sig åt mellan nivåerna av g. 1 betyder bara att vi inte har några ytterligare kovariater i vår modell för varians. Det är möjligt att inkludera en formel som ~ x|g, men det skulle vara felaktigt i det här fallet eftersom vi inte använde x för att generera variansen. En titt på diagrammet visar också att även om variabiliteten i y skiljer sig åt mellan grupperna, så ökar inte variabiliteten när x ökar.
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 märker du längst ner i utmatningen att det skattade residualstandardfelet är 2.0053974, vilket ligger mycket nära det ”sanna” värdet 2. Lägg också märke till att vi i avsnittet ”Variansfunktion” får ett uppskattat värde på 2,58127 för grupp ”m”, vilket ligger mycket nära det ”sanna” värdet 2,5 som vi använde för att generera den olika variansen för grupp ”m”. När du använder varIdent-funktionen för att skatta olika varianser mellan stratanivåer kommer en av nivåerna i allmänhet att sättas till baslinjen, och de andra kommer att skattas som multiplar av standardfelet för residualerna. Eftersom ”f” kommer före ”m” i alfabetisk ordning, har ”f” i detta fall satts till baslinjen, dvs. 1. Det uppskattade standardfelet för gruppen ”f” är \(2,005397 \ gånger 1\). Det uppskattade standardfelet för gruppen ”m” är \(2,005397 \ gånger 2,58127\).
Det är viktigt att notera att detta bara är uppskattningar. För att få en känsla för osäkerheten i dessa uppskattningar kan vi använda funktionen intervals från paketet nlme för att få fram 95 % konfidensintervall. För att minska utdata sparar vi resultatet och visar valda element i objektet.
int <- intervals(vm2)Elementet varStruct innehåller det 95-procentiga konfidensintervallet för parametern som skattades i variansfunktionen. Parametern i det här fallet är multiplikatorn för det kvarvarande standardfelet för den manliga gruppen.
int$varStruct## lower est. upper## m 2.245615 2.58127 2.967095## attr(,"label")## "Variance function:"Elementet sigma innehåller det 95-procentiga konfidensintervallet för det kvarvarande standardfelet.
int$sigma## lower est. upper ## 1.818639 2.005397 2.211335 ## attr(,"label")## "Residual standard error:"Båda intervallen är ganska små och innehåller det ”sanna” värdet som vi använde för att generera data.
varPower
Funktionen varPower gör det möjligt för oss att modellera varians som en potens av det absoluta värdet av en kovariat. För att förklara detta kommer vi återigen att simulera data med denna egenskap. Nedan är det viktigaste att notera sd-argumentet i rnorm-funktionen. Där tar vi standardavvikelsen 2 och multiplicerar den sedan med xs absoluta värde upphöjt till potensen 1,5. Detta liknar hur vi simulerade data för att demonstrera varFixed-funktionen ovan. I det fallet antog vi helt enkelt att exponenten var 0,5. (Kom ihåg att ta kvadratroten av ett tal är likvärdigt med att höja det till potensen 0,5). Här valde vi godtyckligt en potens på 1,5. När vi använder gls med varPower försöker vi återfå det ”sanna” värdet 1,5 utöver 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 att data uppvisar den klassiska formen av att variansen ökar när prediktoren ökar. För att korrekt modellera dessa data med hjälp av gls förser vi den med formeln y ~x och använder argumentet weights med varPower. Lägg märke till att vi anger den ensidiga formeln precis som vi gjorde med funktionen varFixed. I den här modellen får vi dock en uppskattning av 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 uppskattas till 1,52, vilket är mycket nära det ”sanna” värdet 1,5. Det uppskattade residualstandardfelet på 1,8729439 ligger också ganska nära 2. Båda intervallen fångar upp de värden som vi använde för att simulera data. Koefficienterna i modellen är däremot något dåligt uppskattade. Detta är inte förvånande med tanke på hur stor variabilitet som finns i y, särskilt för \(x > 2\).
varExp
Funktionen varExp gör det möjligt för oss att modellera variansen som en exponentiell funktion av en kovariabel. Ännu en gång ska vi förklara denna variansfunktion med hjälp av simulerade data. Den enda förändringen är sd argumentet i rnorm funktionen. Vi har ett fast värde på 2 som vi multiplicerar med x som i sin tur multipliceras med 0,5 och exponeras. Lägg märke till hur snabbt variansen ökar när vi ökar x. Detta beror på den exponentiella tillväxten av 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()
För att arbeta bakåt och återfå dessa värden använder vi varExp-funktionen i weights-argumentet i gls. Den ensidiga formeln ändras inte i detta fall. Den säger att modellera variansen som en funktion av x. Funktionen varExp säger att x har multiplicerats med något värde och exponerats, så gls kommer att försöka skatta det värdet.
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:"Skattningen av ”expon” på 0,4845623 i avsnittet ”variansfunktion” ligger mycket nära det värde som vi har angett, 0,5. På samma sätt ligger det uppskattade residualstandardfelet på 2,1191918 nära det ”sanna” värdet 2. Skattningarna av modellens koefficienter ligger också nära de värden som vi använde för att generera data, men lägg märke till osäkerheten i (Intercept). Hypotesprövningen i sammanfattningen kan inte utesluta ett negativt intercept. Detta är återigen inte förvånande eftersom y har så mycket icke-konstant varians samt det faktum att vi inte har några observationer av x som är lika med 0. Eftersom interceptet är det värde som y antar när x är lika med 0, är vårt uppskattade intercept en extrapolering till en händelse som vi inte har observerat.
varConstPower
Funktionen varConstPower gör det möjligt att modellera varians som en positiv konstant plus en potens av det absoluta värdet av en kovariabel. Det är en munfull, men detta är i princip samma sak som varPower-funktionen, förutom att vi nu lägger till en positiv konstant till den. Följande kod simulerar data för vilka varConstPower-funktionen skulle vara lämplig att använda. Lägg märke till att den är identisk med den kod vi använde för att simulera data för varPower-avsnittet, förutom att vi lägger till 0,7 till x i abs-funktionen. Varför 0,7? Ingen särskild anledning. Det är bara en positiv konstant som vi valde.
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()
Det korrekta sättet att modellera dessa data är att använda gls med funktionen varConstPower. Den ensidiga formeln är densamma som tidigare. Sammanfattningen returnerar tre skattningar för variansmodellen: konstanten, effekten och det kvarvarande standardfelet. Minns att de ”sanna” värdena är 0,7, 1,5 respektive 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:"Intervallen fångar upp de ”sanna” värdena, även om intervallet för konstanten är ganska brett. Om vi inte kände till det sanna värdet skulle det verka som om konstanten rimligen kunde vara 2 eller till och med 4.
varComb
Slutligt gör varComb-funktionen det möjligt för oss att modellera kombinationer av två eller flera variansmodeller genom att multiplicera samman motsvarande variansfunktioner. Det är uppenbart att detta kan rymma några mycket komplexa variansmodeller. Vi ska simulera några grundläggande data som skulle vara lämpliga att använda med varComb-funktionen.
Vad vi gör nedan är att kombinera två olika variansprocesser:
- En som låter standardavvikelsen skilja sig åt mellan könen (”f” = \(2\), ”m” = \(2 \ gånger 2.5\))
- En som låter standardavvikelsen öka när
xökar, därxmultipliceras med 1,5 och exponeras
För att underlätta visualiseringen av data har vi begränsat x till ett intervall från 1 till 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()
Plotten visar ökande varians när x ökar, men också skillnader i varians mellan könen. Variansen för y för grupp ”m” är mycket större än variansen för y i grupp ”f”, särskilt när x är större än 1,5.
För att på ett korrekt sätt modellera den datagenereringsprocess som vi specificerade ovan och för att försöka återskapa de sanna värdena använder vi varComb-funktionen som en omslagsform kring ytterligare två variansfunktioner: varIdent och varExp. Varför dessa två? Därför att vi har olika varianser mellan könen och att variansen ökar exponentiellt som en funktion av 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 residualSammanfattningsavsnittet innehåller två avsnitt för modellering av varians. I det första uppskattas multiplikatorn för gruppen ”m” till 2,437, vilket ligger mycket nära det verkliga värdet 2,5. Exponentialparametern uppskattas till 1,54, vilket ligger extremt nära det sanna värdet 1,5 som vi använde när vi genererade data. Slutligen uppskattas det återstående standardfelet till cirka 1,79, vilket ligger nära det sanna värdet 2. Att ringa intervals(vm6) visar mycket snäva konfidensintervall. Prefixen A och B i elementet varStruct är bara etiketter för de två olika variansmodellerna.
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:"På grund av den stora exponentiella variabiliteten är tyvärr skattningarna av modellkoefficienterna bedrövligt dåliga.
Vad är poängen?
Så varför simulerade vi alla dessa falska data och försökte sedan återskapa ”sanna” värden? Vad är det för nytta med det? Rättvisa frågor. Svaret är att det hjälper oss att förstå vilka antaganden vi gör när vi specificerar och anpassar en modell. När vi specificerar och anpassar följande modell…
m <- lm(y ~ x1 + x2)…säger vi att vi tror att y är ungefär lika med en vägd summa av x1 och x2 (plus ett intercept) med fel som slumpmässiga dragningar från en normalfördelning med medelvärdet 0 och en viss fast standardavvikelse. På samma sätt, om vi specificerar och anpassar följande modell …
m <- gls(y ~ x1 + x2, data = d, weights = varPower(form = ~x1)) … säger vi att vi tror att y är ungefär lika med en vägd summa av x1 och x2 (plus ett intercept) med fel som slumpmässiga dragningar från en normalfördelning med medelvärdet 0 och en standardavvikelse som växer som en multipel av x1 upphöjd med en viss styrka.
Om vi kan simulera data som lämpar sig för dessa modeller har vi en bättre förståelse för de antaganden vi gör när vi använder dessa modeller. Förhoppningsvis har du nu en bättre förståelse för vad du kan göra för att modellera varians med hjälp av nlme-paketet.
För frågor eller förtydliganden om den här artikeln kan du kontakta UVA Library StatLab: [email protected]
Se hela samlingen av artiklar från UVA Library StatLab.
Clay Ford
Statistisk forskningskonsult
University of Virginia Library
7 april 2020
-
1. Paketet
nlmeär kanske mer känt för sinlme-funktion, som används för att anpassa modeller med blandade effekter (dvs. modeller med både fasta och slumpmässiga effekter). I det här blogginlägget demonstreras variansfunktioner med hjälp avgls, som inte anpassar slumpmässiga effekter. Allt vi presenterar i det här blogginlägget kan dock användas medlme.
.