Análisis de los supuestos
Luis Fernando Apáez Álvarez
7/11/2021
Ejercicio
Algunos grupos madereros, silvicultores y biólogos están interesados en poder predecir el volumen de un árbol con sólo conocer su diámetro. Para ello se cuentan con 70 observaciones en la base de datos shortleaf donde se tienen dos variables: el diámetro en pulgadas (Diam) y el volumen en pies cúbicos (Vol). Verifique si puede ajustarse un modelo de regresión lineal simple e interprete los coeficientes.
Primero importamos los datos.
#Cargamos los datos y vemos rapidamente su estructura
datos<- read.table("shortleaf.txt", header = T)
head(datos)## Diam Vol
## 1 4.4 2.0
## 2 4.6 2.2
## 3 5.0 3.0
## 4 5.1 4.3
## 5 5.1 3.0
## 6 5.2 2.9Notamos que tenemos dos variables de tipo numérico y 70 observaciones
str(datos)## 'data.frame': 70 obs. of 2 variables:
## $ Diam: num 4.4 4.6 5 5.1 5.1 5.2 5.2 5.5 5.5 5.6 ...
## $ Vol : num 2 2.2 3 4.3 3 2.9 3.5 3.4 5 7.2 ...Notamos que la variable independiente es Diam y la dependiente es Vol.
Supesto de linealidad
Comenzamos graficando un scatterplot:
par(mar=c(4, 5, 1, 1))
plot(datos$Diam, datos$Vol, xlab="Diámetro (pulgadas)", ylab="Volumen (pies cúbicos)",col="blue")
o alternativamente
# Cargamos la librería necesaria
library(ggplot2)
#graficamos
ggplot(datos, aes(x=Diam, y=Vol)) +
geom_point(colour=4) 
Aparentemente podemos ver que no se cumple el supuesto de linealidad, pero ajustaremos la regresión lineal sobre estos datos para hacerlo más notorio
#Ajustamos el modelo
fit1 <- lm(Vol~Diam, data=datos)
# Mostramos los valores de b0 y b1
fit1$coefficients## (Intercept) Diam
## -41.568090 6.836718Después graficamos la línea de ajuste en el diagrama de dispersión:
par(mar=c(4, 5, 1, 1))
plot(datos$Diam, datos$Vol, xlab="Diámetro (pulgadas)", ylab="Volumen (pies cúbicos)",col="blue")
# Agregamos la línea
abline(fit1$coefficients, col="red")
Más aún podemos utilizar la gráfica que nos brinda R para analizar el supuesto de la linealidad
plot(fit1, 1)
lo cual nos corroborá que este supuesto no se cumple pues nos gustaría que la línea roja estuviera muy cerca de la recta en 0. También, podemos checar rápidamente el supuesto de la homocedasticidad
par(mar=c(4, 5, 1, 1))
plot(fit1, 3) 
el cual también no se cumple. Por lo anterior, debemos de encontrar un mejor modelo transformando nuestros datos.
Luego, utilizaremos la transformación logarítmica para hacer que los supuestos se cumplan, donde se sabe que la transformación que funciona es esa pues ya se probó (para los fines de esta práctica usamos directo dicha transformación, pero para saber qué transformación usar están las transformaciones de tipo Box-Cox o Box-Tidwell).
#Creamos nuestras variables transformadas y las guardamos en el dataframe.
datos$lnDiam=log(datos$Diam)
datos$lnVol=log(datos$Vol)
#Ajustamos un modelo pero con los datos transformados.
fit2 <- lm(lnVol ~ lnDiam, data = datos)
fit2$coefficients## (Intercept) lnDiam
## -2.871790 2.564416Después, comparamos las gráficas de los dos modelos obtenidos, es decir, la gráfica con los datos no transformados versus la gráfica con los datos transformados
par(mfrow=c(1,2))
par(mar=c(4, 5, 1, 1))
# Gráfica datos no transformados
plot(datos$Diam, datos$Vol, main="Datos en escala original",xlab="", ylab="",col="blue")
abline(fit1$coefficients, col="green")
# Gráfica datos transformados
plot(datos$lnDiam, datos$lnVol, main="Datos transformados",xlab="", ylab="", col="brown")
abline(fit2$coefficients, col="red")
Checamos ahora los supuestos de linealidad y homocedasticidad con el modelo ajustado con los datos transformados
# Linealidad
plot(fit2, 1)
luego
# homocedasticidad
plot(fit2, 3)
concluyendo que es mejor trabajar con los datos transformados.
Descripción del summary
Considerando que
summary(fit2)##
## Call:
## lm(formula = lnVol ~ lnDiam, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3323 -0.1131 0.0267 0.1177 0.4280
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.8718 0.1216 -23.63 <2e-16 ***
## lnDiam 2.5644 0.0512 50.09 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1703 on 68 degrees of freedom
## Multiple R-squared: 0.9736, Adjusted R-squared: 0.9732
## F-statistic: 2509 on 1 and 68 DF, p-value: < 2.2e-16tenemos que
- Coefficients: encontramos los valores estimados de \(\hat{\beta_{0}}\) y \(\hat{\beta_{1}}\), además de \(\sqrt{\hat{V}(\hat{\beta_{0}})}\) y \(\sqrt{\hat{V}(\hat{\beta_{1}})}\) en la columna Std.error.
-
Coefficients: también, encontramos los valores de las estadísticas \(t\) asociadas a las pruebas de hipótesis
\[ H_{0}: \beta_{0} \ \ vs\ \ H_{a}: \beta_{0}\neq0 \]
y
\[ H_{0}: \beta_{1} \ \ vs\ \ H_{a}: \beta_{1}\neq0 \] y los_p-values_ correspondientes.
- En el último apartado nos arroja los valores de \(\hat{\sigma}\), el coeficiente de determinación \(R^{2}\) y los resultados asociados a la prueba \(F\) de la tabla anova.
Interpretación de los coeficientes
Para \(\hat{\beta_{1}}\): Por cada unidad log(pulgada) adicional al diámetro de un árbol, el número promedio de log(pies cúbicos) de volumen de un árbol aumenta 2.5644.
Para \(\hat{\beta_{0}}\): Para 0 unidades de log(pulgadas) del diámetro de un árbol en promedio se tienen -2.87 log(pies cúbicos) de volumen.
Sin embargo, a veces será necesario regresar a nuestros datos a su escala original pues la transformación no siempre es fácil de entender. Para ello, partiendo de \(ln(Vol)=-2.87+5.56\cdot ln(Diam)\) tenemos que \(Vol=e^{\hat{\beta_{0}}}\cdot Diam^{\hat{\beta_{1}}}\). Es importante recordar que aunque regresamos a la escala original NO necesariamente estamos modelando la \(\mathbb{E}[y;x]\), de hecho, si el supuesto de normalidad se cumpliera estaríamos modelando la \(Med(Vol; Diam)\), ya que la mediana es invariante bajo transformaciones biyectivas.
Notamos que tenemos un modelo multiplicativo, por lo que nos conviene tomar el cociente en lugar de la diferencia para llevar a cabo la interpretación . Es decir, si consideramos un aumento proporcional en \(\alpha\) unidades al diámetro tenemos
\[ \frac{Med(Y;\alpha Diam)}{Med(Y;Diam)}=\alpha^{2.56}=\alpha^{\hat{\beta_{1}}} \]
Si \(\alpha=2\), por ejemplo, entonces \(\alpha^{2.56}\approx 5.92\) por que se puede interpretar en la escala original que: al aumentar dos veces el diámetro de un árbol la mediana estimada del volumen del volumen cambia en un factor de 5.92. Por ejemplo, la mediana del volumen de un árbol de 20 pulgadas de diámetro, se estima que sea 5.92 veces más grande que la mediana del volumen de un árbol de 10 pulgadas de diámetro.
Finalmente, podemos graficar nuestros datos no transformados y ajustarle una curva asociada a la recta ajustada:
# Gráfico datos no transformados
par(mar=c(4, 5, 1, 1))
plot(datos$Diam, datos$Vol, xlab="Diametro (pulgadas)", ylab="Volumen (pies cúbicos)",
col="blue")
# función para crear la curva
f<- function(diam)
{exp(fit2$coefficients[1])*diam**{fit2$coefficients[2]}}
# Graficamos la curva
curve(f, from=min(datos$Diam), to = max(datos$Diam), col="red", add=T)
donde fit2$coefficients[1] es \(\hat{\beta_{0}}\) y fit2$coefficients[2] es \(\hat{\beta_{1}}\).