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.

Regresión por polinomios locales.

En el paquete KernSmooth existen variedad de funciones para estimar modelos de regresión no paramétricos utilizando polinomios locales, así que me parece interesante exponer el uso de la función locpoly usando como ejemplo un ejercicio propuesto en una clase de modelización no paramétrica. El ejercicio fue tomado de los apuntes de la clase Modelos No Paramétricos del Máster en Estadística e Investigación Operativa impartido en la Universidad Politécnica de Cataluña.

Las instrucciones del ejercicio son las siguientes:
Haz la regresión de lgSpeed frente a Yr usando locpoly. Utiliza distintos valores de los parámetros bandwith, degree (grado del polinomio local ajustado) y drv (derivada estimada).

¿En qué años se produjo un aumento más rápido en la velocidad de los aviones fabricados? (Si no lo ves claro, prueba con bandwith = 7, degree = 1, drv = 1, o con bandwith = 10, degree = 2, drv = 1).

Para contestar la pregunta planteada se tiene que calcular la primera derivada de los polinomios locales y dibujarla para ver el comportamiento de las curvas y así determinar los periodos de mayor crecimiento.

La solución en R es:

# install.packages(c('KernSmooth', 'sm'))
library(KernSmooth)
library(sm)

Para obtener la despción de los datos sólo hace falta hacer lo siguiente:

provide.data(aircraft)

## Data file being loaded

head(aircraft)  # asi lucen los datos

##   Yr Period Power Span Length Weight Speed Range
## 1 14      1  82.0 12.8   7.60   1070   105   400
## 2 14      1  82.0 11.0   9.00    830   145   402
## 3 14      1 223.6 17.9  10.35   2200   135   500
## 4 15      1 164.0 14.5   9.80   1946   138   500
## 5 15      1 119.0 12.9   7.90   1190   140   400
## 6 15      1  74.5  7.5   6.30    653   177   350

Un poco de estimación….

lgSpeed <- log(Speed)
fit1 <- locpoly(Yr, lgSpeed, degree = 1, bandwidth = 7, drv = 1)
fit2 <- locpoly(Yr, lgSpeed, degree = 2, bandwidth = 10, drv = 1)

Graficando

plot(fit1, col=3, lwd=2, 
     type="l",  bty="l",
     main="Derivadas de orden 1 \n de los polinomios locales estimados",

     font.main=1, las=1,
     xlab="Year", xaxs="i",
     ylab="Derivadas de los polinomios locales estimados")

lines(fit2, col=6, lwd=2)

legend("topright", c("grado=1, bandwidth=7", "grado=2, bandwidth=10"),
       col=c(3,6), lty=1, fill= c(3,6), border= c(3,6),
       bty="n", cex=.9)

text(30, 0.023, "Periodo de mayor \n crecimiento en la \n velocidad de \n los aviones", cex=.95)

El primer polinomio sitúa el máximo nivel en incremento de velocidad de los aviones en el año 1934 mientras que el segundo polinomio lo sitúa en el año 1933, véase Figura. En términos de la variable Period esto corresponde al período 1.

1900 + floor(fit1$x[which.max(fit1$y)])  # año de mayor crecimiento según el primer polinomio

## [1] 1934

1900 + floor(fit2$x[which.max(fit2$y)])  # año de mayor crecimiento según el segundo polinomio

## [1] 1933

Unir varios data.frames en un sólo paso: 'merge' y 'Reduce'

Para unir data.frames en R existe la función merge que con todas sus opciones hace que esto sea una tarea fácil, sin embargo, se vuelve una tarea aburrida y repetitiva cuando se tienen muchos data.frames para unir, puesto que la función merge solo permite la unión de dos a la vez, véase la definición de los argumentos y la ayuda de dicha función. Pese a este inconveniente, la función merge puede ser utilizada para unir tantos data.frames como queramos si la combinamos apropiadamente con la función Reduce, el único requisito que se ha de tener en cuenta es que el nombre de la variable identificadora sea el mismo en todos los data.frames a ser unidos. Vemos unos ejemplos.

Consideremos los siguientes data.frames:

A <- data.frame(id = c("A", "B", "C", "D"), age = c(24, 25, 17, 19), height = c(1.8, 
    1.9, 1.75, 1.65))
A

##   id age height
## 1  A  24   1.80
## 2  B  25   1.90
## 3  C  17   1.75
## 4  D  19   1.65

B <- data.frame(gender = c("M", "M", "F", "F"), id = c("A", "B", "C", "D"))
B

##   gender id
## 1      M  A
## 2      M  B
## 3      F  C
## 4      F  D

C <- data.frame(id = c("A", "B", "C", "D"), math = c(6.5, 8.9, 7.4, 9.2), science = c(7.2, 
    8.4, 6.5, 8.7))
C

##   id math science
## 1  A  6.5     7.2
## 2  B  8.9     8.4
## 3  C  7.4     6.5
## 4  D  9.2     8.7

D <- data.frame(id = c("A", "B", "C", "D"), eyes = c("blue", "brown", "green", 
    "black"))
C

##   id math science
## 1  A  6.5     7.2
## 2  B  8.9     8.4
## 3  C  7.4     6.5
## 4  D  9.2     8.7

# Nótese que todos ellos tienen como variable identificadora id

## Uniendo los dataframes con merge
AB <- merge(A, B)  # une A con B
ABC <- merge(AB, C)  # a la unión de A y B le agrega C
ABCD <- merge(ABC, D)  # a la unión de A, B y C le agrega D
ABCD  # resultado final

##   id age height gender math science  eyes
## 1  A  24   1.80      M  6.5     7.2  blue
## 2  B  25   1.90      M  8.9     8.4 brown
## 3  C  17   1.75      F  7.4     6.5 green
## 4  D  19   1.65      F  9.2     8.7 black

Como se ha visto, se tiene que ir uniendo dos data.frames en cada paso. ¿Hay alguna manera de hacerlo en un sólo paso? La respuesta es Sí y sólo requiere combinar la función merge con Reduce y tener en cuenta que la variable identificadora en cada data.frame tenga el mismo nombre (id en este ejemplo), no hace falta que esta variable ocupe la misma posición en cada data.frame, sólo se requiere que sea nombrada de la misma manera. El ejemplo anterior se puede hacer tan sólo con una simple linea de comandos:

Reduce(merge, list(A, B, C, D))

##   id age height gender math science  eyes
## 1  A  24   1.80      M  6.5     7.2  blue
## 2  B  25   1.90      M  8.9     8.4 brown
## 3  C  17   1.75      F  7.4     6.5 green
## 4  D  19   1.65      F  9.2     8.7 black

Como se ve en ambos casos se obtiene el mismo resultado, sin embargo en el segundo, se ha hecho todo de una vez.

Otro ejemplo:

authors <- data.frame(surname = I(c("Tukey", "Venables", "Tierney", "Ripley", 
    "McNeil")), nationality = c("US", "Australia", "US", "UK", "Australia"), 
    deceased = c("yes", rep("no", 4)))

books <- data.frame(name = I(c("Tukey", "Venables", "Tierney", "Ripley", "Ripley", 
    "McNeil", "R Core")), title = c("Exploratory Data Analysis", "Modern Applied Statistics ...", 
    "LISP-STAT", "Spatial Statistics", "Stochastic Simulation", "Interactive Data Analysis", 
    "An Introduction to R"), other.author = c(NA, "Ripley", NA, NA, NA, NA, 
    "Venables & Smith"))

colnames(authors)[1] <- "name"  # cambiando el ID de authors para que sea igual para todos

edition <- data.frame(name = authors[, 1], edition = c(4, 2, 3, 1, 2))  # invento
year <- data.frame(name = authors[, 1], year = 2000:2004)  # invento

Los data.frames authors y books son los ejemplos que se encuentran en la ayuda de merge',editionyyear fueron inventos míos para ilustrar el ejemplo.

Uniendo los data.frames usando merge, dos a la vez.

m1 <- merge(authors, books)
m2 <- merge(m1, edition)
(m3 <- merge(m2, year))  # El resultado final es:

##       name nationality deceased                         title other.author
## 1   McNeil   Australia       no     Interactive Data Analysis         <NA>
## 2   Ripley          UK       no            Spatial Statistics         <NA>
## 3   Ripley          UK       no         Stochastic Simulation         <NA>
## 4  Tierney          US       no                     LISP-STAT         <NA>
## 5    Tukey          US      yes     Exploratory Data Analysis         <NA>
## 6 Venables   Australia       no Modern Applied Statistics ...       Ripley
##   edition year
## 1       2 2004
## 2       1 2003
## 3       1 2003
## 4       3 2002
## 5       4 2000
## 6       2 2001

Uniéndolos todos a la vez:

Reduce(merge, list(authors, books, edition, year))

##       name nationality deceased                         title other.author
## 1   McNeil   Australia       no     Interactive Data Analysis         <NA>
## 2   Ripley          UK       no            Spatial Statistics         <NA>
## 3   Ripley          UK       no         Stochastic Simulation         <NA>
## 4  Tierney          US       no                     LISP-STAT         <NA>
## 5    Tukey          US      yes     Exploratory Data Analysis         <NA>
## 6 Venables   Australia       no Modern Applied Statistics ...       Ripley
##   edition year
## 1       2 2004
## 2       1 2003
## 3       1 2003
## 4       3 2002
## 5       4 2000
## 6       2 2001

Como se ha visto, si se combina merge con Reduce se pueden unir tantos data.frames en un sólo paso, sin embargo, con el enfoque tradicional de usar sólo merge se tiene que repetir la operación \( K-1 \) veces donde \( K \) es el número de data.frames.