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, ...) :
...等
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")
