Como ilustración vamos a usar los datos del ejemplo 2.1 del libro de Montgomery, Peck and Vining (2003). En el ejemplo 2.1 los autores ajustan un modelo de regresión lineal simple para explicar la Resistencia de una soldadura en función de la Edad de la soldadura.
df<-read_csv("propelent.csv", col_names=TRUE, show_col_types = FALSE)
pander(df)| observacion | resistencia | edad |
|---|---|---|
| 1 | 2159 | 15.5 |
| 2 | 1678 | 23.75 |
| 3 | 2316 | 8 |
| 4 | 2061 | 17 |
| 5 | 2208 | 5.5 |
| 6 | 1708 | 19 |
| 7 | 1785 | 24 |
| 8 | 2575 | 2.5 |
| 9 | 2358 | 7.5 |
| 10 | 2257 | 11 |
| 11 | 2165 | 13 |
| 12 | 2400 | 3.75 |
| 13 | 1780 | 25 |
| 14 | 2337 | 9.75 |
| 15 | 1765 | 22 |
| 16 | 2054 | 18 |
| 17 | 2414 | 6 |
| 18 | 2200 | 12.5 |
| 19 | 2654 | 2 |
| 20 | 1754 | 21.5 |
El objetivo es obtener:
IC del 95% para \(E(Y\mid x-0)\) cuando \(X_0=13\) semanas.
IC del 95% para \(E(Y\mid x_0)\) cuando \(x_0=2\) semanas.
IC del 90% para \(\hat{Y_0}\) cuando \(x_0=10\) semanas.
Crear el diagrama de dispersión agregando las líneas de los IC de 95% para \(E(Y\mid x_0)\) y \(\hat{Y_0}\)
mod1<- lm(resistencia~edad, data = df)
pander(mod1)| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 2628 | 44.18 | 59.47 | 4.064e-22 |
| edad | -37.15 | 2.889 | -12.86 | 1.643e-10 |
\[ \hat{Resistencia}= \hat{\beta_0} - \hat{\beta_1}\,edad +Ui \]
\[ \hat{Resistencia}=2628-37.15 \, edad \]
Teniendo en cuenta que en nuestro modelo el valor resultante en \(\hat{\beta_1}=-37.15\) nuestro modelo tendrá un signo menos ya que la relación del mismo es indirecta.
Para obtener el IC del 95% para \(E(Y\mid x_0)\) para \(x_0=13\) y \(x_0=2\) se usa el siguiente código.
Donde con un Intervalo de confianza del 95% o un \(\alpha=0,5\) nuestra \(edad=13\) y \(edad=2\)
#Nuestro nuevo dataframe contiene solo los datos de la edad que deseamos estudiar que son, que pasa cuando edad=13 y edad=2 por lo tanto estas son las unicas variables que incluiremos en el nuevo df (nwdf)
nwdf<- data.frame(edad=c(13, 2))
predictdf<-predict(object=mod1, newdata = nwdf, interval = "confidence", level = 0.95)
pander(predictdf)| fit | lwr | upr |
|---|---|---|
| 2145 | 2100 | 2190 |
| 2554 | 2471 | 2636 |
Donde fit nos indica las predicciones ajustadas para cada fila de nwdf (nuevo data frame) lwr y upr limite inferior y superior del IC al 95% para cada prediccón.
results<- data.frame(
edad=nwdf$edad,
fit=predictdf[, "fit"], #la sintaxis [, "fit"] es para extraer columnas especificas del objeto predictions
lwr=predictdf[, "lwr"],
upr=predictdf[, "upr"]
)#Graficar usando ggplot2
library(ggplot2)
ggplot(results, aes(x=edad, y=fit))+
geom_point()+
geom_errorbar(aes(ymin=lwr, ymax=upr), width=0.1)+ #Intervalo de confianza
labs(x="Edad", y="Resistencia")+
ggtitle("Predicción de Resistecia vs Edad con intervalo de Confianza") 
Para obtener el IC del 90% para \(\hat{Y_0}\) cuando \(x_0=10\) semanas se usa el siguiente código:
nwdata2<- data.frame(edad=10)
predictdf2<-predict(object = mod1, newdata=nwdata2, interval = "prediction", level=0.90)
results2<-data.frame(
edad=nwdata2$edad,
fit=predictdf2[, "fit"], #la sintaxis [, "fit"] es para extraer columnas especificas del objeto predictions
lwr=predictdf2[, "lwr"],
upr=predictdf2[, "upr"]
)#Graficar usando ggplot2
library(ggplot2)
ggplot(results2, aes(x=edad, y=fit))+
geom_point()+
geom_errorbar(aes(ymin=lwr, ymax=upr), width=0.1)+ #Intervalo de confianza
labs(x="Edad", y="Resistencia")+
ggtitle("Predicción de Resistecia vs Edad con intervalo de Confianza") 
Observe que en el primer caso se usó interval=“confidence” mientras que en el segundo se usó interval=“prediction”
Ahora vamos a obtener todos los IC \(\hat{Y_0}\) y los vamos a almacenar en el objeto future_y que luego luego vamos a agregar al marco de datos original.
future_y<-predict(object=mod1, interval="prediction", level=0.95)Warning in predict.lm(object = mod1, interval = "prediction", level = 0.95): predictions on current data refer to _future_ responsesnewdat<- cbind(df, future_y)Con el codigo mostrado a continuacion se construye el diagrama de dispersión y se agrega la linea de regresión (en azu) y los IC para \(E(Y \mid x_0)\) (en rosado) por medio de geom_smooth. Los IC para \(\hat{Y_0}\) (en rojo) se agregan por medio de geom_line.
library(ggplot2)
ggplot(newdat, aes(x=edad, y=resistencia))+
geom_point()+
geom_line(aes(y=lwr), color="red", linetype="dashed")+
geom_line(aes(y=upr), color="red", linetype="dashed")+
geom_smooth(method=lm, formula=y~x, se=TRUE, level=0.95, col="blue", fill="pink2")+
theme_light()
De la figura anterior se observa claramente que los IC para \(\hat{Y_0}\) son siempre más anchos que los IC para \(E(Y\mid x_0)\)
Los puntos: indican los datos tomados del modelo (base de datos inicial) que se ubican en el eje x= edad y su relacion con el eje y=resistencia .
Linea azúl: Indica la regresión ajustada por el modelo entre “edad” y “resistencia”.
Banda rosa: Indica el intervalo de confianza del 95%, muestra el rango en el se espera caiga la media de “edad” y “resistencia”.
Linea discontinua: Indica el intervalo de predicción del 95% como vemos el rango desde la linea de regresión es mas amplia, es decir la presición del intervalo es menor que el del intervalo de confianza.
Relacion negativa: La grafica de regresión es negativa, lo que indica una relación inversa entre “edad” y “resistencia”, es decir que entre mas aumente la edad (en semanas) la resistencia va a disminuir.
\[ \hat{Resistencia}= \hat{\beta_0} - \hat{\beta_1}\,edad \]
\[ \hat{Resistencia}=2628-37.15 \, edad \]
Suponiendo que Edad=0 (semanas) la resistencia o estado inicial de la soldatura será de \(2638\)
\[ 2638-37.15(0)=2638 \]
Suponiendo que Edad=4 (semanas)
\[ 2638-37.15(4)=2415 \]
Suponiendo que Edad=15 (semanas)
\[ 2638-37.15(15)=2080.75 \]
Suponiendo que Edad=57 (semanas)
\[ 2638-37.15(57)=520.45 \]
Podemos observar como la resistencia de la soldadura disminuye a medida que aumenta la variable dependiente, comprobando así la relación negativa de nuestra regresión.