Línea de regresión y desviaciones en R: Ejemplo para evitar el uso de 'for'.

Recientemente en una página de Facebook sobre usuarios de R en Nicaragua me di cuenta de la existencia de un sencillo manual de estadística aplicada con utilización de R el cual fue escrito por un profesor de la Universidad Tecnológica Nacional de Argentina, este documento se encuentra disponible  en CRAN (The Comprehensive R Archive Network). Este es un documento muy sencillo y fácil de leer y proporciona muchos ejemplos sencillo de cómo usar R, sin embargo, la complejidad de algunas de sus implementaciones en R está sobredimensionada y hacen que sea difícil de entender cómo implementar el ejemplo en R.

Tras una rápida lectura al manual pude detectar un ejemplo cuya implementación en R me resultó ser más complicada de lo que realmente debería ser. En este post reproduzco el ejemplo del manual y  luego propongo una alternativa más natural y eliminando, principalmente, el uso del loop y eliminando asignaciones e indexaciones innecesarias.

El ejemplo que se desarrollará corresponde al de la página 43 del documento que puede encontrar aquí. Se trata simplemente de obtener la siguiente gráfica en R. (clic sobre la imagen para verla mejor)

A continuación se reproducen los códigos del manual (sólo he copiado y pegado los códigos del manual y he realizado una leve manipulación para eliminar variables que en este ejemplo no se usarán).

Edad <- c(22,22,23,24,25,25,26,27,28,29,29,29,29,29, 
          31,31,32,33,34,35,35,35,36,38,39,39,42,42,44,44,45,45, 
          45,47,48,52,59,66,67,69,69) 


SBR <- c(6.3,6.8,7,7.1,5.6,7.4,8.6,5.6,6.9,5.2,6.1,5.6,4.8, 6.2,7.4,
         5.1,6,7.5,6.1,6.8,7.9,6.8,4.9,6,4.9,6.8,5.3,5.8,4, 
         4.4,4.9,6.1,5.8,6.2,2.9,5.6,3.9,4.4,5.7,4.9,4.6) 


PAS.SBR <- data.frame(Edad,SBR)




## Página 43
#-----------
# Gráfico de la linea de regresión y las observaciones. Reproduciendo el ejemplo del manual
plot(PAS.SBR$Edad,PAS.SBR$SBR,
     xlim=c(20,70), 
     main="SBR versus Edad y residuos",
     xlab="Edad(años)",
     ylab="SBR (ms/mmHg)") 
abline(lm.SBR.Edad)  # agregando la linea de regresión


# Con esto obtiene los coeficientes del modelo 
interseccion <- lm.SBR.Edad$coefficients[1] 
pendiente <- lm.SBR.Edad$coefficients[2] 


# asigna que n tiene que ser el tamaño de la muestra para introducirla en el 'for'
n <- length(PAS.SBR$Edad) 


# a partir de este loop pinta cada línea punteada que indica las desviaciones de cada 
# observación respecto a la línea de regresión
for(i in 1:n){ 
  xx <- PAS.SBR$Edad[i] 
  y1 <- interseccion+pendiente*xx 
  y2 <- y1+lm.SBR.Edad$residuals[i] 
  lines(c(xx,xx),c(y1,y2),lty="dotted") 
} 




La alternativa


## Una alternativa que evita el for es utilizar la función 'segments' 
attach(PAS.SBR)


# Esto es lo nuevo: no hace falta un 'for' para pintar las desviaciones
xx <- Edad # vector de valores del regresor
y1 <- predict(lm.SBR.Edad) # valores predichos por el modelo 
y2 <- SBR # vector de valores de la variable respuesta
segments(x0=xx, y0=y1, y1=y2, lty=3)

Los principales cambios que se han hecho con esta propuesta son:

1.  xx <- PAS.SBR$Edad[i] se puede sacar del 'for' y hacer la asignación simplemente así: xx <- Edad
2. y1 <- interseccion+pendiente*xx fue reemplazado por y1 <- predict(lm.SBR.Edad), porque si lo que se quiere obtener son los valores predichos por el modelo, la función 'predict' lo hace por nosotros.
3.  y2 <- y1+lm.SBR.Edad$residuals[i], esto es simplemente el valor predicho más el error asociado a esa predicción por teoría ya se sabe que esto es simplemente el valor de la observación, de manera que reemplazo todo eso por el valor de la observación así y2 es simplemente y2 <- SBR ¿Mejor, no?
4. Todo el 'for' fue reemplazado por la función 'segments'. 

Comparando la estructura de mi propuesta con la del manual se ve la ventaja de la función 'segments' pues esta elimina la necesidad del 'for', por otro lado, se puede reducir aún más el código, para esto se tendrían que eliminar las asignaciones innecesarias, así el código final que sustituiría al 'for' es: 

segments(x0=Edad, y0=predict(lm.SBR.Edad), y1=SBR, lty=3)

Y así 8 líneas de código se reducen sólo a 1 simple línea.


Resultado final
El gráfico entonces se obtiene así:
  plot(Edad, SBR,
       xlim=c(20,70), 
       main="SBR versus Edad y residuos",
       xlab="Edad (años)",
       ylab="SBR (ms/mmHg)") 
abline(lm.SBR.Edad) # hasta aquí es igual que el manual
segments(x0=Edad, y0=predict(lm.SBR.Edad), y1=SBR, lty=3) # Esto elimina al for y a todas las asignaciones innecesarias.