# ------------------------------------------------------------------ # Analyse anchor regression simulations: # RMSE and SD increase in CMIP split experiment after distributional change # ------------------------------------------------------------------ # Sebastian Sippel # 31.05.2021 library(raster) library(ncdf4) library(hydroGOF) library(fields) library(matrixStats) library(bigmemory) library(weights) ## 0. Make figures/plots: ## ------------------------------------------------------------------------------------------ library(RColorBrewer) library(ggplot2) set1_all = brewer.pal(n = 8, name = "Set1")[-6] set3 = brewer.pal(n = 12, name = "Set3") set2 = brewer.pal(n = 8, name = "Set2") Dark2 = brewer.pal(n = 8, name = "Dark2") set1 = brewer.pal(n = 8, name = "Dark2")[c(8, 7, 2, 1)] set1 = brewer.pal(n = 8, name = "Set1")[c(1,4,5,2)] library(ggsci) library(scales) # show_col(pal_nejm(alpha = 1)(9)) train.col = pal_nejm(alpha = 1)(8)[2] train.col_ = pal_nejm(alpha = 0.5)(8)[2] test.col = pal_nejm(alpha = 1)(8)[3] test.col_ = pal_nejm(alpha = 0.9)(8)[3] ## Define test.ix and train.ix: train.ix = which(names(cv.scen_cs) %in% c(cv.scen_cs_INFO$train.modcl.un)) test.ix = which((names(cv.scen_cs) %in% c(cv.scen_cs_INFO$test.modcl.un))) train.piC.ix = which(names(cv.piControl_cs) %in% c(cv.piControl_cs_INFO$train.modcl.un)) test.piC.ix = which(names(cv.piControl_cs) %in% c(unique(substr(cv.piControl_cs_INFO$test.mod.phys, 1, 3)))) ### FIGURES 2C 2D: SD / RMSE INCREASE ## ------------------------------------------------------------------------------------------ # Fig 02 c/d: RMSE AND SD PERCENT INCREASE IN CMIP-SPLIT EXPERIMENT: { for (t.ix in 1:length(trend.lengths)) { AGMT = sapply(X = cv.scen_cs, FUN=function(x) x$AGMT$trends[[t.ix]]$test.rmse.2040) EOF = sapply(X = cv.scen_cs, FUN=function(x) x$EOF$trends[[t.ix]]$test.rmse.2040) ridge = sapply(X = cv.scen_cs, FUN=function(x) x$g1$trends[[t.ix]]$test.rmse.2040[ridge.ix[[t.ix+1]]]) anchor0.5 = sapply(X = cv.scen_cs, FUN=function(x) x[[pareto_opti_cs[[t.ix+1]][11,2]]]$trends[[t.ix]]$test.rmse.2040[pareto_opti_cs[[t.ix+1]][11,1]]) anchor0.1 = sapply(X = cv.scen_cs, FUN=function(x) x[[pareto_opti_cs[[t.ix+1]][3,2]]]$trends[[t.ix]]$test.rmse.2040[pareto_opti_cs[[t.ix+1]][3,1]]) AGMT.piC = sapply(X = cv.piControl_cs, FUN=function(x) x$AGMT$trends[[t.ix]]$test.sd) EOF.piC = sapply(X = cv.piControl_cs, FUN=function(x) x$EOF$trends[[t.ix]]$test.sd) ridge.piC = sapply(X = cv.piControl_cs, FUN=function(x) x$g1$trends[[t.ix]]$test.sd[ridge.ix[[t.ix+1]]]) anchor0.5.piC = sapply(X = cv.piControl_cs, FUN=function(x) x[[pareto_opti_cs[[t.ix+1]][11,2]]]$trends[[t.ix]]$test.sd[pareto_opti_cs[[t.ix+1]][11,1]]) anchor0.1.piC = sapply(X = cv.piControl_cs, FUN=function(x) x[[pareto_opti_cs[[t.ix+1]][3,2]]]$trends[[t.ix]]$test.sd[pareto_opti_cs[[t.ix+1]][3,1]]) # Index for : test.scen.ix = which(cv.scen.dist_cs$M.mod$y40$modcl %in% cv.scen_cs_INFO$test.modcl.un) train.scen.ix = which((cv.scen.dist_cs$M.mod$y40$modcl %in% cv.scen_cs_INFO$train.modcl.un)) # train.piC.ix; test.piC.ix test.mod.piC.ix = which( cv.piC.dist_cs$M.mod$y40$train.test == "test" ) # unique(cv.piC.dist_cs$M.mod$y40$modcl[which( cv.piC.dist_cs$M.mod$y40$train.test == "test")]) train.mod.piC.ix = which( (cv.piC.dist_cs$M.mod$y40$train.test == "train") ) # rmse increase: pdf(file = paste("_fig02/_fig_02c_rmse_change_cs_", trend.lengths[t.ix], "y.pdf", sep=""), width = 5.5, height=3.25) par(mfrow=c(1,1), mar = c(4, 5, 0.5, 0)) { plot(c(1,1), type="n", xlim = c(0, 0.01*trend.lengths[t.ix]), ylim = c(0,9), bty="n", yaxt = "n", ylab = "Detection metrics", # xlab = bquote("Root mean square error of 40-year trend estimates [°C "*"("*.(trend.lengths[t.ix])*"y)"^-1*"]") xlab = bquote("Forced trend reconstruction RMSE [°C per 40 years]")) axis(side = 2, at = c(8.75, 6.75, 4.75, 2.75, 0.75), labels = c("GMT", "MWP", "Ridge", "Anchor \n (0.5/0.5)", "Anchor \n (0.1/0.9)"), tick=F, las=2, cex = 1, line = -2) # GMT: points(x = cv.scen.dist_cs$Y.mod[[t.ix+1]]$AGMT[train.scen.ix] * trend.lengths[t.ix], y = rep(9, length(train.scen.ix)), col = train.col, pch = "|", cex = 0.8) points(x = cv.scen.dist_cs$Y.mod[[t.ix+1]]$AGMT[test.scen.ix] * trend.lengths[t.ix], y = rep(8.5, length(test.scen.ix)), col = test.col, pch = "|", cex = 0.8) points(x = mean(AGMT[train.ix]) * trend.lengths[t.ix], y = 8.75, col = train.col, pch = 21, cex = 1, bg = train.col_) points(x = mean(AGMT[test.ix]) * trend.lengths[t.ix], y = 8.75, col = test.col, pch = 21, cex = 1, bg = test.col_) arrows(x0 = mean(AGMT[train.ix]) * trend.lengths[t.ix], y0 = 8.75, x1 = mean(AGMT[test.ix]) * trend.lengths[t.ix], y1 = 8.75, length=0.1, col = "grey25") text(x = mean(AGMT[test.ix]) * trend.lengths[t.ix], y = 9, labels = paste("+", round((mean(AGMT[test.ix]) / mean(AGMT[train.ix]) - 1)*100, 0), "%", sep=""), pos = 4, cex = 1) # EOF: points(x = cv.scen.dist_cs$Y.mod[[t.ix+1]]$EOF[train.scen.ix] * trend.lengths[t.ix], y = rep(7, length(train.scen.ix)), col = train.col, pch = "|", cex = 0.8) points(x = cv.scen.dist_cs$Y.mod[[t.ix+1]]$EOF[test.scen.ix] * trend.lengths[t.ix], y = rep(6.5, length(test.scen.ix)), col = test.col, pch = "|", cex = 0.8) # points(x = EOF[train.cs.ix] * trend.lengths[t.ix], y = rep(7, length(train.cs.ix)), col = train.col, pch = "|", cex = 0.8) # points(x = EOF[test.cs.ix] * trend.lengths[t.ix], y = rep(6.5, length(test.cs.ix)), col = test.col, pch = "|", cex = 0.8) points(x = mean(EOF[train.ix]) * trend.lengths[t.ix], y = 6.75, col = train.col, pch = 22, cex = 1, bg = train.col_) points(x = mean(EOF[test.ix]) * trend.lengths[t.ix], y = 6.75, col = test.col, pch = 22, cex = 1, bg = test.col_) arrows(x0 = mean(EOF[train.ix]) * trend.lengths[t.ix], y0 = 6.75, x1 = mean(EOF[test.ix]) * trend.lengths[t.ix], y1 = 6.75, length=0.1, col = "grey25") text(x = mean(EOF[test.ix]) * trend.lengths[t.ix], y = 7, labels = paste("+", round((mean(EOF[test.ix]) / mean(EOF[train.ix]) - 1)*100, 0), "%", sep=""), pos = 4, cex = 1) # Ridge: points(x = cv.scen.dist_cs$Y.mod[[t.ix+1]]$ridge[train.scen.ix] * trend.lengths[t.ix], y = rep(5, length(train.scen.ix)), col = train.col, pch = "|", cex = 0.8) points(x = cv.scen.dist_cs$Y.mod[[t.ix+1]]$ridge[test.scen.ix] * trend.lengths[t.ix], y = rep(4.5, length(test.scen.ix)), col = test.col, pch = "|", cex = 0.8) # points(x = ridge[train.cs.ix] * trend.lengths[t.ix], y = rep(5, length(train.cs.ix)), col = train.col, pch = "|", cex = 0.8) # points(x = ridge[test.cs.ix] * trend.lengths[t.ix], y = rep(4.5, length(test.cs.ix)), col = test.col, pch = "|", cex = 0.8) points(x = mean(ridge[train.ix]) * trend.lengths[t.ix], y = 4.75, col = train.col, pch = 24, cex = 1, bg = train.col_) points(x = mean(ridge[test.ix]) * trend.lengths[t.ix], y = 4.75, col = test.col, pch = 24, cex = 1, bg = test.col_) arrows(x0 = mean(ridge[train.ix]) * trend.lengths[t.ix], y0 = 4.75, x1 = mean(ridge[test.ix]) * trend.lengths[t.ix], y1 = 4.75, length=0.1, col = "grey25") text(x = mean(ridge[test.ix]) * trend.lengths[t.ix], y = 5, labels = paste("+", round((mean(ridge[test.ix]) / mean(ridge[train.ix]) - 1)*100, 0), "%", sep=""), pos = 4, cex = 1) # anchor0.5: points(x = cv.scen.dist_cs$Y.mod[[t.ix+1]]$opti_0.5[train.scen.ix] * trend.lengths[t.ix], y = rep(3, length(train.scen.ix)), col = train.col, pch = "|", cex = 0.8) points(x = cv.scen.dist_cs$Y.mod[[t.ix+1]]$opti_0.5[test.scen.ix] * trend.lengths[t.ix], y = rep(2.5, length(test.scen.ix)), col = test.col, pch = "|", cex = 0.8) # points(x = anchor0.5[train.cs.ix] * trend.lengths[t.ix], y = rep(3, length(train.cs.ix)), col = train.col, pch = "|", cex = 0.8) # points(x = anchor0.5[test.cs.ix] * trend.lengths[t.ix], y = rep(2.5, length(test.cs.ix)), col = test.col, pch = "|", cex = 0.8) points(x = mean(anchor0.5[train.ix]) * trend.lengths[t.ix], y = 2.75, col = train.col, pch = 23, cex = 1, bg = train.col_) points(x = mean(anchor0.5[test.ix]) * trend.lengths[t.ix], y = 2.75, col = test.col, pch = 23, cex = 1, bg = test.col_) arrows(x0 = mean(anchor0.5[train.ix]) * trend.lengths[t.ix], y0 = 2.75, x1 = mean(anchor0.5[test.ix]) * trend.lengths[t.ix], y1 = 2.75, length=0.1, col = "grey25") text(x = mean(anchor0.5[test.ix]) * trend.lengths[t.ix], y = 3, labels = paste("+", round((mean(anchor0.5[test.ix]) / mean(anchor0.5[train.ix]) - 1)*100, 0), "%", sep=""), pos = 4, cex = 1) # anchor 0.1: points(x = cv.scen.dist_cs$Y.mod[[t.ix+1]]$opti_0.1[train.scen.ix] * trend.lengths[t.ix], y = rep(1, length(train.scen.ix)), col = train.col, pch = "|", cex = 0.8) points(x = cv.scen.dist_cs$Y.mod[[t.ix+1]]$opti_0.1[test.scen.ix] * trend.lengths[t.ix], y = rep(0.5, length(test.scen.ix)), col = test.col, pch = "|", cex = 0.8) # points(x = anchor0.1[train.cs.ix] * trend.lengths[t.ix], y = rep(1, length(train.cs.ix)), col = train.col, pch = "|", cex = 0.8) # points(x = anchor0.1[test.cs.ix] * trend.lengths[t.ix], y = rep(0.5, length(test.cs.ix)), col = test.col, pch = "|", cex = 0.8) points(x = mean(anchor0.1[train.ix]) * trend.lengths[t.ix], y = 0.75, col = train.col, pch = 23, cex = 1.5, bg = train.col_) points(x = mean(anchor0.1[test.ix]) * trend.lengths[t.ix], y = 0.75, col = test.col, pch = 23, cex = 1.5, bg = test.col_) arrows(x0 = mean(anchor0.1[train.ix]) * trend.lengths[t.ix], y0 = 0.75, x1 = mean(anchor0.1[test.ix]) * trend.lengths[t.ix], y1 = 0.75, length=0.1, col = "grey25") text(x = mean(anchor0.1[test.ix]) * trend.lengths[t.ix], y = 1, labels = paste("+", round((mean(anchor0.1[test.ix]) / mean(anchor0.1[train.ix]) - 1)*100, 0), "%", sep=""), pos = 4, cex = 1) dev.off() } pdf(file = paste("_fig02/_fig_02d_sd_change_cs_", trend.lengths[t.ix],"y.pdf", sep=""), width = 5.5, height=3.25) par(mfrow=c(1,1), mar = c(4, 5, 0.5, 0)) { plot(c(1,1), type="n", xlim = c(0, 0.0125*trend.lengths[t.ix]), ylim = c(0,9), bty="n", yaxt = "n", ylab = "Detection metrics", xlab = "") axis(side = 2, at = c(8.75, 6.75, 4.75, 2.75, 0.75), labels = c("GMT", "MWP", "Ridge", "Anchor \n (0.5/0.5)", "Anchor \n (0.1/0.9)"), tick=F, las=2, cex = 1, line = -2) # mtext(text = bquote("Standard deviation of 40-year trends [°C "*"("*.(trend.lengths[t.ix])*"y)"^-1*"]"), side = 1, line = 2.2) mtext(text = bquote("Standard deviation of 40-year trends [°C per 40 years]"), side = 1, line = 2.2) mtext(text = "in pre-industrial control simulations", side = 1, line = 3) # legend("topright", "pre-industrial") # AGMT: points(x = cv.piC.dist_cs$Y.mod[[t.ix+1]]$AGMT[train.mod.piC.ix] * trend.lengths[t.ix], y = rep(9, length(train.mod.piC.ix)), col = train.col, pch = "|", cex = 0.8) points(x = cv.piC.dist_cs$Y.mod[[t.ix+1]]$AGMT[test.mod.piC.ix] * trend.lengths[t.ix], y = rep(8.5, length(test.mod.piC.ix)), col = test.col, pch = "|", cex = 0.8) # points(x = AGMT.piC[test.piC.ix] * trend.lengths[t.ix], y = rep(8.5, length(test.piC.ix)), col = test.col, pch = "|", cex = 0.8) points(x = mean(AGMT.piC[train.piC.ix]) * trend.lengths[t.ix], y = 8.75, col = train.col, pch = 21, cex = 1, bg = train.col_) points(x = mean(AGMT.piC[test.piC.ix]) * trend.lengths[t.ix], y = 8.75, col = test.col, pch = 21, cex = 1, bg = test.col_) arrows(x0 = mean(AGMT.piC[train.piC.ix]) * trend.lengths[t.ix], y0 = 8.75, x1 = mean(AGMT.piC[test.piC.ix]) * trend.lengths[t.ix], y1 = 8.75, length=0.1, col = "grey25") text(x = mean(AGMT.piC[test.piC.ix]) * trend.lengths[t.ix], y = 9, labels = paste("+", round((mean(AGMT.piC[test.piC.ix]) / mean(AGMT.piC[train.piC.ix]) - 1)*100, 0), "%", sep=""), pos = 4, cex = 1) # EOF: points(x = cv.piC.dist_cs$Y.mod[[t.ix+1]]$EOF[train.mod.piC.ix] * trend.lengths[t.ix], y = rep(7, length(train.mod.piC.ix)), col = train.col, pch = "|", cex = 0.8) points(x = cv.piC.dist_cs$Y.mod[[t.ix+1]]$EOF[test.mod.piC.ix] * trend.lengths[t.ix], y = rep(6.5, length(test.mod.piC.ix)), col = test.col, pch = "|", cex = 0.8) # points(x = EOF.piC[train.piC.ix] * trend.lengths[t.ix], y = rep(7, length(train.piC.ix)), col = train.col, pch = "|", cex = 0.8) # points(x = EOF.piC[test.piC.ix] * trend.lengths[t.ix], y = rep(6.5, length(test.piC.ix)), col = test.col, pch = "|", cex = 0.8) points(x = mean(EOF.piC[train.piC.ix]) * trend.lengths[t.ix], y = 6.75, col = train.col, pch = 22, cex = 1, bg = train.col_) points(x = mean(EOF.piC[test.piC.ix]) * trend.lengths[t.ix], y = 6.75, col = test.col, pch = 22, cex = 1, bg = test.col_) arrows(x0 = mean(EOF.piC[train.piC.ix]) * trend.lengths[t.ix], y0 = 6.75, x1 = mean(EOF.piC[test.piC.ix]) * trend.lengths[t.ix], y1 = 6.75, length=0.1, col = "grey25") text(x = mean(EOF.piC[test.piC.ix]) * trend.lengths[t.ix], y = 7, labels = paste("+", round((mean(EOF.piC[test.piC.ix]) / mean(EOF.piC[train.piC.ix]) - 1)*100, 0), "%", sep=""), pos = 4, cex = 1) # ridge.piC: points(x = cv.piC.dist_cs$Y.mod[[t.ix+1]]$ridge[train.mod.piC.ix] * trend.lengths[t.ix], y = rep(5, length(train.mod.piC.ix)), col = train.col, pch = "|", cex = 0.8) points(x = cv.piC.dist_cs$Y.mod[[t.ix+1]]$ridge[test.mod.piC.ix] * trend.lengths[t.ix], y = rep(4.5, length(test.mod.piC.ix)), col = test.col, pch = "|", cex = 0.8) # points(x = ridge.piC[train.piC.ix] * trend.lengths[t.ix], y = rep(5, length(train.piC.ix)), col = train.col, pch = "|", cex = 0.8) # points(x = ridge.piC[test.piC.ix] * trend.lengths[t.ix], y = rep(4.5, length(test.piC.ix)), col = test.col, pch = "|", cex = 0.8) points(x = mean(ridge.piC[train.piC.ix]) * trend.lengths[t.ix], y = 4.75, col = train.col, pch = 24, cex = 1, bg = train.col_) points(x = mean(ridge.piC[test.piC.ix]) * trend.lengths[t.ix], y = 4.75, col = test.col, pch = 24, cex = 1, bg = test.col_) arrows(x0 = mean(ridge.piC[train.piC.ix]) * trend.lengths[t.ix], y0 = 4.75, x1 = mean(ridge.piC[test.piC.ix]) * trend.lengths[t.ix], y1 = 4.75, length=0.1, col = "grey25") text(x = mean(ridge.piC[test.piC.ix]) * trend.lengths[t.ix], y = 5, labels = paste("+", round((mean(ridge.piC[test.piC.ix]) / mean(ridge.piC[train.piC.ix]) - 1)*100, 0), "%", sep=""), pos = 4, cex = 1) # anchor0.5.piC: points(x = cv.piC.dist_cs$Y.mod[[t.ix+1]]$opti_0.5[train.mod.piC.ix] * trend.lengths[t.ix], y = rep(3, length(train.mod.piC.ix)), col = train.col, pch = "|", cex = 0.8) points(x = cv.piC.dist_cs$Y.mod[[t.ix+1]]$opti_0.5[test.mod.piC.ix] * trend.lengths[t.ix], y = rep(2.5, length(test.mod.piC.ix)), col = test.col, pch = "|", cex = 0.8) # points(x = anchor0.5.piC[train.piC.ix] * trend.lengths[t.ix], y = rep(3, length(train.piC.ix)), col = train.col, pch = "|", cex = 0.8) # points(x = anchor0.5.piC[test.piC.ix] * trend.lengths[t.ix], y = rep(2.5, length(test.piC.ix)), col = test.col, pch = "|", cex = 0.8) points(x = mean(anchor0.5.piC[train.piC.ix]) * trend.lengths[t.ix], y = 2.75, col = train.col, pch = 23, cex = 1, bg = train.col_) points(x = mean(anchor0.5.piC[test.piC.ix]) * trend.lengths[t.ix], y = 2.75, col = test.col, pch = 23, cex = 1, bg = test.col_) arrows(x0 = mean(anchor0.5.piC[train.piC.ix]) * trend.lengths[t.ix], y0 = 2.75, x1 = mean(anchor0.5.piC[test.piC.ix]) * trend.lengths[t.ix], y1 = 2.75, length=0.1, col = "grey25") text(x = mean(anchor0.5.piC[test.piC.ix]) * trend.lengths[t.ix] + 0.01, y = 3, labels = paste("+", round((mean(anchor0.5.piC[test.piC.ix]) / mean(anchor0.5.piC[train.piC.ix]) - 1)*100, 0), "%", sep=""), pos = 4, cex = 1) # anchor 0.1: points(x = cv.piC.dist_cs$Y.mod[[t.ix+1]]$opti_0.1[train.mod.piC.ix] * trend.lengths[t.ix], y = rep(1, length(train.mod.piC.ix)), col = train.col, pch = "|", cex = 0.8) points(x = cv.piC.dist_cs$Y.mod[[t.ix+1]]$opti_0.1[test.mod.piC.ix] * trend.lengths[t.ix], y = rep(0.5, length(test.mod.piC.ix)), col = test.col, pch = "|", cex = 0.8) # points(x = anchor0.1.piC[train.piC.ix] * trend.lengths[t.ix], y = rep(1, length(train.piC.ix)), col = train.col, pch = "|", cex = 0.8) # points(x = anchor0.1.piC[test.piC.ix] * trend.lengths[t.ix], y = rep(0.5, length(test.piC.ix)), col = test.col, pch = "|", cex = 0.8) points(x = mean(anchor0.1.piC[train.piC.ix]) * trend.lengths[t.ix], y = 0.75, col = train.col, pch = 23, cex = 1, bg = train.col_) points(x = mean(anchor0.1.piC[test.piC.ix]) * trend.lengths[t.ix], y = 0.75, col = test.col, pch = 23, cex = 1, bg = test.col_) arrows(x0 = mean(anchor0.1.piC[train.piC.ix]) * trend.lengths[t.ix], y0 = 0.75, x1 = mean(anchor0.1.piC[test.piC.ix]) * trend.lengths[t.ix], y1 = 0.75, length=0.1, col = "grey25") text(x = mean(anchor0.1.piC[test.piC.ix]) * trend.lengths[t.ix] + 0.01, y = 1, labels = paste("+", round((mean(anchor0.1.piC[test.piC.ix]) / mean(anchor0.1.piC[train.piC.ix]) - 1)*100, 0), "%", sep=""), pos = 4, cex = 1) # legend("bottomright", c("Training models", "Test models"), col = c(train.col, test.col), pch = "|", inset = 0.02) legend("bottomright", c("Training models", "Test models"), text.col = c(train.col, test.col), inset = 0.02, bty="n") dev.off() } } }