我使用 HH 库中的 ci.plot 为简单的线性回归问题制作了置信区间图
observations2 <- subset(observations, year > 1994 & year < 2049)
model.seaice <- lm(seaice ~ NHtemp, data=observations2)
ci.plot(model.seaice, conf.level=0.95)
而且效果非常好。但是,我想添加到情节中。首先,我想加入积分,以便我可以轻松地及时遵循他们的顺序。我试过了:
plot(observations2[,'NHtemp'], observations2[,'seaice'], add="TRUE", type="b")
但它不喜欢它:
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
任何想法如何向 ci.plot 添加功能?
更新
我最终只是运行:
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)
该图看起来像这样(注意虚线是预测间隔而不是 conf.int。这是我现在意识到我真正想要的):
我仍然不能做的是显示颜色图的图例(在 Matlab 中,这将是命令:colorbar)。这个想法是为了突出数据中的轻微滞后。从 1995 年到 2048 年,颜色从蓝色变为粉红色。
有人知道如何显示颜色条吗?
一些数据
根据此处的要求,是要使用的数据的子集:
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")