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()
