Introducción

En esta actividad, aprenderemos a ajustar modelos de regresión lineal utilizando un conjunto de datos que contiene información sobre la altura total (HTO) y el diámetro a la altura del pecho (DAP) de árboles, entre otras variables. Usaremos R para explorar los datos, ajustar modelos y evaluar su desempeño.

Configuración Inicial

Primero, carguemos las bibliotecas necesarias para el análisis:

library(readr)
library(ggplot2)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ stringr   1.5.1
## ✔ forcats   1.0.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Carga y Exploración de Datos

Carguemos el conjunto de datos desde el archivo proporcionado:

Datos_base <- read_table("C:/Users/gvega/Downloads/DBase_RegresionHDap (1).txt")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   TRAT = col_character(),
##   ARB = col_double(),
##   MNU = col_double(),
##   DAP = col_double(),
##   HTO = col_double(),
##   Edad = col_double()
## )

Exploremos la estructura de los datos:

names(Datos_base)
## [1] "TRAT" "ARB"  "MNU"  "DAP"  "HTO"  "Edad"
summary(Datos_base)
##      TRAT                ARB             MNU             DAP       
##  Length:4417        Min.   : 1.00   Min.   :1.000   Min.   : 5.30  
##  Class :character   1st Qu.: 8.00   1st Qu.:2.000   1st Qu.:12.10  
##  Mode  :character   Median :15.00   Median :4.000   Median :14.10  
##                     Mean   :15.32   Mean   :3.502   Mean   :14.19  
##                     3rd Qu.:23.00   3rd Qu.:5.000   3rd Qu.:16.20  
##                     Max.   :32.00   Max.   :6.000   Max.   :25.30  
##       HTO             Edad       
##  Min.   : 8.20   Min.   : 5.090  
##  1st Qu.:15.10   1st Qu.: 6.130  
##  Median :17.50   Median : 8.270  
##  Mean   :17.49   Mean   : 7.727  
##  3rd Qu.:19.70   3rd Qu.: 9.390  
##  Max.   :68.00   Max.   :10.330

Análisis Visual

Grafiquemos la relación entre HTO y DAP para inspeccionar los datos:

GRAFICO <- ggplot(Datos_base, aes(x = DAP, y = HTO)) +
  geom_point()
GRAFICO  # Podemos observar valores atípicos (outliers)

Filtrado de Datos

Notamos valores atípicos en la variable HTO. Filtrémoslos para mejorar el ajuste de los modelos:

Datos <- Datos_base %>%
  filter(HTO < 35)

GRAFICO <- ggplot(Datos, aes(x = DAP, y = HTO)) +
  geom_point()
GRAFICO

Ajuste de Modelos

A continuación, ajustaremos diferentes modelos de regresión lineal y analizaremos sus resultados.

Modelo 1: HTO ~ B0 + B1*DAP

Ajustemos un modelo lineal simple con intercepto:

Modelo1 <- lm(HTO ~ DAP, Datos)
summary(Modelo1)
## 
## Call:
## lm(formula = HTO ~ DAP, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.8433 -1.2939 -0.0445  1.2642  6.2964 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.449758   0.134077   33.19   <2e-16 ***
## DAP         0.915634   0.009248   99.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.832 on 4410 degrees of freedom
## Multiple R-squared:  0.6897, Adjusted R-squared:  0.6896 
## F-statistic:  9802 on 1 and 4410 DF,  p-value: < 2.2e-16

Modelo 2: HTO ~ B1*DAP (sin intercepto)

Probemos un modelo sin intercepto:

Modelo2 <- lm(HTO ~ 0 + DAP, Datos)
summary(Modelo2)
## 
## Call:
## lm(formula = HTO ~ 0 + DAP, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.6889 -1.2411  0.1771  1.4598  7.6936 
## 
## Coefficients:
##     Estimate Std. Error t value Pr(>|t|)    
## DAP 1.216004   0.002126   571.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.047 on 4411 degrees of freedom
## Multiple R-squared:  0.9867, Adjusted R-squared:  0.9867 
## F-statistic: 3.271e+05 on 1 and 4411 DF,  p-value: < 2.2e-16

Modelo 3: HTO ~ DAP + Edad

Incluyamos la variable Edad como predictor adicional:

Modelo3 <- lm(HTO ~ DAP + Edad, Datos)
summary(Modelo3)
## 
## Call:
## lm(formula = HTO ~ DAP + Edad, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7570 -0.9230 -0.1261  0.9117  5.3560 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.811321   0.109214   16.59   <2e-16 ***
## DAP         0.675964   0.007975   84.76   <2e-16 ***
## Edad        0.781817   0.013128   59.55   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.364 on 4409 degrees of freedom
## Multiple R-squared:  0.828,  Adjusted R-squared:  0.828 
## F-statistic: 1.061e+04 on 2 and 4409 DF,  p-value: < 2.2e-16

Transformaciones de Variables

Algunos modelos requieren transformaciones de las variables. Veamos dos ejemplos.

Modelo 4: LN(HTO) ~ B0 + LN(DAP)

Transformemos las variables a su logaritmo natural:

Datos$LNDAP <- log(Datos$DAP)
Datos$LNHTO <- log(Datos$HTO)
Modelo4 <- lm(LNHTO ~ LNDAP, Datos)
summary(Modelo4)
## 
## Call:
## lm(formula = LNHTO ~ LNDAP, data = Datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.59089 -0.07101  0.00331  0.07530  0.42469 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.903904   0.020275   44.58   <2e-16 ***
## LNDAP       0.736437   0.007685   95.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1111 on 4410 degrees of freedom
## Multiple R-squared:  0.6756, Adjusted R-squared:  0.6755 
## F-statistic:  9184 on 1 and 4410 DF,  p-value: < 2.2e-16

Modelo 5: HTO - 1.3 ~ B1*DAP

Restemos 1.3 a HTO para ajustar este modelo:

Datos$HTO1.3 <- Datos$HTO - 1.3
Modelo5 <- lm(HTO1.3 ~ DAP, Datos)
summary(Modelo5)
## 
## Call:
## lm(formula = HTO1.3 ~ DAP, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.8433 -1.2939 -0.0445  1.2642  6.2964 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.149758   0.134077   23.49   <2e-16 ***
## DAP         0.915634   0.009248   99.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.832 on 4410 degrees of freedom
## Multiple R-squared:  0.6897, Adjusted R-squared:  0.6896 
## F-statistic:  9802 on 1 and 4410 DF,  p-value: < 2.2e-16

Visualización de Ajustes

Visualicemos el ajuste del Modelo 1 sobre los datos:

M1_B0 <- Modelo1$coefficients[1]
M1_B1 <- Modelo1$coefficients[2]

GRAFICO_modelo1 <- ggplot(Datos, aes(x = DAP, y = HTO)) +
  geom_point() +
  geom_line(aes(x = DAP, y = M1_B0 + M1_B1 * DAP), color = "red")
GRAFICO_modelo1

Criterios de Evaluación

Evaluemos el Modelo 1 con el índice de Furnival y el criterio de información de Akaike (AIC).

Índice de Furnival

M1_Furnival <- sigma(Modelo1)
M1_Furnival
## [1] 1.831692

Criterio de Información de Akaike (AIC)

Calculemos el AIC manualmente:

error_m1 <- Modelo1$residuals
cm_m1 <- sum(error_m1^2)
L1 <- length(error_m1)
rmse_m1 <- sqrt(cm_m1 / (L1 - 2))  # Ajuste por 2 parámetros
AIC_m1 <- L1 * log(cm_m1 / L1) + 2 * 2  # 2 parámetros: B0 y B1
AIC_m1
## [1] 5342.636

Ejercicio

Como práctica, calcula el índice de Akaike (AIC) para los Modelos 2, 3, 4 y 5 siguiendo el mismo procedimiento que usamos para el Modelo 1.