Regresión lineal e intervalo de confianza

Una cuestión habitual en bioinformática es la de buscar e interpretar posibles correlaciones entre variables de un conjunto de datos. Cuando la correlación entre dos variables es significativa a menudo es interesante modelar la relación entre ellas por medio de una función lineal de regresión. Como todo modelo, este tipo de funciones tienen asociado un grado de error o incertidumbre, que podrá o no ser tolerable según el uso que queramos darle.

En este ejemplo mostramos gráficamente el error de una función de regresión, partiendo de un conjunto de datos, que guardaremos en un fichero llamado datos_regresion.tab, con datos separados con tabuladores:

days y2000 y2010
1 33.61 20.58
5 31.96 23.70
9 40.42 23.30
15 35.43 33.12
18 46.13 31.13
22 39.19 36.30
25 42.72 33.97
29 46.72 31.26
31 42.62 32.33
36 45.09 34.10 
39 55.05 38.17
43 46.89 37.46
46 54.98 40.08
51 51.83 40.68
54 54.64 40.68
57 49.87 39.47
60 51.33 39.57

De las tres variables vamos a trabajar con las dos primeras, days, que cuenta los días desde el inicio de lexperimento, y y2000, que se corresponde con las mediciones de la variableX en el año 2000. El siguiente código en R muestra como calcular la correlación, el coeficiente de determinación (y su P-valor) y el intervalo de confianza al 95% de la recta de regresión:

1:  # prog4.R , inspirado por ejemplos de   
2:  # http://wwwmaths.anu.edu.au/~johnm/r-book/2edn/scripts/exploit.R  
3:    
4:  # función: grafica el intervalo de confianza 0.95 en torno a la regresión  
5:  curvasIC <- function(form=varX$y2010~varX$days,data=varX,lty=2,  
6:     col=3,lwd=2,newdata=data.frame(days=varX$days))  
7:  {   
8:    varX.lm <- lm(form,data=data)   
9:    predmod <- predict(varX.lm,newdata=newdata,interval="confidence")   
10:    lines(spline(newdata$days,predmod[,"fit"]),lwd=3)   
11:    lines(spline(newdata$days,predmod[,"lwr"]),lty=lty,lwd=lwd,col=col)   
12:    lines(spline(newdata$days,predmod[,"upr"]),lty=lty,lwd=lwd,col=col)   
13:  }  
14:    
15:  varX = read.table("datos_regresion.tab",header=T)  
16:    
17:  ## nos centramos en la serie 2000   
18:  lm2000 = lm(varX$y2000~varX$days) # calcula regresión lineal  
19:  sum = summary(lm2000)  
20:    
21:  ## imprime diagrama XY   
22:  title = sprintf("year 2000 (R^2=%.2f)",sum$r.squared)  
23:  plot(varX$days,varX$y2000,ylab="variableX [arbitrary units]",  
24:     xlab="(days from experiment start)",main=title )  
25:    
26:  ## imprime regresión, intervalos de confianza y función ajustada  
27:  curvasIC(form=varX$y2000~varX$days)  
28:  text(40,35,sprintf("variableX (days) = %1.2f + %1.2f days",  
29:     lm2000$coefficients[1],lm2000$coefficients[2]))  
30:  text(40,32,sprintf("P-value = %1.3g",sum$coefficients[2,4]))  
31:    
32:  ## crea copia en formato PDF  
33:  dev.copy(pdf, file="regresion2000.pdf")  
34:  dev.off()