Modelos lineales generalizados

Ajuste de un Modelo Lineal usando las Funciones de R: lm() y glm()

logo

Especialización en Estadística Aplicada

Roberto Trespalacios

Modelos lineales con generalizados R

Recuerde que un modelo lineal de la forma \( y_i = \boldsymbol{x}'_i \boldsymbol{\beta} + \boldsymbol{\varepsilon}_i \) tiene una estructura de la media

\[ E(y_i|\boldsymbol{x}_i) =\boldsymbol{x}'_i\boldsymbol{\beta} \]

Es decir, debemos especificar que variables y términos (por ejemplo, interacciones) serán incluidos en la media del modelo.
De otro lado, debemos especificar una estructura de varianza para los residuales ya que asumimos un modelo lineal clásico (independencia entre los errores):
\( \varepsilon_i \stackrel{i.i.d}{\sim} N(0, \sigma^2) \). Más adelante nos referiremos a la estructura de correlación de los errores en modelos más complejos.
Estas dos componentes estructura de media + estructura de los errores más los supuestos distribución definen el modelo.

Ejemplo: Nacimientos en PR en el 2011

En esta sección vamos a utilizar una muestra aleatoria de 982 nacimientos en Puerto Rico en el 2011 (newbornpr.csv). Vamos a determinar las variables relevantes para explicar el peso de un recién nacido. La base de datos newbornpr2011.csv esta disponible en la página del curso. Una descripción de las variables esta disponible en la misma página.
Vamos a ver la estructura de la base de datos. La salida de R nos indica que tenemos tanto variables cuantitativas como cualitativas. La variable respuesta en este problema es el peso (gramos) del recién nacido (DBWT).

# Estructura de la base de datos
library(readr)
newbornpr<- read_csv("~/Desktop/libro_glm/datos/newbornpr2011.csv")
dim(newbornpr) # dimension

[1] 982  17

names(newbornpr) # variables

 [1] "X"        "MAGER"    "MRACEREC" "MAR"      "MEDUC"    "PRECARE" 
 [7] "WTGAIN"   "CIG_REC"  "PGWT"     "DWGT"     "SEX"      "DBWT"    
[13] "RF_DIAB"  "RF_GEST"  "RF_PHYP"  "RF_GHYP"  "RF_ECLAM"

head(newbornpr, n=4)

# A tibble: 4 x 17
      X MAGER MRACEREC   MAR MEDUC PRECARE WTGAIN CIG_REC  PGWT  DWGT
  <int> <int>    <int> <int> <int>   <int>  <int> <chr>   <int> <int>
1 23767    34        1     1     6       3     24 N         125   149
2  3882    29        1     1     6       3     29 N         120   149
3 18949    22        1     2     3       3     16 N         171   187
4  6553    27        1     1     5       3      6 N         216   222
# ... with 7 more variables: SEX <chr>, DBWT <int>, RF_DIAB <chr>,
#   RF_GEST <chr>, RF_PHYP <chr>, RF_GHYP <chr>, RF_ECLAM <chr>

Ajuste de un modelo lineal usando las funciones de R: lm() y glm()

El primer paso par ajustar un modelo es tener un data.frame.
En nuestro ejemplo usamos newbornpr.
Antes de ajustar el modelo es necesario especificar la estructura de la media usando las fórmulas en R.
Para esto es necesario entender la sintáxis de las fórmulas en R y el manejo de factores y sus respectivas codificaciones.

Ajuste de un modelo lineal usando las funciones de R: lm() y glm()

El primer paso par ajustar un modelo es tener un data.frame. En nuestro ejemplo usamos newbornpr Antes de ajustar el modelo es necesario especificar la estructura de la media usando las fórmulas en R. Para esto es necesario entender la sintáxis de las fórmulas en R y el manejo de factores y sus respectivas codificaciones.

# Sintaxis de formulas en R
y ~ x1              # Regresion simple
y ~ 0 + x1          # No intercepto
y ~ -1 + x1         # No intercepto
y ~ f1 + x1         # Factor + Continua
y ~ f1 + x1 + f1:x1 # Efectos principales e interaccion factor*continua
y ~ f1 + f2 + f1:f2 # ANOVA a dos vias
y ~ f1*f2           # ANOVA a dos vias
y ~ f1 + f2 + x1 + f1:x1 + f2:x1 + f1:f2 # interacciones dobles
y ~ (x1 + f1 + f2)^2  # interacciones dobles
y ~ f1 + f1:f3      # f3 anidado en f1
y ~ f1 + f1%in%f3   # f3 anidado en f1
y ~ f1/f3           # f3 anidado en f1

Ajuste de un modelo lineal usando las funciones de R: lm() y glm()

Además de especificar fórmulas con los factores y variables numéricas directamente también es posible usar transformaciones en las fórmulas.

# Transformaciones en las formulas
log(y) ~ sqrt(x1) # Transformacion logaritmica y raiz cuadrada
y ~ ordered(x1, breaks) # factor ordenado
y ~ poly(x1, 3) # polinomio de grado 3
y ~ poly(x1, x2, 2) # superficie cuadratica bivariada
y ~ bs(x1, df=3) # B-splines en x1 (nodos + orden)
y ~ I(x1 + 100/x2) # regresion simple sobre (x1+100/x2)

En cuanto a los factores supongamos que tenemos una variable cualitativa con tres niveles llamada f.
- Primero definimos la variable como un factor
- luego debemos especificar el constraste para ese factor.

# Factores y contrastes
f1 = factor(f) # Defina f como un factor
contrasts(f1) <- contr.treatment(3) # Constraste de grupo de referencia (nivel 1). Contraste por defecto en R

Contrastes de factores usando las funciones de R: lm() y glm().

Para entender como funciona contr.treatment veamos la codificación de los niveles.
- Note que el nivel 1 no aparece en la codificación de la columna ya que es el nivel de referencia.
- Esto se puede modificar usando la opción base. Por ejemplo: contr.treatment(3).
De manera similar, existen otros contrastes en R.
El mensaje más importante en esta parte es que el uso de un constraste en particular determina la interpretación de los coeficientes en el modelo.

Ejemplo

Para el nivel de referencia 1

#Nivel 1 = Referencia
contr.treatment(3) # toma el nivel 1 como referencia por defecto

Para el nivel de referencia 3

# Nivel 3 = Referencia
contr.treatment(3, base = 3)

Sintaxis en `R` para extraer información de un modelo ajustado con las funciones lm() y gls()

Objeto	Función `lm()`	Función `gls()`
Resumen	`sum.lm <- summary(lm.fit)`	`sum.gls <- summary(lm.gls)`
\( \hat{\beta} \)	`coef(lm.fit)`	`coef(gls.fit)`
\( \hat{\beta} \), \( \hat{se}(\hat{\beta}) \), \( t-test \)	`vcov(lm.fit)`	`vcov(gls.fit)`
IC 95\% para \( \beta \)	`confint(lm.fit)`	`confint(gls.fit), intervals(gls.fit, which="coef")`
ML	`logLik(lm.fit)`	`logLik(gls.fit, REML=FALSE)`
REML	`logLik(lm.fit, REML=TRUE)`	`logLik(gls.fit, REML=TRUE)`
AIC	`AIC(lm.fit)`	`AIC(gls.fit)`
BIC	`BIC(lm.fit)`	`BIC(gls.fit)`
Valores ajustados	`fitted(lm.fit)`	`fitted(gls.fit)`
Residuales crudos	`residuals(lm.fit, type="response")`	`residuals(lm.fit, type="response")`
Predicciones	`predicted(lm.fit, newdata)`	`predicted(gls.fit, newdata)`
Matriz de diseño	`model.matrix(lm.fit)`

Sintaxis en `R` parametrización de variables categóricas de un modelo ajustado con las funciones lm() y gls()

La parametrización de variables categóricas en R se puede llevar a acabo usando los siguientes comandos:

f = as.factor(f) # f es el factor
m = lm(y ~ f, contrasts = list(f = "contr.SAS")) # último nivel de referencia

El modelo homocedástico (varianza constante) se puede ajustar con la función lm() y gls(). Sin embargo, el modelo heterocedástico (varianza no constante) se puede ajustar únicamente con gls().
La tabla anterior muestra un resumen de las funciones que se usan para extraer los resultados de un objeto resultante de ajustar un modelo con lm() y gls().
Por ahora asumamos que ajustamos un modelo usando las dos funciones y que llamamos lm.fit y gls.fit a los objetos resultantes del ajuste de los modelos.

Ejemplo de los visto hasta ahora

Los siguientes son dos modelos que ilustran lo discutido hasta ahora.

# Ejemplo de un modelo ajustado con las funciones
fit.ml <- lm(y ~ x1 + x2, data = mydata)
fit.gls <- gls(y ~ x1 + x2, data = mydata)

Modelo con varianza homogénea

Vamos a analizar la relación entre algunas de las variables predictoras y el peso del recién nacido (DBWT).
El peso de la madre en el nacimiento es la variable que muestra la correlación más grande con el peso del recién nacido, aunque dicha correlación se puede considerar baja.
Otro punto importante en este análisis gráfico es que el peso de la madre antes del embarazo se correlaciona altamente con el peso al momento del parto (PGWT, DWGT).
Esto se conoce como multicolinealidad en las variables predictoras.
Las consecuencias de multicolinealidad en regresión pueden ser considerables, particularmente los errores se pueden sobre-estimar y las estimaciones de los coeficientes pueden ser muy diferentes de aquellos que se obtendrían en modelos univariados.
De otro lado, para determinar qué factores o variables cuantitativas podrían entrar al modelo se pueden construir diagramas de caja para la respuesta por el factor de interés.

Ejemplo: modelo con varianza homogénea (Datos Newbornpr)

A continuación encontramos las correlaciones entre las variables MAGER, PRECARE, WTGAIN, PGWT y DWGT.

library("GGally")
library("ggplot2")
p <- ggpairs(newbornpr[, c(2, 6, 7, 9, 10)])+
      theme(text=element_text(size = 20))

plot of chunk unnamed-chunk-10

Ejemplo: modelo con varianza homogénea (Datos Newbornpr)

La función give.n calcula promedios y tamaños de muestra. Figura: peso del recién nacido por sexo.

give.n <- function(x) { # fun. cal. el prom. y tam. de mues.
return(c(y = mean(x), label = length(x)))}
ggplot(newbornpr, aes(SEX, DBWT, fill = SEX)) +
  geom_boxplot() +  xlab("SEXO") +
  stat_summary(fun.data = give.n ,geom = "text", colour = "blue")+
  theme(text=element_text(size = 26))

plot of chunk unnamed-chunk-11

A continuación en la figura vemos el peso del recién nacido por raza de la madre.

newbornpr$f.race <- factor(newbornpr$MRACEREC, labels = c("Blanca", "Negra"))
ggplot(newbornpr, aes(f.race, DBWT, fill = f.race)) +
  geom_boxplot() +  xlab("Raza de la madre") +
  stat_summary(fun.data = give.n, geom = "text", colour = "blue")+
  theme(text=element_text(size = 26))

plot of chunk unnamed-chunk-12

Ejemplo: modelo con varianza homogénea (Datos Newbornpr)

Recuerde que la validez de la inferencia en un modelo lineal depende de la verificación de supuestos, particularmente normalidad y varianza constante.
Antes de proceder a interpretar los valores \( p \) (\( p-valor \)) debemos dar un vistazo a los gráficos de residuales.
A pesar de que se pueden llevar a cabo pruebas formales en este caso vamos a asumir que los supuestos son razonables para proceder con el análisis.

Ejemplo: modelo con varianza homogénea (Datos Newbornpr)

A continuación en la figura, vemos el peso del recién nacido por nivel educativo de la madre.

newbornpr$f.educ <- factor(newbornpr$MEDUC, labels = c("<=8", "9-12",
"High School", "Some college", "Assoc", "Bach", "Master", "Doc", "Unkn"))

p1 <- ggplot(newbornpr, aes(f.educ, DBWT, fill = factor(f.educ))) +
        geom_boxplot() +
        stat_summary(fun.data = give.n, geom = "text",colour = "blue") +
        xlab("Nivel Educativo") +
        theme(axis.text.x = element_text(angle = 90))+
        theme(text=element_text(size = 26))

p1

plot of chunk unnamed-chunk-14

Modelo con varianza homogénea (Datos Newbornpr). Modelado con lm()

# Modelo lineal usando lm() - Nacimientos en PR

# Modelo lineal incluyendo las variables MAGER, PRECARE, WTGAIN, PGWT, DWGT
lm.fit1 <- lm(DBWT ~ MAGER + PRECARE + WTGAIN + PGWT + DWGT, data = newbornpr)
summary(lm.fit1)


Call:
lm(formula = DBWT ~ MAGER + PRECARE + WTGAIN + PGWT + DWGT, data = newbornpr)

Residuals:
     Min       1Q   Median       3Q      Max 
-2706.84  -272.07    33.26   320.37  2441.90 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2128.7056   112.0223  19.003  < 2e-16 ***
MAGER          6.9438     2.9738   2.335  0.01975 *  
PRECARE       17.1212    12.6458   1.354  0.17608    
WTGAIN         0.3164     3.0170   0.105  0.91649    
PGWT          -6.1503     3.0240  -2.034  0.04224 *  
DWGT           9.1383     3.0135   3.032  0.00249 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 529.7 on 976 degrees of freedom
Multiple R-squared:  0.0921,    Adjusted R-squared:  0.08745 
F-statistic:  19.8 on 5 and 976 DF,  p-value: < 2.2e-16

# Modelo lineal usando lm() -- excluyendo PGWT
lm.fit2 <- lm(DBWT ~ MAGER + PRECARE + WTGAIN + DWGT, data = newbornpr)
summary(lm.fit2)


Call:
lm(formula = DBWT ~ MAGER + PRECARE + WTGAIN + DWGT, data = newbornpr)

Residuals:
     Min       1Q   Median       3Q      Max 
-2713.30  -266.73    32.62   321.35  2456.30 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2131.2748   112.1948  18.996  < 2e-16 ***
MAGER          6.6718     2.9756   2.242   0.0252 *  
PRECARE       17.3177    12.6658   1.367   0.1718    
WTGAIN         5.8961     1.2573   4.690 3.13e-06 ***
DWGT           3.0865     0.4772   6.468 1.56e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 530.6 on 977 degrees of freedom
Multiple R-squared:  0.08825,   Adjusted R-squared:  0.08452 
F-statistic: 23.64 on 4 and 977 DF,  p-value: < 2.2e-16

Ejemplo: modelo con varianza homogénea (Datos Newbornpr). Matriz de covarianza

Matriz de covarianza

# matriz de covarianza del modelo  lm.fit2
vcov(lm.fit2)

            (Intercept)        MAGER       PRECARE      WTGAIN
(Intercept)  12587.6754 -195.8187000 -622.05524288 -30.0975047
MAGER         -195.8187    8.8539512    5.75349030   0.3810544
PRECARE       -622.0552    5.7534903  160.42140798   0.6066499
WTGAIN         -30.0975    0.3810544    0.60664987   1.5807573
DWGT           -27.4225   -0.3139394   -0.06674166  -0.1390530
                    DWGT
(Intercept) -27.42250350
MAGER        -0.31393937
PRECARE      -0.06674166
WTGAIN       -0.13905301
DWGT          0.22769494

Intervalos de confianza

 # Intervalos del 95% de confianza del modelo  lm.fit2
confint(lm.fit2)

                   2.5 %      97.5 %
(Intercept) 1911.1042993 2351.445382
MAGER          0.8325422   12.510987
PRECARE       -7.5374860   42.172953
WTGAIN         3.4288019    8.363369
DWGT           2.1500933    4.022901

Modelo con varianza homogénea (Datos Newbornpr). Bondad de ajuste

A continuación veremos los diferentes criterios de bondad de ajuste del modelo lm.fit2.

logLik(lm.fit2) # loglikelihood (log maxima verosimilitud)

'log Lik.' -7551.884 (df=6)

logLik(lm.fit2, REML = T) # loglikelihood (log maxima verosimilitud) REML

'log Lik.' -7541.415 (df=6)

AIC(lm.fit2) # Akaike's Information criterion

[1] 15115.77

BIC(lm.fit2) # Bayesian Information Criterion

[1] 15145.11

head(model.matrix(lm.fit2)) # Construccion de la matriz de diseno

  (Intercept) MAGER PRECARE WTGAIN DWGT
1           1    34       3     24  149
2           1    29       3     29  149
3           1    22       3     16  187
4           1    27       3      6  222
5           1    34       1     40  193
6           1    22       3     33  148

Ejemplo: modelo con varianza homogénea (Datos Newbornpr). Residuales

Gráfico de residuales del modelo lm.fit2.

par(mfrow = c(2, 2), oma = c(0, 0, 4, 0))
plot(lm.fit2)

plot of chunk unnamed-chunk-20

Ejemplo: modelo con varianza homogénea (Datos Newbornpr). Residuales

Note que podemos llevar a cabo la prueba de hipótesis entre el modelo 1 y 2 usando la función anova() para determinar si los modelos son significativamente diferentes.
Existen otras formas de llevar a cabo pruebas de hipótesis para varios coeficientes simúltaneamente usando la función I().
Por ejemplo, si queremos probar que \( H_0: \beta_1 = \beta_2 \) entonces debe ajustar un modelo con el comando lm(y~I(x1+x2)) y compararlo con el modelo lm(y~x1+x2).
Otros tipos de pruebas de hipótesis se pueden llevar a cabo usando la función anterior.

# Prueba F entre los modelos lm.fit2 y lm.fit1
anova(lm.fit2, lm.fit1)

Analysis of Variance Table

Model 1: DBWT ~ MAGER + PRECARE + WTGAIN + DWGT
Model 2: DBWT ~ MAGER + PRECARE + WTGAIN + PGWT + DWGT
  Res.Df       RSS Df Sum of Sq      F  Pr(>F)  
1    977 275014126                              
2    976 273853491  1   1160635 4.1364 0.04224 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ajustemos el modelo anterior usando la función gls().

Modelo lineal usando gls() - Natos en PR Modelo lineal excluyendo PGWT

#Modelo lineal usando gls() - Natos en PR Modelo lineal excluyendo PGWT
library(nlme)
gls.fit <- gls(DBWT ~ MAGER + PRECARE + WTGAIN + DWGT, data = newbornpr)
summary(gls.fit)

Generalized least squares fit by REML
  Model: DBWT ~ MAGER + PRECARE + WTGAIN + DWGT 
  Data: newbornpr 
       AIC      BIC    logLik
  15094.83 15124.14 -7541.415

Coefficients:
                Value Std.Error   t-value p-value
(Intercept) 2131.2748 112.19481 18.996198  0.0000
MAGER          6.6718   2.97556  2.242189  0.0252
PRECARE       17.3177  12.66576  1.367288  0.1718
WTGAIN         5.8961   1.25728  4.689550  0.0000
DWGT           3.0865   0.47717  6.468285  0.0000

 Correlation: 
        (Intr) MAGER  PRECAR WTGAIN
MAGER   -0.587                     
PRECARE -0.438  0.153              
WTGAIN  -0.213  0.102  0.038       
DWGT    -0.512 -0.221 -0.011 -0.232

Standardized residuals:
        Min          Q1         Med          Q3         Max 
-5.11408102 -0.50274128  0.06147393  0.60569308  4.62968757 

Residual standard error: 530.5548 
Degrees of freedom: 982 total; 977 residual

Ajustemos el modelo anterior usando la función gls().

Matriz de covarianza

# matriz de covarianza del modelo  lm.fit2
vcov(gls.fit)

            (Intercept)        MAGER       PRECARE      WTGAIN
(Intercept)  12587.6754 -195.8187000 -622.05524288 -30.0975047
MAGER         -195.8187    8.8539512    5.75349030   0.3810544
PRECARE       -622.0552    5.7534903  160.42140798   0.6066499
WTGAIN         -30.0975    0.3810544    0.60664987   1.5807573
DWGT           -27.4225   -0.3139394   -0.06674166  -0.1390530
                    DWGT
(Intercept) -27.42250350
MAGER        -0.31393937
PRECARE      -0.06674166
WTGAIN       -0.13905301
DWGT          0.22769494

Intervalos de confianza

 # Intervalos del 95% de confianza del modelo  lm.fit2
confint(gls.fit)

                  2.5 %      97.5 %
(Intercept) 1911.377054 2351.172628
MAGER          0.839776   12.503753
PRECARE       -7.506695   42.142162
WTGAIN         3.431858    8.360312
DWGT           2.151253    4.021741

Ajustemos el modelo anterior usando la función gls().

Los criterios de bondad de ajuste del modelo gls.fit

logLik(gls.fit) # loglikelihood (log maxima verosimilitud)

'log Lik.' -7541.415 (df=6)

AIC(gls.fit) # Akaike's Information criterion

[1] 15094.83

BIC(gls.fit) # Bayesian Information Criterion

[1] 15124.14

Ejemplo: modelo con varianza homogénea (Datos Newbornpr). Residuales

Gráfico de residuales del modelo gls.fit.

par(mfrow = c(2, 2), oma = c(0, 0, 1, 0))
plot(gls.fit)

plot of chunk unnamed-chunk-26

Comentarios sobre pruebas de hipótesis sobre los modelos anidados

Podemos llevar a cabo una prueba de razón de verosimilitud entre dos modelos anidados usando la función anova() y gls().
Es muy importante que se ajusten los modelos usando el método ML ya que las medias de los modelos son diferentes.
Si los modelos anidados tienen medias iguales pero son anidados por los parámetros de varianza (por ejemplo, varianza constante vs. varianza no constante) entonces se puede usar el método REML y la prueba de razón de verosimilitud, de lo contrario no es adecuado usar la prueba LRT.
En este caso como \( p < 0.05 \) entonces rechazamos la hipótesis nula (modelo reducido) y nos quedamos con el modelo que incluye el término adicional.
Note que para propósitos de interpretación de los coeficientes este modelo quizás no sea conveniente debido a la multicolinedad. Además, los errores estándar del modelo completo tienden a ser mayores a los del modelo reducido como consecuencia de la multicolinedad.
Así que,a pesar del resultado de la prueba LRT la mejor estrategía es adoptar el modelo sin la variable PGWT.

Prueba LRT para modelos ajustados con gls()

# LR Test Modelo con PGWT
gls.fit1 <- gls(DBWT ~ MAGER + PRECARE + WTGAIN + PGWT + DWGT,
                method = "ML", data = newbornpr)

# Modelo con PGWT
gls.fit2 <- gls(DBWT ~ MAGER + PRECARE + WTGAIN + DWGT, method = "ML",
                data = newbornpr)

anova(gls.fit1, gls.fit2)

         Model df      AIC      BIC    logLik   Test L.Ratio p-value
gls.fit1     1  7 15113.61 15147.84 -7549.807                       
gls.fit2     2  6 15115.77 15145.10 -7551.884 1 vs 2 4.15308  0.0416

Recuerde que!

REML vs ML

ML: estimación de máxima verosimilitud
- Desarrolla una ecuación de probabilidad para la probabilidad de datos dados el modelo.
- Elige los valores de los parámetros que maximizan este probabilidad.
- La estimación de ML está ligeramente sesgada debido a la articulación estimación de múltiples parámetros.
REML: máxima verosimilitud restringida
- REML proporciona una solución para eliminar el sesgo
Al comparar dos modelos, ambos necesitan estar en forma usando el mismo procedimiento de estimación.

Ejercicios: comparación de métodos de estimación

Ejercicio 2: comparación de métodos de estimación

Ajuste el modelo lineal anterior (gls.fit2) usando la función gls() y el método ML. Compare las estimaciones de este modelo (ML) con las obtenidas con la función lm() y gls() usando REML.

Ejercicio 3: modelo lineal con factores

Construya cada modelo pedido, y luego responda cada inciso.

Ajuste el modelo lineal anterior (gls.fit2) pero agregue el sexo como una de las variables predictoras. Determine si esta nueva variable es significativa.
Ajuste el modelo lineal anterior (gls.fit2) pero agregue el sexo, la raza y su interacción como variables predictoras. Determine si la interacción es significativa.
Ajuste el modelo del punto 1 pero cambie el nivel de referencia de mujeres a hombres. Use el término relevel(SEX,ref=2) en lugar de SEX.
Finalmente, para el modelo en gls.fit2 agregue el nivel educativo y determine si esta variable es significativa usando la prueba LRT.

Modelo con varianza heterogénea

Recuerde que en un modelo con varianza heterogénea asumimos que

\[ \mu_i = E(y_i|\boldsymbol{x}_i ) = \boldsymbol{x}'_i \boldsymbol{\beta} \] y

\[ Var(\varepsilon_i |\boldsymbol{x}_i) = \sigma_i^2 \]

donde \( \sigma_i^2 \) podría tomar alguna forma funcional.

Para ajustar un modelo con errores heterocedásticos o varianza no constante (heterogénea) la función lm() ya no se puede utilizar.
En estos casos necesitamos recurrir a la función gls().

Primero vamos a definir lo que se conoce como la función de varianza, la cual denotaremos con \( \lambda(\boldsymbol{\sigma}, \mu_i; \boldsymbol{v}_i) \) siguiendo la notación del libro de referencia.
En la definición anterior \( \boldsymbol{\sigma} \) es un vector de parámetros común a todas las observaciones, \( \mu_i \) es la media para la observación \( i \) y \( \boldsymbol{v_i} \) es un vector de covariables relacionadas a la varianza.
Es decir, la función de varianza puede depender tanto de la media como de covariables adicionales y puede tener parámetros adicionales. Usando esta función de varianza podemos definir

\[ Var(\varepsilon_i) = \sigma^2\lambda^2(\boldsymbol{\sigma}, \mu_i; \boldsymbol{v}_i) \]

Modelo con varianza heterogénea

Note que la función \( \lambda(\cdot) \) es llamada la función de varianza aunque en realidad es la desviación estándar.
La tabla funciones de varianza en la función gls() de R anterior, muestra las diferentes funciones de varianza disponibles en la función gls() de R.
En esta tabla la covariable \( v_i \) se puede reemplazar por la media \( \mu_i \) en el caso que la varianza dependa de la media.
Cuando la varianza depende de pesos conocidos, como en muestreo, se usa el algoritmo WLS para estimación.
Si la varianza depende de la media entonces se usa el algoritmo IRLS.
En caso que la media dependa de una covariable y parámetros desconocidos \( \boldsymbol{\delta} \) entonces se usan ML o REML como métodos de estimación.

Función	\( \lambda_i \)	Descripción
varFixed()	\( w_i \)	Pesos fijos
varPower()	\( \|v_i\|^{\delta_{s_i}} \)	Potenciade una covariable \( v_i \)
varExp()	\( \exp(v_i \delta_{s_i}) \)	Exponencial de una covariable \( v_i \)
varConstPower()	\( \delta_{1,s_i}+ \|v_i\|^{\delta_{2,s_i}} \)	Constante más función potencia, \( \delta_{1,s_i}>0 \)
varIdent()	\( \delta_{s_i} \)	Varianzas diferentes por grupos/estratos: \( \delta_1=1 \); \( \delta_s>0 \), \( s \neq 1 \)

Modelo con varianza heterogénea

Un aspecto importante cuando se ajustan modelos con varianza heterogénea en R, es tener claro la parametrización de \( \lambda(\cdot) \).
Por ejemplo, supongamos que la varianza depende de un factor con dos niveles (e.g., sexo).
En este caso se usa la función varIdent(); R usa la siguiente parametrización: \( Var(\varepsilon_i) = \sigma_{s_i}^2\delta_{s_i}^2 \) pero fija \( \delta_{s_i}^2 = 1 \) para el primer nivel del factor.
Vamos a ilustrar más adelante con una aplicación.

Ejemplo: datos Brain weight

Para esta parte vamos utilizar los datos del tamaño de la cabeza y el cerebro de 237 individuos. Esta es la descripción de las variables: Gender (1=Male, 2=Female), Age Range (1=20-46, 2=46+), Head size (cm\( ^3 \)) y Brain weight (grams).

library(readr)
brainhead <- read_csv("~/Desktop/libro_glm/datos/brainhead.csv",
                      quote = "\"", col_names =FALSE)

# dando nombres a las variables
colnames(brainhead) = c("gender", "age", "head", "brain")
head(brainhead)

# A tibble: 6 x 4
  gender   age  head brain
   <int> <int> <int> <int>
1      1     1  4512  1530
2      1     1  3738  1297
3      1     1  4261  1335
4      1     1  3777  1282
5      1     1  4177  1590
6      1     1  3585  1300

str(brainhead)

Classes 'tbl_df', 'tbl' and 'data.frame':   237 obs. of  4 variables:
 $ gender: int  1 1 1 1 1 1 1 1 1 1 ...
 $ age   : int  1 1 1 1 1 1 1 1 1 1 ...
 $ head  : int  4512 3738 4261 3777 4177 3585 3785 3559 3613 3982 ...
 $ brain : int  1530 1297 1335 1282 1590 1300 1400 1255 1355 1375 ...
 - attr(*, "spec")=List of 2
  ..$ cols   :List of 4
  .. ..$ X1: list()
  .. .. ..- attr(*, "class")= chr  "collector_integer" "collector"
  .. ..$ X2: list()
  .. .. ..- attr(*, "class")= chr  "collector_integer" "collector"
  .. ..$ X3: list()
  .. .. ..- attr(*, "class")= chr  "collector_integer" "collector"
  .. ..$ X4: list()
  .. .. ..- attr(*, "class")= chr  "collector_integer" "collector"
  ..$ default: list()
  .. ..- attr(*, "class")= chr  "collector_guess" "collector"
  ..- attr(*, "class")= chr "col_spec"

Ejemplo: datos Brain weight (descripción de los datos)

Veamos algunas gráficas que muestre la relación entre las variables.

La figura abajo muestra que el tamaño del cerebro tiende a aumentar a medida que aumenta el tamaño de la cabeza y dicha relación se puede modelar linealmente.
Además, es claro que hay más variabilidad en el tamaño del cerebro a medida que crece el tamaño de la cabeza.

ggplot(data=brainhead, aes(x = head, y = brain))+
   geom_point()+
   theme(text=element_text(size = 26))

plot of chunk unnamed-chunk-32

Ejemplo: datos Brain weight (modelación)

Vamos a comenzar ajustando un modelo lineal con varianza homogénea.
Con base en los gráficos de residuales podemos ver que hay indicios de que los residuales tienen varianza no constante que podría depender de la media o el tamaño de la cabeza.
Cuando queremos detectar varianza no constante, es importante trabajar con los residuales estudentizados(ya que estos eliminan el efecto de los puntos de leverage), es decir \[ r_i=\frac{\hat{\varepsilon}_i}{\hat{\sigma}_{R}{\sqrt{1-h_{ii}}}}, \]
Es decir, al graficar los residuales crudos se podría concluir erróneamente, que tienen varianza no constante ya que la varianza de \( r_i \) depende de los valores de apalancamiento (leverage points \( h_{ii} \)).

# Modelo lineal con varianza homogénea
lm.brain = lm(brain ~ head, data = brainhead)
summary(lm.brain)


Call:
lm(formula = brain ~ head, data = brainhead)

Residuals:
     Min       1Q   Median       3Q      Max 
-171.935  -48.636   -2.397   46.832  242.578 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 361.88685   48.05983    7.53 1.08e-12 ***
head          0.25382    0.01316   19.29  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 73.84 on 235 degrees of freedom
Multiple R-squared:  0.6129,    Adjusted R-squared:  0.6112 
F-statistic: 372.1 on 1 and 235 DF,  p-value: < 2.2e-16

Ejemplo: datos Brain weight (modelación)

# Figura
library("ggplot2")
library("gridExtra")
brain.res=data.frame(head=brainhead$head, stud.resid=rstudent(lm.brain),
            gender=factor(brainhead$gender), age=factor(brainhead$age))

p1 <- ggplot(data = brain.res) +
        stat_qq(aes(sample = stud.resid)) +
        geom_abline(color="blue")+
        labs(title = "Plot de Normalidad",
             x = "Teóricos", y = "Resid. Student")+
        theme(text=element_text(size = 18))

p2 <- ggplot(brain.res, aes(x=gender, y=stud.resid, fill=gender))+
        geom_boxplot() +
        scale_fill_discrete(name="Gender",
                         breaks=c("1", "2"),
                         labels=c("Male", "Female"))+
        labs(title = "Boxplot de Resid. Stud vs. gender",
             x = "Gender", y = "Resid. Student")+
        theme(text=element_text(size = 18))

p3 <- ggplot(brain.res, aes(x=age, y=stud.resid, fill=age)) +
        geom_boxplot() +
        scale_fill_discrete(name="Age Group",
                         breaks=c("1", "2"),
                         labels=c("20 < age < 46", "age > 46"))+
        labs(title = "Boxplot de Resid. Stud vs. age",
             x = "Age Group", y = "Resid. Student")+
        theme(text=element_text(size = 18))

p4 <- ggplot(brain.res, aes(x=head, y=stud.resid)) +
        geom_point() +
        geom_hline(yintercept = 0, color="blue")+
        labs(title = "Plot de Resid. Stud vs. head",
             x = "Head", y = "Resid. Student")+
        theme(text=element_text(size = 18))

grid.arrange(p1, p2, p3, p4, ncol=2)

plot of chunk unnamed-chunk-35

Ejemplo: Ajuste de varianza heterogenea

Para propósitos de ilustración vamos a ajustar un modelo con varianza heterogénea en los residuales que dependen de la edad.

Para este modelo M1 tenemos:
- \( Var(\varepsilon_i|edad = 1) = \sigma^2 \delta_1^2 = (73.19265^2)(1) = 5357.164 \) y
- \( Var(\varepsilon_i|edad = 2) = \sigma^2 \delta_2^2 = (73.19265^2)(0.9930937) = 5320.166 \)
Ahora vamos a determinar si adoptamos el modelo más sencillo con varianza constante y el modelo cuya varianza depende de la edad usando la prueba LRT. Note que en este caso el \( p-valor \) es 0.9409, así que la evidencia sugiere que el modelo reducido (varianza constante) es preferido.

Modelo con varianza homogénea (M0)

library(nlme)
# Modelo lineal con varianza homogénea (M0)
gls.brain = gls(brain ~ as.factor(gender) + as.factor(age) + head,
                data = brainhead)
summary(gls.brain)

Generalized least squares fit by REML
  Model: brain ~ as.factor(gender) + as.factor(age) + head 
  Data: brainhead 
       AIC      BIC    logLik
  2700.609 2717.865 -1345.305

Coefficients:
                      Value Std.Error   t-value p-value
(Intercept)        453.0072  60.03404  7.545839  0.0000
as.factor(gender)2 -22.7488  11.29986 -2.013194  0.0452
as.factor(age)2    -22.1110   9.68811 -2.282287  0.0234
head                 0.2347   0.01539 15.248955  0.0000

 Correlation: 
                   (Intr) as.fctr(gn)2 as.factr(g)2
as.factor(gender)2 -0.589                          
as.factor(age)2    -0.265  0.167                   
head               -0.990  0.528        0.177      

Standardized residuals:
        Min          Q1         Med          Q3         Max 
-2.20443534 -0.67950464 -0.02642223  0.59539209  3.39128977 

Residual standard error: 72.92205 
Degrees of freedom: 237 total; 233 residual

Modelo con varianza heterogénea por edad como factor (M1)

# Modelo lineal con varianza heterogénea por edad como factor (M1)
gls.brain2 = gls(brain ~ factor(gender) + factor(age) + head,
                 weights = varIdent(form = ~1|factor(age)),
                 data = brainhead)
summary(gls.brain2)

Generalized least squares fit by REML
  Model: brain ~ factor(gender) + factor(age) + head 
  Data: brainhead 
       AIC     BIC    logLik
  2702.604 2723.31 -1345.302

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | factor(age) 
 Parameter estimates:
        1         2 
1.0000000 0.9930937 

Coefficients:
                   Value Std.Error   t-value p-value
(Intercept)     452.4986  60.03207  7.537615  0.0000
factor(gender)2 -22.6329  11.29977 -2.002955  0.0463
factor(age)2    -22.0913   9.69273 -2.279166  0.0236
head              0.2349   0.01539 15.258179  0.0000

 Correlation: 
                (Intr) fctr(gn)2 factr(g)2
factor(gender)2 -0.589                    
factor(age)2    -0.265  0.167             
head            -0.990  0.528     0.177   

Standardized residuals:
        Min          Q1         Med          Q3         Max 
-2.21195105 -0.67757343 -0.02631317  0.59305409  3.40240572 

Residual standard error: 73.19265 
Degrees of freedom: 237 total; 233 residual

Compararión de los modelos M0 y M1 usando la prueba LRT (anova)

# prueba LRT
anova(gls.brain, gls.brain2)

           Model df      AIC      BIC    logLik   Test     L.Ratio p-value
gls.brain      1  5 2700.609 2717.865 -1345.305                           
gls.brain2     2  6 2702.604 2723.310 -1345.302 1 vs 2 0.005501798  0.9409

El gráfico de residuales sugiere que la varianza de los residuales quizás dependa del tamaño de la cabeza.
Así que procedemos a ajustar un modelo con la función de varianza potencia para determinar si la varianza de los errores aumenta con el tamaño de la cabeza, es decir,

\[ Var(\varepsilon) = \sigma^2 |head|^\delta = 1.969301^2|head|^{0.440503} \]

Note que en este caso la prueba de LRT sugiere que el modelo con varianza heterogénea para los errores es preferido (\( p-valor=0.2681 \)).

Modelo con varianza igual a la potencia de una covariable=head (M3)

# Modelo lineal con varianza igual a la potencia de una covariable=head (M3)
gls.brain3 = gls(brain ~ factor(gender) + factor(age) + head,
                 weights = varPower(form = ~head), data = brainhead)
s <- summary(gls.brain3)

Generalized least squares fit by REML
  Model: brain ~ factor(gender) + factor(age) + head 
  Data: brainhead 
       AIC      BIC    logLik
  2701.383 2722.089 -1344.691

Variance function:
 Structure: Power of variance covariate
 Formula: ~head 
 Parameter estimates:
   power 
0.440503 

Coefficients:
                   Value Std.Error   t-value p-value
(Intercept)     458.2897  59.99305  7.639046  0.0000
factor(gender)2 -23.4738  11.22086 -2.091978  0.0375
factor(age)2    -22.5717   9.65548 -2.337708  0.0202
head              0.2334   0.01547 15.087823  0.0000

 Correlation: 
                (Intr) fctr(gn)2 factr(g)2
factor(gender)2 -0.594                    
factor(age)2    -0.268  0.169             
head            -0.990  0.532     0.180   

Standardized residuals:
        Min          Q1         Med          Q3         Max 
-2.06109501 -0.67689163 -0.03424125  0.60693993  3.72251851 

Residual standard error: 1.969301 
Degrees of freedom: 237 total; 233 residual

Comparemos los modelos M0 y M3, usando la prueba LRT (anova)

# prueba LRT
anova(gls.brain, gls.brain3)

           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
gls.brain      1  5 2700.609 2717.865 -1345.305                        
gls.brain3     2  6 2701.383 2722.089 -1344.691 1 vs 2 1.226684  0.2681

Intervalo de confianza para los parámetros en la varianza

# Intervalo de confianza para los parámetros en la varianza
intervals(gls.brain3, which = "var-cov")$varStruct

           lower     est.    upper
power -0.3154461 0.440503 1.196452
attr(,"label")
[1] "Variance function:"

Gráficos de varianza para los modelos lm.brain, y M3

Ahora veamos los gráficos de varianza para los modelos lm.brain, y M3.

library("ggplot2")
library("gridExtra")
lm.brain.res=data.frame(head=lm.brain$fitted,
                        std.resid=rstandard(lm.brain))
p1 <- ggplot(lm.brain.res, aes(x=head, y=std.resid)) +
        geom_point() +
        geom_hline(yintercept = 0, color="blue")+
        labs(title = "Varianza homogénea",
             x = "Fitted", y = "Resid. standarized")+
        theme(text=element_text(size = 26))

gls.brain3.res=data.frame(fitted=gls.brain3$fitted,
                          std.resid=resid(gls.brain3, type = "pearson"))
p2 <- ggplot(gls.brain3.res, aes(x=fitted, y=std.resid)) +
        geom_point() +
        geom_hline(yintercept = 0, color="blue")+
        labs(title = "Varianza heterogénea",
             x = "Fitted", y = "Resid. standarized")+
        theme(text=element_text(size = 26))

grid.arrange(p1, p2, ncol=2)

plot of chunk unnamed-chunk-44

Gráfica del cuadrado de los residuales del modelo versus el tamaño de la cabeza y su fun. de varianza

Vamos graficar el cuadrado de los residuales del modelo versus el tamaño de la cabeza y la respectiva función de varianza que obtuvimos en el modelo anterior figura del modelo gls.brain3.
Note que esta gráfica intenta usar el hecho que \( Var(\varepsilon) = E(\varepsilon^2) \); así, un gráfico de los residuales al cuadrado versus la covariable podría indicar que función de varianza ajustar.

Gráfica del cuadrado de los residuales del modelo versus el tamaño de la cabeza y su fun. de varianza

head.sort = sort(brainhead$head)
pot = gls.brain3$modelStruct$varStruct
d.brain <- data.frame(x =brainhead$head, y=resid(gls.brain3)^2)
d.brain.l <- data.frame(hs=head.sort, 
      gls.b.h.s.p = (gls.brain3$sigma^2)*(head.sort^pot)^2)

p1 <- ggplot(data=d.brain, aes(x=x, y=y)) +
        geom_point() +
        geom_line(data=d.brain.l, 
                  aes(x=hs, y=gls.b.h.s.p), 
                  color="blue")+
        labs(title = "Función de varianza ajustada",
             x = "Head", y = expression(paste(Res.^2)))+
        theme(text=element_text(size = 26))

p1

plot of chunk unnamed-chunk-46

Comentarios finales

Una forma de trabajar el problema de varianza no constante consiste en transformar la variable respuesta usando por ejemplo:
- logarítmica
- exponencial
- transformaciones de Box-Cox
- o cualquier otra función que tenga como resultado un modelo con varianza homogénea para los errores.
Cuando existen muchas variables predictoras es muy complicado modelar la función de varianza directamente.
- puede de difícil establecer la fuente de varianza no constante.
Tranformar los datos es una estrategia comúnmente usada y no tiene muchos inconvenientes si el modelo se ajusta principalmente para predecir.
Si el modelo se ajusta para interpretar los coeficientes de las variables explicativas en términos de cambios promedio en la respuesta entonces la estrategia de transformación se debe manejar cuidadosamente para no hacer intepretaciones incorrectas en la escala original.
- Todos los coeficientes del modelo ajustado con la variable transformada se deben interpretar en la escala transformada. La interpretación de coeficientes en la escala original requiere de un esfuerzo adicional de parte del analista.

Ejercicio: Salarios de profesores de la librería (car)

Descripción El salario académico de nueve meses de 2008-09 para Profesores Auxiliares (397 observaciones sobre las siguientes 6 variables.), Profesores Asociados y Profesores en una universidad de los EE.UU. Los datos fueron recopilados como parte del esfuerzo continuo de la administración de la universidad para monitorear las diferencias salariales entre los miembros de la facultad masculina y femenina.

rango: factor con niveles AssocProf AsstProf Prof
disciplina: factor con niveles A (departamentos “teóricos”) o B (departamentos “aplicados”).
yrs.since.phd: años desde PhD.
yrs.service: años de servicio.
sexo: factor con niveles Mujer Hombre
salario: salario de nueve meses, en dólares.

Ejercicio: Salarios de profesores

Usar las variables salario como respuesta y las variables yrs.since.phd yrs.since como variables predictoras continuas y sexo (o disciplina) como variable predictora categórica.

Hacer un estudio de los datos (boxplot, diagrama de puntos)
Ajuste el modelo salario.lm: en el cual verificará los residuales (en particular observe si hay homocedasticidad)
Si no hay homocedasticidad, entonces haga un estudio para encontrar la mejor función que estructure la varianza de los residuales.