Taller1
John Jairo Bedoya Saenz
2023-02-22
Importar datos biomasa
#install.packages("devtools") # solo la primera vez
#devtools::install_github("dgonxalex80/paqueteMOD", force =TRUE)
library(paqueteMOD)
data(arboles)
head(arboles,10)## 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
## 7 7 21.44 6.6 8.3
## 8 8 13.81 5.3 7.3
## 9 9 11.88 4.9 6.7
## 10 10 16.62 5.9 7.1attach(arboles)Exploración univariada de datos
summarytools::descr(arboles[,2:4])## Descriptive Statistics
## arboles
## N: 90
##
## altura diametro peso
## ----------------- -------- ---------- --------
## Mean 6.63 5.45 18.77
## Std.Dev 1.80 1.45 8.16
## Min 3.30 2.20 5.98
## Q1 5.20 4.50 13.61
## Median 6.45 5.40 17.48
## Q3 7.90 6.50 22.82
## Max 11.30 8.80 47.87
## MAD 1.93 1.48 6.35
## IQR 2.65 1.97 9.16
## CV 0.27 0.27 0.43
## Skewness 0.38 -0.07 1.11
## SE.Skewness 0.25 0.25 0.25
## Kurtosis -0.44 -0.36 1.55
## N.Valid 90.00 90.00 90.00
## Pct.Valid 100.00 100.00 100.00Normalidad de las variables
# Diametro
shapiro.test(arboles$diametro)##
## Shapiro-Wilk normality test
##
## data: arboles$diametro
## W = 0.99109, p-value = 0.8079# Altura
shapiro.test(arboles$altura)##
## Shapiro-Wilk normality test
##
## data: arboles$altura
## W = 0.97909, p-value = 0.156H0: La distribución es normal H1: La distribución no es normal alpha: 0.05
Se rechaza H0 si el p-valor de la prueba es <= a 0.05
Diametro: p-value = 0.8079 ; No se rechaza la hipotesis nula Altura: p-value = 0.156 ; No se rechaza la hipotesis nula
require(ggplot2)
require(ggpubr)
bp1=ggplot(arboles,aes(x=peso))+geom_boxplot(fill="gray", color=1)+ggtitle("Gráfico de Cajas - Peso")+theme_bw()
bp2=ggplot(arboles,aes(x=diametro))+geom_boxplot(fill="lightblue", color=1)+ggtitle("Gráfico de Cajas - Diametro")+theme_bw()
bp3=ggplot(arboles,aes(x=altura))+geom_boxplot(fill="lightgreen", color=1)+ggtitle("Gráfico de Cajas - Altura")+theme_bw()
arboles$diametro_grupo=cut(arboles$diametro,breaks = c(2,4,6,8,10))
arboles$altura_grupo=cut(arboles$altura,breaks = c(3,4,5,6,7,8,9,10,11,12))
g1=ggplot(arboles,aes(x=diametro_grupo))+geom_bar(fill="lightblue", color=2)+geom_text(stat='count', aes(x=diametro_grupo,label=..count..),vjust=0.8, col="black")+theme_bw()+ggtitle("Gráfico de Barras - Diametro")
g2=ggplot(arboles,aes(x=altura_grupo))+geom_bar(fill="lightgreen", color=4)+geom_text(stat='count', aes(x=altura_grupo,label=..count..),vjust=0.8, col="black")+theme_bw()+ggtitle("Gráfico de Barras - Altura")
D1 = ggplot(arboles, aes(x=diametro))+geom_density()+ggtitle("Gráfico de Densidad - Diametro")
D2 = ggplot(arboles, aes(x=altura))+geom_density()+ggtitle("Gráfico de Densidad - Altura")
D3 = ggplot(arboles, aes(x=peso))+geom_density()+ggtitle("Gráfico de Densidad - Peso")
ggarrange(bp2, g1, D1, bp3, g2, D2, bp1, D3, ncol = 2, nrow = 2)## $`1`
##
## $`2`
##
## attr(,"class")
## [1] "list" "ggarrange"Análisis univariado de datos - variables cuantitativas
Altura: El promedio de altura de los arboles es de 6.63 metros, con una desviación estandar de 1.8 metros, evidenciando que los metros de altura se desvian de la media aproximadamente en 1.8 metros, indicando que la media no es una medida certera en este caso. Es de notar que el 50% de los arboles tienen alturas por debajo de los 6.45 metros. Esto indica que la grafica es simétrica.
Diametro: El promedio del diametro de los arboles es de 5.45 metros, con una desviación estandar de 1.45 metros, evidenciando que los metros de diametro se desvian de la media aproximadamente en 1.45 metros, indicando que la media no es una medida certera en este caso. Es de notar que el 50% de los arboles tienen diametros inferiores a los 5.4 metros. Esto indica que la grafica es simétrica.
Análisis bivariable de datos
require(table1)
y <- table1::table1(~peso|diametro_grupo, data = arboles)
y| (2,4] (N=15) |
(4,6] (N=45) |
(6,8] (N=27) |
(8,10] (N=3) |
Overall (N=90) |
|
|---|---|---|---|---|---|
| peso | |||||
| Mean (SD) | 9.52 (2.68) | 16.2 (3.36) | 25.4 (4.97) | 43.7 (5.28) | 18.8 (8.16) |
| Median [Min, Max] | 8.99 [5.98, 14.9] | 15.9 [9.89, 24.5] | 23.8 [18.7, 34.4] | 45.4 [37.7, 47.9] | 17.5 [5.98, 47.9] |
y <- table1::table1(~peso|altura_grupo, data = arboles)
y| (3,4] (N=5) |
(4,5] (N=15) |
(5,6] (N=20) |
(6,7] (N=12) |
(7,8] (N=19) |
(8,9] (N=9) |
(9,10] (N=6) |
(10,11] (N=3) |
(11,12] (N=1) |
Overall (N=90) |
|
|---|---|---|---|---|---|---|---|---|---|---|
| peso | ||||||||||
| Mean (SD) | 7.07 (0.706) | 12.2 (2.47) | 15.9 (2.81) | 16.2 (4.68) | 20.6 (3.95) | 24.9 (6.20) | 31.9 (7.50) | 34.0 (3.49) | 47.9 (NA) | 18.8 (8.16) |
| Median [Min, Max] | 7.03 [5.98, 7.87] | 13.2 [6.98, 15.4] | 15.9 [8.73, 21.1] | 17.1 [9.89, 22.8] | 20.4 [13.8, 27.5] | 25.2 [16.2, 34.4] | 31.5 [23.8, 45.4] | 33.4 [30.8, 37.7] | 47.9 [47.9, 47.9] | 17.5 [5.98, 47.9] |
Peso con respecto al diametro: Se evidencia que de las 15 especies de arboles que tienen diametro entre los 2 y los 4 metros tienen en promedio un peso de 9.52 kg con una desviacion estandar de 2.68 kg, las 45 especies de arboles que tienen diametro dentre los 4 y 6 metros tienen en promedio un peso de 16.2 kg con una desvest de 3.36 kg, las 27 especies de arboles que tienen diametro dentre los 6 y 8 metros tienen en promedio un peso de 25.4 kg con una desvest de 4.97 kg, y las 3 especies de arboles que tienen diametro dentre los 8 y 10 metros tienen en promedio un peso de 43.7 kg con una desvest de 5.28 kg.
Peso con respecto a la altura: Se evidencia que de las 5 especies de arboles que tienen altura entre los 3 y los 4 metros tienen en promedio un peso de 7.07 kg con una desviacion estandar de 0.706 kg, 15 especies entre 4 y 5 metros de altura pesan en promedio 12.2 kg con desvest de 2.47kg, 20 especies entre 5 y 6 metros de altura pesan en promedio 15.9 kg con desvest de 2.81 kg, 12 especies entre 6 y 7 metros de altura pesan en promedio 16.2 kg con desvest de 4.68 kg, 19 especies entre 7 y 8 metros de altura pesan en promedio 20.6 kg con desvest de 3.95 kg, 9 especies entre 8 y 9 metros de altura pesan en promedio 24.9 kg con desvest de 6.2 kg, 6 especies entre 9 y 10 metros de altura pesan en promedio 31.9 kg con desvest de 7.5 kg, 3 especies entre 10 y 11 metros de altura pesan en promedio 34 kg con desvest de 3.49 kg, y 1 especie entre 11 y 12 metros que pesa 47.9 kg.
Correlacion de variables
g1=ggplot(arboles,aes(x=diametro,y=peso))+geom_point()+theme_bw()+geom_smooth()
g2=ggplot(arboles,aes(x=altura,y=peso))+geom_point()+theme_bw()+geom_smooth()
ggarrange(g1,g2,ncol = 2, nrow = 1)
De acuerdo con el grafico de dispresión, donde se grafican las variables diametro (eje x) y peso (eje y), se logra percibir que la dispersion de los datos tiene una tendencia directa o creciente (relacion lineal positiva).
Coeficiente de correlacion
library(kableExtra)
G=matrix(c(cor.test(arboles$diametro,arboles$peso)$estimate,cor.test(arboles$altura,arboles$peso)$estimate), ncol=1)
rownames(G)<- c("Correlacion peso/diametro","Correlacion peso/altura")
G %>%
kbl(booktabs = T,) %>%
kable_classic_2(full_width = F)| Correlacion peso/diametro | 0.9081230 |
| Correlacion peso/altura | 0.8582009 |
Diametro: El resultado obtenido de coeficiente de correlación r = 0.908123 indica que entre el diametro y el peso del arbol existe una asociación lineal positiva fuerte, es decir, que entre mas sea el diametro de la especie, mayor será el peso.
Altura: El resultado obtenido de coeficiente de correlación r = 0.8582009 indica que entre la altura y el peso del arbol existe una asociación lineal positiva fuerte, es decir, que entre mas sea la altura de la especie, mayor será el peso.
Pruebas de hipótesis para verificar que las correlaciones entre las variables: peso, diametro y altura son diferentes de cero
H0: La correlación entre las variables peso, diametro y altura son iguales a cero
H1: La correlación entre las variables peso, diametro y altura son diferentes de cero
alpha = 0.05
NO Se rechaza H0 si el p-valor de la prueba es > a 0.05
Se rechaza H0 si el p-valor de la prueba es <= a 0.05
Diametro: p-value < 2.2e-16 Se rechaza la hipotesis nula, el diametro es diferente de 0, si hay correlación entre el diametro y el peso
Altura: p-value < 2.2e-16 Se rechaza la hipotesis nula, la altura es diferente de 0, si hay correlación entre la altura y el peso
Estimación del modelo por Minimos Cuadrados Ordinarios - MCO
X1 = diametro (variable independiente)
X2 = altura (variable independiente)
Y = peso (variable dependiente)
Ecuación de la recta de regresión: Y = B0 + B1.X1 + B2.X2 + e
# Diametro
library("stargazer")
mod1=lm(arboles$peso~ arboles$diametro)
stargazer(mod1, type = 'html')| Dependent variable: | |
| peso | |
| diametro | 5.103*** |
| (0.251) | |
| Constant | -9.020*** |
| (1.413) | |
| Observations | 90 |
| R2 | 0.825 |
| Adjusted R2 | 0.823 |
| Residual Std. Error | 3.435 (df = 88) |
| F Statistic | 413.961*** (df = 1; 88) |
| Note: | p<0.1; p<0.05; p<0.01 |
peso = -9.0203 + 5.1026 * diametro
B0 = -9.0203 para este valor de B0, nos indica que en el caso hipotetico que un arbol tenga 0 metros de diametro, este pesaría -9.0203 Kg aproximadamente
B1 = 5.1026 para el valor de B1, nos indica que por cada metro que incremente el diametro del arbol, se espera un aumento en 5,1026 kg en el peso del arbol.
# Altura
mod2=lm(arboles$peso~arboles$altura)
stargazer(mod2, type = 'html')| Dependent variable: | |
| peso | |
| altura | 3.891*** |
| (0.248) | |
| Constant | -7.046*** |
| (1.705) | |
| Observations | 90 |
| R2 | 0.737 |
| Adjusted R2 | 0.734 |
| Residual Std. Error | 4.211 (df = 88) |
| F Statistic | 245.977*** (df = 1; 88) |
| Note: | p<0.1; p<0.05; p<0.01 |
peso = -7.0456 + 3.8906 * altura
B0 = -7.0456 para este valor de B0, nos indica que en el caso hipotetico que un arbol tenga 0 metros de altura, este pesaría -7.0456 Kg aproximadamente
B1 = 3.8906 para el valor de B1, nos indica que por cada metro que incremente la altura del arbol, se espera un aumento en 3.8906 kg en el peso del arbol.
# Altura + Diametro
mod3=lm(arboles$peso~arboles$altura+arboles$diametro) # variable de respuesta en función de la altura y el diametro
stargazer(mod2, type = 'html')| Dependent variable: | |
| peso | |
| altura | 3.891*** |
| (0.248) | |
| Constant | -7.046*** |
| (1.705) | |
| Observations | 90 |
| R2 | 0.737 |
| Adjusted R2 | 0.734 |
| Residual Std. Error | 4.211 (df = 88) |
| F Statistic | 245.977*** (df = 1; 88) |
| Note: | p<0.1; p<0.05; p<0.01 |
\(y = \beta_0 + \beta_1X1 + \beta_2X2\) peso = -9.1205 + 0.3132 * altura + 4.7395 * diametro
\(\beta_0\) = -9.1205 para este valor de B0, nos indica que en el caso hipotetico que un arbol tenga 0 metros de altura y diametro, este pesaría -9.1205 Kg aproximadamente
\(\beta_1\) = 0.3132 para el valor de B1, nos indica que por cada metro que incremente la altura del arbol, se espera un aumento en 0.3132 kg en el peso del arbol.
\(\beta_2\) = 4.7395 para el valor de B2, nos indica que por cada metro que incremente el diametro del arbol, se espera un aumento en 4.7395 kg en el peso del arbol.
Se observa entonces que para el modelo, la altura no representa mucha significacia dado que su aporte en el peso del arbol por metro de altura no aporta valor en su peso, la variable entonces estadisticamente no es relavante en el modelo.
De acuerdo con los 3 modelos (\(peso~diametro\) “Adjusted R-squared: 0.8227”, \(peso~altura\) “Adjusted R-squared: 0.7335” y \(peso~diametro+altura\) “Adjusted R-squared: 0.8213”) obsevamos que el indicador de \(R^2\) para para el modelo 1 (mod1) tiene un ajuste de un 82.27% del peso o nos explica el peso en un 82.27%, obteniendo un \(R^2\) por encima del modelo 2 y modelo 3, sin embargo el modelo 3 nos explica el peso en un 82.13%, muy cercano al modelo 1, no obstante se continuará el ejercicio con el modelo 1.
Validación de los supuestos sobre los errores
par(mfrow = c(2, 2))
plot(mod1)
Mediante el uso de las graficas de los errores y algunos test estadisticos se evalúan los supuestos de los residuales para el modelo \(Peso=\beta_0 +\beta_1 Diametro\), los cuales son:
\(*\) Normalidad de los errores
\(*\) Homocedasticidad de los errores
\(*\) Ausencia de autocorrelación de los errores
\(*\) Media cero de los errores.
Supuesto de normalidad de los errores
shapiro.test(mod1$residuals)##
## Shapiro-Wilk normality test
##
## data: mod1$residuals
## W = 0.95356, p-value = 0.002793H0: los errores tienen una distribución normal
Se rechaza H0 si el p-valor de la prueba es <= a 0.05
p-value = 0.002793 se rechaza la hipotesis nula, los errores no tienen una distribución normal.
Supuesto de varianza constante (homocedasticidad)
library(lmtest)
gqtest(mod1)##
## Goldfeld-Quandt test
##
## data: mod1
## GQ = 2.0131, df1 = 43, df2 = 43, p-value = 0.01196
## alternative hypothesis: variance increases from segment 1 to 2H0: la varianza de los erroes es constante
p-value = 0.01196 se rechaza la hipotesis nula
Supuesto no autocorrelacion
En cuanto a la no autocorrelación de los errores, la violacion de este supuesto se da principalmente en muestras de serie de tiempo y no corte transversal, dado que no es posible tener autocorrelacion seriales en los errores si las observaciones se realizan todas en el mismo momento. Por lo tanto se asume que este supuesto se cumple y los errores no estarian correlacionados
Adicionalmente analizaremos la linealidad del modelo
ggplot(arboles, aes(diametro, peso)) +
geom_point() +
geom_smooth(method = "lm", level = 0.95, formula = y ~ x)
Dado que se puede evidenciar un problema en la linealidad, normalidad y homocedasticidad del modelo, lo que puede estar siguiriendo el analisis es que la forma funcional de la regresion no es la adecuada, y para esto realizaron una transformacion, de manera que podamos validad los supuestos de los residuales y poder realizar inferencias estadisticas sobre los coeficientes y los estadisticos respectivos.
Para poder encontrar la mejor transformacion realizaremos un \(Boxcox\):
##Transformacion del modelo
require(MASS)
bc=boxcox(lm(arboles$peso ~ arboles$diametro), lambda = -1:1)
(lambda <- bc$x[which.max(bc$y)])## [1] 0.07070707Con aproximado de a\(\lambda = 0\) al \(95\%\) de confianza, la transformacion propuesta por la metodologia de boxcox seria \(log(Peso)=\beta_0+\beta_1Diametro\).
Modelo log-lin
mod4=lm(log(arboles$peso)~arboles$diametro)
stargazer(mod4, type = 'html')| Dependent variable: | |
| peso) | |
| diametro | 0.278*** |
| (0.011) | |
| Constant | 1.328*** |
| (0.060) | |
| Observations | 90 |
| R2 | 0.887 |
| Adjusted R2 | 0.885 |
| Residual Std. Error | 0.145 (df = 88) |
| F Statistic | 687.562*** (df = 1; 88) |
| Note: | p<0.1; p<0.05; p<0.01 |
Con esta nueva especificacion del modelo, realizaremos una validación de los supuestos sobre los errores, para saber si los coeficientes son consistentes, insesgados y eficientes.
Validacion de los supuestos
par(mfrow = c(2, 2))
plot(mod4)
Para el modelo 4, donde se evalúan los residuales vs valores ajustados: Se observa completa aleatoriedad de los residuales, los datos no presentan algún comportamiento. Grafico de normalidad: se observa que los datos se encuentran cercanos a la linea de la distribución normal
Supuesto de normalidad de los errores
shapiro.test(mod4$residuals)##
## Shapiro-Wilk normality test
##
## data: mod4$residuals
## W = 0.98394, p-value = 0.3338H0: los errores tienen una distribución normal
Se rechaza H0 si el p-valor de la prueba es <= a 0.05
p-value = 0.3338 NO se rechaza la hipotesis nula, los errores tienen una distribución normal.
Supuesto de varianza constante
gqtest(mod4)##
## Goldfeld-Quandt test
##
## data: mod4
## GQ = 1.1538, df1 = 43, df2 = 43, p-value = 0.3206
## alternative hypothesis: variance increases from segment 1 to 2H0: la varianza de los erroes es constante p-value = 0.3206 NO se rechaza la hipotesis nula de homocedasticidad
A pesar de que hemos validados los supuestos de los errores y la regresion planteada es consistente, insesgada y eficiente, compararemos mediante el Criterio de información de Akaike y el \(R^2 Ajustado\) cual de las transformaciones mas populares en la literatura nos generan un mejor ajuste en los datos.
Otras tranformaciones:
#install.packages("stargazer")
mod5=lm(arboles$peso~log(arboles$diametro))
mod6=lm(log(arboles$peso)~log(arboles$diametro))
stargazer(mod1, mod4, mod5, mod6, type="html", title="Comparación de modelos", df=FALSE)| Dependent variable: | ||||
| peso | peso) | peso | peso) | |
| (1) | (2) | (3) | (4) | |
| diametro | 5.103*** | 0.278*** | ||
| (0.251) | (0.011) | |||
| diametro) | 23.369*** | 1.344*** | ||
| (1.564) | (0.058) | |||
| Constant | -9.020*** | 1.328*** | -19.909*** | 0.618*** |
| (1.413) | (0.060) | (2.629) | (0.098) | |
| Observations | 90 | 90 | 90 | 90 |
| R2 | 0.825 | 0.887 | 0.717 | 0.858 |
| Adjusted R2 | 0.823 | 0.885 | 0.714 | 0.857 |
| Residual Std. Error | 3.435 | 0.145 | 4.362 | 0.162 |
| F Statistic | 413.961*** | 687.562*** | 223.224*** | 532.232*** |
| Note: | p<0.1; p<0.05; p<0.01 | |||
Al comparar los valores de \(R^2\), el modelo 4 (mod4 = Modelo log-lin) presenta el mayor porcentaje de explicación de la variabilidad del peso.
Indicador AIC
AIC(mod1, mod4, mod5, mod6)## df AIC
## mod1 3 481.50130
## mod4 3 -87.82198
## mod5 3 524.52212
## mod6 3 -67.70739Criterio de información de Akaike - AIC: El modelo que tiene el mínimo AIC es el modelo 4 (mod4)
Dado que con ambos estadisticos se puede comprobar que el modelo \(log(Peso)=\beta_0+\beta_1Diametro\) es el mas adecuado, y ademas sus residuales se comportan de acuerdo a los supuestos de una regresion MCO, realizaremos un ejercicio de pronostico de la variable dependiente.
Predicción de Y0
library(kableExtra)
p=exp(predict(mod4,list(diametro=7),interval = "prediction",level = 0.95))
G=matrix(c(p[1], p[2], p[3]), ncol=3)
colnames(G)<- c("Media","Inferior","Superior")
rownames(G)<- c("Peso estimado")
G %>%
kbl(booktabs = T,) %>%
kable_classic_2(full_width = F)| Media | Inferior | Superior | |
|---|---|---|---|
| Peso estimado | 13.94957 | 16.48349 | 14.34307 |