1. Descripción del problema y caraterización de los datos.

Uno de los propósitos a nivel mundial es el de estimar la captura de dióxido de carbono, uno principales gases causantes del calentamiento global y por tanto del cambio climático.

Un grupo de investigadores están interesados en poder construir modelos que permitan la valoración de estos beneficios, a través de información recolectada sobre características de los árboles de una región húmeda, que les permita una estimación de la biomasa y así facilitar la toma de decisiones y la generación de políticas públicas.

Como no es posible obtener el valor del peso del árbol sin cortar este mismo, se plantea la opción de estimarla a través de variables como la altura y el diámetro del tronco, información requerida para la estimación del valor de la biomasa.

Para ayudar en las estimaciones al grupo de cientificos, se utiliza una base de datos denominada arboles que contiene las variables: id, altura, diametro y peso. Para conocer los datos y sacar conclusiones iniciales generamos un resumen estadístico de la base de datos de árboles 

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.00

El problema que se intenta resolver es usar la información de la base de datos de árboles con base a las variables altura, diametro para predecir el peso sin usar tecnicas como la tala de los arboles

Para dar solución, se usará un modelo de regresión lineal simple, el cuál se basa en la relación entre dos variables, una variable dependiente y una variable independiente, en este caso la variable dependiente es el peso y la variable independiente es el diametro del arbol. Se usa las variable diametro como dato de entrada y variable independiente para calcular el peso, esto es porque la obtención del diametro es menos invasivo y costoso que calcular la altura del arbol.

En el siguiente grupo de gráficas se observa la distribución de cada conjunto de datos que vamos a utilizar, el diametro y peso.

En la siguiente gráfica, observamos la Aproximación lineal de los puntos de las variables diametro vs peso, y se concluye que la distribución de punto, a simple vista, no cumple un comportamiento lineal. De igual forma, más delante a través de supuestos, tendremos pruebas concluyente para reconocer el comportamiento de los puntos. 

library(ggplot2)
ggplot(arboles, aes(x = diametro, y = peso)) + geom_point() + geom_smooth(method = "lm")

2. Normalidad de de la variables

2.1. Normalidad de de la variable diametro.

H0: Los datos del diámetro tienen una distribución normal.

shapiro.test(arboles$diametro)
## 
##  Shapiro-Wilk normality test
## 
## data:  arboles$diametro
## W = 0.99109, p-value = 0.8079

Los datos de la variable diametro tiene un comportamiento de distribución normal debido a que p-value > 5% de significancia. Por lo tanto, no se rechaza la hipotesis H0.

2.2. Normalidad de de la variable peso.

H0: Los datos del Peso tienen una distribución normal.

H1: Los datos del Peso no tienen una distribución normal.

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

Los datos de la variable peso no tienen un comportamiento de distribución normal debido a que p-value < 5% de significancia. Por lo tanto, se rechaza la hipotesis H0 y acepta la hipotesis H1.

2.3. Inferencia sobre la correlación entre la variable diametro y peso. 

cor.test(arboles$diametro, arboles$peso)
## 
##  Pearson's product-moment correlation
## 
## data:  arboles$diametro and arboles$peso
## 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

La variables diametro y peso se correlacionan directamente.

3. Estimación del modelo por MCO. 

modelo = lm(peso~diametro, arboles)
summary(modelo)
## 
## 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

En el cálculo de los valores de Beta0 y Beta1 para la función lineal se obtiene para Beta0 = -9.0203 y Beta1 = 5.1026 para un función lineal y = -9.0203 + 5.1026x

3.1. Validación de los supuestos sobre los errores. 

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

Obtención de los errores: 

e = modelo$residuals

3.2. Supuesto de normalidad de los errores.

H0: los errores tienen distribución normal.

shapiro.test(e)
## 
##  Shapiro-Wilk normality test
## 
## data:  e
## W = 0.95356, p-value = 0.002793

Los residuales NO cumple los criterios de distribution normal debido a que p-value < 5% de significancia. Por lo tanto, se rechaza la hipótesis H0.

3.3. Supuesto de independencia de errores.

H0: Los errores sin independientes (no están relacionados).

dwtest(modelo)
## 
##  Durbin-Watson test
## 
## data:  modelo
## DW = 1.0719, p-value = 1.035e-06
## alternative hypothesis: true autocorrelation is greater than 0

Los errores son dependientes debido a que p-value < 5% de significancia. Por lo tanto, se rechaza la hipótesis H0 y acepta la hipótesis alternativa.De igual forma, este supuesto no aplica para este conjunto de datos debido a que no existe una relación temporal(tiempo) en la recolección de los mismos.

3.4. Supuesto de varianza constante.

H0: la varianza de los errores es constante.

gqtest(modelo)
## 
##  Goldfeld-Quandt test
## 
## data:  modelo
## GQ = 2.0131, df1 = 43, df2 = 43, p-value = 0.01196
## alternative hypothesis: variance increases from segment 1 to 2

Los errores NO cumple con el supuesto de varianza constante debido a que p-value < 5% de significancia. Por lo tanto, se rechaza la hipótesis H0.

4. Tranformaciones.

4.1. Evaluación de linealidad.

Para realizar la transformación de la regresión lineal aplicamos la técnica de Boxcox para hallar el valor de λ, dónde de acuerdo al valor λ se establece una transformación que permita ajustar de mejor forma la regresión. 

par(mfrow = c(1,2))
boxcox(lm(arboles$peso ~ arboles$diametro), lambda = -2:2)
bc<-boxcox(lm(arboles$peso ~ arboles$diametro), lambda = 0:1)

Hallamos el valor de lambda para establecer la tranformación. 

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

El valor calculado de λ = 0, por lo cual, la transformación aplicar sobre la regresión es log(peso) dentro de la fórmula de la regresión para obtener log(peso)~diametro. 

data(arboles)
modelo_transformado=lm(log(peso) ~ diametro, arboles)
summary(modelo_transformado)
## 
## Call:
## lm(formula = log(peso) ~ diametro, data = arboles)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27395 -0.10180 -0.00328  0.10073  0.33742 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.32798    0.05977   22.22   <2e-16 ***
## diametro     0.27818    0.01061   26.22   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1453 on 88 degrees of freedom
## Multiple R-squared:  0.8865, Adjusted R-squared:  0.8852 
## F-statistic: 687.6 on 1 and 88 DF,  p-value: < 2.2e-16

El valor calculado de λ = 0, por lo cuál, se transforma la variable peso dentro de la fórmula de la regresión para obtener log(peso)~diametro, por lo cual, tendremos un función lineal ajustada con un nuevo Beta0 = 1.32798 y Beta1 = 0.27818.

En el siguiente grupo de gráficas se observa la distribución de cada uno del conjunto de datos que vamos a utilizar, el diametro y peso después de realizar la transformación y el ajuste de datos.

Graficamos la regresión con los datos ajustados. 

library(ggplot2)
ggplot(arboles, aes(x = diametro, y = log(peso))) + geom_point() + geom_smooth(method = "lm")

En la grafica se dibujan los datos en la escala logarítmica transformada log(y) según con el λ calculada por la función Boxcox y al agregar la línea ajustada, la relación parece lineal, y ya podemos ver que la variación sobre la línea ajustada parece constante. 

library(rlang)
# f_ajustado_inv = exp(ajust_b0 + ajust_b1 * x) = exp(ajust_y)
f_ajustada_inv <- function(x) exp(1.32798+0.27818*x)
ggplot(arboles, aes(x = diametro, y = peso)) + geom_point(colour = "black") + 
geom_function(fun = f_ajustada_inv, colour = "blue", size = 1) + geom_smooth(size=0)

Al graficar los datos en la escala original aplicamos exp(ajust_y) (función inversa) a la regresión ajustada observando una relación exponencial. Sin embargo, este sigue siendo un modelo lineal, ya que la nueva respuesta transformada, log(y), sigue siendo una combinación lineal de los predictores.

4.2. Validación de los supuestos sobre los errores.

Graficamos los supuestos sobre lo errores. 

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

Obtención de los errores: 

e = modelo_transformado$residuals

4.3. Supuesto de normalidad de los errores

H0: los errores tienen distribución normal.

shapiro.test(e)
## 
##  Shapiro-Wilk normality test
## 
## data:  e
## W = 0.98394, p-value = 0.3338

Los residuales SI cumplen los criterios de distribution normal debido a que p-value > 5% de significancia. Por lo tanto, se acepta la hipótesis H0.

4.4. Supuesto de independencia de errores

H0: Los errores son independientes (no están relacionados).

library(lmtest)
dwtest(modelo_transformado)
## 
##  Durbin-Watson test
## 
## data:  modelo_transformado
## DW = 0.67803, p-value = 1.716e-13
## alternative hypothesis: true autocorrelation is greater than 0

Los errores son dependientes debido a que p-value < 5% de significancia. Por lo tanto, se rechaza la hipótesis H0 y se acepta la hipótesis alternativa.De igual forma, este supuesto no aplica para este conjunto de datos debido a que no existe una relación temporal(tiempo) en la recolección de los mismos.

4.5. Supuesto de varianza constante.

H0: la varianza de los errores es constante.

gqtest(modelo_transformado)
## 
##  Goldfeld-Quandt test
## 
## data:  modelo_transformado
## GQ = 1.1538, df1 = 43, df2 = 43, p-value = 0.3206
## alternative hypothesis: variance increases from segment 1 to 2

Los errores SI cumple con el supuesto de varianza constante debido que p-value > 5% de significancia. Por lo tanto, se acepta la hipótesis H0.

5. Predicción de Y

Para la predicción Y tenemos la funnción lineal ajustada ajust_y = 1.32798+0.27818x. Para la prueba escogemos el valor de x = 6.0 y obtenemos los siguiente datos. 

datos = filter(arboles, arboles$diametro == 6.0)
datos
## # A tibble: 5 × 4
##      id  peso diametro altura
##   <int> <dbl>    <dbl>  <dbl>
## 1    46  19.0        6    6  
## 2    50  20.1        6    6.6
## 3    62  24.5        6    7.2
## 4    66  22.8        6    7.3
## 5    67  23.1        6    7.5

Realizamos la predicción Y tenemos la función lineal ajustada ajust_y = 1.32798+0.27818x. Para la prueba escogemos el valor de x(diametro) = 6.0 y obtenemos: 

predict(modelo_transformado, data.frame(diametro=6.0), 
                     interval = "predict", level = 0.95)
##        fit      lwr      upr
## 1 2.997089 2.706495 3.287682

Los valores anteriores es con ajust_y, para llevarlos a escala original aplicamos exp(ajust_y) su inversa. 

exp(predict(modelo_transformado, data.frame(diametro=6.0), 
                     interval = "predict", level = 0.95))
##        fit      lwr      upr
## 1 20.02714 14.97669 26.78072

6. CONCLUSIÓN GENERAL.

Los análisis de regresión y de correlación lineal se constituyen como un métodos que nos permite conocer las relaciones y significación entre variables dentro conjunto de datos. Lo anterior, es de suma importancia porque se presentan variables de respuesta e independencias que interactúan para caracterizar un proceso en particular y por ende; analizar, predecir valores de la variable dependiente y examinar el grado de relacionan de dichas variables.

Las técnicas de análisis de la estadística inferencial permiten identificar para un conjunto de datos los modelos que más se ajusten a la caracterización de los mismos, es decir, determinar con cuál función de distribución de probabilidad se relacionan los datos para poder ajustar y usar adecuadamente modelos como en el caso de estudio de regresión líneal.

Identificar adecuadamente las variables independientes de las dependientes es el fundamento del ejercicio, pues este le permite al modelo de regresión lineal tener una mejor entrada para su formulación. Es importante, el análisis de la aleatoriedad en los datos, el Supuesto de varianza constante, Supuesto de independencia de errores, con el fin de verificar las condiciones del modelo inicial, y contrarrestar los valores p con el modelo una vez se realiza la transformación según el lambda.