Taller1
Pontificia Universidad Javeriana - Cali
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(peso~diametro, data = arboles)
stargazer(mod1, type = 'text')=============================================== 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 |
|---|
| 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(peso~diametro+ altura, data = arboles)# variable de respuesta en función de la altura y el diametro
stargazer(mod3, type = 'text')=============================================== Dependent
variable:
————————— peso
———————————————– diametro 4.739***
(0.713)
altura 0.313
(0.575)
Constant -9.121***
(1.430)
| Observations 90 R2 0.825 Adjusted R2 0.821 Residual Std. Error 3.449 (df = 87) F Statistic 205.474*** (df = 2; 87) =============================================== Note: p<0.1; p<0.05; p<0.01 |
|---|
| 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:
require(stargazer)
mod5 <- lm(peso~log(diametro), data = arboles)
mod6 <- lm(log(peso)~log(diametro), data = arboles)
stargazer(mod1, mod4, mod5, mod6, type="text", title="Comparación de modelos", df=FALSE)Comparación de modelos
Dependent variable:
-------------------------------------------
peso log(peso) peso log(peso)
(1) (2) (3) (4) | diametro 5.103*** 0.278*** (0.251) (0.011) |
| log(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, newdata = data.frame(diametro=9),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 | 46.13825 | 34.18439 | 62.27224 |
p1<-predict(mod1, newdata = data.frame(diametro=9),interval = "prediction",level = 0.95)
G=matrix(c(p1[1], p1[2], p1[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 | 36.90294 | 29.81417 | 43.99171 |