Taller 1. Regresi贸n lineal

#install.packages("devtools") # solo la primera vez
#devtools::install_github("dgonxalex80/paqueteMOD", force =TRUE)
library(paqueteMOD)
data(arboles)
head(arboles)
##   id  peso diametro altura
## 1  1 13.73      4.7    5.0
## 2  2 14.58      5.3    5.6
## 3  3 15.88      4.8    5.8
## 4  4  8.99      3.2    4.3
## 5  5  6.99      2.2    3.3
## 6  6 19.34      6.3    7.9

1. An谩lisis descriptivo

En las gr谩ficas A y B se puede identificar que el peso de los 谩rboles tiene una mayor concentraci贸n en los pesos de 5 a 25 toneladas, esta distribuci贸n tiene una asimetr铆a positiva que indica que hay pocos arboles con pesos superiores a 25 toneladas. El promedio del peso de los 谩rboles es de 18.8 toneladas y la mediana que corresponde a 17.5 toneladas, por otra parte, la desviaci贸n est谩ndar es de 8.16 toneladas.

En cuanto a los di谩metros, se puede observar en los gr谩ficos C y D que esta variable tiene una distribuci贸n normal. La mayor铆a de los 谩rboles tiene un di谩metro de 4.5 a 6.5 metros, con un promedio de di谩metro de 5.45 metros y una desviaci贸n est谩ndar de 1.45 metros. Con respecto a la altura, se identifica en los grafico E y F que la altura de los arboles tiene una mayor concentraci贸n en los valores de 5 y 7.5 metros, esta distribuci贸n tambi茅n tiene una asimetr铆a positiva, con pocos arboles que tienen altura superior a los 10 metros. El promedio de altura de los arboles es de 6.63 metros y desviaci贸n est谩ndar de 1.80 metros.

#install.packages("ggplot2")
#install.packages("ggpubr")
require(ggplot2)
## Loading required package: ggplot2
require(ggpubr)
## Loading required package: ggpubr
g1=ggplot(arboles,aes(x=peso))+geom_histogram(bins=30)+theme_bw()
g3=ggplot(arboles, aes(x=peso))+geom_boxplot(width=0.5)+theme_bw()
g2=ggplot(arboles,aes(x=diametro))+geom_histogram(bins=30)+theme_bw()
g4=ggplot(arboles, aes(x=diametro))+geom_boxplot(width=0.5)+theme_bw()
g5=ggplot(arboles,aes(x=altura))+geom_histogram(bins=30)+theme_bw()
g6=ggplot(arboles, aes(x=altura))+geom_boxplot(width=0.5)+theme_bw()
ggarrange(g1, g3, g2, g4, g5, g6, labels = c("A", "B", "C", "D", "E", "F"),ncol = 2, nrow = 3)

y <- table1::table1(~peso+diametro+altura, data = arboles)
y
Overall
(N=90)
peso
Mean (SD) 18.8 (8.16)
Median [Min, Max] 17.5 [5.98, 47.9]
diametro
Mean (SD) 5.45 (1.45)
Median [Min, Max] 5.40 [2.20, 8.80]
altura
Mean (SD) 6.63 (1.80)
Median [Min, Max] 6.45 [3.30, 11.3]

2. Pruebas de normalidad

Se realizan pruebas de normalidad sobre las variables Y y X, donde se obtuvo que los datos observados para la variable Y no siguen una distribuci贸n normal, ya que se rechaza la hip贸tesis nula con nivel de confianza del 95%. Para las variables predictoras no se rechaza la hip贸tesis nula y se puede concluir que los datos observados del di谩metro y la altura siguen una distribuci贸n normal.

shapiro.test(arboles$peso)
## 
##  Shapiro-Wilk normality test
## 
## data:  arboles$peso
## W = 0.92778, p-value = 9.028e-05
shapiro.test(arboles$diametro)
## 
##  Shapiro-Wilk normality test
## 
## data:  arboles$diametro
## W = 0.99109, p-value = 0.8079
shapiro.test(arboles$altura)
## 
##  Shapiro-Wilk normality test
## 
## data:  arboles$altura
## W = 0.97909, p-value = 0.156

3. Gr谩ficos de dispersi贸n

En el gr谩fico de dispersi贸n del peso de los 谩rboles vs el di谩metro se puede identificar una relaci贸n positiva; sin embargo, se observa un crecimiento exponencial del peso de los 谩rboles ante el aumento del di谩metro, con lo cual es muy posible que debamos realizar ciertos ajustes a las variables para poder calcular un modelo de regresi贸n lineal que cumpla con los supuestos de los errores. En cuanto al grafico de dispersi贸n del peso de los arboles vs la altura tambi茅n se identifica una relaci贸n positiva; sin embargo, se observa un cambio de nivel en la relaci贸n del peso y la altura, en los valores mas altos de la variable predictora altura.

g3=ggplot(arboles,aes(y=peso,x=diametro))+geom_point()+theme_bw()+geom_smooth()
g3
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

g3=ggplot(arboles,aes(y=peso,x=altura))+geom_point()+theme_bw()+geom_smooth()
g3
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

4. Coeficientes de correlaci贸n

De acuerdo con los coeficientes de correlaci贸n que corresponden a 90.8% y 85.8%, se puede concluir que el peso de los 谩rboles y el di谩metro, y el peso de los 谩rboles y la altura tienen una relaci贸n lineal positiva fuerte.

cor.test(arboles$peso, arboles$diametro)
## 
##  Pearson's product-moment correlation
## 
## data:  arboles$peso and arboles$diametro
## t = 20.346, df = 88, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8634081 0.9386817
## sample estimates:
##      cor 
## 0.908123
cor.test(arboles$peso, arboles$altura)
## 
##  Pearson's product-moment correlation
## 
## data:  arboles$peso and arboles$altura
## t = 15.684, df = 88, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7918402 0.9045332
## sample estimates:
##       cor 
## 0.8582009

5. Modelo de regresi贸n lineal

De acuerdo con el P-valor de 饾憦2, no se rechaza la hip贸tesis nula con un 95% de confianza, lo que nos indica que la variable altura no es significativa, y que el peso de los 谩rboles no depende de esta variable. Por esto, se calcula un segundo modelo sin incluir esta variable.

modelo1=lm(peso ~ diametro+altura, arboles)
summary(modelo1)
## 
## Call:
## lm(formula = peso ~ diametro + altura, data = arboles)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3083 -2.5121  0.1608  2.0088 11.7446 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -9.1205     1.4305  -6.376 8.44e-09 ***
## diametro      4.7395     0.7128   6.649 2.49e-09 ***
## altura        0.3132     0.5751   0.544    0.587    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.449 on 87 degrees of freedom
## Multiple R-squared:  0.8253, Adjusted R-squared:  0.8213 
## F-statistic: 205.5 on 2 and 87 DF,  p-value: < 2.2e-16

6. Modelo de regresi贸n lineal: con las variables significativas

De acuerdo con el P-valor, se rechaza la hip贸tesis nula con un 95% de confianza, lo que nos indica que la variable di谩metro es significativa, y que el peso de los 谩rboles si depende de esta variable.

Se tiene que 饾憥= -9.0203. Este resultado indica que: Con un di谩metro de 0, el peso esperado de los 谩rboles seria de -9.0203 toneladas.

Se tiene que 饾憦1= 5.1026. Este resultado indica que: Por cada di谩metro adicional del 谩rbol se espera un aumento en 5.1026 toneladas en el peso del 谩rbol.

De acuerdo con el coeficiente de determinaci贸n ajustado, se puede concluir que el modelo explica en 82.27% la variabilidad del peso de los arboles. Este coeficiente indica que el modelo se ajusta bien a los datos, y que es un modelo confiable y tendr铆a un buen desempe帽o para realizar predicciones.

modelo2=lm(peso ~ diametro, arboles)
summary(modelo2)
## 
## Call:
## lm(formula = peso ~ diametro, data = arboles)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3775 -2.6594  0.0237  1.8758 11.9876 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -9.0203     1.4129  -6.384 7.86e-09 ***
## diametro      5.1026     0.2508  20.346  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.435 on 88 degrees of freedom
## Multiple R-squared:  0.8247, Adjusted R-squared:  0.8227 
## F-statistic:   414 on 1 and 88 DF,  p-value: < 2.2e-16

6. Validaci贸n de los supuestos de los errores

- Media cero: En la primera gr谩fica (residuales vs ajustados) se observa que los errores no tienen un comportamiento aleatorio con media cero.

- Varianza constante: En esta misma gr谩fica se identifica un aumento de la varianza, por lo tanto, se puede inferir que los errores no cumplen con la propiedad de varianza constante. Adicionalmente, de acuerdo con las diferentes pruebas que eval煤a si la varianza constante, se puede concluir que los errores del modelo de regresi贸n lineal no tienen varianza constante ya que en todos los casos se rechaza la hip贸tesis nula.

- Independencia: Con base en la prueba estadistica de Durbin-Watson, se identifica que los errores del modelo de regresi贸n lineal no cumplen con la propiedad de independencia ya que se rechaza la hip贸tesis nula. Sin embargo, se conoce de antemano que los datos se tomaron con corte transversal y no dependen del tiempo, por esta raz贸n, los resultados de esta prueba no se consideran para evaluar el supuesto de independencia.

- Normalidad: De acuerdo con el segundo gr谩fico (Q-Q Plot), se identifica que los errores en su mayor铆a siguen una distribuci贸n normal; sin embargo, para algunos datos extremos se evidencia que se alejan de la l铆nea de normalidad. En base a la prueba de shapiro-wilk, se puede concluir que los errores del modelo de regresi贸n lineal no siguen una distribuci贸n normal ya que se rechaza la hip贸tesis nula con un nivel de confianza del 95%.

En base a la validaci贸n de supuestos descrita anteriormente, se decide hacer una transformaci贸n a las variables, con el fin de ajustar el modelo de regresi贸n lineal y as铆 tener un mejor rendimiento en la validaci贸n de los supuestos de los errores.

A. Normalidad

par(mfrow = c(2, 2))
plot(modelo2)

e= modelo2$residuals
# Ho: los errores tienen distribuci贸n normal
shapiro.test(e)
## 
##  Shapiro-Wilk normality test
## 
## data:  e
## W = 0.95356, p-value = 0.002793

B. Independencia

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
# Ho: los errores son independientes (no estan relacionados)
dwtest(modelo2)
## 
##  Durbin-Watson test
## 
## data:  modelo2
## DW = 1.0719, p-value = 1.035e-06
## alternative hypothesis: true autocorrelation is greater than 0

C. Varianza constante

# Ho: la varianza de los erroes es constante
gqtest(modelo2)
## 
##  Goldfeld-Quandt test
## 
## data:  modelo2
## GQ = 2.0131, df1 = 43, df2 = 43, p-value = 0.01196
## alternative hypothesis: variance increases from segment 1 to 2
bptest(modelo2)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo2
## BP = 11.76, df = 1, p-value = 0.0006052
#install.packages("glmtoolbox")
library(glmtoolbox)
vdtest(modelo2)
## 
##              Score test for varying dispersion parameter
## 
##           Statistic =  17.09096 
##  degrees of freedom =  1 
##             p-value =  3.5632e-05

7. Transformaci贸n Box-Cox

A partir del m茅todo de box-cox se calculo el 位 optimo, para as铆 identificar cual era la transformaci贸n que se deb铆a realizar sobre la variable de respuesta Y. Se obtuvo un valor de 位 de 0.0202 y en base a esto debemos transformar la variable respuesta con la funci贸n logaritmo.

#install.packages("paqueteMOD")
library(MASS)
bc=boxcox(lm(peso ~ diametro+altura, arboles), lambda = -2:3)

#(lambda <- bc$x[which.max(bc$y)])

8. Transformaci贸n de la variable Y

Una vez realizada la transformaci贸n de la variable de respuesta (logaritmo), se realizo la prueba de shapiro-wilk sobre esta variable, en la cual no se rechaz贸 la hip贸tesis nula y se concluye que los datos observados transformados siguen una distribuci贸n normal. Tambi茅n se realizaron nuevamente los gr谩ficos de dispersi贸n, y en estos se ve una mejora en la relaci贸n lineal de las variables peso y di谩metro.

arboles2=arboles
arboles2$peso_transformado=log(arboles$peso)
shapiro.test(arboles2$peso_transformado)
## 
##  Shapiro-Wilk normality test
## 
## data:  arboles2$peso_transformado
## W = 0.99033, p-value = 0.7545
g3=ggplot(arboles2,aes(y=peso_transformado,x=diametro))+geom_point()+theme_bw()+geom_smooth()
g3
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

g3=ggplot(arboles2,aes(y=peso_transformado,x=altura))+geom_point()+theme_bw()+geom_smooth()
g3
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

9. Modelo de regresi贸n lineal: variable Y transformada

De acuerdo con los P-valores, se rechazan las hip贸tesis nulas con un 95% de confianza, lo que nos indica que las variables di谩metro y altura son significativas, y que el peso de los 谩rboles si depende de estas variables.

Se tiene que 饾憥= 1.34378. Este resultado indica que: Con un di谩metro y altura de 0, el peso esperado de los 谩rboles seria de 3.8335 toneladas.

Se tiene que 饾憦1= 0.33539. Este resultado indica que: Por cada metro adicional en el di谩metro del 谩rbol se espera un aumento promedio de 1.3984 toneladas en el peso del 谩rbol.

Se tiene que 饾憦2= -0.04934. Este resultado indica que: Por cada metro adicional en la altura del 谩rbol se espera un aumento promedio de 0.9518 toneladas en el peso del 谩rbol.

De acuerdo con el coeficiente de determinaci贸n ajustado, se puede concluir que el modelo explica en 88.94% la variabilidad del peso de los arboles. Este coeficiente indica que el modelo se ajusta bien a los datos, y que es un modelo confiable y tendr铆a un buen desempe帽o para realizar predicciones.

modelo3=lm(arboles2$peso_transformado ~ arboles2$diametro+ arboles2$altura)
summary(modelo3)
## 
## Call:
## lm(formula = arboles2$peso_transformado ~ arboles2$diametro + 
##     arboles2$altura)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.2722 -0.1168 -0.0100  0.1018  0.3198 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        1.34378    0.05917  22.710   <2e-16 ***
## arboles2$diametro  0.33539    0.02949  11.374   <2e-16 ***
## arboles2$altura   -0.04934    0.02379  -2.074   0.0411 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1427 on 87 degrees of freedom
## Multiple R-squared:  0.8919, Adjusted R-squared:  0.8894 
## F-statistic: 358.8 on 2 and 87 DF,  p-value: < 2.2e-16

10. Validaci贸n de los supuestos de los errores

- Media cero: En la primera gr谩fica (residuales vs ajustados) se observa que los errores tienen un comportamiento aleatorio con media cero.

- Varianza constante: En esta misma gr谩fica se identifica que la varianza de los errores no aumenta, por lo tanto, se puede inferir que los errores cumplen con la propiedad de varianza constante. Adicionalmente, de acuerdo con las diferentes pruebas que eval煤a si la varianza es constante, se puede concluir que los errores del modelo de regresi贸n lineal tienen varianza constante ya que en todos los casos no se rechaza la hip贸tesis nula con un nivel de confianza del 95%.

- Independencia: Con base en la prueba estadistica de Durbin-Watson, se identifica que los errores del modelo de regresi贸n lineal con la variable de respuesta transformada no cumplen con la propiedad de independencia ya que se rechaza la hip贸tesis nula. Sin embargo, se conoce de antemano que los datos se tomaron con corte transversal y no dependen del tiempo, por esta raz贸n, los resultados de esta prueba no se consideran para evaluar el supuesto de independencia.

- Normalidad: De acuerdo con el segundo gr谩fico (Q-Q Plot), se identifica que los errores siguen una distribuci贸n normal; solo para algunos datos extremos se evidencia que se alejan de la l铆nea de normalidad. En base a la prueba de shapiro-wilk, se puede concluir que los errores del modelo de regresi贸n lineal siguen una distribuci贸n normal ya que no se rechaza la hip贸tesis nula con un nivel de confianza del 95%.

A. Normalidad

par(mfrow = c(2, 2))
plot(modelo3)

e= modelo3$residuals
# Ho: los errores tienen distribuci贸n normal
shapiro.test(e)
## 
##  Shapiro-Wilk normality test
## 
## data:  e
## W = 0.97487, p-value = 0.07832

B. Independencia

library(lmtest)
# Ho: los errores son independientes (no estan relacionados)
dwtest(modelo3)
## 
##  Durbin-Watson test
## 
## data:  modelo3
## DW = 0.83921, p-value = 2.663e-10
## alternative hypothesis: true autocorrelation is greater than 0

C. Varianza de los errores

# Ho: la varianza de los erroes es constante
gqtest(modelo3)
## 
##  Goldfeld-Quandt test
## 
## data:  modelo3
## GQ = 1.0399, df1 = 42, df2 = 42, p-value = 0.4499
## alternative hypothesis: variance increases from segment 1 to 2
bptest(modelo3)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo3
## BP = 2.4797, df = 2, p-value = 0.2894
library(glmtoolbox)
vdtest(modelo3)
## 
##              Score test for varying dispersion parameter
## 
##           Statistic =  1.63094 
##  degrees of freedom =  2 
##             p-value =  0.44243