Este material hace parte del proyecto de investigación “EVALUACIÓN DE UN MODELO PARA EL APRENDIZAJE ADAPTATIVO SOBRE LA TÉCNICA DE REGRESIÓN LINEAL SIMPLE” de la Escuela de Ciencias Básicas y Aplicadas de la Universidad de La Salle, Bogotá, Colombia.

Si ya has revisado el material publicado en en el link https://tatiana-pamela-jimenez-valderrama.shinyapps.io/RegresionNivel2/, sigue ahora las siguientes instrucciones para realizar el ejercicio propuesto como ejemplo (Ejercicio1).

Ingreso de datos

Pureza=c(86.91, 89.85, 90.28, 86.34, 92.58, 87.33, 86.29, 91.86, 95.61, 89.86,
         96.73, 99.42, 98.66, 96.07, 93.65, 87.31, 95.00, 96.85, 85.20, 90.56)
Hidrocarburos=c(1.02, 1.11, 1.43, 1.11, 1.01, 0.95, 1.11, 0.87, 1.43, 1.02,
                1.46, 1.55, 1.55, 1.55, 1.40, 1.15, 1.01, 0.99, 0.95, 0.98)
datos=cbind.data.frame(Pureza,Hidrocarburos)
n=dim(datos)[1]

Gráfico de dispersión

library(ggplot2)
ggplot(datos, aes(x = Hidrocarburos, y = Pureza)) +
  geom_point(color = "blue", size = 3) +             # puntos azules
  labs(
    title = "Relación entre Hidrocarburos y Pureza",
    x = "Hidrocarburos (%)",
    y = "Pureza (%)"
  ) +
  theme_minimal()

No olvides ir interpretando cada elemento del análisis: ¿Qué puedes observar en el gráfico de dispersión sobre la posible relación entre el porcentaje de hidrocarburos y el porcentaje de pureza?

Coeficientes de correlación y de determinación

Coeficiente de correlación

(r=cor(datos$Hidrocarburos,datos$Pureza))
## [1] 0.6237968

No olvides ir interpretando cada elemento del análisis: A partir de este valor estimado del coeficiente de correlación, ¿que puedes decir sobre la posible relación entre las dos variables?

Intervalo de confianza sobre el coeficiente de correlación

inferencia_rho=cor.test(datos$Pureza,datos$Hidrocarburos)
(intervalo_rho=inferencia_rho$conf.int)
## [1] 0.2503961 0.8356439
## attr(,"conf.level")
## [1] 0.95

No olvides ir interpretando cada elemento del análisis: A partir del intervalo construido, al 95% de confianza, ¿que puedes decir sobre la posible relación entre las dos variables?

Prueba de hipótesis sobre el coeficiente de correlación

\[Ho: \rho=0\]

\[Ha: \rho \ne 0\]

(p.valor_rho=inferencia_rho$p.value)
## [1] 0.003291122

No olvides ir interpretando cada elemento del análisis: A partir del p_valor calcualdo, ¿que puedes decir sobre la posible relación lineal entre las dos variables?

Coeficiente de determinación

(R2=r^2)
## [1] 0.3891224

No olvides ir interpretando cada elemento del análisis: A partir del \(R^2\) estimado, ¿que puedes decir sobre la viabilidad de usar un modelo de regresión lineal simple para expresar la relación entre las dos variables?

Ecuación de ajuste

Valores de bo y b1

Estimamos los valores para \(\beta_0\) y para \(\beta_1\)

(ecuacion=lm(datos$Pureza~datos$Hidrocarburos))
## 
## Call:
## lm(formula = datos$Pureza ~ datos$Hidrocarburos)
## 
## Coefficients:
##         (Intercept)  datos$Hidrocarburos  
##               77.86                11.80

Recuerda: Escribir la ecuación resultante e interpretar los valores obtenidos.

valor de error estándar del modelo

(MS_res=sum((ecuacion$residuals)^2)/(n-2))
## [1] 12.93524

Inferencia sobre los parámetros \(\beta_0\) y \(\beta_1\)

Intervalos de confianza

confint(ecuacion, level = 0.95) #El valor de la confianza la puedes cambiar según el ejercicio a resolver 
##                         2.5 %   97.5 %
## (Intercept)         69.041747 86.68482
## datos$Hidrocarburos  4.479066 19.12299

Recuerda: Interpreta cada uno de los intervalos que acabas de obtener. Ten encuenta que, en el ejercicio se está suponiendo que el valor real de \(\beta_0\) es 70.

Pruebas de hipótesis

Con las siguientes instrucciones encontrarás la información nececaria para las pruebas de hipótesis correspondientes a los parámetros del modelo:

\(Ho: \beta_0 = 0\) vs \(Ha: \beta_0 \ne 0\)

\(Ho: \beta_1=0\) vs \(Ha: \beta_1 \ne 0\)

resumen_ecuacion=summary(ecuacion) #en la pestaña help puedes consultar los resultados que este resumen guarda.
(resumen_ecuacion$coefficients)
##                     Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)         77.86328   4.198888 18.543786 3.537382e-13
## datos$Hidrocarburos 11.80103   3.485119  3.386119 3.291122e-03

Recuerda: interpreta el p.valor, que se encuentra en la columna \((Pr(>|t|))\).

Ahora bien, con este procedimiento rápido, solo obtenemos información para la prueba de hipótesis donde se somete a prueba el valor 0 para \(\beta_0\). Sin emabargo, en el ejercicio nos pide que sometamos a prueba el valor 70.

\(Ho: \beta_0 = 70\) vs \(Ha: \beta_0 \ne 70\)

Para esto debemos realizar algunos cálculos paso a paso:

# Extraer coeficientes
intercepto_est <- coef(ecuacion)[1]
se_intercepto <- summary(ecuacion)$coefficients[1, 2] #Error estándar de la estimación de beta_o

# Hipótesis nula: intercepto = 70
intercepto_esperado <- 70

# Calcular estadístico t
t <- (intercepto_est - intercepto_esperado) / se_intercepto

# Grados de libertad
gl <- df.residual(ecuacion)

# Calcular p-value (prueba bilateral)
p_value <- 2 * pt(-abs(t), df = gl)

# Mostrar resultado
cat("El valor.p, es de", p_value, "\n")
## El valor.p, es de 0.07744474

¿Qué puedes concluir a partir de este valor.p sobre el posible valor de 70 para \(\beta_0\)

validación del modelo

Para la validación del modelo trabajamos con los residuales, por ello nuestro primer paso es guardar los residuales:

datos$residuales=ecuacion$residuals

Simetría

A nivel descriptivo

Se puede realizar un gráfico descriptivo para observar el comportamiento de los residuales y la posibilidad de que ellos provengan de una población con distribución normal.

ggplot(datos, aes(x = "", y = residuales)) +
  geom_boxplot(fill = "lightgreen", color = "black") +
  labs(
    title = "Boxplot de Pureza",
    x = "",
    y = "Residuales"
  ) +
  theme_minimal()

También podemos hacer uso del valor del coeficiente de asimetría de los residuales:

library(moments)
skewness(datos$residuales)
## [1] 0.265499

¿Será este valor lo suficientemente cercano a cero para considerar que los residuales provienen de una población con distribución simétrica?

A nivel inferencial

Para tomar una decisión sobre si la distribución es simétrica llevaremos acabo la siguiente prueba de hipótesis, llamada Test de Agostino:

\(Ho:G1=0\) (distribución simétrica)

\(Ha: G1 \ne 0\) (distribució asimétrica)

(prueba_simetria=agostino.test(datos$residuales))
## 
##  D'Agostino skewness test
## 
## data:  datos$residuales
## skew = 0.26550, z = 0.59266, p-value = 0.5534
## alternative hypothesis: data have a skewness

A partir del valor_p, ¿cuál es tu conclusión sobre el valor del coeficiente de asimetría?

Curtosis

A nivel descriptivo

Podemos construir un histograma y calcular el coeficiente de curtosis de los residuales

# Calcular la media y desviación estándar de los residuales
media <- mean(datos$residuales, na.rm = TRUE)
sd <- sd(datos$residuales, na.rm = TRUE)

# Crear el histograma con la curva normal superpuesta
ggplot(datos, aes(x = residuales)) +
  geom_histogram(aes(y = ..density..), 
                 bins = 5, 
                 fill = "skyblue", 
                 color = "black") +
  stat_function(fun = dnorm, 
                args = list(mean = media, sd = sd), 
                color = "red", 
                size = 1) +
  labs(title = "Histograma de Residuales con Curva Normal",
       x = "Residuales",
       y = "Densidad") +
  theme_minimal()
## 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.
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

¿Cuál es tu opinión sobre la curtosis del gráfico anterior?

Calculamos ahora el coeficiente de curtosis de los residuales

(kurtosis(datos$residuales))
## [1] 2.190697

¿Crees que este valor es lo suficientemente cercano a 3 para considerar una disribución mesocúrtica?

A nivel inferencial

Vamos a llevar a cabo la prueba de Anscombe sobre el coeficiente de curtosis

\(Ho: G2=3\) (distribución mesocúrtica)

\(Ha: G2 \ne 3\)

(anscombe.test(datos$residuales))
## 
##  Anscombe-Glynn kurtosis test
## 
## data:  datos$residuales
## kurt = 2.19070, z = -0.67672, p-value = 0.4986
## alternative hypothesis: kurtosis is not equal to 3

¿Cuál es tu conclusión sobre el tipo de curtosis de la distribución de los errores?

Normalidad

A nivel descriptivo

Utilizaremos el gráico de probabilidad normal

ggplot(datos, aes(sample = residuales)) +
  stat_qq(color = "blue", size = 2) +
  stat_qq_line(color = "red", linetype = "dashed") +
  labs(title = "Gráfico de Probabilidad Normal (Q-Q plot)",
       x = "Cuantiles teóricos",
       y = "Cuantiles de los residuales") +
  theme_minimal()

¿Qué puedes concluir de este gráfico?

A nivel inferencial

Utilizaremos la prueba de jarque-bera

\[Ho: X \sim N(0,\sigma^2)\]

\[Ha: X \nsim N(0,\sigma^2)\]

(jarque.test(datos$residuales))
## 
##  Jarque-Bera Normality Test
## 
## data:  datos$residuales
## JB = 0.78078, p-value = 0.6768
## alternative hypothesis: greater

¿Cuál es la decisió estadística sobre el ajuste de la distribución normal de los errores?

Independencia

A nivel descriptivo

Construyamos un gráfico con los residuales para revisar si se observa algún patrón de comportamiento entre ellos. Si no se observa ninguno, podemos suponer que los errores son independientes.

plot(datos$residuales,main="Residuales",xlab="Orden de los residuales",ylab="Residuales",col="darkblue")

¿Encuentras algún patrón de comportamiento en los datos? ¿Pueden suponerse independientes?

A nivel inferencial

\(Ho: DW=2\) (Errores independientes)

\(Ha: DW \ne 2\) (Errores autocorrelacionados)

library(lmtest)
dwtest(ecuacion, alternative = "two.sided") #Se debe colocar el nombre del modelo lineal
## 
##  Durbin-Watson test
## 
## data:  ecuacion
## DW = 1.9108, p-value = 0.7331
## alternative hypothesis: true autocorrelation is not 0

A partir del valor_p, ¿qué se puede concluir sobre la independencia de los errores?

Homogeneidad

En este caso debemos mostrar que los errores provienen de distribuciones con igual varianza

\(Ho: \sigma_1^2=\sigma_2^2= \dots =\sigma_n^2\) (Las varianzas son iguales)

\(Ha: \sigma_i^2 > \sigma_j^2\) (por lo menos hay una que es mayor a otras)

A nivel descriptivo

Construiremos un gráfico de dispersión entre los residuales y los valores de \(X\). Si no se observa algún patrón de comportamiento se puede suponer que los errores provienen de poblaciones con igual varianza. Pero, sí observa algún patrón de comportamiento no es posible suponer la igualdad de varianzas.

plot(datos$residuales,datos$Hidrocarburos,main = "Gráfico para mostrar homogeneidad", ylab = "Residuales",
     xlab = "Hidrocarburo (%)", col="darkgreen")

¿Qué puedes concluir sobre la homogeneidad de los errores?

A nivel inferencial

Vamos a utilizarla prueba de Bruch-Pagan

\(Ho: \sigma_1^2=\sigma_2^2= \dots =\sigma_n^2\) (Las varianzas son iguales)

\(Ha: \sigma_i^2 > \sigma_j^2\) (por lo menos hay una que es mayor a otras)

(bptest(ecuacion))
## 
##  studentized Breusch-Pagan test
## 
## data:  ecuacion
## BP = 1.9169, df = 1, p-value = 0.1662

¿Qué puedes decir sobre el cumplimiento de esta condición?

Importante: Ahora, construye tu análisis completo sobre la relación entre estas dos variables y la adecuancia de la ecuación de ajuste.

  1. Montgomery, D. C., Peck, E. A., & Vining, G. Introducción al análisis de regresión lineal. Tercera edición. Editorial CECSA., México. 2005. Doi: 10.4272.↩︎