Importación de librerias

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)

Importación de data de arboles

options(repos = list(CRAN="http://cran.rstudio.com/"))
data= read.csv("C:/Users/franc/OneDrive/Escritorio/Maestria_Ciencia_Datos/Metodos_Estadisticos/Actividad Arboles/arboles.txt", sep=";")

dataRda = load("C:/Users/franc/OneDrive/Escritorio/Maestria_Ciencia_Datos/Metodos_Estadisticos/Actividad Arboles/arboles.rda")

Exploración de datos

glimpse(data)

## Rows: 90
## Columns: 4
## $ id       <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18…
## $ peso     <dbl> 13.73, 14.58, 15.88, 8.99, 6.99, 19.34, 21.44, 13.81, 11.88, …
## $ diametro <dbl> 4.7, 5.3, 4.8, 3.2, 2.2, 6.3, 6.6, 5.3, 4.9, 5.9, 2.5, 4.7, 2…
## $ altura   <dbl> 5.0, 5.6, 5.8, 4.3, 3.3, 7.9, 8.3, 7.3, 6.7, 7.1, 3.7, 5.7, 3…

str(data)

## 'data.frame':    90 obs. of  4 variables:
##  $ id      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ peso    : num  13.73 14.58 15.88 8.99 6.99 ...
##  $ diametro: num  4.7 5.3 4.8 3.2 2.2 6.3 6.6 5.3 4.9 5.9 ...
##  $ altura  : num  5 5.6 5.8 4.3 3.3 7.9 8.3 7.3 6.7 7.1 ...

Se puede ver que el corpus de datos cuenta con 4 variables (id,peso,diametro y altura), y con 90 observaciones. En este caso la variable de interes (y) es peso, y las variables predictoras (xi) serán diametro y altura. Adicionalmente, se tiene que las 3 varaibes son de naturaleza númerica.

1. Análisis de Modelos de Regresión Simple

A continuación se realizará el análisis de la variable de interes en función de cada una de las variables predictoras, verificando los supuestos que se deben cumplir para asegurar un correcto desempeño del modelo de regresión lineal

pairs(data)

Se puede ver en la matriz de graficos de dispersión que existe una posible relación positiva entre el peso y el diametro, y entre el peso y la altura. Dichas relaciones se analizarán a continuación, primero en dos análisis bi-variados, y posteriormente, en un análisis multivariado

1.1 análisis de modelo 1: Peso vs Diametro

A continuación se realizará el análisis del modelo de regresión lineal 1, el cual cruzará a la variable objetivo “Peso” contra la variable predictora “Diametro”.

1.1.1 Desarrollo del modelo 1:

# Modelo LM 1
mod1 = lm(data$peso ~ data$diametro)
summary(mod1)

## 
## Call:
## lm(formula = data$peso ~ data$diametro)
## 
## 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 ***
## data$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

plot(data$diametro,data$peso, xlab="Diametro", ylab="Peso")
abline(mod1)

Por otro lado, el modelo arroja que la recta que minimiza la sumatoria de los cuadrados de los errores es y = -9.02 + 5.1x , con un coeficiente de determinación del 82,4%, lo cual indica la bondad de ajuste de la recta de regresión lineal al comportamiento de los datos.

Adicionalmente, se realizará el análisis de correlación entre las variables del modelo 1, con el objetivo de determinar la fuerza de la relación entre el Peso y el Diametro mediante el coeficiente del Coeficiente de Producto-Momento de Pearson.

cor.test(data$peso, data$diametro)

## 
##  Pearson's product-moment correlation
## 
## data:  data$peso and data$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

Se puede ver que el coeficiente de pearson para Peso vs Diametro es del 90%, lo cual indica que las dos variables estan altamente correlacionadas de manera positiva.

1.1.2 Validación de supuestos

Para asegurarse que existe una relación lineal entre el peso y el diametro, se deben validar los siguientes supuestos:

La relación entre la variable Y y X es lineal

ggplot(data, aes(x=diametro, y=peso)) + geom_point() + geom_smooth()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

en primer lugar, se puede ver que el modelo de regresión lineal entre el peso y el diametro muestra graficamente una relación lineal entre las dos variable, por ende, se valida el primer supuesto de aplicabilidad del modelo.

La variable predictora X no es una variable aleatoria, y la varianza en la variable objetivo Y se debe a los errores existentes entre las predicciones del modelo y los valores reales de Y.

## Ploteo de los valores predichos por el modelo RL vs los errores de cada preducción
plot(mod1$fitted.values,mod1$residuals)

##Creación de dataframe con los valores predichos por el modelo y sus respectivos errores

res1 = data.frame(mod1$fitted.values, mod1$residuals )
ggplot(res1, aes(x=mod1$fitted.values, y=mod1$residuals))+geom_point() + geom_smooth()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Se puede ver que al graficarse los errores del modelo 1 los errores tienden a distribuirse de manera homogenea al rededor del 0, lo cual podría indicar que la varianza del modelo se debe al factor aleatorio de los errores y que estos no estan incorporando información adicional que podrian determinar el comportamiento de la variable objetivo Y.

Los errores se distribuyen como una variable aleatoria con distribución normal y media = 0, de tal modo, que los errores no esten incluyendo información adicional que ayude a describir la variablidad de la variable objetivo.

# Media de los errores
media_error = mean(res1$mod1.residuals)
desv1 = sd(res1$mod1.residuals)
media_error

## [1] 2.553898e-17

desv1

## [1] 3.415494

Se puede ver que la media de los errores del modelo 1 (Peso vs Diametro) es 0,00000000000000000255 y una desviación estandar de 3,4, lo cual cumple con el terce de los criterios de validación.

#Grafica de histograma de los errores del modelo 1
ggplot(res1, aes(x=mod1$residuals)) + geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

qqnorm(mod1$residuals)
qqline(mod1$residuals, col = "blue")

Al analizar el comportamiento del histograma, se puede ver que los errores siguen aparentemente una distribuación normal, adicionalmente, el qqplot muestra que los errores tienden a coplarse a la linea de normalidad, con un leve comportamiento de alejarse en las colas del grafico, lo cual puede significar que hay presencia de outliers, pero que en general si se comportan como una variable aleatoria con distribución normal.

Adicionalmente, se realizará un test de Shapiro Wilks para corroborar la normalidad

shapiro.test(mod1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  mod1$residuals
## W = 0.95356, p-value = 0.002793

Se puede ver en el resultado de la prueba de Shapiro-Wilks que el P-Value es de 0.002, el cual es inferior al nivel de significancia de la prueba 5%, por ende no se rechaza la hipotesis nula de la prueba, la cual es que los datos de los residuos del modelo 1 se distribuye normal.

Supuesto de independencia de errores

Para validar la independencia de los errores se utilizará la prueba de Durbin-Watson, la cual tiene por hipotesis nula que los datos no estan autocorrelacionados.

#Prueba Durbin-Watson
library(lmtest)

## Warning: package 'lmtest' was built under R version 4.2.2

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

dwtest(mod1)

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

El resultado de la prueba Durbin-Watson arroja que los errores del modelo 1 no estan autocorrelacionados, ya que el P-Value de la prueba es de 0.00000135, el cual es inferior a la significancia de la prueba del 5%.

Supuesto de varianza constante en los errores

Para validar el supuesto varianza constante en los errores se utilizará la prueba de Goldfeld-Quandt, la cual tiene por hipotesis nula que hay homocedasticidad en los errores, lo cual quiere decir que la varianza en los datos es homogenea de manera trasversal, lo cual dice que no existe una dispersión heterogenea de los residuos (heterocedasticidad)

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 2

Se puede ver que el P-Value de la prueba es de 0.011, el cual es inferior al nivel de significancia de la prueba del 5%, por ende, se puede decir que hay homocedasticidad en los errores del modelo 1, lo cual verifica el supuesto correspondiente

1.2 análisis de modelo 1: Peso vs Altura

A continuación se realizará el análisis del modelo de regresión lineal 1, el cual cruzará a la variable objetivo “Peso” contra la variable predictora “Diametro”.

1.2.1 Desarrollo del modelo 2:

# Modelo LM 2
mod2 = lm(data$peso ~ data$altura)
summary(mod2)

## 
## Call:
## lm(formula = data$peso ~ data$altura)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.228 -1.969  0.572  2.377 15.106 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.0456     1.7046  -4.133 8.14e-05 ***
## data$altura   3.8906     0.2481  15.684  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.211 on 88 degrees of freedom
## Multiple R-squared:  0.7365, Adjusted R-squared:  0.7335 
## F-statistic:   246 on 1 and 88 DF,  p-value: < 2.2e-16

plot(data$altura,data$peso, xlab="altura", ylab="Peso")
abline(mod2)

Por otro lado, el modelo arroja que la recta que minimiza la sumatoria de los cuadrados de los errores es y = -7.04 + 3.8x , con un coeficiente de determinación del 73.6%, lo cual indica la bondad de ajuste de la recta de regresión lineal al comportamiento de los datos.

1.2.2 Validación de supuestos

Para asegurarse que existe una relación lineal entre el peso y el diametro, se deben validar los siguientes supuestos:

La relación entre la variable Y y X es lineal

ggplot(data, aes(x=altura, y=peso)) + geom_point() + geom_smooth()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

en primer lugar, se puede ver que el modelo de regresión lineal entre el peso y la altura muestra graficamente una relación lineal entre las dos variable, por ende, se valida el primer supuesto de aplicabilidad del modelo.

La variable predictora X no es una variable aleatoria, y la varianza en la variable objetivo Y se debe a los errores existentes entre las predicciones del modelo y los valores reales de Y.

## Ploteo de los valores predichos por el modelo RL vs los errores de cada preducción
plot(mod2$fitted.values,mod2$residuals)

##Creación de dataframe con los valores predichos por el modelo y sus respectivos errores

res2 = data.frame(mod2$fitted.values, mod2$residuals )
ggplot(res2, aes(x=mod2$fitted.values, y=mod2$residuals))+geom_point() + geom_smooth()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Se puede ver que al graficarse los errores del modelo 2 los errores tambien tienden a distribuirse de manera homogenea al rededor del 0, lo cual podría indicar que la varianza del modelo se debe al factor aleatorio de los errores y que estos no estan incorporando información adicional que podrian determinar el comportamiento de la variable objetivo Y.

Los errores se distribuyen como una variable aleatoria con distribución normal y media = 0, de tal modo, que los errores no esten incluyendo información adicional que ayude a describir la variablidad de la variable objetivo.

# Media de los errores
media_error2 = mean(res2$mod2.residuals)
desv2 = sd(res2$mod2.residuals)
media_error2

## [1] -1.28331e-16

desv2

## [1] 4.18726

Se puede ver que la media de los errores del modelo 2 (Peso vs altura) es 0,0000000000000000128 y una desviación estandar de 4.1, lo cual cumple con el terce de los criterios de validación.

#Grafica de histograma de los errores del modelo 1
ggplot(res2, aes(x=mod2$residuals)) + geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

qqnorm(mod2$residuals)
qqline(mod2$residuals, col = "red")

Al analizar el comportamiento del histograma de los errores del modelo 2, se puede ver que los errores siguen aparentemente una distribuación normal, adicionalmente, el qqplot muestra que los errores tienden a coplarse a la linea de normalidad, con un leve comportamiento de alejarse en las colas del grafico, lo cual puede significar que hay presencia de outliers, pero que en general si se comportan como una variable aleatoria con distribución normal.

Adicionalmente, se realizará un test de Shapiro Wilks para corroborar la normalidad

shapiro.test(mod2$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  mod2$residuals
## W = 0.95439, p-value = 0.003157

Se puede ver en el resultado de la prueba de Shapiro-Wilks para el modelo 2 que el P-Value es de 0.003, el cual es inferior al nivel de significancia de la prueba 5%, por ende no se rechaza la hipotesis nula de la prueba, la cual es que los datos de los residuos del modelo 2 se distribuye normal.

Supuesto de independencia de errores

Para validar la independencia de los errores se utilizará la prueba de Durbin-Watson, la cual tiene por hipotesis nula que los datos no estan autocorrelacionados.

#Prueba Durbin-Watson
library(lmtest)
dwtest(mod2)

## 
##  Durbin-Watson test
## 
## data:  mod2
## DW = 0.90216, p-value = 4.234e-09
## alternative hypothesis: true autocorrelation is greater than 0

El resultado de la prueba Durbin-Watson arroja que los errores del modelo 2 no estan autocorrelacionados, ya que el P-Value de la prueba es de 0.0000000042, el cual es inferior a la significancia de la prueba del 5%.

Supuesto de varianza constante en los errores

gqtest(mod2)

## 
##  Goldfeld-Quandt test
## 
## data:  mod2
## GQ = 1.1017, df1 = 43, df2 = 43, p-value = 0.3761
## alternative hypothesis: variance increases from segment 1 to 2

Se puede ver que el P-Value de la prueba es de 0.2, el cual es superior al nivel de significancia de la prueba del 5%, por ende, se puede decir que si hay heterocedasticidad en los errores del modelo 2, lo cual rechaza el supuesto correspondiente

2. Selección del modelo de regresión lineal

Teniendo en cuenta los análisis previos, se puede ver que:

El R2 de modelo 1 (Peso vs Diametro) es del 82.4%, mientras que el R2 del modelo 2 (Peso vs Altura) es de 73.6%, lo cual indica que el modelo 1 explica en un mayor grado el comportamiento de la variable objetivo (Peso).
Los modelos 1 y 2, cumplen con los supuestos de 1) linealidad, 2) errores del modelo como variables aleatorios con distribución normal (mean = 0, st= sigma) y 3) independencia de los errores.
Solo el modelo 1 cumple con el supuesto de homocedasticidad de la varianza en los errores, mientras que el modelo 2 presenta heterocedasticidad de la varianza.

Por ende, se selecciona al modelo 1 como el modelo adecuado para determinar la biomasa de los arboles, el cual es:

y = -9.02 + 5.1x

3. Conclusiones y recomendaciones

Con referencia al modelo propuesto, el diametro de los arboles es la variable que mejor ayuda a pronosticar el peso (biomasa) de los arboles, teniendo en cuenta que el modelo de regresión lineal de estas dos variables explica el 82% de los datos, lo cual muestra que el modelo es robusto.
Se puede determinar, según la validación de los supuesto correspondientes, que existe una relación lienal entre el diametro de los arboles y la cantidad de biomasa de estos.
Para la selección de un modelo predictor no solo es util la comparación del R2 (coeficiente de determinación), si no que se deben validar los supuestos del modelo, teniendo en cuenta que la validación de estos podrá asegurar que el modelo lineal refleje de manera consistente una herramienta de predicción para las variables incluidas.

Actividad 1

Francisco Pizarro

2023-02-25