Regresión lineal simple

En esta sección trabajaremos con los datos contenidos en el data.frame trees. La siguiente sucesión de órdenes nos permiten cargarlos, conocer su contenido y ponerlo en ruta de búsqueda.
data(trees)
str(trees)
?trees
attach(trees)
Vemos que consiste en una serie de datos de altura, diámetro y volumen de 31 árboles. Como en todo análisis de regresión simple, intentaremos analizar las relaciones existentes entre las variables con la intención de predecir valores de un de ellas en función de otra. Antes de comenzar, haremos un análisis preliminar de las correlaciones bivariadas. Con
cor(trees)
obtendremos la siguiente matriz de correlaciones

 		   Girth 		   Height 		  Volume


Girth 		  1.0000000 		 0.5192801 		 0.9671194


Height 		 0.5192801 		 1.0000000 		 0.5982497


Volume 		 0.9671194 		 0.5982497 		 1.0000000

Con
pairs(trees)
obtendremos la Figura .

Figura: Relación entre las variables de trees

Tanto en la matriz como en la figura se puede apreciar que hay una relación lineal muy fuerte entre el volumen del árbol y su diámetro. También se observa que la relación entre el volumen y la altura no lo es tanto. Finalmente, hay también relación entre el diámetro y la altura.

En R, para efectuar el análisis de regresión linear se utiliza la función lm(). Esta genera un tipo especial de objeto que puede ser argumento de muchas funciones que podremos utilizar para obtener información sobre el modelo de regresión calculado. Por ese motivo, es conveniente asignar el resultado de lm() a algún objeto con su correspondiente nombre. En nuestro caso, para hacer la regresión simple entre volumen y diámetro teclearemos
lm.vol.dia<-lm(Volume~Girth)
lo que almacena en el objeto lm.vol.dia todos los resultados del análisis de regresión.
Como puede observarse, lm() utiliza el sistema de fórmulas (con la tilde) para expresar la variable dependiente y la independiente.
La función summary() aplicada al objeto de regresión, nos va a proporcionar la mayor parte de la información que nos interesa en un principio como vemos a continuación
summary(lm.vol.dia)
Call:
lm(formula = Volume  Girth)
Residuals:
Min 1Q Median 3Q Max
-8.0654 -3.1067 0.1520 3.4948 9.5868
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -36.9435 3.3651 -10.98 7.62e-12 ***
Girth 5.0659 0.2474 20.48 < 2e-16 ***
--
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 4.252 on 29 degrees of freedom
Multiple R-Squared: 0.9353, Adjusted R-squared: 0.9331
F-statistic: 419.4 on 1 and 29 DF, p-value: 0
Tanto el valor como el grado de significación los coeficientes de la recta de regresión aparecen en el apartado coefficients. En la penúltima linea aparecen los valores de $r$ y $r^2$; y, en la última, la probabilidad de que el modelo obtenido sea inadecuado.
En este caso, el coeficiente de determinación ($r^2$) nos indica que el 93% de la variación en el volumen del árbol puede explicarse por variaciones en el diámetro del mismo. Por otro lado, el modelo de regresión obtenido es apropiado puesto que la probabilidad de que no lo sea es cero.

Para hacer predicciones del volumen de una serie de árboles a partir de su diámetro debemos crear un data.frame con una columna que almacene los diámetros de los mismos. Por ejemplo, si tenemos tres ejemplares, A,B,C con diámetros 11.5,10.6,18.7, para obtener los volúmenes predichos habría que proceder como sigue.
muestra<-data.frame(Girth=c(11.5,10.6,18.7),

row.names=c(``A'';''B'';''C''))
predict(lm.vol.dia,muestra)
Con la primera instrucción creamos el data.frame muestra que contiene una columna con el nombre de la variable independiente (en este caso el diámetro) que alberga los valores de la muestra. Además se indica la etiqueta de cada uno de los registros (A,B,C).
Con la segunda instrucción pedimos que calcule los valores predichos según el modelo guardado en lm.vol.dia. El resultado es

A 		   B 		    C   


21.31389 		16.75462 		 57.78806

Si, además del valor predicho, deseamos conocer el intervalo de confianza a un nivel determinado (digamos del 95%), deberemos teclear
predict(lm.vol.dia,muestra,level=.95,interval=``confidence'')
el resultado será

 		 fit 		    lwr 		   upr                                      


A 		21.31389		 19.51889 		23.10889                         


B 		16.75462 		14.69672 		18.81252                            


C 		57.78806		 54.61832 		60.95779 

donde podemos conocer el valor inferior (lwr) y superior (upr) para ese nivel de confianza. Para modificar el nivel de confianza, sólo hay que introducir otro valor en el argumento level de la función predict().

Finalmente, vamos a ver cómo incorporar al diagrama de dispersión algunos elementos importantes en análisis de regresión.
El primero es la recta de regresión. Para ello vamos a utilizar la función abline() que, pasada por consola tras haber generado un gráfico, le incorpora una recta con los parámetros que se le indiquen como argumentos. Veamos cómo hacerlo si queremos estudiar la relación entre el volumen y la altura.
plot(Height,Volume)
abline(lm(Volume~Height))
El resultado se presenta en la Figura .

Figura: Inclusión de la recta de regresión con abline()

En muchas ocasiones sería conveniente identificar ciertos puntos del diagrama de dispersión que destacan por alguna circunstancia (en general, valores anómalos). Para hacer esto podemos utlizar la función identify() de la siguiente manera
plot(Height,Volume)
identify(Height,Volume)
Tras haber realizado el gráfico con plot(), tecleamos la función identify() aplicada al mismo. Esto es así porque la función interactúa con el dispositivo gráfico. Tras teclearla, deberemos utilizar el ratón para marcar con el botón izquierdo los puntos que deseamos etiquetar. Por defecto se añadirá al gráfico un número índice, pero esto puede modificarse si existe alguna otra etiqueta. Cuando hayamos concluido, deberemos hacer clic con el botón derecho sobre el gráfico para cerrar la función. En ese momento aparacerán sobre la consola los índices de los puntos seleccionados. En la Figura  puede apreciarse cómo se ha seleccionado y etiquetado algunos de los casos.

Figura: Identificación de casos con identify()