Eén van de basisveronderstellingen van lineaire modellering is constante, of homogene, variantie. Wat betekent dat nu precies? Laten we wat gegevens simuleren die aan deze voorwaarde voldoen om het concept te illustreren.
Hieronder maken we een gesorteerde vector van getallen variërend van 1 tot 10 genaamd x, en vervolgens maken we een vector van getallen genaamd y die een functie is van x. Als we x uitzetten tegen y, krijgen we een rechte lijn met een intercept van 1,2 en een helling van 2,1.
x <- seq(1,10, length.out = 100)y <- 1.2 + 2.1 * xplot(x, y)
Nu gaan we wat “ruis” toevoegen aan onze gegevens, zodat y niet volledig wordt bepaald door x. Dat kunnen we doen door willekeurig waarden te trekken uit een theoretische normale verdeling met gemiddelde 0 en een bepaalde variantie, en die vervolgens toe te voegen aan de formule die y genereert. De functie rnorm in R stelt ons in staat dit eenvoudig te doen. Hieronder trekken we 100 willekeurige waarden uit een normale verdeling met gemiddelde 0 en standaardafwijking 2 en slaan die op als een vector genaamd noise. (Bedenk dat de standaardafwijking eenvoudigweg de vierkantswortel van de variantie is.) Vervolgens genereren we y met de toegevoegde ruis. Met de functie set.seed(222) kunt u dezelfde “willekeurige” gegevens krijgen, voor het geval u wilt meekijken. Tenslotte tekenen we de gegevens opnieuw.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2)y <- 1.2 + 2.1 * x + noiseplot(x, y)
Nu hebben we gegevens die een combinatie zijn van een lineaire, deterministische component (y = 1,2 + 2,1x)) en willekeurige ruis getrokken uit een verdeling met een waarde van N(0, 2)\). Dit zijn de basisveronderstellingen die we maken over onze gegevens wanneer we een traditioneel lineair model toepassen. Hieronder gebruiken we de functie lm om de “ware” waarden te “herstellen” die we hebben gebruikt om de gegevens te genereren, waarvan er drie zijn:
- Het intercept: 1,2
- De helling: 2.1
- De standaardafwijking: 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-16De schattingen voor (Intercept) en x (d.w.z. de helling) in de sectie Coëfficiënten, 1,27 en 2,09, liggen vrij dicht bij respectievelijk 1,2 en 2,1. De residuele standaardafwijking van 1,926 ligt ook vrij dicht bij de constante waarde van 2. We produceerden een “goed” model omdat we wisten hoe y werd gegenereerd en we de lm functie het “juiste” model gaven om te passen. Hoewel we dit in het echt niet kunnen doen, is het een nuttige oefening om ons te helpen de veronderstellingen van lineaire modellering te begrijpen.
Nu, wat als de variantie niet constant was? Wat als we de standaardafwijking van 2 vermenigvuldigen met de vierkantswortel van x? Zoals we in de onderstaande plot zien, neemt de verticale spreiding van de punten toe naarmate x toeneemt.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2 * sqrt(x))y <- 1.2 + 2.1 * x + noiseplot(x, y)
We vermenigvuldigen 2 met sqrt(x) omdat we de standaardafwijking specificeren. Als we variantie konden specificeren, zouden we 4 hebben vermenigvuldigd met slechts x.
Als we hetzelfde model fitten met lm, krijgen we een Residuele Standaardfout van 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-16We weten dat dit niet klopt, omdat we de gegevens hebben gesimuleerd. Er was geen constante standaardafwijking toen we de “ruis” creëerden. Elke willekeurige waarde werd getrokken uit een andere normale verdeling, elk met gemiddelde 0 maar een standaardafwijking die varieerde volgens x. Dit betekent dat onze aanname van constante variantie geschonden is. Hoe zouden we dit in het echt kunnen vaststellen?
De meest gebruikelijke manier is het plotten van residuen versus de geteste waarden. Dit is eenvoudig te doen in R. Roep gewoon plot aan op het model object. Dit genereert vier verschillende plots om de traditionele modelaannames te beoordelen. Zie deze blog post voor meer informatie. De plots waarin we geïnteresseerd zijn, zijn de 1e en 3e plot, die we kunnen specificeren met het which argument.
plot(m2, which = 1)
plot(m2, which = 3)
In de eerste plot zien we de variabiliteit rond 0 toenemen naarmate we verder naar rechts gaan met grotere fitted values. In de derde grafiek zien we ook een toenemende variabiliteit naarmate we verder naar rechts gaan, hoewel de residuen deze keer gestandaardiseerd zijn en getransformeerd naar de vierkantswortelschaal. Het positieve traject van de gladde rode lijn wijst op een toename van de variantie.
Dus nu we bevestigd hebben dat onze veronderstelling van niet-constante variantie ongeldig is, wat kunnen we doen? Eén benadering is de gegevens log te transformeren. Dit werkt soms als uw responsvariabele positief is en sterk scheefgetrokken. Maar dat is hier niet echt het geval. y is slechts in geringe mate scheefgetrokken. (Roep hist() op y aan om te verifiëren.) Bovendien weten we dat onze gegevens niet op een logschaal liggen.
Een andere benadering is het modelleren van de niet-constante variantie. Dat is het onderwerp van deze blog post.
Om dit te doen, zullen we functies gebruiken uit het nlme pakket dat is inbegrepen bij de basisinstallatie van R. Het werkpaard van de functie is gls, dat staat voor “gegeneraliseerde kleinste kwadraten”. We gebruiken deze functie net zoals we de lm functie zouden gebruiken, behalve dat we ook het weights argument gebruiken samen met een handvol variantie functies. Laten we dit met een voorbeeld demonstreren.
Hieronder passen we het “juiste” model op onze gegevens die een niet-constante variantie vertoonden. We laden het nlme pakket zodat we de gls functie1 kunnen gebruiken. We specificeren de modelsyntaxis zoals voorheen, y ~ x. Dan gebruiken we het weights-argument om de variantiefunctie op te geven, in dit geval varFixed, ook onderdeel van het nlme-pakket. Dit zegt dat onze variantiefunctie geen parameters heeft en een enkele covariaat, x, wat precies is hoe we de gegevens hebben gegenereerd. De syntaxis, ~x, is een eenzijdige formule die kan worden gelezen als “modelvariantie als functie van 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 residualHet resultaat geeft een residuele standaardafwijking van 1,925 die heel dicht bij 2 ligt, de “echte” waarde die we hebben gebruikt om de gegevens te genereren. De (Intercept) en helling waarden, 1.37 en 2.09, liggen ook zeer dicht bij 1.2 en 2.1.
Wederom kunnen we plot aanroepen op het model object. In dit geval wordt er slechts één plot gegenereerd: gestandaardiseerde residuen versus de geteste waarden:
plot(vm1)
Deze plot ziet er goed uit. Zolang we onze variantie modelleren als een functie van x, past het model niet systematisch te veel of te weinig (in tegenstelling tot wanneer we lm gebruikten om het model te passen en een constante variantie veronderstelden)
De functie varFixed creëert gewichten voor elke waarneming, gesymboliseerd als \(w_i\). Het idee is dat hoe hoger het gewicht van een gegeven waarneming is, des te kleiner de variantie van de normale verdeling is, waaruit de ruiscomponent is getrokken. In ons voorbeeld, als x toeneemt, neemt de variantie toe. Daarom moeten hogere gewichten worden geassocieerd met kleinere x waarden. Dit kan worden bereikt door de reciproke van x te nemen, d.w.z. \(w_i = 1/x). Dus als x = 2, dan is w = 1/2. Als x = 10, dan is w = 1/10. Grotere x krijgen kleinere gewichten.
Eindelijk, om er zeker van te zijn dat grotere gewichten gepaard gaan met kleinere varianties, delen we de constante variantie door het gewicht. Wiskundig uitgedrukt,
Of door de vierkantswortel te nemen om de standaardafwijking uit te drukken,
Dus hoe groter de noemer (d.w.z. hoe groter het gewicht en dus hoe kleiner de x), des te kleiner de variantie en des te nauwkeuriger de waargenomen y.
Toevallig kunnen we lm ook gebruiken om waarnemingen te wegen. Ook deze functie heeft een weights argument, net als de gls functie. Het enige verschil is dat we explicieter moeten zijn over hoe we de gewichten uitdrukken. In dit geval moeten we de reciproke van x opgeven. Merk op dat het resultaat bijna identiek is aan wat we krijgen met gls en 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-16De kracht van het nlme pakket is dat het een verscheidenheid aan variantie functies toestaat. De varFixed functie die we zojuist hebben geïllustreerd is de eenvoudigste en iets dat kan worden gedaan met lm. De andere variantiefuncties zijn:
varIdentvarPowervarExpvarConstPowervarComb
Laten we elk van deze onderzoeken met behulp van gesimuleerde gegevens.
varIdent
De varIdent functie stelt ons in staat om verschillende varianties per stratum te modelleren. Om te zien hoe het werkt, zullen we eerst gegevens met deze eigenschap simuleren. Hieronder gebruiken we set.seed(11) voor het geval iemand dezelfde willekeurige gegevens wil simuleren. Daarna stellen we n gelijk aan 400, het aantal waarnemingen dat we zullen simuleren. x wordt op dezelfde manier gegenereerd als voorheen. In dit voorbeeld nemen we een extra voorspeller op, g, die kan worden beschouwd als geslacht. We nemen een willekeurige steekproef van 400 waarden van “m” en “f”. Vervolgens simuleren we een vector van twee standaardafwijkingen en slaan die op als msd. Als g == "m", dan is de standaardafwijking (2 maal 2.5). Anders is het 2 keer 2,5 voor g == "f". We gebruiken deze vector in de volgende regel om y te genereren. Merk op waar msd in de rnorm functie is gestopt. Merk ook op hoe we y genereren. We voegen een wisselwerking toe tussen x en y. Bij g == "f" zijn de intercept en de helling 1,2 en 2,1. Wanneer g == "m", zijn de intercept en de helling respectievelijk (1,2 + 2,8) en (2,1 + 2,8). Ten slotte plaatsen we onze gesimuleerde gegevens in een dataframe en plotten we ze met 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()
Let op de verschillende varianties voor elk niveau van g. De variantie voor “m” is veel groter dan de variantie voor “f”. Het heeft veel meer spreiding dan “f”. We hebben de gegevens op die manier gesimuleerd. We stellen de variantie voor “m” in op 2,5 maal die van “f”.
Nu gebruiken we de gls functie met varIdent in een poging om deze werkelijke waarden terug te vinden. We gebruiken dezelfde manier als voorheen: definieer onze variantie functie in het weights argument. Hierna specificeren we in het form argument dat de formule voor het modelleren van variantie voorwaardelijk is voor g. De uitdrukking ~ 1|g is een eenzijdige formule die zegt dat de variantie verschilt tussen de niveaus van g. De 1 betekent alleen dat we geen extra covariaten in ons model voor variantie hebben. Het is mogelijk om een formule als ~ x|g op te nemen, maar dat zou in dit geval onjuist zijn omdat wij x niet hebben gebruikt bij het genereren van de variantie. Als we naar de plot kijken, zien we ook dat de variabiliteit in y tussen de groepen verschilt, maar dat de variabiliteit niet toeneemt naarmate x toeneemt.
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 residualMerk onderaan de uitvoer op dat de geschatte residuele standaardafwijking 2 is.0053974, wat heel dicht ligt bij de “ware” waarde van 2. Merk ook op dat we in de sectie “Variantiefunctie” een geschatte waarde van 2,58127 krijgen voor groep “m”, wat heel dicht ligt bij de “ware” waarde van 2,5 die we hebben gebruikt om de afwijkende variantie voor groep “m” te genereren. In het algemeen, wanneer u de varIdent functie gebruikt om verschillende varianties tussen strataniveaus te schatten, zal één van de niveaus op de basislijn worden gezet, en zullen de andere worden geschat als veelvouden van de residuele standaardfout. Aangezien in dit geval “f” alfabetisch vóór “m” komt, werd “f” op basislijn, of 1, gezet. De geschatte residuele standaardafwijking voor groep “f” is (2,005397 maal 1). De geschatte standaardafwijking voor groep “m” is 2,005397 maal 2,58127. Het is belangrijk op te merken dat dit slechts schattingen zijn. Om een idee te krijgen van de onzekerheid van deze schattingen, kunnen we de functie intervals uit het pakket nlme gebruiken om 95% betrouwbaarheidsintervallen te krijgen. Om de uitvoer te beperken, slaan we het resultaat op en bekijken we geselecteerde elementen van het object.
int <- intervals(vm2)Het varStruct-element bevat het 95%-betrouwbaarheidsinterval voor de parameter die in de variantiefunctie is geschat. De parameter in dit geval is de residuele standaardfout multiplicator voor de mannelijke groep.
int$varStruct## lower est. upper## m 2.245615 2.58127 2.967095## attr(,"label")## "Variance function:"Het sigma element bevat het 95% betrouwbaarheidsinterval voor de residuele standaardfout.
int$sigma## lower est. upper ## 1.818639 2.005397 2.211335 ## attr(,"label")## "Residual standard error:"Beide intervallen zijn vrij klein en bevatten de “ware” waarde die we hebben gebruikt om de gegevens te genereren.
varPower
De varPower functie stelt ons in staat om variantie te modelleren als een macht van de absolute waarde van een covariaat. Nogmaals, om dit te helpen verklaren, zullen wij gegevens met deze eigenschap simuleren. Hieronder valt vooral het sd-argument van de rnorm-functie op. Daar nemen we de standaardafwijking van 2 en vermenigvuldigen die met de absolute waarde van x, verheven tot de macht 1,5. Dit is vergelijkbaar met de manier waarop we hierboven gegevens hebben gesimuleerd om de functie varFixed te demonstreren. In dat geval namen we gewoon aan dat de exponent 0,5 was. (Herinneren we eraan dat de vierkantswortel nemen van een getal gelijk is aan het verheffen tot de macht 0,5.) Hier hebben we willekeurig een macht van 1,5 gekozen. Wanneer we gls met varPower gebruiken, proberen we de “ware” waarde van 1,5 te achterhalen, naast 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()
We zien dat de gegevens de klassieke vorm van variantie vertonen die toeneemt naarmate de voorspeller toeneemt. Om deze gegevens correct te modelleren met gls, voorzien wij deze van de formule y ~x en gebruiken wij het argument weights met varPower. Merk op dat we de eenzijdige formule net zo specificeren als we met de varFixed-functie hebben gedaan. In dit model krijgen we echter een schatting voor het vermogen.
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:"Het vermogen wordt geschat op 1,52, hetgeen zeer dicht bij de “ware” waarde van 1,5 ligt. De geschatte residuele standaardafwijking van 1,8729439 ligt ook vrij dicht bij 2. Beide intervallen komen overeen met de waarden die wij hebben gebruikt om de gegevens te simuleren. De coëfficiënten in het model worden daarentegen vrij slecht geschat. Dit is niet verwonderlijk gezien de grote variabiliteit in y, vooral voor x > 2.
varExp
Met de functie varExp kunnen we de variantie modelleren als een exponentiële functie van een covariaat. Ook deze variantiefunctie leggen we uit aan de hand van gesimuleerde gegevens. De enige verandering is in het sd argument van de rnorm functie. We hebben een vaste waarde van 2 die we vermenigvuldigen met x, die zelf weer wordt vermenigvuldigd met 0,5 en geëxponeerd. Merk op hoe snel de variantie toeneemt als we x verhogen. Dit komt door de exponentiële groei van de variantie.
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()
Om terug te werken en deze waarden terug te krijgen gebruiken we de functie varExp in het weights argument van gls. De eenzijdige formule verandert in dit geval niet. Het zegt modelleer de variantie als een functie van x. De functie varExp zegt dat x met een bepaalde waarde is vermenigvuldigd en geëxponeerd, zodat gls zal proberen die waarde te schatten.
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:"De “expon”-schatting van 0,4845623 in het gedeelte “variantiefunctie” ligt zeer dicht bij onze opgegeven waarde van 0,5. Ook de geschatte residuele standaardafwijking van 2,1191918 ligt dicht bij de “echte” waarde van 2. De schattingen van de modelcoëfficiënten liggen ook dicht bij de waarden die we hebben gebruikt om de gegevens te genereren, maar let op de onzekerheid in de (Intercept). De hypothesetest in de samenvatting kan een negatief intercept niet uitsluiten. Ook dit is niet verwonderlijk, aangezien y zoveel niet-constante variantie heeft en wij geen waarnemingen hebben van x gelijk aan 0. Aangezien het intercept de waarde is die y aanneemt wanneer x gelijk is aan 0, is ons geschatte intercept een extrapolatie naar een gebeurtenis die wij niet hebben waargenomen.
varConstPower
De varConstPower functie stelt ons in staat variantie te modelleren als een positieve constante plus een macht van de absolute waarde van een covariaat. Dat is een hele mond vol, maar dit is in feite hetzelfde als de varPower functie, behalve dat we er nu een positieve constante aan toevoegen. De volgende code simuleert gegevens waarvoor de varConstPower-functie geschikt zou zijn om te gebruiken. Merk op dat het identiek is aan de code die we gebruikten om gegevens te simuleren voor de varPower-sectie, behalve dat we 0,7 toevoegen aan x in de abs-functie. Waarom 0,7? Geen speciale reden. Het is gewoon een positieve constante die we hebben gekozen.
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()
De juiste manier om deze gegevens te modelleren is door gls te gebruiken met de varConstPower functie. De eenzijdige formule is dezelfde als voorheen. De samenvatting geeft drie schattingen voor het variantiemodel: de constante, de macht en de residuele standaardfout. Herinner u dat de “ware” waarden respectievelijk 0,7, 1,5 en 2 zijn.
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:"De intervallen geven de “ware” waarden weer, hoewel het interval voor de constante vrij breed is. Als we de werkelijke waarde niet zouden kennen, lijkt het aannemelijk dat de constante 2 of zelfs 4 zou kunnen zijn.
varComb
Ten slotte stelt de functie varComb ons in staat om combinaties van twee of meer variantiemodellen te modelleren door de overeenkomstige variantiefuncties met elkaar te vermenigvuldigen. Uiteraard kunnen hiermee zeer complexe variantiemodellen worden gemodelleerd. We zullen wat basisgegevens simuleren die geschikt zouden zijn om met de varComb-functie te gebruiken.
Wat we hieronder doen is twee verschillende variantieprocessen combineren:
- Een die standaarddeviatie laat verschillen tussen de geslachten (“f” = 2x, “m” = 2x 2x.5))
- Een die de standaardafwijking laat toenemen als
xtoeneemt, waarbijxmet 1,5 wordt vermenigvuldigd en geëxponeerd
Om de gegevens te helpen visualiseren, hebben we x beperkt tot een bereik van 1 tot 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()
De grafiek toont een toenemende variantie naarmate x toeneemt, maar ook verschillen in variantie tussen de geslachten. De variantie van y voor groep “m” is veel groter dan de variantie van y in groep “f”, vooral wanneer x groter is dan 1,5.
Om het hierboven gespecificeerde proces voor het genereren van gegevens correct te modelleren en te proberen de werkelijke waarden te achterhalen, gebruiken we de functie varComb als een omhulsel rond nog twee variantiefuncties: varIdent en varExp. Waarom deze twee? Omdat we verschillende varianties tussen de geslachten hebben, en omdat de variantie exponentieel toeneemt als functie van 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 residualDe samenvatting bevat twee secties voor het modelleren van de variantie. In de eerste wordt de multiplicator voor de “m”-groep op 2,437 geschat, hetgeen zeer dicht bij de werkelijke waarde van 2,5 ligt. De exponentiële parameter wordt geschat op 1,54, hetgeen zeer dicht ligt bij de werkelijke waarde van 1,5 die wij bij het genereren van de gegevens hebben gebruikt. Tenslotte wordt de residuele standaardafwijking geschat op ongeveer 1,79, dicht bij de werkelijke waarde van 2. Het oproepen van intervals(vm6) vertoont zeer nauwe betrouwbaarheidsintervallen. De voorvoegsels A en B in het varStruct element zijn slechts labels voor de twee verschillende variantiemodellen.
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:"Gevolge de grote exponentiële variabiliteit zijn de schattingen van de modelcoëfficiënten helaas bedroevend slecht.
Wat is het nut?
Waarom hebben we dan al die nepgegevens gesimuleerd en vervolgens geprobeerd de “ware” waarden te achterhalen? Wat heeft dat voor zin? Eerlijke vragen. Het antwoord is dat het ons helpt te begrijpen welke aannames we maken als we een model specificeren en passen. Als we het volgende model specificeren en fitten…
m <- lm(y ~ x1 + x2)… zeggen we dat we denken dat y ongeveer gelijk is aan een gewogen som van x1 en x2 (plus een intercept) met fouten als willekeurige trekkingen uit een normale verdeling met een gemiddelde van 0 en een vaste standaardafwijking. Evenzo, als we het volgende model specificeren en passen…
m <- gls(y ~ x1 + x2, data = d, weights = varPower(form = ~x1))…zeggen we dat we denken dat y ongeveer gelijk is aan een gewogen som van x1 en x2 (plus een intercept) met fouten als willekeurige trekkingen uit een normale verdeling met een gemiddelde van 0 en een standaardafwijking die groeit als een veelvoud van x1, vermeerderd met een bepaalde macht.
Als we gegevens kunnen simuleren die geschikt zijn voor die modellen, dan hebben we een beter begrip van de veronderstellingen die we maken als we die modellen gebruiken. Hopelijk hebt u nu een beter begrip van wat u kunt doen om variantie te modelleren met behulp van het nlme-pakket.
Voor vragen of verduidelijkingen over dit artikel kunt u contact opnemen met het UVA Library StatLab: [email protected]
Bekijk de hele collectie UVA Library StatLab-artikelen.
Clay Ford
Statistical Research Consultant
University of Virginia Library
April 7, 2020
-
1. Het
nlme-pakket is misschien beter bekend om zijnlme-functie, die wordt gebruikt om mixed-effect-modellen te fitten (dat wil zeggen, modellen met zowel vaste als willekeurige effecten). Deze blog post demonstreert variantie functies met behulp vangls, die geen random effecten fit. Alles wat we in deze blogpost presenteren, kan echter worden gebruikt metlme.