Regresión lineal multiple

#install.packages("ggplot2")
#install.packages("dplyr")
#install.packages("tidyr")  
#install.packages("corrplot")
#install.packages("car")
#install.packages("MASS")
#install.packages("lmtest")  # Instala el paquete
library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

library(MASS)
library(car)

## Loading required package: carData

library(corrplot)

## corrplot 0.95 loaded

library(tidyr)
library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following object is masked from 'package:MASS':
## 
##     select

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

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

library(dplyr)

Contexto del ejercicio

Se cree que la energía eléctrica consumida cada mes por una planta química está relacionada con la temperatura ambiental promedio (x1), el número de días en el mes (x2), la pureza promedio del producto (x3) y las toneladas de producto producido (x4). Los datos históricos del año pasado se encuentran a continuación.

12-29. Considera los datos de consumo de energía eléctrica.

Encuentra intervalos de confianza del 95% para 1, 2, 3 y 4.
Encuentra un intervalo de confianza del 95% para la media de Y cuando x₁ = 75, x₂ = 24, x₃ = 90, y x₄ = 98.
Encuentra un intervalo de predicción del 95% para el consumo de energía cuando x₁ = 75, x₂ = 24, x₃ = 90, y x₄ = ## Cargar los datos

data <- data.frame(
  y = c(240, 236, 270, 274, 301, 316, 300, 296, 267, 276, 288, 261),
  x1 = c(25, 31, 45, 60, 65, 72, 80, 84, 75, 60, 50, 38),
  x2 = c(24, 21, 24, 25, 25, 26, 25, 25, 24, 25, 25, 23),
  x3 = c(91, 90, 88, 87, 91, 94, 87, 86, 88, 91, 90, 89),
  x4 = c(100, 95, 110, 88, 94, 99, 97, 96, 110, 105, 100, 98)
)

Descriptivo de los datos

summary(data)

##        y               x1              x2              x3       
##  Min.   :236.0   Min.   :25.00   Min.   :21.00   Min.   :86.00  
##  1st Qu.:265.5   1st Qu.:43.25   1st Qu.:24.00   1st Qu.:87.75  
##  Median :275.0   Median :60.00   Median :25.00   Median :89.50  
##  Mean   :277.1   Mean   :57.08   Mean   :24.33   Mean   :89.33  
##  3rd Qu.:297.0   3rd Qu.:72.75   3rd Qu.:25.00   3rd Qu.:91.00  
##  Max.   :316.0   Max.   :84.00   Max.   :26.00   Max.   :94.00  
##        x4        
##  Min.   : 88.00  
##  1st Qu.: 95.75  
##  Median : 98.50  
##  Mean   : 99.33  
##  3rd Qu.:101.25  
##  Max.   :110.00

str(data)

## 'data.frame':    12 obs. of  5 variables:
##  $ y : num  240 236 270 274 301 316 300 296 267 276 ...
##  $ x1: num  25 31 45 60 65 72 80 84 75 60 ...
##  $ x2: num  24 21 24 25 25 26 25 25 24 25 ...
##  $ x3: num  91 90 88 87 91 94 87 86 88 91 ...
##  $ x4: num  100 95 110 88 94 99 97 96 110 105 ...

# Calcular y mostrar la matriz de correlación
cor_matrix <- cor(data[, -1])
print(cor_matrix)

##             x1          x2          x3          x4
## x1  1.00000000  0.66045595 -0.28756637 -0.02355867
## x2  0.66045595  1.00000000  0.11273901 -0.02532776
## x3 -0.28756637  0.11273901  1.00000000  0.07891362
## x4 -0.02355867 -0.02532776  0.07891362  1.00000000

# Visualizar la matriz de correlación
corrplot(cor_matrix, method = "circle", type = "upper", tl.cex = 0.8)

length(data)

## [1] 5

# Tamaño de cada columna
column_sizes <- sapply(data, length)

# Mostrar el tamaño de cada columna
column_sizes

##  y x1 x2 x3 x4 
## 12 12 12 12 12

pairs(data, main = "Gráfico de Pares", pch = 19)

Visualización de los datos

library(ggplot2)

# Crear un bucle para graficar todas las variables independientes contra y
for (var in colnames(data)[-1]) {  # Excluye la primera columna (y)
  p <- ggplot(data = data, aes_string(x = var, y = "y")) +
    geom_point(color = "firebrick", size = 2) +
    labs(
      title = paste("Diagrama de dispersión: y vs", var),
      x = var,
      y = "y"
    ) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5))
  
  print(p) # Muestra cada gráfico en el panel
}

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#Mostar la posible línea de relación
# Supongamos que tienes varias columnas en tu data.frame
for (var in colnames(data)[-1]) { # Excluye la primera columna (y)
  p <- ggplot(data, aes_string(x = var, y = "y")) +
    geom_point(color = "blue", size = 2) +
    geom_smooth(method = "lm", se = FALSE, color = "red") +
    labs(title = paste("Dispersión: y vs", var),
         x = var, y = "y") +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5))
  
  print(p) # Imprime cada gráfico
}

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

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

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

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

Tratamiento de datos en caso de ser necesario

# Limpieza de datos: eliminar filas con NA, imputar medias si es necesario, filtrar outliers y eliminar duplicados
data_clean <- data %>%
  filter_all(all_vars(!is.na(.))) %>%  # Eliminar filas con valores faltantes
  mutate(across(where(is.numeric), ~ ifelse(is.na(.), mean(., na.rm = TRUE), .))) %>%  # Imputar NA con la media
  filter(if_any(where(is.numeric), 
                ~ . >= (quantile(., 0.25, na.rm = TRUE) - 1.5 * IQR(.)) & 
                   . <= (quantile(., 0.75, na.rm = TRUE) + 1.5 * IQR(.)))) %>%  # Eliminar outliers
  distinct()  # Eliminar duplicados

# Ver el resultado
data_clean

##      y x1 x2 x3  x4
## 1  240 25 24 91 100
## 2  236 31 21 90  95
## 3  270 45 24 88 110
## 4  274 60 25 87  88
## 5  301 65 25 91  94
## 6  316 72 26 94  99
## 7  300 80 25 87  97
## 8  296 84 25 86  96
## 9  267 75 24 88 110
## 10 276 60 25 91 105
## 11 288 50 25 90 100
## 12 261 38 23 89  98

Comprobar supuestos para la elección correcta del coeficiente de correlación

Prueba de normalidad

data_long <- data %>%
  pivot_longer(cols = everything(), names_to = "variable", values_to = "valor")

ggplot(data_long, aes(x = valor)) +
  geom_histogram(bins = 10, fill = "skyblue", color = "black", alpha = 0.7) +
  facet_wrap(~ variable, scales = "free") +
  labs(title = "Distribución de Variables", x = "Valor", y = "Frecuencia") +
  theme_minimal()

# Realizar la prueba de normalidad para cada columna numérica del dataframe
results <- lapply(names(data)[sapply(data, is.numeric)], function(var_name) {
  x <- data[[var_name]]  # Selecciona la columna según el nombre
  test_result <- shapiro.test(x)
  
  # Crear el mensaje con el nombre de la variable y el p-valor
  if (test_result$p.value > 0.05) {
    return(paste("La variable", var_name, "sigue una distribución normal. p-valor:", round(test_result$p.value, 4)))
  } else {
    return(paste("La variable", var_name, "NO sigue una distribución normal. p-valor:", round(test_result$p.value, 4)))
  }
})

# Mostrar los resultados
results

## [[1]]
## [1] "La variable y sigue una distribución normal. p-valor: 0.8258"
## 
## [[2]]
## [1] "La variable x1 sigue una distribución normal. p-valor: 0.7322"
## 
## [[3]]
## [1] "La variable x2 NO sigue una distribución normal. p-valor: 0.0122"
## 
## [[4]]
## [1] "La variable x3 sigue una distribución normal. p-valor: 0.6468"
## 
## [[5]]
## [1] "La variable x4 sigue una distribución normal. p-valor: 0.4748"

#forma 1 a 1 pa sacar la normalidad
shapiro.test(data$y)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$y
## W = 0.96302, p-value = 0.8258

shapiro.test(data$x1)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$x1
## W = 0.95645, p-value = 0.7322

shapiro.test(data$x2)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$x2
## W = 0.80987, p-value = 0.01216

shapiro.test(data$x3)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$x3
## W = 0.95067, p-value = 0.6468

shapiro.test(data$x4)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$x4
## W = 0.93818, p-value = 0.4748

Prueba de linealidad

library(ggplot2)

# Tu dataframe
data <- data.frame(
  y = c(240, 236, 270, 274, 301, 316, 300, 296, 267, 276, 288, 261),
  x1 = c(25, 31, 45, 60, 65, 72, 80, 84, 75, 60, 50, 38),
  x2 = c(24, 21, 24, 25, 25, 26, 25, 25, 24, 25, 25, 23),
  x3 = c(91, 90, 88, 87, 91, 94, 87, 86, 88, 91, 90, 89),
  x4 = c(100, 95, 110, 88, 94, 99, 97, 96, 110, 105, 100, 98)
)

# Crear gráficos de dispersión para cada variable independiente con respecto a 'y'
for (var in names(data)[2:5]) {  # Las columnas de x1 a x4
  # Crear el gráfico de dispersión
  p <- ggplot(data = data, aes_string(x = var, y = "y")) +
    geom_point(aes(color = y)) +  # Puntos de dispersión, con color dependiendo de 'y'
    scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +  # Colores
    geom_hline(yintercept = 0) +  # Línea horizontal en y = 0
    geom_segment(aes(xend = !!sym(var), yend = 0), alpha = 0.2) +  # Líneas de los puntos a 0
    labs(title = paste("Distribución de los residuos:", var, "vs y"),
         x = var,
         y = "y") +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
  
  # Mostrar el gráfico con print
  print(p)
}

Calcular coeficientes de correlación

determine_method <- function(x, y) {
  test_x <- shapiro.test(x)$p.value
  test_y <- shapiro.test(y)$p.value
  
  if (test_x > 0.05 & test_y > 0.05) {
    return("pearson")
  } else {
    return("spearman")
  }
}

# Calcular la correlación entre 'y' y todas las demás variables
correlations <- sapply(names(data)[-1], function(var) {
  method <- determine_method(data[[var]], data$y)
  cor(data[[var]], data$y, method = method)
})

# Mostrar resultados
correlations <- data.frame(Variable = names(correlations),
                           Correlacion = correlations)
print(correlations)

##    Variable Correlacion
## x1       x1  0.80253849
## x2       x2  0.91427065
## x3       x3  0.09285061
## x4       x4 -0.13266048

# Función para comprobar la normalidad con Kolmogorov-Smirnov
#es_normal <- function(x) {
  #test <- ks.test(x, "pnorm", mean(x), sd(x))  # Kolmogorov-Smirnov para normalidad
  #return(test$p.value > 0.05)  # Si el p-valor es mayor que 0.05, asumimos normalidad
#} 

es_normal <- function(x) {
  test <- shapiro.test(x)  # Shapiro-Wilk para normalidad
  return(test$p.value > 0.05)  # Si el p-valor es mayor que 0.05, asumimos normalidad
}

# Función para calcular la correlación según la normalidad
calcular_correlacion <- function(data) {
  n <- ncol(data)
  correlacion <- matrix(NA, nrow = n, ncol = n)
  colnames(correlacion) <- rownames(correlacion) <- colnames(data)
  
  for (i in 1:n) {
    for (j in i:n) {
      # Verificamos si ambas variables son normales
      if (es_normal(data[[i]]) && es_normal(data[[j]])) {
        correlacion[i, j] <- cor(data[[i]], data[[j]], method = "pearson")
      } else {
        correlacion[i, j] <- cor(data[[i]], data[[j]], method = "spearman")
      }
      correlacion[j, i] <- correlacion[i, j]  # Matriz simétrica
    }
  }
  
  return(correlacion)
}

# Aplicar la función a tu dataframe
correlacion <- calcular_correlacion(data)

# Graficar la matriz de correlación
corrplot(correlacion, method = "circle", type = "full", 
         col = colorRampPalette(c("red", "white", "blue"))(200), 
         tl.col = "black", tl.srt = 45, 
         addCoef.col = "black", diag = TRUE)

Modelo de Regresión

Elección del modelo adecuado

Aquí se crea un modelo de regresión lineal múltiple usando las variables independientes x1, x2, x3, etc., para predecir la variable dependiente y. Esto se hace con la función lm(). Luego, usas summary() para ver el resumen del modelo, que te da información sobre los coeficientes, errores estándar, valores p, y otros detalles estadísticos.

modelo_inical = lm(y ~ x1 + x2 + x3 + x4, data)
summary(modelo_inical)

## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.098  -9.778   1.767   6.798  13.016 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -123.1312   157.2561  -0.783    0.459  
## x1             0.7573     0.2791   2.713    0.030 *
## x2             7.5188     4.0101   1.875    0.103  
## x3             2.4831     1.8094   1.372    0.212  
## x4            -0.4811     0.5552  -0.867    0.415  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.79 on 7 degrees of freedom
## Multiple R-squared:  0.852,  Adjusted R-squared:  0.7675 
## F-statistic: 10.08 on 4 and 7 DF,  p-value: 0.00496

Verificación de multicolinealidad con VIF

El VIF (Variance Inflation Factor - Factor de inflación de la varianza) es una medida para identificar multicolinealidad en las variables predictoras de un modelo. Si el VIF de una variable es mayor que 5 o 10, eso indica que esa variable está altamente correlacionada con otras y podría ser redundante, lo que afecta la precisión del modelo.

vif(modelo_inical)

##       x1       x2       x3       x4 
## 2.322836 2.160759 1.335410 1.008733

Como ningún VIF es mayor que 5 o 10, eso quiere decir que no existe una colienealidad alta entre ninguna de las variables.

Ahora veremos el AIC del modelo inicial.

Cálculo del AIC del modelo inicial

El AIC (Akaike Information Criterion - Criterio de información de Akaike) es una métrica usada para comparar modelos. Mientras menor sea el AIC, mejor es el modelo, ya que balancea el ajuste del modelo y la complejidad (número de variables). Un AIC más bajo indica un modelo más eficiente en términos de la cantidad de variables que usa para ajustar los datos. El valor de k sería 5, ya que el modelo tiene 5 parámetros (intercepto + 4 variables independientes)

AIC(modelo_inical, k = 5)

## [1] 116.7936

# Manual
log_verosimilitud <- logLik(modelo_inical)  # Esto da la log-verosimilitud
k <- length(coef(modelo_inical))  # El número de parámetros (coeficientes), incluyendo el intercepto

# AIC manual
AIC_manual <- -2 * log_verosimilitud + 2 * k
AIC_manual

## 'log Lik.' 96.79364 (df=6)

Método Stepwise Both

Modelo de selección de variables utilizando el método Stepwise Both, que forma parte de la técnica de regresión paso a paso (Stepwise Regression). Esta técnica se utiliza para construir un modelo ajustado seleccionando de manera automática las variables más relevantes (predictoras) y descartando las que no aportan mucho al modelo.

#Este código crea un modelo vacío donde y es la variable dependiente y 1 representa el intercepto. Este modelo no tiene ninguna variable predictora (independiente), solo el intercepto. Es un modelo de referencia para el proceso de selección de variables. En otras palabras, el modelo comienza solo con el intercepto y luego se van agregando o eliminando variables.
empty.model <- lm(y ~ 1, data = data)

#Se define meta de variables que se van a considerar para entrar en el modelo. Es decir, se crea una fórmula que define las variables x1, x2, x3, x4 como las posibles variables predictoras para y.

#Este es el conjunto de variables que el algoritmo stepAIC considerará al realizar el proceso de selección paso a paso.
meta <- formula(y ~ x1 + x2 + x3 + x4)

Esta línea utiliza la función stepAIC() de R, que realiza la selección paso a paso de las variables más significativas según el AIC (Criterio de Información de Akaike). El proceso de selección se puede hacer de dos maneras:

direction = “both”: Esto significa que el proceso de selección incluye tanto la adición de nuevas variables al modelo como la eliminación de las variables menos significativas. El algoritmo evaluará si agregar o eliminar variables mejora el modelo, según el AIC.
trace = TRUE: Esta opción permite visualizar los pasos intermedios del proceso, mostrando qué variables se agregan o eliminan durante el proceso de selección.

La función stepAIC seleccionará las variables más significativas para el modelo, es decir, aquellas que contribuyen a mejorar el ajuste del modelo y que tienen un impacto significativo en la variable dependiente y.

Cuanto más bajo sea el AIC, mejor será el modelo.

El valor de representa el modelo base en cada paso, es decir, el modelo actual que no tiene ninguna variable adicional en ese paso específico de la selección. En cada iteración, el proceso de selección evalúa qué variable añadir o eliminar, y en cada uno de esos pasos, el se refiere al modelo antes de agregar o quitar la variable que se está considerando en ese paso.

La columna Df (Degrees of Freedom - Grados de libertad) muestra los grados de libertad asociados con el cambio que se hace en el modelo en cada paso del proceso de selección. Si en todos los pasos que se ve que Df = 1, significa que solo se está agregando o eliminando una variable a la vez en el modelo.
Sum of Sq (Suma de los cuadrados): Esta columna muestra la cantidad de variación explicada por el modelo al agregar o quitar una variable. Es la diferencia entre el RSS (Residual Sum of Squares) antes y después de la inclusión o eliminación de una variable.

Por ejemplo, si al agregar una variable, la suma de los cuadrados de los residuos baja considerablemente, significa que esa variable explica una cantidad significativa de la variabilidad en los datos. Suma de cuadrados (Sum of Sq) muestra cuánto mejora (o empeora) el modelo al agregar la variable en cuestión.

RSS (Residual Sum of Squares - Suma de los cuadrados residuales): Esta columna muestra la suma de los cuadrados de los residuos, es decir, la cantidad de variabilidad en los datos que no es explicada por el modelo.

Cuanto menor sea el RSS, mejor es el modelo, ya que indica que el modelo ajusta más de cerca los datos. El RSS cambia en cada paso según si se agregan o se eliminan variables del modelo.

# Selección Stepwise Both (stepAIC)
metodoboth <- stepAIC(empty.model, trace=TRUE, direction="both", scope=meta)

## Start:  AIC=77.67
## y ~ 1
## 
##        Df Sum of Sq    RSS    AIC
## + x2    1    4495.0 2077.9 65.851
## + x1    1    4233.4 2339.5 67.273
## <none>              6572.9 77.670
## + x4    1     115.7 6457.2 79.457
## + x3    1      56.7 6516.2 79.566
## 
## Step:  AIC=65.85
## y ~ x2
## 
##        Df Sum of Sq    RSS    AIC
## + x1    1     766.2 1311.7 62.330
## <none>              2077.9 65.851
## + x4    1      82.1 1995.8 67.367
## + x3    1       0.0 2077.9 67.851
## - x2    1    4495.0 6572.9 77.670
## 
## Step:  AIC=62.33
## y ~ x2 + x1
## 
##        Df Sum of Sq    RSS    AIC
## + x3    1    234.88 1076.8 61.962
## <none>              1311.7 62.330
## + x4    1     77.59 1234.1 63.598
## - x1    1    766.22 2077.9 65.851
## - x2    1   1027.82 2339.5 67.273
## 
## Step:  AIC=61.96
## y ~ x2 + x1 + x3
## 
##        Df Sum of Sq     RSS    AIC
## <none>              1076.80 61.962
## - x3    1    234.88 1311.69 62.330
## + x4    1    104.34  972.46 62.739
## - x2    1    512.21 1589.01 64.632
## - x1    1   1001.11 2077.91 67.851

Una vez que stepAIC ha realizado la selección de variables, metodoboth$anova muestra el Análisis de la Varianza (ANOVA) del modelo final. El ANOVA permitirá observar cómo cambian los valores de F y el p-valor a medida que las variables se agregan o eliminan del modelo. Esto puede ayudar a ver qué variables aportan significativamente a la variabilidad explicada en la variable dependiente y.

metodoboth$anova

## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## y ~ 1
## 
## Final Model:
## y ~ x2 + x1 + x3
## 
## 
##   Step Df  Deviance Resid. Df Resid. Dev      AIC
## 1                          11   6572.917 77.66968
## 2 + x2  1 4495.0060        10   2077.911 65.85054
## 3 + x1  1  766.2227         9   1311.688 62.32996
## 4 + x3  1  234.8840         8   1076.804 61.96215

Dado a que se obtiene un modelo con un AIC más bajo que el modelo inicial, se procede a usar este como el nuevo modelo.

summary(metodoboth) #y ~ x2 + x1 + x3

## 
## Call:
## lm(formula = y ~ x2 + x1 + x3, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.810  -6.770   3.257   7.978   9.531 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -162.1350   148.3154  -1.093   0.3061  
## x2             7.6906     3.9424   1.951   0.0869 .
## x1             0.7487     0.2745   2.727   0.0260 *
## x3             2.3434     1.7739   1.321   0.2230  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.6 on 8 degrees of freedom
## Multiple R-squared:  0.8362, Adjusted R-squared:  0.7747 
## F-statistic: 13.61 on 3 and 8 DF,  p-value: 0.001652

AIC(metodoboth, k = 4) #k es la cantidad de parametros (Los Beta_i)

## [1] 108.0167

nuevo_modelo <- metodoboth
summary(nuevo_modelo)

## 
## Call:
## lm(formula = y ~ x2 + x1 + x3, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.810  -6.770   3.257   7.978   9.531 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -162.1350   148.3154  -1.093   0.3061  
## x2             7.6906     3.9424   1.951   0.0869 .
## x1             0.7487     0.2745   2.727   0.0260 *
## x3             2.3434     1.7739   1.321   0.2230  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.6 on 8 degrees of freedom
## Multiple R-squared:  0.8362, Adjusted R-squared:  0.7747 
## F-statistic: 13.61 on 3 and 8 DF,  p-value: 0.001652

Verificar supuestos

qqnorm(): Esta función genera un gráfico Q-Q (Quantile-Quantile) que compara la distribución de los residuos del modelo (residuals(nuevo_modelo)) con una distribución normal estándar. El propósito de este gráfico es evaluar si los residuos siguen una distribución normal, que es una de las suposiciones clave de los modelos de regresión.
qqline(): Añade una línea recta al gráfico Q-Q que representa los cuantiles esperados bajo una distribución normal. Si los puntos del gráfico siguen esta línea recta, eso sugiere que los residuos tienen una distribución normal.

qqnorm(residuals(nuevo_modelo))  # Crear gráfico Q-Q
qqline(residuals(nuevo_modelo))  # Añadir la línea de referencia

Este gráfico muestra la relación entre los residuos y las predicciones del modelo. Aquí, las predicciones (prediccion) son los valores ajustados por el modelo (es decir, los valores predichos), y los residuos son la diferencia entre los valores observados y los valores predichos.

geom_point(): Dibuja los puntos, donde cada punto representa un valor de predicción frente al residuo correspondiente. El color de los puntos se ajusta según el valor del residuo, con la escala de colores que va de azul a rojo.
geom_hline(yintercept = 0): Añade una línea horizontal en 0, lo que ayuda a visualizar si los residuos están distribuidos alrededor de 0, como debería ser.
geom_segment(): Dibuja líneas desde los puntos de la predicción hasta la línea de 0, mostrando cuánto se desvían los residuos para cada observación.
labs() y theme(): Ajusta la apariencia del gráfico, como el título y los ejes.

# La función lm() calcula y almacena los valores predichos por el modelo y los
# residuos.
data$prediccion <- nuevo_modelo$fitted.values
data$residuos   <- nuevo_modelo$residuals


#Distribución de los residuos 
ggplot(data = data, aes(x = prediccion, y = residuos)) +
  geom_point(aes(color = residuos)) +
  scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
  geom_hline(yintercept = 0) +
  geom_segment(aes(xend = prediccion, yend = 0), alpha = 0.2) +
  labs(title = "Distribución de los residuos", x = "predicción modelo",
       y = "residuo") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")

Prueba de normalidad

shapiro.test(data$residuos)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$residuos
## W = 0.85199, p-value = 0.03886

# Realizar la prueba de normalidad de Shapiro-Wilk
resultado_shapiro <- shapiro.test(data$residuos)

# Mostrar los resultados y conclusión
print(resultado_shapiro)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$residuos
## W = 0.85199, p-value = 0.03886

# Evaluar si los residuos siguen una distribución normal
if (resultado_shapiro$p.value > 0.05) {
  cat("\nConclusión: Los residuos siguen una distribución normal.\n")
} else {
  cat("\nConclusión: Los residuos NO siguen una distribución normal.\n")
}

## 
## Conclusión: Los residuos NO siguen una distribución normal.

Prueba de Homocedasticidad

El comando bptest(modelo) en R realiza el test de Breusch-Pagan para detectar heterocedasticidad en un modelo de regresión lineal.

Heterocedasticidad: Es una condición en la que la varianza de los residuos no es constante a lo largo de las observaciones. Esto viola una de las principales suposiciones de los modelos de regresión lineal clásica (OLS), que asume que los residuos tienen varianza constante (homocedasticidad).

Si hay heterocedasticidad, los coeficientes del modelo aún pueden ser estimados, pero:

Las pruebas estadísticas (como t-test y F-test) pueden no ser confiables.
Los intervalos de confianza y los valores p pueden ser incorrectos.

Test de Breusch-Pagan:

Hipótesis:
- H0H0: Homocedasticidad (la varianza de los residuos es constante).
- HaHa: Heterocedasticidad (la varianza de los residuos no es constante).
Procedimiento:
- Calcula los residuos del modelo de regresión.
- Regresa los cuadrados de estos residuos contra las variables explicativas o alguna función de ellas.
- Evalúa si estas variables explicativas están asociadas con los residuos cuadrados.
Interpretación:
- Si el valor pp es pequeño (menor que el nivel de significancia, típicamente 0.05), se rechaza H0H0, lo que indica que hay evidencia de heterocedasticidad.
- Si el valor pp es grande (mayor que 0.05), no se rechaza H0H0, lo que sugiere que no hay evidencia de heterocedasticidad.

bptest(nuevo_modelo)

## 
##  studentized Breusch-Pagan test
## 
## data:  nuevo_modelo
## BP = 0.14025, df = 3, p-value = 0.9866

# Realizar la prueba de Breusch-Pagan
resultado <- bptest(nuevo_modelo)

# Mostrar los resultados y conclusión
print(resultado)

## 
##  studentized Breusch-Pagan test
## 
## data:  nuevo_modelo
## BP = 0.14025, df = 3, p-value = 0.9866

if (resultado$p.value > 0.05) {
  cat("\nConclusión: Se cumple la homocedasticidad (varianza constante).\n")
} else {
  cat("\nConclusión: No se cumple la homocedasticidad (hay heterocedasticidad).\n")
}

## 
## Conclusión: Se cumple la homocedasticidad (varianza constante).

BP: Estadístico de la prueba Breusch-Pagan.
dfdf: Grados de libertad (número de variables explicativas incluidas).
p−valuep−value: Si es mayor a 0.05, no hay evidencia de heterocedasticidad.

ggplot(data = data, aes(x = prediccion, y = residuos)) +
  geom_point(aes(color = residuos)) +
  scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
  geom_segment(aes(xend = prediccion, yend = 0), alpha = 0.2) +
  geom_smooth(se = FALSE, color = "firebrick") +
  labs(title = "Distribución de los residuos", x = "predicción modelo",
       y = "residuo") +
  geom_hline(yintercept = 0) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")

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

Autocorrelacion De Los Residuos

dw_resultado <- dwtest(nuevo_modelo) 
print(dw_resultado)

## 
##  Durbin-Watson test
## 
## data:  nuevo_modelo
## DW = 1.4787, p-value = 0.1336
## alternative hypothesis: true autocorrelation is greater than 0

ggplot(data = data, aes(x = seq_along(residuos), y = residuos)) +
  geom_point(aes(color = residuos)) +
  scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
  geom_line(size = 0.3) +
  labs(title = "Distribución de los residuos", x = "index", y = "residuo") +
  geom_hline(yintercept = 0) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Se cumplen los supuestos.

Outliers y Alto Leverage

Prueba De Bonferroni

outlierTest(nuevo_modelo)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 9 -2.212231           0.062591      0.75109

El mensaje indica que no se detectaron residuos studentizados que sean significativos después de aplicar la corrección de Bonferroni. En detalle:

“No Studentized residuals with Bonferroni p < 0.05”: Ningún residuo estandarizado tiene un valor p (ajustado por Bonferroni) menor que 0.05, lo que significa que no hay evidencia estadística de outliers en el modelo.

“Largest |rstudent|:” Indica el residuo estandarizado más grande (en valor absoluto) encontrado en el modelo, incluso si no es significativo. Si no aparece un número después de esta línea, significa que el cálculo no encontró valores destacados.

# Definir la función para testear outliers con un formato más amigable
test_outliers <- function(modelo) {
  # Realizar el test de outliers
  resultado <- outlierTest(modelo)
  
  # Verificar si el resultado tiene outliers
  if (is.null(resultado)) {
    cat("No se han encontrado outliers significativos según la prueba de Bonferroni.\n")
  } else {
    cat("Prueba de outliers (con corrección de Bonferroni):\n")
    
    # Imprimir la lista de outliers
    cat("Outliers detectados:\n")
    print(resultado)
    
    # Explicación del proceso de Bonferroni
    cat("\nExplicación: La prueba de Bonferroni ajusta el valor p para múltiples pruebas para controlar el error tipo I.\n")
    cat("Esto significa que, con una corrección de Bonferroni, solo se consideran significativos los residuos con un valor p inferior a un umbral más estricto.\n")
    cat("El umbral p ajustado es: alpha/n, donde 'alpha' es el nivel de significancia (usualmente 0.05) y 'n' es el número de observaciones.\n")
    
    # Mostrar el umbral de significancia ajustado
    umbral_bonferroni <- 0.05 / length(resultado$obs)
    cat("Umbral p ajustado (Bonferroni) para determinar outliers: ", round(umbral_bonferroni, 4), "\n")
  }
}

# Ejemplo de uso (asumiendo que tienes un modelo llamado 'nuevo_modelo')
test_outliers(nuevo_modelo)

## Prueba de outliers (con corrección de Bonferroni):
## Outliers detectados:
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 9 -2.212231           0.062591      0.75109
## 
## Explicación: La prueba de Bonferroni ajusta el valor p para múltiples pruebas para controlar el error tipo I.
## Esto significa que, con una corrección de Bonferroni, solo se consideran significativos los residuos con un valor p inferior a un umbral más estricto.
## El umbral p ajustado es: alpha/n, donde 'alpha' es el nivel de significancia (usualmente 0.05) y 'n' es el número de observaciones.
## Umbral p ajustado (Bonferroni) para determinar outliers:  Inf

# Valores de leverage
hatvalues_model <- hatvalues(nuevo_modelo)
data$Leverage <- hatvalues_model
n <- nrow(data)
p <- length(coef(nuevo_modelo))
umbral_leverage <- 2 * (p / n)
plot(hatvalues_model, main = "Valores de Leverage", xlab = "Observaciones", ylab = "Leverage")
abline(h = umbral_leverage, col = "red", lty = 2)

# Distancia de Cook
cooks_dist <- cooks.distance(nuevo_modelo)
plot(cooks_dist, main = "Distancia de Cook", xlab = "Observaciones", ylab = "Distancia de Cook")
abline(h = 4 / n, col = "red", lty = 2)

data$Cooks_Distance <- cooks_dist

# Identificar observaciones influyentes
influyentes <- which(cooks_dist > (4 / n))
cat("Observaciones influyentes:", influyentes, "\n")

## Observaciones influyentes: 1

# Prueba de Bonferroni para Outliers
bonferroni_test <- outlierTest(nuevo_modelo, cutoff=Inf)
print("Resultados del Test de Bonferroni para Outliers:")

## [1] "Resultados del Test de Bonferroni para Outliers:"

print(bonferroni_test)

##      rstudent unadjusted p-value Bonferroni p
## 9  -2.2122309           0.062591      0.75109
## 1  -1.9395623           0.093595           NA
## 10 -1.1773927           0.277510           NA
## 11  0.9035969           0.396240           NA
## 12  0.8776307           0.409250           NA
## 5   0.8311617           0.433310           NA
## 3   0.7190465           0.495410           NA
## 7   0.5844037           0.577290           NA
## 6   0.5571649           0.594770           NA
## 4  -0.4783796           0.646960           NA

# Gráfico de diagnóstico usando influenceIndexPlot
influenceIndexPlot(nuevo_modelo, vars=c("Cook", "hat", "Bonf"), id.n=4, las=1)

## Warning in plot.window(...): "id.n" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "id.n" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "id.n" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "id.n" is not a
## graphical parameter

## Warning in box(...): "id.n" is not a graphical parameter

## Warning in title(...): "id.n" is not a graphical parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.n" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.n" is not a
## graphical parameter

## Warning in plot.window(...): "id.n" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "id.n" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "id.n" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "id.n" is not a
## graphical parameter

## Warning in box(...): "id.n" is not a graphical parameter

## Warning in title(...): "id.n" is not a graphical parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.n" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.n" is not a
## graphical parameter

## Warning in plot.window(...): "id.n" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "id.n" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "id.n" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "id.n" is not a
## graphical parameter

## Warning in box(...): "id.n" is not a graphical parameter

## Warning in title(...): "id.n" is not a graphical parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.n" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.n" is not a
## graphical parameter

## Manejo de Outliers Usando M- estimación

modelo_robusto <- rlm(y ~ x1 + x2+ x3, data = data)
AIC(modelo_robusto)

## [1] 98.02932

AIC(nuevo_modelo)

## [1] 98.01668

Predicción

# Para los nuevos datos de predicción (nuevo conjunto de predictores)
nuevos_datos <- data.frame(x1 = 5, x2 = 10, x3= 15)

# Predicción con intervalos de confianza
pred_conf <- predict(modelo_robusto, newdata = nuevos_datos, interval = "confidence", level = 0.95)

# Predicción con intervalos de predicción
pred_pred <- predict(modelo_robusto, newdata = nuevos_datos, interval = "prediction", level = 0.95)

## Warning in predict.lm(modelo_robusto, newdata = nuevos_datos, interval = "prediction", : Assuming constant prediction variance even though model fit is weighted

# Ver los resultados
print(pred_conf)

##        fit       lwr      upr
## 1 -45.0246 -298.4615 208.4123

print(pred_pred)

##        fit       lwr      upr
## 1 -45.0246 -299.5828 209.5336

# Predicciones del modelo
predicciones <- predict(modelo_robusto)

# Gráfico
plot(data$y, predicciones, 
     xlab = "Valores Reales", 
     ylab = "Predicciones", 
     main = "Predicciones vs Valores Reales",
     pch = 19, col = "blue")
abline(0, 1, col = "red", lwd = 2)  # Línea y = x para referencia

Regresión lineal multiple

Eliana R

2024-11-29

Contexto del ejercicio

Descriptivo de los datos

Visualización de los datos

Tratamiento de datos en caso de ser necesario

Comprobar supuestos para la elección correcta del coeficiente de correlación

Prueba de normalidad

Prueba de linealidad

Calcular coeficientes de correlación

Modelo de Regresión

Elección del modelo adecuado

Verificación de multicolinealidad con VIF

Cálculo del AIC del modelo inicial

Método Stepwise Both

Verificar supuestos

Prueba de normalidad

Prueba de Homocedasticidad

Test de Breusch-Pagan:

Autocorrelacion De Los Residuos

Outliers y Alto Leverage

Prueba De Bonferroni

Predicción