Regresión múltiple

Los datos trees: presentación y análisis gráfico

Vamos a utilizar el conjunto de datos trees (ya incluido en R) que contiene medidas tomadas en 31 cerezos. Para cada uno de ellos se mide el volumen de la madera que se obtuvo de ellos, la altura y el diámetro (medido a cierta altura del suelo). Los nombres de las variables son:

  • Girth: Diámetro en pulgadas.
  • Height: Altura en pies.
  • Volume: Volumen de madera en pies cúbicos.

El objetivo es explicar el volumen de madera obtenido en función de las otras dos variables. Para tener una primera idea de cómo son las relaciones entre las variables, podemos representar una matriz de diagramas de dispersión:

# put histograms on the diagonal
panel.hist <- function(x, ...)
{
    usr <- par("usr")
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}

pairs(trees,pch=16,cex=0.5,oma=c(2,2,5,2),
      diag.panel = panel.hist,bg = "light blue",
      main = 'Matriz de diagramas de dispersión (trees)')

Ejercicio 1: ¿Cuáles de las hipótesis habituales del modelo de regresión lineal no parece que se vayan a cumplir?

Compara el gráfico anterior con el que se obtiene cuando consideramos los logaritmos de las tres variables en lugar de las variables originales:

pairs(log(trees), pch=16,cex=0.5,oma=c(2,2,5,2),
      diag.panel = panel.hist,bg = "light blue", 
      main = 'Matriz de diagramas de dispersión (log(trees))',
      labels=c("log(Girth)","log(Height)","log(Volume)"))

Ajuste del modelo

Si \(h\) es la altura, \(D\) el diámetro del árbol y \(V\) es el volumen, debemos esperar \(V\approx \pi h D^2/4\). La relación anterior sugiere ajustar el modelo de regresión múltiple: \[ \log V = \beta_0 + \beta_1\log h + \beta_2\log D + \epsilon \] Ajustamos el modelo y obtenemos un resumen de los resultados:

reg = lm(log(Volume) ~ log(Height) + log(Girth), data = trees)
summary(reg)
## 
## Call:
## lm(formula = log(Volume) ~ log(Height) + log(Girth), data = trees)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.168561 -0.048488  0.002431  0.063637  0.129223 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.63162    0.79979  -8.292 5.06e-09 ***
## log(Height)  1.11712    0.20444   5.464 7.81e-06 ***
## log(Girth)   1.98265    0.07501  26.432  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08139 on 28 degrees of freedom
## Multiple R-squared:  0.9777, Adjusted R-squared:  0.9761 
## F-statistic: 613.2 on 2 and 28 DF,  p-value: < 2.2e-16

Ejercicio 2: ¿Qué valores es de esperar que tomen \(\beta_1\) y \(\beta_2\)? ¿Se parecen los estimadores obtenidos a estos valores? ¿Son los regresores individualmente significativos?

La matriz de covarianzas estimadas del vector \(\hat{\boldsymbol\beta}=(\hat\beta_0,\hat\beta_1,\hat\beta_2)'\) (es decir, \(s^2_R(\mathbb X'\mathbb X)^{-1}\), donde \(\mathbb X\) es la matriz de diseño y \(s^2_R\) es la varianza residual) se obtiene con:

vcov(reg)
##             (Intercept)  log(Height)   log(Girth)
## (Intercept)  0.63966361 -0.160062124  0.020793539
## log(Height) -0.16006212  0.041794512 -0.008130512
## log(Girth)   0.02079354 -0.008130512  0.005626592

Con el comando confint obtenemos intervalos de confianza para los parámetros \(\beta_0\), \(\beta_1\) y \(\beta_2\), al nivel de confianza del 95%:

confint(reg)
##                 2.5 %    97.5 %
## (Intercept) -8.269912 -4.993322
## log(Height)  0.698353  1.535894
## log(Girth)   1.828998  2.136302

Ejercicio 3: Calcula intervalos de confianza de nivel 0.9 para \(\beta_1\) y \(\beta_2\).

Para calcular la tabla ANOVA del contraste global de la regresión, primero ajustamos el modelo sin regresores (sólo con término independiente) y luego comparamos la variabilidad incremental del modelo reducido y el completo:

reg0 = lm(log(Volume) ~ 1, data = trees)
anova(reg0,reg)
## Analysis of Variance Table
## 
## Model 1: log(Volume) ~ 1
## Model 2: log(Volume) ~ log(Height) + log(Girth)
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     30 8.3087                                  
## 2     28 0.1855  2    8.1232 613.19 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ejercicio 4: Contrasta la hipótesis \(H_0: \beta_1=\beta_2=0\). ¿Es el modelo globalmente significativo? ¿Es coherente con lo obtenido en los contrastes individuales?

Ejercicio 5: Contrasta \(H_0: \beta_2 = 2\beta_1\) mediante el método de incremento relativo de la variabilidad:

Z = log(trees$Height) +2 * log(trees$Girth)
regH0 = lm(log(trees$Volume) ~ Z)
regH1 = lm(log(trees$Volume) ~ log(trees$Height) + log(trees$Girth))
anova(regH0,regH1)

o, equivalentemente, usando la función I:

regR = lm(log(Volume) ~ I(log(Height) +2* log(Girth)), data = trees)
anova(regR,reg)

Diagnóstico del modelo

La diagnosis de un modelo estadístico (como el de regresión lineal que utilizamos en esta práctica) tiene como objetivo comprobar si se cumplen las hipótesis del modelo (en este caso, independencia, linealidad, homocedasticidad y normalidad). Si los datos no satisfacen las hipótesis del modelo, entonces las conclusiones de inferencia estadística podrían ser incorrectas. Las hipótesis de linealidad, homocedasticidad y normalidad equivalen a que \(\epsilon_i\sim N(0,\sigma^2)\), \(i=1,\ldots,n\). Como se considera que los residuos \(e_i\) son los valores observados de las perturbaciones \(\epsilon_i\), lo que se hace es comprobar si los residuos son compatibles con esa distribución. En la práctica se llevan a cabo comprobaciones gráficas sencillas.

Para comprobar si los residuos (estandarizados) tienen media aproximadamente 0 (hipótesis de linealidad) y si la variabilidad no cambia con los regresores (hipótesis de homocedasticidad), representamos los residuos \(e_i\) frente a los valores previstos \(\hat y_i\):

# residuos frente a valores ajustados
residuos.estandarizados <- rstandard(reg)
valores.ajustados <- fitted(reg)
plot(valores.ajustados, residuos.estandarizados)
abline(h=0, lty=2)  # añade recta y=0

Podemos comparar el anterior gráfico con el gráfico de residuos frente a valores previstos cuando no transformamos las variables con logaritmos (que ya sabemos que NO es el modelo correcto):

regNT <- lm(Volume ~ Height + Girth, data = trees)
plot(regNT$fitted.values,regNT$residuals)
abline(h=0, lty=2)
title("Obsérvese la falta de linealidad")

regNT2 <- lm(Volume ~ Height, data = trees)
plot(regNT2$fitted.values,regNT2$residuals)
abline(h=0, lty=2)
title("Obsérvese la heterocedasticidad")

Para comprobar si los residuos son aproximadamente normales (hipótesis de normalidad):

qqnorm(residuos.estandarizados)
qqline(residuos.estandarizados)

hist(residuos.estandarizados,freq=F)
v = seq(-3,3,0.01)
d = dnorm(v)
lines(v,d,lwd=2,col="red")

Hay que tener cuidado con los contrastes de bondad de ajuste a la distribución normal, pues habitualmente se supone independencia de las observaciones de la muestra, pero en este caso los residuos no son independientes (porque \(\sum_{i=1}^n e_i=0\)). Se puede encontrar más información en Peña (2010).

Sobre las transformaciones de variables

Es habitual transformar alguna o todas las variables del conjunto de datos para “linealizar” las relaciones entre ellas. Es frecuente que la falta de linealidad en la relación entre las variables vaya acompañada de heterocedasticidad y ambos fallos en las suposiciones del modelo de regresión lineal se resuelvan transformando adecuadamente las variables. Las transformaciones que suelen funcionar bien son las que “acercan” los datos transformados al modelo normal multivariante, puesto que en este modelo las únicas relaciones de dependencia entre las variables son las lineales. Así que una manera de elegir la transformación adecuada para una variable individual es que la variable transformada tenga un histograma simétrico (parecido al de una normal univariante).

Las transformaciones suelen ser funciones sencillas. Si la variable tiene asimetría hacia la derecha, la transformación debe “contraer” los valores altos y expandir los bajos. En este caso se usan funciones crecientes y cóncavas, como el logaritmo (neperiano) o las potencias con exponente menor que 1, como tomar la raíz cuadrada de la variable. Si la variable tiene asimetría hacia la izquierda (fenómeno mucho menos corriente que la asimetría hacia la derecha), entonces se aplicará una transformación dada por funciones crecientes y convexas como la exponencial o potencias con exponente mayor que 1, como elevar al cuadrado.

Respecto a qué variables transformar, conviene ir probando a transformar una a una cada variable (si se considera necesario) y evaluando el efecto (si mejora o no el ajuste del modelo a los datos, si se corrige la falta de linealidad o heterocedasticidad). Como regla general (ver más detalles en Peña 2010), si hay heterocedasticidad y falta de linealidad, se suele transformar por lo menos la variable respuesta. Si no hay heterocedasticidad, pero sí falta de linealidad, basta transformar regresores.

Ejemplo de falta de linealidad y heterocedasticidad:

Ejemplo de falta de linealidad sin heterocedasticidad:

Referencias

Peña, D. (2010). Regresión y diseño de experimentos. Alianza Editorial.