Utilizaremos los datos “trees”. Se trata de un conjunto de datos que consta de 31 observaciones sobre mediciones de la circunferencia, altura y volumen de la madera sobre cerezos negros talados.
Descripción del conjunto de datos:
Girth: Diámetro en pulgadas del árbol.
Height: Altura en pies del árbol.
Volume: Volumen de la madera en pies cúbicos.
Se desea estimar un modelo de regresión dque estime el volumen (Volume) en función de la circunferencia del árbol (Girth) y la altura del árbol (Height); para ello:
1 - Realice un análisis exploratorio y visualice la relación entre las tres variables.
2 - Realice un modelo de regresión múltiple por medio de la función lm() solo con efectos principales. Interprete la bondad de ajuste del modelo. Calcule la tasa de error del modelo.
3 - Interprete los coeficientes del modelo.
4 - Calcule los intervalos de confianza para los coeficientes.
5 - Calcule la importancia relativa de los predictores.
6 - Evalué los supuestos del modelo.
7 - Realice un modelo de regresión con efectos de interacción entre los 2 predictores. Evalúe posibles problemas de multicolinealidad y en tal caso centre los predictores.
8 - Interprete los coeficientes del modelo.
EN CADA UNO DE LOS PASOS INTERPRETE LOS RESULTADOS
Cargamos y revisamos los datos:
data("trees")
head(trees)
## Girth Height Volume
## 1 8.3 70 10.3
## 2 8.6 65 10.3
## 3 8.8 63 10.2
## 4 10.5 72 16.4
## 5 10.7 81 18.8
## 6 10.8 83 19.7
str(trees)
## 'data.frame': 31 obs. of 3 variables:
## $ Girth : num 8.3 8.6 8.8 10.5 10.7 10.8 11 11 11.1 11.2 ...
## $ Height: num 70 65 63 72 81 83 66 75 80 75 ...
## $ Volume: num 10.3 10.3 10.2 16.4 18.8 19.7 15.6 18.2 22.6 19.9 ...
library(GGally)
ggpairs(trees[, c("Girth", "Height", "Volume")])
En la última fila de gráficos se observa la variable volume en el eje y y los predictores en el eje x. Existe una relación fuerte, positiva y muy lineal entre el Volume y el Girth, mientras que la relación entre el volume y el Height es fuerte, aunque más cerca del rango de una relación moderada, positiva y no menos lineal.
model_trees <- lm(Volume ~ Girth + Height, data = trees)
summary(model_trees)
##
## Call:
## lm(formula = Volume ~ Girth + Height, data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.4065 -2.6493 -0.2876 2.2003 8.4847
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -57.9877 8.6382 -6.713 2.75e-07 ***
## Girth 4.7082 0.2643 17.816 < 2e-16 ***
## Height 0.3393 0.1302 2.607 0.0145 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.882 on 28 degrees of freedom
## Multiple R-squared: 0.948, Adjusted R-squared: 0.9442
## F-statistic: 255 on 2 and 28 DF, p-value: < 2.2e-16
Al interpretar la bondad de ajuste según los resultados de la prueba F global, el modelo de regresión predice el volumen de la madera del árbol, en función del diámetro y la altura, significativamente bien (F (2, 28) = 255, p−value < .001). El modelo funciona bien para los datos disponibles.
sigma(model_trees)
## [1] 3.881832
sigma(model_trees)/mean(trees$Volume)
## [1] 0.1286612
Al continuar con la evaluación de la calidad general del modelo, el error estándar residual es de 3.88 pies cúbicos, con una tasa error de 12.87%.
summary(model_trees)$adj.r.squared
## [1] 0.9442322
El coeficiente de determinación ajustado es R2 = 94.42, por lo que el modelo explica el 94.42% de la variabilidad de la respuesta, mediante los datos de los demas atributos del arbol (Girth y Height).
library(rsq)
rsq.partial(model_trees)
## $adjustment
## [1] FALSE
##
## $variable
## [1] "Girth" "Height"
##
## $partial.rsq
## [1] 0.9189376 0.1952712
Cuando vemos el porcentaje de variación de la respuesta que es explicada por cada predictor, de forma independiente, el Girth explica un 91.9% de variación, mientras que el Height explica solo un 19.52%.
summary(model_trees)$coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -57.9876589 8.6382259 -6.712913 2.749507e-07
## Girth 4.7081605 0.2642646 17.816084 8.223304e-17
## Height 0.3392512 0.1301512 2.606594 1.449097e-02
El intercepto vale -57.9876589, en este caso no tiene sentido interpretar su valor ya que es imposible estimar el volumen del árbol en pies cúbicos, para un árbol con 0 pulgadas de diámetro y 0 pies de altura.
El coeficiente de regresión para el predictor Girth es de 4.7081605. Este valor representa el cambio en volumen medio en pies cúbicos del árbol al aumentar 1 pulgada de diámetro, cuando la altura media permanece constante. Según la prueba T, obtenemos un coeficiente significativamente distinto de 0 con un p-value < 0.001.
El coeficiente de regresión para el predictor Height vale 0.3392512. Este valor representa el cambio en volumen medio en pies cúbicos del árbol al aumentar 1 pie de altura, cuando el diámetro medio permanece constante. El resultado de la prueba T nos indica que el coeficiente es significativamente distinto de 0 con un p-value < 0.05.
confint(model_trees)
## 2.5 % 97.5 %
## (Intercept) -75.68226247 -40.2930554
## Girth 4.16683899 5.2494820
## Height 0.07264863 0.6058538
Con un 95% de confianza podemos decir que, en la población, el intercepto tomará un valor entre [-75.68226247, -40.2930554], el coeficiente para Girth estará comprendido entre [4.16683899, 5.2494820] y el de Height entre [0.07264863, 0.6058538].
library(relaimpo)
crlm_trees <- calc.relimp(model_trees,
type = c("lmg"),
rela = TRUE)
crlm_trees
## Response variable: Volume
## Total response variance: 270.2028
## Analysis based on 31 observations
##
## 2 Regressors:
## Girth Height
## Proportion of variance explained by model: 94.8%
## Metrics are normalized to sum to 100% (rela=TRUE).
##
## Relative importance metrics:
##
## lmg
## Girth 0.804561
## Height 0.195439
##
## Average coefficients for different model sizes:
##
## 1X 2Xs
## Girth 5.065856 4.7081605
## Height 1.543350 0.3392512
plot(crlm_trees)
Según el método LMG, la variable Girth (80.46%) tiene mucha más importancia relativa en el modelo en comparación con Height (19.54%).
par(mfrow = c(2,2))
plot(model_trees)
1 - Los residuos presentan un patrón sesgado y heteroscedástico, por lo que no se cumple el supuesto de linealidad y homocedasticidad.
2 - Al ver la distribución de los residuos en el gráfico cuantil-cuantil normal, vemos que los residuos tienen una distribución normal, por lo que se cumple el supuesto de normalidad.
3 - La varianza de los residuos no es constante. La variación de los residuos cambia a lo largo de los valores ajustados y no se reparten de manera aleatoria alrededor de la linea, por lo que observa un patrón que implica que no se cumple el supuesto de homocedasticidad.
4 - Cuando observamos el gráfico, se puede identificar un valor influyente que pasa la distancia de Cook arriba a la derecha, correspondiente a la observación numero 31, por lo que tampoco se cumple este supuesto.
model_trees2 <- lm(Volume ~ Girth*Height, data = trees)
library(car)
vif(model_trees2)
## Girth Height Girth:Height
## 148.66145 15.93884 210.97302
Al analizar los factores de inflación de varianza, los VIF son > 5, lo que implica un problema de multicolinealidad entre los predictores. Esto ocurre porque las variables no están centradas y hemos incluido un término de interacción.
Girth2 <- scale(trees$Girth, center = TRUE, scale = FALSE)
Height2 <- scale(trees$Height, center = TRUE, scale = FALSE)
model_trees3 <- lm(trees$Volume ~ Girth2*Height2)
vif(model_trees3)
## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif
## Girth2 Height2 Girth2:Height2
## 1.513180 1.487784 1.126849
Luego de centrar las variables, los VIF son menores que 5, por lo que se determina que ya no hay problemas de multicolinealidad.
Anova(model_trees3)
## Anova Table (Type II tests)
##
## Response: trees$Volume
## Sum Sq Df F value Pr(>F)
## Girth2 4783.0 1 651.965 < 2.2e-16 ***
## Height2 102.4 1 13.956 0.0008867 ***
## Girth2:Height2 223.8 1 30.512 7.484e-06 ***
## Residuals 198.1 27
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Al evaluar nuevamente la significación, vemos que la interacción es significativa (F (1, 27) = 30.512, p < .001) y por tanto debe incluirse en el modelo final.
summary(model_trees3)
##
## Call:
## lm(formula = trees$Volume ~ Girth2 * Height2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.5821 -1.0673 0.3026 1.5641 4.6649
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.81791 0.54466 52.910 < 2e-16 ***
## Girth2 4.37789 0.19384 22.585 < 2e-16 ***
## Height2 0.48687 0.09466 5.143 2.07e-05 ***
## Girth2:Height2 0.13465 0.02438 5.524 7.48e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.709 on 27 degrees of freedom
## Multiple R-squared: 0.9756, Adjusted R-squared: 0.9728
## F-statistic: 359.3 on 3 and 27 DF, p-value: < 2.2e-16
Al ver el resumen final, se determinó que el coeficiente de interacción Girth2:Height2 es significativo con un valor de 0.13465 y un p-value < 0.001.
Cuando Height2 = 0, un aumento de 1 pulgada en el Girth2 incrementa el volumen en 4.37789 pies cúbicos.
Por cada aumento de 1 pie en la altura respecto al promedio, el efecto del diametro sobre el volumen aumenta en 0.13465 pies cúbicos.
sigma(model_trees3)
## [1] 2.70855
sigma(model_trees3)/mean(trees$Volume)
## [1] 0.08977339
Al analizar nuevamente la calidad general del modelo, el error estándar residual se redujo de 3.88 a 2.70 pies cúbicos, con una reducción en la tasa de error del 12.87% anterior a 8.98%.
summary(model_trees3)$adj.r.squared
## [1] 0.9728491
Asi mismo, el coeficiente de determinación ajustado aumentó de R2 = 94.42 a 97.28, por lo que el modelo final explica el 97.28% de la variabilidad de la respuesta.
El modelo final es:
Volume = 28.81791 + 4.37789(Girth2) + 0.48687(Height2) + 0.13465(Girth2 × Height2)