I've made a confidence interval plot to a simple linear regression problem using ci.plot from the HH library
observations2 <- subset(observations, year > 1994 & year < 2049)
model.seaice <- lm(seaice ~ NHtemp, data=observations2)
ci.plot(model.seaice, conf.level=0.95)
And it works really well. However, I want to add to the plot. Firstly, I'd like to join the points so I can easily follow their order in time. I tried:
plot(observations2[,'NHtemp'], observations2[,'seaice'], add="TRUE", type="b")
but it doesn't like it:
Warning messages:
1: In plot.window(...) : "add" is not a graphical parameter
2: In plot.xy(xy, type, ...) : "add" is not a graphical parameter
3: In axis(side = side, at = at, labels = labels, ...) :
...etc
Any ideas how to add features to a ci.plot?
UPDATE
I've ended up just running with:
model.seaice <- lm(seaice ~ NHtemp, data=observations2)
plot(observations2[,'NHtemp'],observations2[,'seaice'],type="l",
+ xlab=expression(paste(Delta, ~degree~C)),
+ ylab=expression(min.~~sea~~ice~~(italic(millions~~sqr.~~km))),
+ main="Temperature anomaly and minimum sea ice extent")
+ points(observations2[,'NHtemp'],observations2[,'seaice'],pch=21,
+ col=topo.colors(dim(observations2)[1]),bg=topo.colors(dim(observations2)[1]))
newdata = data.frame(NHtemp=seq(0,2,0.1))
conf.vals <- predict(model.seaice, newdata, interval="prediction")
lines(newdata[,'NHtemp'],conf.vals[,'fit'],col="red")
lines(newdata[,'NHtemp'],conf.vals[,'lwr'],col="red",lty=2)
lines(newdata[,'NHtemp'],conf.vals[,'upr'],col="red",lty=2)
The figure looks like this (note the dashed lines are the prediction interval rather than the conf. int. which is what I now realise I actually wanted):
What I still can't do is show a legend of the colormap (in Matlab this would be the command: colorbar). The idea is to highlight the slight hysteresis in the data. The color progresses from blue to pink as the years progress from 1995 to 2048.
Anyone know how to show the colorbar?
SOME DATA
As requested here is a subset of the data to play with:
observations2 <- structure(list(year = c(1995, 1996, 1997, 1998, 1999, 2000, 2001,
2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009), NHtemp = c(0.314632441576551,
0.186366426772796, 0.476974553197026, 0.60389558449144, 0.808296103505494,
0.960909839541133, 0.755078064632368, 0.698608832025563, 0.795763098867946,
0.929551971074, 0.978281142325487, 1.11821925610747, 1.23718650073483,
1.19223388599054, 0.966671466003422), SHtemp = c(0.210929776358558,
0.139817106972881, 0.21783040846463, 0.43431186700407, 0.370475749149358,
0.255331083239238, 0.385363664392244, 0.403967011006417, 0.473766310099707,
0.440515909854607, 0.69104778063082, 0.716519830684257, 0.678544576020567,
0.590446589601271, 0.430452631622458), global = c(0.262781108967546,
0.163091766872832, 0.347402480830821, 0.519103725747748, 0.589385926327418,
0.608120461390177, 0.570220864512299, 0.551287921515984, 0.634764704483819,
0.685033940464296, 0.834664461478146, 0.917369543395855, 0.957865538377691,
0.891340237795898, 0.698562048812933), seaice = c(-0.79633165604255,
-0.0484464471858699, -0.821451994832936, -0.761853932432746,
-1.21404527404276, -1.80462439656434, -2.33021932073994, -1.62220305218969,
-1.75453917196037, -2.32457807586636, -2.58044546640576, -2.93648428237656,
-2.24364361423478, -2.63830194428149, -2.53740259589355)), .Names = c("year",
"NHtemp", "SHtemp", "global", "seaice"), row.names = 137:151, class = "data.frame")