Análisis de Regresión Lineal Simple: Corrosión del Hierro

1 Universidad Rafael Landívar

1.1 Proyecto - Programación estadística con R

2 Instalación y carga de paquetes

# Instalar paquetes (solo ejecutar una vez)
install.packages("tidyverse")
install.packages("broom")
install.packages("pubh")
install.packages("sjPlot")
install.packages("lmtest")

# Cargar librerías
library(ggplot2)
library(tidyverse)
library(broom)
library(pubh)
library(sjPlot)
library(lmtest)

3 Importación de datos

Creamos el banco de datos con información sobre contenido de hierro y pérdida de peso por corrosión:

hierro <- c(0.01, 0.48, 0.71, 0.95, 1.19, 0.01, 0.48, 1.44, 0.71, 
            1.96, 0.01, 1.44, 1.96)
peso <- c(127.6, 124, 110.8, 103.9, 101.5, 130.1, 122, 92.3, 113.1, 
          83.7, 128, 91.4, 86.2)
corrosion <- data.frame(hierro, peso)

# Mostrar los primeros datos
head(corrosion)

hierro	peso
0.01	128
0.48	124
0.71	111
0.95	104
1.19	102
0.01	130

3.1 Visualización de los datos

ggplot(corrosion, aes(x = hierro, y = peso)) +
  geom_point(size = 3, color = "steelblue") +
  labs(x = "Contenido en hierro", 
       y = "Pérdida de peso",
       title = "Relación entre contenido de hierro y pérdida de peso") +
  theme_minimal()

Diagrama de dispersión entre contenido de hierro y pérdida de peso

4 Recta de regresión

Representamos la recta de regresión lineal ajustada a los datos:

ggplot(corrosion, aes(x = hierro, y = peso)) + 
  geom_point(size = 3, color = "steelblue") + 
  stat_smooth(method = "lm", se = FALSE, color = "red") + 
  labs(x = "Contenido en hierro", 
       y = "Pérdida de peso",
       title = "Recta de regresión lineal") +
  theme_bw()

Recta de regresión lineal ajustada

4.1 Ajuste del modelo

fit <- lm(peso ~ hierro, data = corrosion)

4.2 Coeficientes del modelo

tidy(fit)

term	estimate	std.error	statistic	p.value
(Intercept)	130	1.4	92.5	2.93e-17
hierro	-24	1.28	-18.8	1.06e-09

glm_coef(fit)

Parameter	Coefficient	Pr(>|t|)
(Intercept)	129.79 (126.7, 132.87)	< 0.001
hierro	-24.02 (-26.84, -21.2)	< 0.001

4.3 Gráfico del ajuste

plot_model(fit, "pred", terms = "hierro", 
           ci.lvl = NA, 
           show.data = TRUE, 
           axis.title = c("Contenido en hierro", "Peso"),
           title = "Modelo de regresión lineal simple")

Modelo de regresión lineal simple

5 Análisis de la varianza

5.1 Medidas de bondad de ajuste

glance(fit)

r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual	nobs
0.97	0.967	3.06	352	1.06e-09	1	-31.9	69.8	71.5	103	11	13

5.2 Tabla ANOVA

anova(fit)

Df	Sum Sq	Mean Sq	F value	Pr(>F)
1	3.29e+03	3.29e+03	352	1.06e-09
11	103	9.35

5.3 Interpretación

Se concluye que el modelo obtenido resulta útil para explicar la mayor parte de la variabilidad existente en la respuesta a partir de la variabilidad explicada por dicho modelo.

Como R² = 0.9697, significa que el 97% de la variabilidad de la pérdida de peso es explicada por la recta ajustada.

5.4 Resumen del modelo

summary(fit)

## 
## Call:
## lm(formula = peso ~ hierro, data = corrosion)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7980 -1.9464  0.2971  0.9924  5.7429 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  129.787      1.403   92.52  < 2e-16 ***
## hierro       -24.020      1.280  -18.77 1.06e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058 on 11 degrees of freedom
## Multiple R-squared:  0.9697, Adjusted R-squared:  0.967 
## F-statistic: 352.3 on 1 and 11 DF,  p-value: 1.055e-09

6 Diagnóstico del modelo de RLS

6.1 Residuos y medidas de diagnóstico

diagnostico <- fortify(fit)

6.2 Gráfico de residuos estandarizados vs ajustados

ggplot(diagnostico, aes(x = .fitted, y = .stdresid)) + 
  geom_point(size = 3, color = "steelblue") + 
  stat_smooth(method = "lm", se = FALSE, color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  labs(x = "Valores ajustados",
       y = "Residuos estandarizados",
       title = "Residuos vs Valores ajustados") +
  theme_bw()

Residuos estandarizados vs valores ajustados

6.3 Test de homocedasticidad (Breusch-Pagan)

bptest(fit)

## 
##  studentized Breusch-Pagan test
## 
## data:  fit
## BP = 0.024539, df = 1, p-value = 0.8755

Interpretación: Dado que el p-valor del contraste es superior a 0.05, no podemos rechazar la hipótesis nula del test de homocedasticidad, y por tanto podemos concluir que se verifica la hipótesis de homogeneidad o varianza constante.

6.4 Normalidad de los residuos

6.4.1 Gráfico Q-Q Plot

ggplot(diagnostico, aes(sample = .stdresid)) + 
  stat_qq(color = "steelblue", size = 3) + 
  geom_abline(color = "red") +
  labs(title = "Gráfico Q-Q de residuos estandarizados",
       x = "Cuantiles teóricos",
       y = "Cuantiles muestrales") +
  theme_bw()

Gráfico Q-Q de residuos estandarizados

6.4.2 Test de Shapiro-Wilk

shapiro.test(diagnostico$.stdresid)

## 
##  Shapiro-Wilk normality test
## 
## data:  diagnostico$.stdresid
## W = 0.92905, p-value = 0.3312

Interpretación: Dado que el p-valor del contraste es superior a 0.05, no podemos rechazar la hipótesis nula del contraste de normalidad, y por tanto podemos considerar que los residuos se distribuyen normalmente.

6.5 Independencia de los residuos

6.5.1 Gráfico de autocorrelación

acf(diagnostico$.stdresid, main = "Función de autocorrelación")

Función de autocorrelación de los residuos

6.5.2 Test de Durbin-Watson

dwtest(fit, alternative = "two.sided")

## 
##  Durbin-Watson test
## 
## data:  fit
## DW = 2.5348, p-value = 0.2952
## alternative hypothesis: true autocorrelation is not 0

Interpretación: Dado que el p-valor del contraste es superior a 0.05, no podemos rechazar la hipótesis nula del contraste de incorrelación, y por tanto podemos considerar que los residuos se distribuyen de forma independiente.

7 Conclusión

Se puede concluir que se verifican todas las hipótesis del modelo de regresión lineal simple:

1. Linealidad: Se verifica mediante la recta horizontal en el gráfico de residuos.
2. Homocedasticidad: Los residuos se comportan de forma aleatoria sin agrupaciones ni tendencias (confirmado por el test de Breusch-Pagan, p-valor = 0.8755).
3. Normalidad: Los residuos se distribuyen normalmente (confirmado por el test de Shapiro-Wilk, p-valor = 0.3312 y el gráfico Q-Q).
4. Independencia: Los residuos son independientes (confirmado por el test de Durbin-Watson, p-valor = 0.2952).

7.1 Ecuación de regresión estimada

La ecuación de regresión lineal simple obtenida es:

\[\hat{y} = 129.787 + -24.02x\]

Donde:

• \(\hat{y}\) es la pérdida de peso estimada
• \(x\) es el contenido en hierro

7.2 Tabla con los residuos estandarizados

residuos <- data.frame(
  Observacion = 1:nrow(diagnostico),
  Hierro = corrosion$hierro,
  Peso_Observado = corrosion$peso,
  Peso_Ajustado = round(diagnostico$.fitted, 2),
  Residuo_Estandarizado = round(diagnostico$.stdresid, 3)
)
knitr::kable(residuos, 
             caption = "Residuos estandarizados del modelo",
             col.names = c("Obs.", "Hierro", "Peso Obs.", "Peso Ajust.", "Resid. Std."))

Residuos estandarizados del modelo

Obs.	Hierro	Peso Obs.	Peso Ajust.	Resid. Std.
1	0.01	127.6	129.55	-0.715
2	0.48	124.0	118.26	1.984
3	0.71	110.8	112.73	-0.659
4	0.95	103.9	106.97	-1.045
5	1.19	101.5	101.20	0.102
6	0.01	130.1	129.55	0.203
7	0.48	122.0	118.26	1.293
8	1.44	92.3	95.20	-1.018
9	0.71	113.1	112.73	0.125
10	1.96	83.7	82.71	0.384
11	0.01	128.0	129.55	-0.568
12	1.44	91.4	95.20	-1.334
13	1.96	86.2	82.71	1.350

---

Nota: Hay que tener mucho cuidado con la interpretación de los gráficos de diagnóstico cuando el tamaño de muestra es pequeño, ya que resulta difícil apreciar tendencias o agrupaciones en los residuos. Por este motivo es importante complementar con los tests de diagnóstico formales.