Eine der Grundannahmen der linearen Modellierung ist konstante oder homogene Varianz. Was bedeutet das genau? Simulieren wir einige Daten, die diese Bedingung erfüllen, um das Konzept zu veranschaulichen.
Nachfolgend erstellen wir einen sortierten Zahlenvektor von 1 bis 10 mit der Bezeichnung x und erstellen dann einen Zahlenvektor mit der Bezeichnung y, der eine Funktion von x ist. Wenn wir x gegen y auftragen, erhalten wir eine gerade Linie mit einem Achsenabschnitt von 1,2 und einer Steigung von 2,1.
x <- seq(1,10, length.out = 100)y <- 1.2 + 2.1 * xplot(x, y)
Nun fügen wir etwas „Rauschen“ zu unseren Daten hinzu, damit y nicht vollständig von x bestimmt wird. Dazu ziehen wir nach dem Zufallsprinzip Werte aus einer theoretischen Normalverteilung mit Mittelwert 0 und einer bestimmten Varianz und fügen sie dann in die Formel ein, die y erzeugt. Mit der Funktion rnorm in R lässt sich dies leicht bewerkstelligen. Im Folgenden ziehen wir 100 Zufallswerte aus einer Normalverteilung mit Mittelwert 0 und Standardabweichung 2 und speichern sie als Vektor mit dem Namen noise. (Die Standardabweichung ist einfach die Quadratwurzel der Varianz.) Dann erzeugen wir y mit dem hinzugefügten Rauschen. Mit der Funktion set.seed(222) können Sie dieselben „zufälligen“ Daten erhalten, falls Sie dies nachvollziehen möchten. Schließlich stellen wir die Daten erneut dar.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2)y <- 1.2 + 2.1 * x + noiseplot(x, y)
Jetzt haben wir Daten, die eine Kombination aus einer linearen, deterministischen Komponente (\(y = 1,2 + 2,1x\)) und zufälligem Rauschen aus einer \(N(0, 2)\)-Verteilung sind. Dies sind die Grundannahmen, die wir über unsere Daten machen, wenn wir ein traditionelles lineares Modell anpassen. Nachfolgend verwenden wir die Funktion lm, um die „wahren“ Werte wiederherzustellen, die wir zur Generierung der Daten verwendet haben, von denen es drei gibt:
- Der Achsenabschnitt: 1.2
- Die Steigung: 2.1
- Die Standardabweichung: 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-16Die (Intercept) und x (d.h. die Steigung) Schätzungen im Abschnitt Koeffizienten, 1,27 und 2,09, liegen recht nahe bei 1,2 bzw. 2,1. Der Reststandardfehler von 1,926 liegt ebenfalls recht nahe am konstanten Wert von 2. Wir haben ein „gutes“ Modell erstellt, weil wir wussten, wie y generiert wurde und der Funktion lm das „richtige“ Modell zur Anpassung gegeben haben. Auch wenn dies im wirklichen Leben nicht möglich ist, so ist es doch eine nützliche Übung, die uns hilft, die Annahmen der linearen Modellierung zu verstehen.
Was wäre nun, wenn die Varianz nicht konstant wäre? Was wäre, wenn wir die Standardabweichung von 2 mit der Quadratwurzel von x multiplizieren würden? Wie wir in der folgenden Grafik sehen, nimmt die vertikale Streuung der Punkte mit zunehmender x zu.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2 * sqrt(x))y <- 1.2 + 2.1 * x + noiseplot(x, y)
Wir multiplizieren 2 mit sqrt(x), weil wir die Standardabweichung angeben. Wenn wir die Varianz angeben könnten, hätten wir 4 nur mit x multipliziert.
Wenn wir dasselbe Modell mit lm anpassen, erhalten wir einen Reststandardfehler von 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-16Wir wissen, dass das nicht richtig ist, weil wir die Daten simuliert haben. Es gab keine konstante Standardabweichung, als wir das „Rauschen“ erzeugten. Jeder Zufallswert wurde aus einer anderen Normalverteilung gezogen, jeweils mit einem Mittelwert von 0, aber einer Standardabweichung, die entsprechend x variierte. Das bedeutet, dass unsere Annahme einer konstanten Varianz verletzt ist. Wie würden wir dies im wirklichen Leben feststellen?
Die gebräuchlichste Methode ist das Auftragen der Residuen gegen die angepassten Werte. Dies ist in R leicht zu machen. Rufen Sie einfach plot für das Modellobjekt auf. Dies erzeugt vier verschiedene Diagramme, um die traditionellen Modellierungsannahmen zu bewerten. Weitere Informationen finden Sie in diesem Blogbeitrag. Die Diagramme, die uns interessieren, sind das erste und das dritte Diagramm, die wir mit dem Argument which angeben können.
plot(m2, which = 1)
plot(m2, which = 3)
Im ersten Diagramm sehen wir, dass die Variabilität um 0 herum zunimmt, wenn wir uns mit größeren angepassten Werten weiter nach rechts bewegen. Im dritten Diagramm sehen wir ebenfalls eine zunehmende Variabilität, wenn wir uns nach rechts bewegen, obwohl dieses Mal die Residuen standardisiert und auf die Quadratwurzelskala transformiert wurden. Der positive Verlauf der glatten roten Linie deutet auf eine Zunahme der Varianz hin.
Nachdem wir nun bestätigt haben, dass unsere Annahme einer nicht konstanten Varianz ungültig ist, was können wir tun? Ein Ansatz ist die Log-Transformation der Daten. Das funktioniert manchmal, wenn die Antwortvariable positiv und stark verzerrt ist. Aber das ist hier nicht der Fall. y ist nur leicht verzerrt. (Rufen Sie hist() zu y auf, um das zu überprüfen.) Außerdem wissen wir, dass unsere Daten nicht auf einer logarithmischen Skala liegen.
Ein anderer Ansatz ist die Modellierung der nicht konstanten Varianz. Das ist das Thema dieses Blogbeitrags.
Dazu verwenden wir Funktionen aus dem Paket nlme, das in der Basisinstallation von R enthalten ist. Die wichtigste Funktion ist gls, die für „generalisierte kleinste Quadrate“ steht. Wir verwenden sie genauso wie die Funktion lm, mit dem Unterschied, dass wir auch das Argument weights zusammen mit einer Handvoll Varianzfunktionen verwenden. Wir wollen das an einem Beispiel demonstrieren.
Nachfolgend passen wir das „richtige“ Modell an unsere Daten an, die eine nicht konstante Varianz aufweisen. Wir laden das nlme-Paket, damit wir die gls-Funktion1 verwenden können. Wir geben die Modellsyntax wie zuvor an, y ~ x. Dann verwenden wir das Argument weights, um die Varianzfunktion anzugeben, in diesem Fall varFixed, ebenfalls Teil des Pakets nlme. Dies bedeutet, dass unsere Varianzfunktion keine Parameter und eine einzige Kovariate, x, hat, was genau der Art und Weise entspricht, wie wir die Daten generiert haben. Die Syntax ~x ist eine einseitige Formel, die als „Modellvarianz als Funktion von x“ gelesen werden kann.
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 residualDas Ergebnis ergibt einen Reststandardfehler von 1,925, der sehr nahe an 2 liegt, dem „wahren“ Wert, den wir zur Generierung der Daten verwendet haben. Die Werte für (Achsenabschnitt) und Steigung, 1,37 und 2,09, liegen ebenfalls sehr nahe bei 1,2 und 2,1.
Wir können erneut plot für das Modellobjekt aufrufen. In diesem Fall wird nur eine Darstellung erzeugt: standardisierte Residuen gegen angepasste Werte:
plot(vm1)
Diese Darstellung sieht gut aus. Solange wir unsere Varianz als eine Funktion von x modellieren, über- oder unterschreitet das angepasste Modell weder systematisch (anders als bei der Verwendung von lm zur Anpassung des Modells und der Annahme einer konstanten Varianz)
Die Funktion varFixed erstellt Gewichte für jede Beobachtung, symbolisiert als \(w_i\). Je höher die Gewichtung einer Beobachtung ist, desto geringer ist die Varianz der Normalverteilung, aus der ihre Rauschkomponente gezogen wurde. In unserem Beispiel nimmt die Varianz mit x zu. Daher sollten höhere Gewichte mit kleineren x-Werten verbunden sein. Dies kann erreicht werden, indem man den Kehrwert von x nimmt, d. h. \(w_i = 1/x\). Wenn also \(x = 2\), \(w = 1/2\). Wenn \(x = 10\), \(w = 1/10\). Größere x erhalten kleinere Gewichte.
Um schließlich sicherzustellen, dass größere Gewichte mit kleineren Varianzen verbunden sind, teilen wir die konstante Varianz durch das Gewicht. Mathematisch ausgedrückt,
\
oder durch Ziehen der Quadratwurzel, um die Standardabweichung auszudrücken,
\
Je größer also der Nenner (d.h. je größer die Gewichtung und damit kleiner die x), desto kleiner die Varianz und desto präziser die beobachtete y.
Inzidentell können wir auch lm zur Gewichtung von Beobachtungen verwenden. Auch sie hat ein weights Argument wie die gls Funktion. Der einzige Unterschied besteht darin, dass wir die Gewichtung expliziter ausdrücken müssen. In diesem Fall müssen wir den Kehrwert von x angeben. Das Ergebnis ist fast identisch mit dem, das wir mit gls und varFixed erhalten.
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-16Die Stärke des Pakets nlme besteht darin, dass es eine Vielzahl von Varianzfunktionen ermöglicht. Die varFixed-Funktion, die wir gerade illustriert haben, ist die einfachste und kann mit lm ausgeführt werden. Die anderen Varianzfunktionen sind:
varIdentvarPowervarExpvarConstPowervarComb
Lassen Sie uns jede dieser Funktionen anhand von simulierten Daten untersuchen.
varIdent
Die varIdent Funktion erlaubt uns, verschiedene Varianzen pro Schicht zu modellieren. Um zu sehen, wie sie funktioniert, simulieren wir zunächst Daten mit dieser Eigenschaft. Im Folgenden verwenden wir set.seed(11), falls jemand die gleichen Zufallsdaten simulieren möchte. Dann setzen wir n auf 400, die Anzahl der Beobachtungen, die wir simulieren werden. x wird auf die gleiche Weise erzeugt wie zuvor. In diesem Beispiel fügen wir einen zusätzlichen Prädiktor namens g hinzu, den man sich als Geschlecht vorstellen kann. Es werden 400 Werte von „m“ und „f“ zufällig ausgewählt. Als nächstes simulieren wir einen Vektor mit zwei Standardabweichungen und speichern ihn als msd. Wenn g == "m", dann ist die Standardabweichung \(2 \mal 2,5\). Ansonsten ist es \(2\) für g == "f". Wir verwenden diesen Vektor in der nächsten Zeile, um y zu erzeugen. Beachten Sie, wo msd in die Funktion rnorm eingefügt wird. Beachten Sie auch, wie wir y erzeugen. Wir schließen eine Wechselwirkung zwischen x und y ein. Bei g == "f" sind der Achsenabschnitt und die Steigung 1,2 und 2,1. Bei g == "m" sind der Achsenabschnitt und die Steigung (1,2 + 2,8) bzw. (2,1 + 2,8). Abschließend werden die simulierten Daten in einen Datenrahmen eingefügt und mit ggplot2 dargestellt.
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()
Beachten Sie die unterschiedlichen Varianzen für jedes Niveau von g. Die Varianz für „m“ ist viel größer als die Varianz für „f“. Sie hat eine viel größere Streuung als „f“. Wir haben die Daten auf diese Weise simuliert. Wir setzen die Varianz für „m“ auf das 2,5-fache von „f“.
Nun wollen wir die Funktion gls mit varIdent verwenden, um die wahren Werte zu ermitteln. Wir verwenden die gleiche Methode wie zuvor: Wir definieren unsere Varianzfunktion im Argument weights. Im Folgenden geben wir im Argument form an, dass die Formel für die Modellierung der Varianz von g abhängig ist. Der Ausdruck ~ 1|g ist eine einseitige Formel, die besagt, dass sich die Varianz zwischen den Stufen von g unterscheidet. Die 1 bedeutet lediglich, dass wir keine zusätzlichen Kovariaten in unserem Varianzmodell haben. Es ist möglich, eine Formel wie ~ x|g einzubeziehen, aber das wäre in diesem Fall nicht korrekt, da wir x bei der Erstellung der Varianz nicht verwendet haben. Ein Blick auf das Diagramm zeigt auch, dass die Variabilität in y zwar zwischen den Gruppen unterschiedlich ist, die Variabilität aber nicht zunimmt, wenn x zunimmt.
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 residualZunächst fällt unten in der Ausgabe auf, dass der geschätzte Reststandardfehler 2.0053974 beträgt, was dem „wahren“ Wert von 2 sehr nahe kommt. Beachten Sie auch, dass wir im Abschnitt „Varianzfunktion“ einen geschätzten Wert von 2,58127 für die Gruppe „m“ erhalten, was dem „wahren“ Wert von 2,5 sehr nahe kommt, den wir zur Erzeugung der unterschiedlichen Varianz für die Gruppe „m“ verwendet haben. Im Allgemeinen wird bei der Verwendung der Funktion varIdent zur Schätzung unterschiedlicher Varianzen zwischen Ebenen von Schichten eine der Ebenen auf die Grundlinie gesetzt, und die anderen werden als Vielfache des Reststandardfehlers geschätzt. Da in diesem Fall „f“ alphabetisch vor „m“ steht, wurde „f“ auf die Grundlinie, also 1, gesetzt. Der geschätzte Reststandardfehler für die Gruppe „f“ ist \(2,005397 \mal 1\). Der geschätzte Reststandardfehler für die Gruppe „m“ ist \(2,005397 \mal 2,58127\).
Es ist wichtig zu beachten, dass es sich hierbei nur um Schätzungen handelt. Um ein Gefühl für die Unsicherheit dieser Schätzungen zu bekommen, können wir die Funktion intervals aus dem Paket nlme verwenden, um 95%-Konfidenzintervalle zu erhalten. Um die Ausgabe zu reduzieren, speichern wir das Ergebnis und betrachten ausgewählte Elemente des Objekts.
int <- intervals(vm2)Das Element varStruct enthält das 95%-Konfidenzintervall für den in der Varianzfunktion geschätzten Parameter. Der Parameter ist in diesem Fall der Reststandardfehler-Multiplikator für die männliche Gruppe.
int$varStruct## lower est. upper## m 2.245615 2.58127 2.967095## attr(,"label")## "Variance function:"Das Element sigma enthält das 95%-Konfidenzintervall für den Reststandardfehler.
int$sigma## lower est. upper ## 1.818639 2.005397 2.211335 ## attr(,"label")## "Residual standard error:"Beide Intervalle sind recht klein und enthalten den „wahren“ Wert, den wir zur Generierung der Daten verwendet haben.
varPower
Die Funktion varPower ermöglicht es uns, die Varianz als Potenz des absoluten Werts einer Kovariate zu modellieren. Um dies zu erklären, werden wir erneut Daten mit dieser Eigenschaft simulieren. Im Folgenden ist vor allem das sd-Argument der rnorm-Funktion zu beachten. Hier nehmen wir die Standardabweichung von 2 und multiplizieren sie mit dem absoluten Wert von x, der mit 1,5 potenziert wird. Das ist ähnlich wie bei der Simulation von Daten zur Veranschaulichung der Funktion varFixed oben. In diesem Fall haben wir einfach angenommen, dass der Exponent 0,5 ist. (Zur Erinnerung: Die Quadratwurzel einer Zahl ist gleichbedeutend mit der Potenzierung durch 0,5). Hier haben wir willkürlich eine Potenz von 1,5 gewählt. Wenn wir gls mit varPower verwenden, versuchen wir, den „wahren“ Wert von 1,5 zusätzlich zu 2 zu erhalten.
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()
Wir sehen, dass die Daten die klassische Form der Varianz aufweisen, die mit dem Anstieg des Prädiktors zunimmt. Um diese Daten mit gls korrekt zu modellieren, geben wir ihm die Formel y ~x und verwenden das Argument weights mit varPower. Beachten Sie, dass wir die einseitige Formel genauso angeben wie bei der Funktion varFixed. In diesem Modell erhalten wir jedoch eine Schätzung für die Potenz.
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:"Die Potenz wird auf 1,52 geschätzt, was dem „wahren“ Wert von 1,5 sehr nahe kommt. Der geschätzte Standardfehler der Residuen von 1,8729439 liegt ebenfalls recht nahe bei 2. Beide Intervalle erfassen die Werte, die wir zur Simulation der Daten verwendet haben. Die Koeffizienten des Modells sind dagegen etwas schlecht geschätzt. Dies ist nicht überraschend, wenn man bedenkt, wie groß die Variabilität in y ist, insbesondere für \(x > 2\).
varExp
Die Funktion varExp ermöglicht es uns, die Varianz als Exponentialfunktion einer Kovariate zu modellieren. Auch diese Varianzfunktion wird anhand von simulierten Daten erklärt. Die einzige Änderung betrifft das sd-Argument der rnorm-Funktion. Wir haben einen festen Wert von 2, den wir mit x multiplizieren, das wiederum mit 0,5 multipliziert und potenziert wird. Beachten Sie, wie schnell die Varianz zunimmt, wenn wir x erhöhen. Dies ist auf das exponentielle Wachstum der Varianz zurückzuführen.
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()
Um rückwärts zu arbeiten und diese Werte wiederherzustellen, verwenden wir die Funktion varExp im Argument weights von gls. Die einseitige Formel ändert sich in diesem Fall nicht. Sie besagt, dass die Varianz als eine Funktion von x modelliert wird. Die Funktion varExp besagt, dass x mit einem bestimmten Wert multipliziert und potenziert wurde, so dass gls versuchen wird, diesen Wert zu schätzen.
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:"Der „expon“-Schätzwert von 0,4845623 im Abschnitt „Varianzfunktion“ liegt sehr nahe an unserem angegebenen Wert von 0,5. Ebenso liegt der geschätzte Standardfehler der Residuen von 2,1191918 nahe am „wahren“ Wert von 2. Die Schätzungen der Modellkoeffizienten liegen ebenfalls nahe an den Werten, die wir zur Generierung der Daten verwendet haben, aber beachten Sie die Unsicherheit beim (Intercept). Der Hypothesentest in der Zusammenfassung kann einen negativen Achsenabschnitt nicht ausschließen. Auch dies ist nicht überraschend, da y so viel nicht konstante Varianz hat und wir keine Beobachtungen von x gleich 0 haben. Da der Achsenabschnitt der Wert ist, den y annimmt, wenn x gleich 0 ist, ist unser geschätzter Achsenabschnitt eine Extrapolation auf ein Ereignis, das wir nicht beobachtet haben.
varConstPower
Die Funktion varConstPower ermöglicht es uns, die Varianz als positive Konstante plus eine Potenz des absoluten Werts einer Kovariate zu modellieren. Das ist ein langes Wort, aber es ist im Grunde dasselbe wie die Funktion varPower, nur dass wir jetzt eine positive Konstante hinzufügen. Der folgende Code simuliert Daten, für die die Funktion varConstPower geeignet wäre. Er ist identisch mit dem Code, den wir zur Simulation von Daten für den Abschnitt varPower verwendet haben, außer dass wir in der Funktion abs 0,7 zu x hinzufügen. Warum 0,7? Dafür gibt es keinen besonderen Grund. Es ist einfach eine positive Konstante, die wir gewählt haben.
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()
Der korrekte Weg, diese Daten zu modellieren, ist die Verwendung von gls mit der Funktion varConstPower. Die einseitige Formel ist die gleiche wie zuvor. Die Zusammenfassung liefert drei Schätzungen für das Varianzmodell: die Konstante, die Potenz und den Reststandardfehler. Die „wahren“ Werte sind 0,7, 1,5 bzw. 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:"Die Intervalle erfassen die „wahren“ Werte, wobei das Intervall für die Konstante recht breit ist. Wenn wir den wahren Wert nicht kennen würden, könnte die Konstante plausibelerweise 2 oder sogar 4 sein.
varComb
Schließlich ermöglicht die Funktion varComb die Modellierung von Kombinationen aus zwei oder mehr Varianzmodellen, indem die entsprechenden Varianzfunktionen miteinander multipliziert werden. Auf diese Weise lassen sich natürlich sehr komplexe Varianzmodelle erstellen. Wir werden einige grundlegende Daten simulieren, die für die Verwendung der Funktion varComb geeignet wären.
Was wir im Folgenden tun, ist die Kombination von zwei verschiedenen Varianzprozessen:
- einer, der die Standardabweichung zwischen den Geschlechtern unterscheidet („f“ = \(2\), „m“ = \(2 \mal 2.5\))
- Eine, bei der die Standardabweichung mit zunehmender
xzunimmt, wobeixmit 1,5 multipliziert und potenziert wird
Um die Visualisierung der Daten zu erleichtern, haben wir x auf einen Bereich von 1 bis 3 begrenzt.
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()
Das Diagramm zeigt eine zunehmende Varianz, wenn x zunimmt, aber auch Unterschiede in der Varianz zwischen den Geschlechtern. Die Varianz von y für die Gruppe „m“ ist viel größer als die Varianz von y in der Gruppe „f“, vor allem, wenn x größer als 1,5 ist.
Um den oben spezifizierten Datenerzeugungsprozess korrekt zu modellieren und zu versuchen, die wahren Werte wiederherzustellen, verwenden wir die Funktion varComb als Wrapper um zwei weitere Varianzfunktionen: varIdent und varExp. Warum diese beiden? Weil wir unterschiedliche Varianzen zwischen den Geschlechtern haben und die Varianz exponentiell als Funktion von x ansteigt.
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 residualDer Abschnitt Zusammenfassung enthält zwei Abschnitte zur Modellierung der Varianz. Im ersten wird der Multiplikator für die Gruppe „m“ auf 2,437 geschätzt, was dem wahren Wert von 2,5 sehr nahe kommt. Der Exponentialparameter wird auf 1,54 geschätzt, was sehr nahe an dem wahren Wert von 1,5 liegt, den wir bei der Erstellung der Daten verwendet haben. Schließlich wird der Standardfehler der Residuen auf etwa 1,79 geschätzt, was dem wahren Wert von 2 nahe kommt. Der Aufruf intervals(vm6) zeigt sehr enge Konfidenzintervalle. Die Präfixe A und B im Element varStruct sind nur Bezeichnungen für die beiden verschiedenen Varianzmodelle.
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:"Die Schätzungen der Modellkoeffizienten sind aufgrund der großen exponentiellen Variabilität leider sehr schlecht.
Was soll das bringen?
Warum haben wir also all diese falschen Daten simuliert und dann versucht, „wahre“ Werte zu ermitteln? Wozu ist das gut? Eine berechtigte Frage. Die Antwort ist, dass es uns hilft zu verstehen, welche Annahmen wir machen, wenn wir ein Modell spezifizieren und anpassen. Wenn wir das folgende Modell spezifizieren und anpassen…
m <- lm(y ~ x1 + x2)…sagen wir damit, dass wir glauben, dass y ungefähr gleich einer gewichteten Summe von x1 und x2 (plus einem Achsenabschnitt) ist, wobei die Fehler zufällige Ziehungen aus einer Normalverteilung mit einem Mittelwert von 0 und einer festen Standardabweichung sind. Wenn wir das folgende Modell spezifizieren und anpassen…
m <- gls(y ~ x1 + x2, data = d, weights = varPower(form = ~x1))…sagen wir, dass wir glauben, dass y ungefähr gleich einer gewichteten Summe von x1 und x2 (plus einem Achsenabschnitt) ist, wobei die Fehler zufällige Ziehungen aus einer Normalverteilung mit einem Mittelwert von 0 und einer Standardabweichung sind, die als ein Vielfaches von x1, erhöht um eine Potenz, wächst.
Wenn wir Daten simulieren können, die für diese Modelle geeignet sind, dann haben wir ein besseres Verständnis für die Annahmen, die wir machen, wenn wir diese Modelle verwenden. Hoffentlich haben Sie jetzt ein besseres Verständnis dafür, was Sie tun können, um Varianz mit dem nlme-Paket zu modellieren.
Für Fragen oder Klarstellungen zu diesem Artikel wenden Sie sich bitte an das StatLab der UVA-Bibliothek: [email protected]
Sehen Sie sich die gesamte Sammlung der StatLab-Artikel der UVA-Bibliothek an.
Clay Ford
Statistical Research Consultant
University of Virginia Library
April 7, 2020
-
1. Das
nlme-Paket ist vielleicht besser bekannt für seinelme-Funktion, die zur Anpassung von Modellen mit gemischten Effekten (d.h. Modelle mit sowohl festen als auch zufälligen Effekten) verwendet wird. In diesem Blogbeitrag werden Varianzfunktionen mitglsdemonstriert, das keine zufälligen Effekte anpasst. Alles, was wir in diesem Blogbeitrag vorstellen, kann jedoch mitlmeverwendet werden.