Presentación del problema.

Trabajaremos con los datos de la base de datos airquality de R; se trata de medidas diarias de la calidad del aire y de variables climatológicas tomadas en Nueva Yory entre mayo y septiembre de 1973. En concreto, se consideran las variables

Se pretende establecer un modelo lineal entra ambas variales, con Temp como variable explicativa y Wind como variable respuesta.

Carga de datos.

Teclea en un script y ejecuta las siguientes órdenes

x = airquality$Temp
y = airquality$Wind

Enunciado.

Ejercicio 1

Comprueba que la lectura ha sido correcta con la funcióin summary(): visualiza las primeras lineas de cada vector

summary(x)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  56.00   72.00   79.00   77.88   85.00   97.00 
summary(y)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.700   7.400   9.700   9.958  11.500  20.700

las variables son numéricas y no parece haber datos ausentes.

Ejercicio 2

Visualiza la nube de puntos para determinar de forma visual si tiene sentido utilizar un modelo lineal. 

plot(x, y)

El modelo lineal parece adecuado. Recuerda poner primero la variable independiente y depués la dependiente.

Ejercicio 3

Crea el modelo de regresión lineal y comprueba si se cumplen las hipótesis necesarias.

Recuerda que para crear el modelo lineal hay que poner primero la variable dependiente y después la independiente:

(lmXY = lm(y ~ x))

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
    23.2337      -0.1705

Visualizamos la nube de puntos y la recta de regresión

plot(x, y)
abline(lmXY, lwd = 3, col = "red")

Para comprobar gráficamente las hipótesis:

par(mfrow = c(1, 2))

for(i in 1:2){
  plot(lmXY, which = i)
}

par(mfrow = c(1, 1))

Los gráficos sugieren normalidad y homocedasticidad en los residuos. Aún así, se aconseja hacer los correspondientes contrastes:

# contastes de normalidad
shapiro.test(lmXY$residuals)

    Shapiro-Wilk normality test

data:  lmXY$residuals
W = 0.9856, p-value = 0.1132

El contraste de normalidad es satisfactorio y el de homocedasticidad

# ---- homocedasticityTests
# H0: homocedasticidad en poblaciones normales
library(gvlma)
gvlma(lm(y ~ x), alphalevel = 0.01)

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
    23.2337      -0.1705  


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.01 

Call:
 gvlma(x = lm(y ~ x), alphalevel = 0.01) 

                       Value p-value                Decision
Global Stat        5.1614845 0.27113 Assumptions acceptable.
Skewness           3.2306304 0.07227 Assumptions acceptable.
Kurtosis           0.0005899 0.98062 Assumptions acceptable.
Link Function      0.0417659 0.83807 Assumptions acceptable.
Heteroscedasticity 1.8884982 0.16937 Assumptions acceptable.

también se cumple. Recuerda que hay que fijarse en la última línea, y que en este caso el rótulo puede llevar a equívoco: H0 es “los residuos son homocedásticos” (aunque el rótulo hable de heterocedasticidad). El modelo de regresión lineal tiene sentido.

Ejercicio 4

Calcula los coeficientes de la recta de regresión y sus intervalos de confianza al 95% de confianza Interprétalos

Los coeficientes de la recta son

lmXY$coefficients
(Intercept)           x 
 23.2336881  -0.1704644

Observa que la pendiente es megativa, lo que quiere decir que temperaturas altas tienen asociados vientos menos veloces.

Por otro lado, los intervalos de confianza al 95% son

########################################
# INFERENCIA SOBRE LA RECTA DE REGRESION

# intervalo confianza coeficientes recta
confint(lmXY, level = 0.95)
                 2.5 %    97.5 %
(Intercept) 19.0600210 27.407355
x           -0.2236649 -0.117264

  • El intervalo de confianza para el término independiente, en general, no tiene interés, porque se haya fuera de la zona en la que hemos recogido valores (extrapolación).

  • El intervalo de confianza para la pendiente indica que por cada unidad (grado F) que aumenta Temp se espera que la velocidad del viento disminuya en entre 0.12 y 0.22 millas por hora.

Ejercicio 5

Contrasta las hipótesis \(H_0:\,\beta_1 = -0.12\) frente a \(H_0:\,\beta_1 \neq -0.12\) al nivel de significación 0.05.

Hay que ajustar \(H_0\) en el código de la plantilla

library(smatr)
CH_b1 = slope.test(y, x, 
                   test.value = -0.12,  
                   method = "OLS", alpha = 0.05)

CH_b1$p # P-value
[1] 0.06283671

de donde se decide no rechazar \(H_0\).

Ejercicio 6

Calcula, para las temperaturas Temp = 65 y con Temp = 70 :

  • El valor de Wind predicho por el modelo.
  • El intervalo de confianza para la velocidad media del viento predicha (intervalo de confianza) al nivel de confianza del 95%.
  • El intervalo para el valor de la velocidad del viento predicha (intervalo de predicción) al nivel de confianza del 95%.

La predicción puntual es la que aparece bajo la etiqueta fit. Las etiquetas lwr y upr hacen referencia a lower y upper, los extremos inferior y superior del correspondiente intervalo:

x = airquality$Temp
y = airquality$Wind
lmXY = lm(y ~ x)

 # intervalo de confianza valor predicho 
predict(lmXY, newdata = data.frame(x = c(65, 70)), 
        interval = "confidence", level = 0.95)
       fit      lwr      upr
1 12.15350 11.30402 13.00298
2 11.30118 10.64714 11.95521

Este es el intervalo de confianza para el valor medio predicho

# intervalo de prediccion valor predicho 
predict(lmXY, newdata = data.frame(x = c(65, 70)), 
        interval = "prediction", level = 0.95)
       fit      lwr      upr
1 12.15350 5.887386 18.41961
2 11.30118 5.058555 17.54380

Y este el intervalo de confianza para el valor predicho. Observa que en este caso el intervalo es mucho mayor, ya que los valores siempre están más dispersos que su promedio.