Una de las suposiciones básicas del modelado lineal es la varianza constante, u homogénea. ¿Qué significa eso exactamente? Vamos a simular algunos datos que satisfacen esta condición para ilustrar el concepto.
Abajo creamos un vector ordenado de números que van del 1 al 10 llamado x, y luego creamos un vector de números llamado y que es una función de x. Cuando representamos x frente a y, obtenemos una línea recta con una intercepción de 1,2 y una pendiente de 2,1.
x <- seq(1,10, length.out = 100)y <- 1.2 + 2.1 * xplot(x, y)
Ahora vamos a añadir algo de «ruido» a nuestros datos para que y no esté completamente determinado por x. Podemos hacerlo extrayendo aleatoriamente valores de una distribución Normal teórica con media 0 y alguna varianza establecida, y luego añadiéndolos a la fórmula que genera y. La función rnorm en R nos permite hacer esto fácilmente. A continuación extraemos 100 valores aleatorios de una distribución Normal con media 0 y desviación estándar 2 y los guardamos como un vector llamado noise. (Recordemos que la desviación estándar es simplemente la raíz cuadrada de la varianza). Luego generamos y con el ruido añadido. La función set.seed(222) permite obtener los mismos datos «aleatorios» en caso de que quieras seguirlos. Por último, volvemos a representar los datos.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2)y <- 1.2 + 2.1 * x + noiseplot(x, y)
Ahora tenemos datos que son una combinación de una componente lineal y determinista (\(y = 1,2 + 2,1x\)) y un ruido aleatorio extraído de una distribución \(N(0, 2)\N. Estas son las suposiciones básicas que hacemos sobre nuestros datos cuando ajustamos un modelo lineal tradicional. A continuación utilizamos la función lm para «recuperar» los valores «verdaderos» que utilizamos para generar los datos, de los cuales hay tres:
- El intercepto: 1,2
- La pendiente: 2.1
- La desviación estándar: 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-16Las estimaciones de (Intercept) y x (es decir, la pendiente) en la sección de Coeficientes, 1,27 y 2,09, son bastante cercanas a 1,2 y 2,1, respectivamente. El error estándar residual de 1,926 también está bastante cerca del valor constante de 2. Produjimos un «buen» modelo porque sabíamos cómo se generaba y y le dimos a la función lm el modelo «correcto» para ajustar. Aunque no podemos hacer esto en la vida real, es un ejercicio útil para ayudarnos a entender los supuestos de la modelización lineal.
¿Ahora qué pasaría si la varianza no fuera constante? Qué pasa si multiplicamos la desviación estándar de 2 por la raíz cuadrada de x? Como vemos en el gráfico de abajo, la dispersión vertical de los puntos aumenta a medida que x aumenta.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2 * sqrt(x))y <- 1.2 + 2.1 * x + noiseplot(x, y)
Multiplicamos 2 por sqrt(x) porque estamos especificando la desviación estándar. Si pudiéramos especificar la varianza, habríamos multiplicado 4 por sólo x.
Si ajustamos el mismo modelo utilizando lm, obtenemos un error estándar residual de 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-16Sabemos que no es correcto, porque hemos simulado los datos. No había una desviación estándar constante cuando creamos el «ruido». Cada valor aleatorio se extrajo de una distribución Normal diferente, cada una con media 0 pero con una desviación estándar que variaba según x. Esto significa que se viola nuestra suposición de varianza constante. ¿Cómo podríamos detectar esto en la vida real?
La forma más común es trazar los residuos frente a los valores ajustados. Esto es fácil de hacer en R. Simplemente llame a plot en el objeto modelo. Esto genera cuatro gráficos diferentes para evaluar los supuestos de modelado tradicionales. Consulte esta entrada del blog para obtener más información. Los gráficos que nos interesan son el primero y el tercero, que podemos especificar con el argumento which.
plot(m2, which = 1)
plot(m2, which = 3)
En el primer gráfico vemos que la variabilidad en torno a 0 aumenta a medida que nos desplazamos hacia la derecha con valores ajustados más grandes. En el tercer gráfico también vemos que la variabilidad aumenta a medida que nos desplazamos hacia la derecha, aunque esta vez los residuos han sido estandarizados y transformados a la escala de la raíz cuadrada. La trayectoria positiva de la línea roja suave indica un aumento de la varianza.
Así que ahora que hemos confirmado que nuestra suposición de varianza no constante no es válida, ¿qué podemos hacer? Un enfoque es transformar los datos en logaritmos. Esto a veces funciona si la variable de respuesta es positiva y está muy sesgada. Pero ese no es realmente el caso aquí. y sólo está ligeramente sesgada. (Llame a hist() en y para verificar.) Además sabemos que nuestros datos no están en una escala logarítmica.
Otro enfoque es modelar la varianza no constante. Ese es el tema de esta entrada del blog.
Para hacer esto, vamos a utilizar las funciones del paquete nlme que se incluye con la instalación base de R. La función de caballo de batalla es gls, que significa «mínimos cuadrados generalizados». La utilizamos igual que la función lm, excepto que también utilizamos el argumento weights junto con un puñado de funciones de varianza. Vamos a demostrarlo con un ejemplo.
A continuación ajustamos el modelo «correcto» a nuestros datos que presentaban una varianza no constante. Cargamos el paquete nlme para poder utilizar la función gls1. Especificamos la sintaxis del modelo como antes, y ~ x. Luego usamos el argumento weights para especificar la función de varianza, en este caso varFixed, también parte del paquete nlme. Esto dice que nuestra función de varianza no tiene parámetros y una sola covariable, x, que es precisamente como hemos generado los datos. La sintaxis, ~x, es una fórmula unilateral que puede leerse como «varianza del modelo en función de 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 residualEl resultado produce un error estándar residual de 1,925 que está muy cerca de 2, el valor «verdadero» que utilizamos para generar los datos. Los valores de la (Intercepción) y la pendiente, 1,37 y 2,09, también están muy cerca de 1,2 y 2,1.
Una vez más podemos llamar a plot en el objeto modelo. En este caso sólo se genera un gráfico: residuos estandarizados frente a los valores ajustados:
plot(vm1)
Este gráfico se ve bien. Siempre que modelemos nuestra varianza como una función de x, el modelo ajustado no se sobreajusta ni se infraajusta de ninguna manera sistemática (a diferencia de cuando usamos lm para ajustar el modelo y asumimos una varianza constante.)
La función varFixed crea pesos para cada observación, simbolizados como \(w_i). La idea es que cuanto mayor sea el peso de una determinada observación, menor será la varianza de la distribución Normal de la que se extrajo su componente de ruido. En nuestro ejemplo, a medida que x aumenta la varianza. Por lo tanto, los pesos más altos deben estar asociados con valores más pequeños de x. Esto puede lograrse tomando el recíproco de x, es decir, \(w_i = 1/x\). Así que cuando \(x = 2\), \(w = 1/2\). Cuando \(x = 10\), \(w = 1/10\). Los mayores x obtienen pesos más pequeños.
Por último, para asegurar que los pesos más grandes están asociados con varianzas más pequeñas, dividimos la varianza constante por el peso. Expresado matemáticamente,
\Nde la raíz cuadrada para expresar la desviación estándar,
\Nde la raíz cuadrada para expresar la desviación estándar,
Por lo tanto, cuanto mayor sea el denominador (es decir, cuanto mayor sea el peso y, por lo tanto, menor sea el x), menor será la varianza y más precisa será la y.
Por cierto, también podemos utilizar lm para ponderar las observaciones. También tiene un argumento weights como la función gls. La única diferencia es que tenemos que ser más explícito acerca de cómo expresar los pesos. En este caso, tenemos que especificar el recíproco de x. Observe que el resultado es casi idéntico al que obtenemos utilizando gls y 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-16El poder del paquete nlme es que permite una variedad de funciones de varianza. La función varFixed que acabamos de ilustrar es la más simple y algo que se puede hacer con lm. Las otras funciones de varianza incluyen:
varIdentvarPowervarExpvarConstPowervarComb
Exploremos cada una de ellas utilizando datos simulados.
varIdent
La función varIdent nos permite modelar diferentes varianzas por estrato. Para ver cómo funciona, primero simularemos los datos con esta propiedad. A continuación utilizamos set.seed(11) por si alguien quiere simular los mismos datos aleatorios. Luego ponemos n igual a 400, el número de observaciones que simularemos. x se genera igual que antes. En este ejemplo incluimos un predictor adicional llamado g que puede considerarse como el género. Muestreamos aleatoriamente 400 valores de «m» y «f». A continuación, simulamos un vector de dos desviaciones estándar y lo guardamos como msd. Si g == "m", entonces la desviación estándar es \(2 \times 2,5\). De lo contrario, es \(2\) para g == "f". Utilizamos este vector en la siguiente línea para generar y. Fíjate en que msd se introduce en la función rnorm. También observe cómo generamos y. Incluimos una interacción entre x y y. Cuando g == "f", el intercepto y la pendiente son 1,2 y 2,1. Cuando g == "m", el intercepto y la pendiente son (1,2 + 2,8) y (2,1 + 2,8), respectivamente. Por último, colocamos nuestros datos simulados en un marco de datos y trazamos con 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()
Nota las diferentes varianzas para cada nivel de g. La varianza para «m» es mucho mayor que la varianza para «f». Tiene mucha más dispersión que «f». Simulamos los datos de esta manera. Fijamos la varianza para «m» en 2,5 veces la de «f».
Ahora vamos a utilizar la función gls con varIdent para intentar recuperar estos valores reales. Usamos la misma forma que antes: definir nuestra función de varianza en el argumento weights. A continuación especificamos en el argumento form que la fórmula para modelar la varianza es condicional a g. La expresión ~ 1|g es una fórmula unilateral que dice que la varianza difiere entre los niveles de g. El 1 sólo significa que no tenemos covariables adicionales en nuestro modelo para la varianza. Es posible incluir una fórmula como ~ x|g, pero eso sería incorrecto en este caso ya que no utilizamos x en la generación de la varianza. Mirando el gráfico también muestra que mientras que la variabilidad en y es diferente entre los grupos, la variabilidad no aumenta a medida que x aumenta.
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 residualNote primero en la parte inferior de la salida que el error estándar residual estimado es 2.0053974, muy cerca del valor «verdadero» de 2. Observe también en la sección «Función de varianza» que obtenemos un valor estimado de 2,58127 para el grupo «m», que está muy cerca del valor «verdadero» de 2,5 que utilizamos para generar la varianza diferente para el grupo «m». En general, cuando se utiliza la función varIdent para estimar las diferentes varianzas entre los niveles de los estratos, uno de los niveles se establecerá en la línea de base, y los otros se estimarán como múltiplos del error estándar residual. En este caso, ya que «f» viene antes de «m» alfabéticamente, «f» se estableció en la línea de base, o 1. El error estándar residual estimado para el grupo «f» es \ (2,005397 \ por 1\). El error estándar residual estimado para el grupo «m» es \ (2,005397 \times 2,58127\).
Es importante tener en cuenta que estos son sólo estimaciones. Para tener una idea de la incertidumbre de estas estimaciones, podemos utilizar la función intervals del paquete nlme para obtener intervalos de confianza del 95%. Para reducir la salida, guardamos el resultado y vemos los elementos seleccionados del objeto.
int <- intervals(vm2)El elemento varStruct contiene el intervalo de confianza del 95% para el parámetro estimado en la función de varianza. El parámetro en este caso es el multiplicador del error estándar residual para el grupo de hombres.
int$varStruct## lower est. upper## m 2.245615 2.58127 2.967095## attr(,"label")## "Variance function:"El elemento sigma contiene el intervalo de confianza del 95% para el error estándar residual.
int$sigma## lower est. upper ## 1.818639 2.005397 2.211335 ## attr(,"label")## "Residual standard error:"Ambos intervalos son bastante pequeños y contienen el valor «verdadero» que utilizamos para generar los datos.
varPower
La función varPower nos permite modelar la varianza como una potencia del valor absoluto de una covariable. Una vez más, para ayudar a explicar esto vamos a simular datos con esta propiedad. A continuación lo principal a notar es el argumento sd de la función rnorm. Ahí es donde tomamos la desviación estándar de 2 y luego la multiplicamos por el valor absoluto de x elevado a la potencia de 1,5. Esto es similar a cómo simulamos los datos para demostrar la función varFixed anteriormente. En ese caso, simplemente asumimos que el exponente era 0,5. (Recordemos que sacar la raíz cuadrada de un número es equivalente a elevarlo a la potencia de 0,5). Aquí hemos elegido arbitrariamente una potencia de 1,5. Cuando utilicemos gls con varPower intentaremos recuperar el valor «verdadero» de 1,5 además de 2.
set.seed(4)n <- 1000x <- seq(1,10,length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*abs(x)^1.5)d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()
Vemos que los datos muestran la forma clásica de aumento de la varianza a medida que aumenta el predictor. Para modelar correctamente estos datos utilizando gls, le suministramos la fórmula y ~x y utilizamos el argumento weights con varPower. Observe que especificamos la fórmula unilateral igual que hicimos con la función varFixed. En este modelo, sin embargo, obtendremos una estimación de la potencia.
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:"Se estima que la potencia es de 1,52, que está muy cerca del valor «verdadero» de 1,5. El error estándar residual estimado de 1,8729439 también se acerca bastante a 2. Ambos intervalos recogen los valores que hemos utilizado para simular los datos. Por otro lado, los coeficientes del modelo están algo mal estimados. Esto no es sorprendente dada la cantidad de variabilidad en y, especialmente para \(x > 2\).
varExp
La función varExp nos permite modelar la varianza como una función exponencial de una covariable. Una vez más, explicaremos esta función de varianza utilizando datos simulados. El único cambio está en el argumento sd de la función rnorm. Tenemos un valor fijo de 2 que multiplicamos por x que a su vez se multiplica por 0,5 y se exponentiza. Fíjate en lo rápido que aumenta la varianza a medida que incrementamos x. Esto se debe al crecimiento exponencial de la varianza.
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()
Para trabajar hacia atrás y recuperar estos valores utilizamos la función varExp en el argumento weights de gls. La fórmula unilateral no cambia en este caso. Dice modelar la varianza como una función de x. La función varExp dice que x se ha multiplicado por algún valor y se ha exponenciado, por lo que gls tratará de estimar ese valor.
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:"La estimación «expon» de 0,4845623 en la sección «función de varianza» está muy cerca de nuestro valor especificado de 0,5. Asimismo, el error estándar residual estimado de 2,1191918 se aproxima al valor «verdadero» de 2. Las estimaciones de los coeficientes del modelo también se aproximan a los valores que utilizamos para generar los datos, pero observe la incertidumbre en el (Intercepto). La prueba de hipótesis en el resumen no puede descartar un intercepto negativo. Esto, de nuevo, no es sorprendente ya que y tiene mucha varianza no constante, así como el hecho de que no tenemos observaciones de x igual a 0. Dado que el intercepto es el valor que toma y cuando x es igual a 0, nuestro intercepto estimado es una extrapolación a un evento que no observamos.
varConstPower
La función varConstPower nos permite modelar la varianza como una constante positiva más una potencia del valor absoluto de una covariable. Es un trabalenguas, pero es básicamente lo mismo que la función varPower excepto que ahora le añadimos una constante positiva. El siguiente código simula los datos para los que sería conveniente utilizar la función varConstPower. Observe que es idéntico al código que utilizamos para simular los datos de la sección varPower, excepto que añadimos 0,7 a x en la función abs. ¿Por qué 0,7? No hay ninguna razón especial. Es simplemente una constante positiva que elegimos.
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()
La forma correcta de modelar estos datos es utilizar gls con la función varConstPower. La fórmula unilateral es la misma que antes. El resumen devuelve tres estimaciones para el modelo de varianza: la constante, la potencia y el error estándar residual. Recuerde que los valores «verdaderos» son 0,7, 1,5 y 2, respectivamente.
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:"Los intervalos capturan los valores «verdaderos», aunque el intervalo para la constante es bastante amplio. Si no conociéramos el valor verdadero, parecería que la constante podría ser 2 o incluso 4.
varComb
Por último, la función varComb nos permite modelar combinaciones de dos o más modelos de varianza multiplicando juntas las funciones de varianza correspondientes. Obviamente, esto puede acomodar algunos modelos de varianza muy complejos. Simularemos algunos datos básicos que serían apropiados para usar con la función varComb.
Lo que hacemos a continuación es combinar dos procesos de varianza diferentes:
- Uno que permite que la desviación estándar difiera entre géneros («f» = \(2\), «m» = \(2 \times 2.5\))
- Una que permite que la desviación estándar aumente a medida que
xaumenta, dondexse multiplica por 1,5 y se exponentiza
Para ayudar a visualizar los datos hemos limitado x a un rango de 1 a 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()
El gráfico muestra el aumento de la varianza a medida que x aumenta, pero también las diferencias de varianza entre géneros. La varianza de y para el grupo «m» es mucho mayor que la varianza de y en el grupo «f», especialmente cuando x es mayor que 1,5.
Para modelar correctamente el proceso de generación de datos que especificamos anteriormente e intentar recuperar los valores verdaderos, utilizamos la función varComb como envoltorio de otras dos funciones de varianza: varIdent y varExp. ¿Por qué estas dos? Porque tenemos varianzas diferentes entre géneros, y tenemos varianza que aumenta exponencialmente en función de x.
vm6 <- gls(y ~ x*g, data = d, weights = varComb(varIdent(form = ~ 1|g), varExp(form = ~x)))summary(vm6)## Generalized least squares fit by REML## Model: y ~ x * g ## Data: d ## AIC BIC logLik## 4431.45 4459.32 -2208.725## ## Combination of variance functions: ## Structure: Different standard deviations per stratum## Formula: ~1 | g ## Parameter estimates:## f m ## 1.000000 2.437046 ## Structure: Exponential of variance covariate## Formula: ~x ## Parameter estimates:## expon ## 1.540608 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) -5.996337 6.539552 -0.9169339 0.3597## x 8.829971 4.849272 1.8208860 0.0694## gm -17.572238 16.548540 -1.0618603 0.2889## x:gm 16.932938 12.192793 1.3887661 0.1657## ## Correlation: ## (Intr) x gm ## x -0.972 ## gm -0.395 0.384 ## x:gm 0.387 -0.398 -0.974## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.74441479 -0.67478610 -0.00262221 0.66079254 3.14975288 ## ## Residual standard error: 1.794185 ## Degrees of freedom: 400 total; 396 residualLa sección de resumen contiene dos secciones para modelar la varianza. En la primera se estima que el multiplicador del grupo «m» es de 2,437, lo que se aproxima mucho al valor real de 2,5. El parámetro exponencial se estima en 1,54, muy cercano al valor real de 1,5 que utilizamos al generar los datos. Por último, el error estándar residual se estima en 1,79, cerca del valor real de 2. La llamada intervals(vm6) muestra intervalos de confianza muy ajustados. Los prefijos A y B en el elemento varStruct son sólo etiquetas para los dos modelos de varianza diferentes.
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:"Desgraciadamente, debido a la gran variabilidad exponencial, las estimaciones de los coeficientes del modelo son lamentablemente malas.
¿Qué sentido tiene?
Entonces, ¿por qué simulamos todos estos datos falsos y luego intentamos recuperar los valores «verdaderos»? ¿De qué sirve eso? Es una pregunta justa. La respuesta es que nos ayuda a entender qué suposiciones estamos haciendo cuando especificamos y ajustamos un modelo. Cuando especificamos y ajustamos el siguiente modelo…
m <- lm(y ~ x1 + x2)…estamos diciendo que pensamos que y es aproximadamente igual a una suma ponderada de x1 y x2 (más un intercepto) con errores que son extracciones aleatorias de una distribución Normal con media de 0 y alguna desviación estándar fija. Del mismo modo, si especificamos y ajustamos el siguiente modelo…
m <- gls(y ~ x1 + x2, data = d, weights = varPower(form = ~x1))… estamos diciendo que pensamos que y es aproximadamente igual a una suma ponderada de x1 y x2 (más un intercepto) con errores que son extracciones aleatorias de una distribución Normal con media de 0 y una desviación estándar que crece como un múltiplo de x1 elevado por alguna potencia.
Si podemos simular datos adecuados para esos modelos, entonces tenemos una mejor comprensión de los supuestos que estamos haciendo cuando usamos esos modelos. Esperemos que ahora tenga una mejor comprensión de lo que puede hacer para modelar la varianza utilizando el paquete nlme.
Para preguntas o aclaraciones sobre este artículo, póngase en contacto con el StatLab de la Biblioteca de la UVA: [email protected]
Vea toda la colección de artículos del StatLab de la Biblioteca de la UVA.
Clay Ford
Consultor de Investigación Estadística
Biblioteca de la Universidad de Virginia
7 de abril de 2020
-
1. El paquete
nlmees quizás más conocido por su funciónlme, que se utiliza para ajustar modelos de efectos mixtos (es decir, modelos con efectos fijos y aleatorios). Esta entrada de blog demuestra las funciones de varianza utilizandogls, que no ajusta efectos aleatorios. Sin embargo, todo lo que presentamos en esta entrada del blog se puede utilizar conlme.