######################################################################################################## # # Expanding windows of linear regression # # Kilian Eichenseer # # Linear regressions seperately in the Mesozoic-Cenozoic and in the Palaeozoic # library(nlme) # has to be installed # read the ASI and SCORara data asi_scor <- read.csv("C://Users/keichenseer/OneDrive - University of Plymouth/Manuscripts/Nature_Geoscience/Figshare_Uploads/CSVs/ASI_and_SCORara.csv") asi <- asi_scor$ASI scor <- asi_scor$SCORara stagetab <- readRDS("stagetable.rds") stage_base <- stagetab$stage_base stage_top <- stagetab$stage_top stage_midp <- (stage_base+stage_top)/2 ######################################################################################################## ####### First in the Mesozoic #### Calculate OLS linear models in expanding windows for the Mesozoic start <- 39 end <- 44 pars <- matrix(NA,nrow = 42, ncol = 3) colnames(pars) <- c("slope", "R2", "p") head(pars) for(i in 1:42) { end <- 43 + i ols <- summary(lm(scor[start:end]~asi[start:end])) pars[i,1] <- coef(ols)[2] pars[i,2] <- ols$r.squared pars[i,3] <- ols$coefficients[2,4] } ##### Calculate AR 1 models pars1 <- matrix(NA,nrow = 42, ncol = 6) colnames(pars1) <- c("slope", "R2", "p", "phi", "lratio", "lp") for(i in 1:42) { end <- 43 + i scor1 <- scor[start:end] asi1 <- asi[start:end] index <- start:end glsm <- gls(scor1~asi1, method="ML", correlation=corAR1(form = ~ index)) gls_null <- gls(scor1 ~ 1,correlation=corAR1(form = ~ index), method="ML") m <- glsm ll_gls <- logLik(m) ll_gls_null <- logLik(gls_null) ll_R2 <- 1-exp(-(2/length(index))*(ll_gls[1]-ll_gls_null[1])) pars1[i,1] <- coef(glsm)[2] pars1[i,2] <- ll_R2 pars1[i,3] <- summary(m)$tTable[2,4] phi <- m$modelStruct$corStruct phi <- coef(phi,unconstrained = FALSE) pars1[i,4] <- phi gls_ar0 <- gls(scor1~asi1, method="ML", correlation=NULL) test <- anova(gls_ar0,m) pars1[i,5] <- test$L.Ratio[2] pars1[i,6] <- test$`p-value`[2] } ##### from 11 (38 + 11 = 49) onwards use AR 1 model as it is subst. better than AR 0 pars_m <- rbind(pars[1:10,],pars1[11:42, 1:3]) ######################################################################################################## ####### Now in the Palaeozoic #### Calculate OLS linear models in expanding windows for the Palaeozoic start <- 1 end <- 6 pars <- matrix(NA,nrow = 33, ncol = 3) colnames(pars) <- c("slope", "R2", "p") head(pars) for(i in 1:33) { end <- 5 + i ols <- summary(lm(scor[start:end]~asi[start:end])) pars[i,1] <- coef(ols)[2] pars[i,2] <- ols$r.squared pars[i,3] <- ols$coefficients[2,4] } ##### Calculate AR 1 models pars1 <- matrix(NA,nrow = 33, ncol = 6) colnames(pars1) <- c("slope", "R2", "p", "phi", "lratio", "lp") for(i in 1:33) { end <- 5 + i scor1 <- scor[start:end] asi1 <- asi[start:end] index <- start:end glsm <- gls(scor1~asi1, method="ML", correlation=corAR1(form = ~ index)) gls_null <- gls(scor1 ~ 1,correlation=corAR1(form = ~ index), method="ML") m <- glsm ll_gls <- logLik(m) ll_gls_null <- logLik(gls_null) ll_R2 <- 1-exp(-(2/length(index))*(ll_gls[1]-ll_gls_null[1])) pars1[i,1] <- coef(glsm)[2] pars1[i,2] <- ll_R2 pars1[i,3] <- summary(m)$tTable[2,4] phi <- m$modelStruct$corStruct phi <- coef(phi,unconstrained = FALSE) pars1[i,4] <- phi gls_ar0 <- gls(scor1~asi1, method="ML", correlation=NULL) test <- anova(gls_ar0,m) pars1[i,5] <- test$L.Ratio[2] pars1[i,6] <- test$`p-value`[2] } ##### from 26 onwards use AR 1 model as it is subst. better than AR 0 pars_p <- rbind(pars[1:25,],pars1[26:33, 1:3]) ################################################################################################### ################################################################################################### ##### ##### PLOT - Figure 2b ##### ################################################################################################### ################################################################################################### windows( width = 6, height = 3.8) par(mfrow = c(1,1)) par(mar = c(5.5,4,1,1)) labelsize <- 1 title <- NA ylabel <- "" ylimit <- c(-0.305,0.8) # minimum and maximum of y-axis xlabfactor <- 3 y1 <- ylimit[1] # lower boundary of plot area y2 <- ylimit[1]-0.11*(ylimit[2]-ylimit[1]) # position of x axis ylabpush = -68 plot(0,0,type = "n", axes = FALSE, ylim = ylimit, xlim = c(-stagetab$stage_base[1],0), ylab = NA, xaxs = "i", yaxs = "i", yaxt = "n", xaxt = "n", xlab = NA) periodbase <- c(-485.40, -443.80, -419.20, -358.90, -298.90, -252.17, -201.30, -145.00, -66.00, -23.03) abline(v = periodbase[1:length(periodbase)], lty = 5, col = "grey75", lwd = 1.5) abline(h = 0, col = "grey70") #### Boxes for OLS and AR1 y1 <- -0.21 y2 <- -0.305 polygon(c(-stagetab$stage_base[1],0,0,-stagetab$stage_base[1]), c(y1,y1,y2,y2), border = "black", col = "white", xpd = NA) leftend <- -stagetab$stage_base[c(6,31,39,44,54)] rightend <- -stagetab$stage_top[c(30,38,43,53,85)] polygon(c(-stagetab$stage_base[1], -stagetab$stage_top[5], -stagetab$stage_top[5], -stagetab$stage_base[1]), c(y1-0.002,y1-0.002,y2,y2), border = "black", angle=45, density=18, xpd = NA) polygon(c(leftend[1], rightend[1], rightend[1], leftend[1]), c(y1-0.002,y1-0.002,y2,y2), border = "black", col = "white", xpd = NA) polygon(c(leftend[2], rightend[2], rightend[2], leftend[2]), c(y1-0.002,y1-0.002,y2,y2), border = "black", col = "grey75", xpd = NA) polygon(c(leftend[3], rightend[3], rightend[3], leftend[3]), c(y1-0.002,y1-0.002,y2,y2), border = "black", angle=45, density=18, xpd = NA) polygon(c(leftend[4], rightend[4], rightend[4], leftend[4]), c(y1-0.002,y1-0.002,y2,y2), border = "black", col = "white", xpd = NA) polygon(c(leftend[5], rightend[5], rightend[5], leftend[5]), c(y1-0.002,y1-0.002,y2,y2), border = "black", col = "grey75", xpd = NA) abline(h = c(y1,y2), lwd = 1.2) text((c(leftend[1],leftend[4])+c(rightend[1],rightend[4]))/2, rep((y1+y2)/2,2), labels = rep("OLS",2), xpd = NA, cex = 0.7) text((leftend[c(2,5)]+rightend[c(2,5)])/2, rep((y1+y2)/2,2), labels = c(expression(phi), expression(phi)), xpd = NA, cex = 1) abline(v = c(-stage_base[39]+0.2), lwd = 1.3, col = "#CC0C36") ### Mesozoic xv <- -stagetab$stage_midp[44:85] polygon(c(xv[1],xv,xv[42]),c(0,pars_m[,2],0), col= rgb(0,0.8,1,0.45), border = NA) points(-stagetab$stage_midp[44:85], pars_m[,1], type = "l", col = rgb(0,0,0,1), lwd = 2.3) #points(-stagetab$stage_midp[44:85], pars_m[,1], type = "p", col = rgb(0,0,0,1), pch = 4, cex = 0.75) #points(-stagetab$stage_midp[44:85], pars[,2], type = "l", col = rgb(0,0.2,1,1), lwd = 1.2) psymbol <- rep(NA,47) psymbol[which(pars_m[,3] < 0.01)] <- 18 psymbol[which(pars_m[,3] >= 0.01 & pars_m[,3] < 0.05)] <- 18 psymbol[which(pars_m[,3] >= 0.05 & pars_m[,3] < 0.1)] <- 18 pcol <- rep(NA,47) pcol[which(pars_m[,3] < 0.01)] <- "red" pcol[which(pars_m[,3] >= 0.01 & pars_m[,3] < 0.05)] <- "orange" pcol[which(pars_m[,3] >= 0.05 & pars_m[,3] < 0.1)] <- "#F1F100" points(-stagetab$stage_midp[44:85], pars_m[,3], type = "p", col = pcol, pch = psymbol, cex = 1, lwd = 2.3) #points(-stagetab$stage_midp[44:85], pars1[,4], type = "l", col = rgb(0.7,0.7,0,1), lwd = 1.2) ### Paleozoic xv <- -stagetab$stage_midp[6:38] polygon(c(xv[1],xv,xv[33]),c(0,pars_p[,2],0), col= rgb(0,0.8,1,0.45), border = NA) points(-stagetab$stage_midp[6:38], pars_p[,1], type = "l", col = rgb(0,0,0,1), lwd = 2.3) #points(-stagetab$stage_midp[6:38], pars_p[,1], type = "p", col = rgb(0,0,0,1), pch = 4, cex = 0.75) #points(-stagetab$stage_midp[44:85], pars[,2], type = "l", col = rgb(0,0.2,1,1), lwd = 1.2) psymbol <- rep(NA,33) psymbol[which(pars_p[,3] < 0.01)] <- 18 psymbol[which(pars_p[,3] >= 0.01 & pars_p[,3] < 0.05)] <- 18 psymbol[which(pars_p[,3] >= 0.05 & pars_p[,3] < 0.1)] <- 18 pcol <- rep(NA,33) pcol[which(pars_p[,3] < 0.01)] <- "red" pcol[which(pars_p[,3] >= 0.01 & pars_p[,3] < 0.05)] <- "orange" pcol[which(pars_p[,3] >= 0.05 & pars_p[,3] < 0.1)] <- "#F1F100" points(-stagetab$stage_midp[6:38], pars_p[,3], type = "p", col = pcol, pch = psymbol, cex = 1, lwd = 1.2) #points(-stagetab$stage_midp[44:85], pars1[,4], type = "l", col = rgb(0.7,0.7,0,1), lwd = 1.2) # axes axis(side= 2, at = c(-0.4, -0.2,0,0.2,0.4,0.6, 0.8), labels = rep(NA,6), cex.axis = labelsize) text(rep(-19-stagetab$stage_base[1],5), c(-0.2,0,0.2,0.4,0.6, 0.8),labels = c(-0.2,"0",0.2,0.4,0.6, 0.8), xpd = T, adj = 1) ### Legend x1 <- -95 x2 <- -10 y1 <- 0.40 y2 <- 0.75 y12 <- y2-y1 polygon(c(x1,x2,x2,x1),xpd=TRUE, c(y1-0.095,y1-0.095,y2,y2), col= "white", border= "black", xpd = NA) segments(x1+5,y2-0.125*y12,x1+17,y2-0.125*y12, lwd = 2.3) #points(x1+10,y2-0.125*y12, pch = 4, col = "black", cex = 0.75) polygon(c(x1+5,x1+17,x1+17,x1+5),xpd=TRUE, c(y2-0.30*y12,y2-0.30*y12,y2-0.45*y12,y2-0.45*y12), col= rgb(0,0.8,1,0.45), border= NA) points(x1+11,y2-0.625*y12, pch = 18, col = "#F1F100", cex = 1) points(x1+11,y2-0.875*y12, pch = 18, col = "orange", cex = 1) points(x1+11,y2-1.125*y12, pch = 18, col = "red", cex = 1) text(rep(x1+22,4), y2-c(0.125,0.375,0.625,0.875, 1.125)*y12, labels = c("slope", expression("R"^"2"), "p < 0.1", "p < 0.05", "p < 0.01"), adj = 0, cex = 0.8) abline(v = 0, h = ylimit, lwd = 1.8) text(-56 - stage_base[1],0.3, adj = 0.5, srt = 90, labels = "regression coefficients", xpd = T) ylimit <- c(-0.302,0.8) xlable <- "Time (Ma)" lablesize <- 1 boxfactor <- 0.11 ytop <- max(ylimit) ybase <- min(ylimit) labelsize <- 1 xlabsize <- labelsize ylabel <- "" y1 <- ylimit[1] # lower boundary of plot area y2 <- ylimit[1]-boxfactor*(ylimit[2]-ylimit[1]) # position of x axis period <- c("O", "S", "D", "C", "P", "T", "J", "K", "Pg", "N") periodbase <- -c(485.4,443.8,419.2,358.9,298.9,251.902,201.3,145,66,23.03) periodtop <- -c(443.8,419.2,358.9,298.9,251.902,201.3,145,66,23.03,0) periodmid <- (periodbase + periodtop) / 2 for (i in 1:length(periodbase)) { polygon(c(periodbase[i],periodtop[i],periodtop[i],periodbase[i]),xpd=TRUE, c(y1,y1,y2,y2), col= "white", border= "black") } periodmid <- (periodbase + periodtop) / 2 text(periodmid, mean(c(y1,y2)), labels = period, srt = 0, adj = c(0.5,0.5), xpd = TRUE, cex=lablesize) axis_at <- c(-400,-300,-200,-100,0) lablesize <- 1 xlaxtadjust <- 2.71828^(-lablesize) axis(side=1, at = axis_at, labels = -axis_at, tick = TRUE, line = NA, pos = y2, outer = FALSE, font = NA, lty = "solid", lwd = 1, lwd.ticks = 1, col = NULL, col.ticks = NULL, hadj = NA, padj = -xlaxtadjust, cex.axis=lablesize) text(mean(c(min(periodbase), 0)),y1 -xlabfactor*(y1-y2),labels = xlable, xpd = TRUE) abline(h = ylimit, v = 0)