Regresión Lineal

Datos y paquetes

• CR-1983.csv como clima
• ggplot2
• tmt2insectos.csv como dummy
• bt_model3.csv como btm

Los archivos se encuentran aquí:

https://sites.google.com/site/digitcognem/boxfile_pract

Regresion lineal Simple

El análisis de regresion con R es relativamente simple. Se utiliza el comando lm (de linear model) y se especifica el modelo, sin incluir el intercepto: lm(vardep~varpredictora,data), data es opcional. la funcion lm tiene numerosas opciones que se pueden consultar con help(lm).

Por ejemplo la regresion de tmax como variable dependiente (o variable respuesta) y la temperatura media como variable predictora (= variable independiente) es:

r1<-lm(tmax~tmedia,data=clima)

#ahora ver los resultados

summary(r1)

Call:

lm(formula = tmax ~ tmedia, data = clima)

Residuals:

Min 1Q Median 3Q Max

-6.1041 -1.2196 -0.1837 0.9676 9.9138

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 6.98293 1.03404 6.753 5.76e-11 ***

tmedia 0.96413 0.03826 25.200 < 2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.453 on 363 degrees of freedom

Multiple R-squared: 0.6363, Adjusted R-squared: 0.6353

F-statistic: 635 on 1 and 363 DF, p-value: < 2.2e-16

Se observa que la prueba de F para el modelo tiene un valor de p < 2.2e-16 (<0.05), por lo tanto si existe efecto significativo de la temperatura media en la temperatura máxima.

Para graficar se puede utilizar la funcion plot:

plot(clima$tmedia,clima$tmax,ylim=c(0,50),pch=8,cex=1.2)

Se puede sobreponer la linea de regresión con la funcion abline(r1):

abline(r1)

Para calcular los residuales estandarizados, se calcula lo siguiente:

residuales, son la diferencia entre observados (y) y predichos por el modelo.

obsv<-clima$tmax

pred<- 6.98293+0.96413*clima$tmedia

rsd<-obsv-pred

Los residuales estandarizados son los residuales entre el Error estándar de los residuales (2.453):

rsdn<-rsd/2.453

Lo anterior se puede realizar en una sola operación:

rsdn<-(obsv-pred)/2.453

Ahora se pueden graficar los residuales estandarizados contra los predichos:

plot(pred,rsdn,ylim=c(-5,5),pch=19,cex=1.2)

Una suposición de la regresión lineal es que los residuales tienen una distribución normal, así que un histograma de estos residuales sería:

hist(resdn,main="Histograma Residuales Normalizados",xlab="Residual",ylab="Frecuencia",col="orange")

Se observa que los residuales son relativamente normales, aunque presentan un sesgo hacia la derecha.

Regresion Lineal Múltiple

En este caso, la formulación del modelo es similar, con la inclusión de las diferentes variables predictoras en el modelo. Los datos del archivo tmt2insectos.csv son de un experimento con seis tratamientos que corresponden a diferentes tipos de Bacillus thuringiensis aplicados a una crucífera. La variable respuesta es el número de insectos. El conjunto de datos solamente contiene los primeros dos tratamientos.

head(dummy)

tmt tiempo insectos logInsec logtime

1 0 0 15.75 1.2240148 0.0000000

2 0 2 17.25 1.2612629 0.4771213

3 0 3 3.12 0.6148972 0.6020600

4 0 4 11.87 1.1095785 0.6989700

5 0 6 1.62 0.4183013 0.8450980

6 0 7 1.87 0.4578819 0.9030900

En el caso de la regresión lineal múltiple se pueden formular diferentes modelos y entonces la pregunta es: ¿Cuál modelo es el apropiado? Existen diferentes índices y pruebas para determinar esto. Aquí solamente se empleará el valor de r^2. También se pueden graficar los residuales como en el caso anterior.

El ajuste de varios modelos se ilustra a continuación:

m1<-lm(insectos~tmt+tiempo,data=dummy)

m2logtime<-lm(insectos~tmt+logtime,data=dummy)

m3logins<-lm(logInsec~tmt+tiempo,data=dummy)

m4loglog<-lm(logInsec~tmt+logtime,data=dummy)

El primer modelo es uno lineal, el segundo se utiliza la transformación logarítmica en tiempo; el tercero emplea la transformación logarítmica del conteo de insectos, el cuarto modelo es doble logarítmico.

Para ver los resultados, aplique la función summary(rg), donde rg es la variable que contiene los resultados de lm, es decir, de la regresión, como se hizo previamente.

summary(m3logins)

Call:

lm(formula = logInsec ~ tmt + tiempo, data = dummy)

Residuals:

Min 1Q Median 3Q Max

-0.39003 -0.15875 -0.00799 0.18523 0.43631

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 1.26083 0.12282 10.266 3.53e-08 ***

tmt 0.38344 0.11905 3.221 0.00571 **

tiempo -0.08530 0.01463 -5.829 3.32e-05 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2525 on 15 degrees of freedom

Multiple R-squared: 0.7473, Adjusted R-squared: 0.7136

F-statistic: 22.17 on 2 and 15 DF, p-value: 3.312e-05

Si se revisan los cuatro modelos y, basados en el valor de r^2, el mejor modelo es el tres. NOTA: el valor de r^2 no es el único indice y quizás tampoco el mejor para seleccionar modelos, aquí unicamente se usa con propósitos ilustrativos.

Gráfica de Regresión Múltiple

En este caso, el problema es graficar las dos series de puntos (cuando tmt=0 y cuando tmt=1) y sus respectivas lineas de regresión.

El modelo 3 de regresión es:

Y = 1.26 + 0.3834*tmt -0.0853*tiempo,

Cuando tmt=0:

Y1= 1.26 - 0.0853*tiempo,

Cuando tmt=1:

Y = 1.6 +0.3834(1) - 0.0853*tiempo

Y= 1.64427 - 0.0853*tiempo

Para graficar observados contra predichos (puntos=observados, lineas= modelos) es necesario construir un segundo archivo de datos que contenga los valores predichos por el modelo. Puesto que el modelo es una recta, son solamente necesarios dos puntos para cada tratamiento, dos cuando tmt=0 y dos cuando tmt=1. Estos puntos pueden ser tiempo= 0 y tiempo= 13

Aplicando los modelos previos se tiene:

Para tmt=0 y tiempo = 0:

y= 1.26 - 0.0853*(0) = 1.26

Para tmt= 0 y tiempo = 13:

y= 1.26 - 0.0853*(13) = 0.1511

Para tmt=1 y tiempo = 0:

Y= 1.64427 - 0.0853*(0) = 1.64427

Para tmt= 1 y tiempo = 13:

Y= 1.64427 - 0.0853**(13) = 0.53537

El archivo de datos debe concordar con los nombres de los encabezados del archivo previo (tmt2insectos.csv= dummy) de otra manera ggplot no sabrá que graficar. El contenido de este archivo (bt_model3.csv=btm) es:

tmt, tiempo,logInsec

0,0,1.26

0,13,0.1511

1,0,1.64427

1,13,0.53537

Ahora se procede a generar la información basica con los datos de dummy:

fs<-ggplot(dummy,aes(x=tiempo,y=logInsec,shape=factor(tmt)))

fs+geom_point(size=6)+theme_bw()+geom_line(data=btm,aes(x=tiempo,y=logInsec,group=factor(tmt)))

Observe que en la primera linea se toman los datos del primer frame=dummy, tambien se requiere que la forma de los símbolos: shape=factor(tmt) sea el tipo de tratamiento, se utiliza la función factor(tmt) para convertir los valores numéricos de 0 y 1 de tmt a valores nominales ("0" y "1")

En la segunda linea se solicita graficar los puntos con tamaño 6, del archivo dummy. Mientras que la linea se pide con los datos específicos del frame= btm. la opcion group=factor(tmt) indica a ggplot que considere los datos en grupos por categorías de tmt, convertidos a factor también. Los nombres de las columnas debe concordar, de lo contrario ggplot se confunde y no hace bien las gráfica.