FAPAR correlations under warm and fert effects

This is the R code to the manuscript
Deciduous tundra shrubs shift toward more acquisitive light absorption strategy under climate change treatments
by R. J. Heim, M. Iturrate-Garcia, M. Reji Chacko, S. Karsanaev, T. C. Maximov, M. M. P. D. Heijmans, G. Schaepman-Strub

Data preparation

lops_salp = read.csv("FAPAR_salix.csv")
lops_salp$species = "salix"
lops_betn = read.csv("FAPAR_betula.csv")
lops_betn$species = "betula"
lops = as.data.frame(rbind(lops_salp, lops_betn))
agg_final = read.csv("WE_traits_all.csv")
agg_final = subset(agg_final, (betula == "Y" | salix == "Y"), select = c(blk, plt, ind, betula, cable, heating, ftr, LNC, LPC, CN_ratio, NP_ratio,LA, SLA, LDMC, LA, LMA1, LMA2))
agg_final$thickness = 1/(agg_final$SLA*agg_final$LDMC)
agg_final$LMA = (agg_final$LMA1+agg_final$LMA2)/2
agg_final$species[agg_final$betula == "Y"] = "betula"
agg_final$species[agg_final$betula == "N"] = "salix"
agg_final$Treatment = agg_final$ftr
agg_final$Plot = agg_final$plt
agg_final <- merge(agg_final, lops, by=c("species","Treatment", "Plot"))
agg_final<-as.data.frame(aggregate(cbind(FAPAR, thickness, SLA, LMA,LNC, LPC, LDMC,CN_ratio,NP_ratio, LA) ~ 1 + 
                                     Plot + blk + species + cable + heating + ftr, data=agg_final, FUN=mean))
agg_final$heat_treatment<-as.factor(paste(agg_final$heating, agg_final$cable))
agg_final$fert_treatment<-as.factor(paste(agg_final$ftr))
agg_final$species<-as.factor(paste(agg_final$species))

load packages

library(igraph)
library(brms)
library(rstan)
library(DHARMa)
library(ggplot2)
library(bayesplot)
library(ggExtra)
library(data.table)
library(png)
library(grid)

Model preparation

Models that allow group level variance parameters to covary across traits (1| p |blk)

agg_final$scale_FAPAR<-scale(agg_final$FAPAR)
agg_final$scale_thickness<-scale(agg_final$thickness)
agg_final$scale_SLA<-scale(agg_final$SLA)
agg_final$scale_LDMC<-scale(agg_final$LDMC)
agg_final$scale_LNC<-scale(agg_final$LNC)
agg_final$scale_LPC<-scale(agg_final$LPC)
agg_final$scale_CN_ratio<-scale(agg_final$CN_ratio)
agg_final$scale_NP_ratio<-scale(agg_final$NP_ratio)
agg_final$scale_LA<-scale(agg_final$LA)
agg_final$scale_LMA<-scale(agg_final$LMA)

Model for complete dataset

mod_fapar_cor4<- brm(bf(mvbind(scale_FAPAR,scale_thickness,
                               scale_SLA,scale_LMA,scale_LA,scale_LDMC,scale_LNC,scale_LPC,scale_CN_ratio,
                               scale_NP_ratio) ~ 1 + species+
                          (1| p |blk),
                        (sigma ~ species)) + 
                       set_rescor(TRUE),
                     family = student(),
                     data   = agg_final, 
                     iter   = 4000, 
                     warmup = 2000,
                     chains = 4,
                     control=list(adapt_delta=0.99,max_treedepth=15)
)
mod_fapar_cor4
summary(mod_fapar_cor4)$rescor_pars

#plot(mod_fapar_cor4)
##pp_check(mod_fapar_cor4)
##pp_check(mod_fapar_cor4, type = "ecdf_overlay")

hypothesis test if correlation >0

post<-posterior_samples(mod_fapar_cor4)
hypothesis(post, "rescor__scaleFAPAR__scalethickness > 0") #0.92 

Hypothesis Tests for class :

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(post, "rescor__scaleFAPAR__scaleSLA > 0") #0.51 

Hypothesis Tests for class :

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(post, "rescor__scaleFAPAR__scaleLMA > 0") #0.54

Hypothesis Tests for class :

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(post, "rescor__scaleFAPAR__scaleLA > 0") #0.88 

Hypothesis Tests for class :

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(post, "rescor__scaleFAPAR__scaleLDMC > 0") #0.18     

Hypothesis Tests for class :

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(post, "rescor__scaleFAPAR__scaleLNC > 0") #1

Hypothesis Tests for class :

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(post, "rescor__scaleFAPAR__scaleLPC > 0") #0.98

Hypothesis Tests for class :

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(post, "rescor__scaleFAPAR__scaleCNratio > 0") #0.00

Hypothesis Tests for class :

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(post, "rescor__scaleFAPAR__scaleNPratio > 0") #0.28

Hypothesis Tests for class :

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

###get rescor summary
#all
library(data.table)
ar_mod_fapar_cor4<-as.array(summary(mod_fapar_cor4))
resc_mod_fapar_cor4<-ar_mod_fapar_cor4$rescor_pars
resc_mod_fapar_cor4_F<- resc_mod_fapar_cor4[row.names(resc_mod_fapar_cor4) %like% "FAPAR", ]
resc_mod_fapar_cor4_F$y <- c("Thickness","SLA","LMA","LA","LDMC" ,"LNC","LPC","CN ratio","NP ratio")
resc_mod_fapar_cor4_F$y
resc_mod_fapar_cor4_F$lower<-resc_mod_fapar_cor4_F[,3]
resc_mod_fapar_cor4_F$upper<-resc_mod_fapar_cor4_F[,4]
resc_mod_fapar_cor4_F

resc_mod_fapar_cor4_F$facet<-"Both shrub species"
resc_mod_fapar_cor4_F$y <- factor(resc_mod_fapar_cor4_F$y,levels=c("NP ratio","CN ratio","LPC","LNC","LDMC" ,"LA","LMA","SLA","Thickness"))
plot_cor4<-ggplot(data=resc_mod_fapar_cor4_F, aes(x = Estimate, y = y))+
  labs(y="")+
  labs(x="")+
  xlim(c(-0.9,0.9))+
  facet_grid(~facet)+
  geom_vline(xintercept = 0,col="black") +
  geom_linerange(data=resc_mod_fapar_cor4_F, aes(xmin = lower, xmax=upper,y = y,col=Estimate))+
  geom_point(data=resc_mod_fapar_cor4_F, aes(x = Estimate, y = y,col=Estimate))+
  scale_color_gradient2(low="blue",mid="grey60",high="#C90076")+
  theme_bw()+
  theme(legend.position="n")+
  guides(colour = guide_legend(override.aes = list(alpha=1)))+
  theme(axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        strip.text = element_text(face = "italic"))+
  theme(plot.margin = unit(c(0.1,0.1,0.1, 0.1), "cm"))
plot_cor4

Models for species subsets

bet<-subset(agg_final,species=="betula")
sal<-subset(agg_final,species=="salix")

mod_fapar_bet<- brm(bf(mvbind(scale_FAPAR,scale_thickness,
                              scale_SLA,scale_LMA,scale_LA,scale_LDMC,scale_LNC,scale_LPC,scale_CN_ratio,
                              scale_NP_ratio) ~ 1 +
                         (1| p |blk)) + 
                      set_rescor(TRUE),
                    family = student(),
                    data   = bet, 
                    iter   = 4000, 
                    warmup = 2000,
                    chains = 4,
                    control=list(adapt_delta=0.99,max_treedepth=15)
)
mod_fapar_bet
#plot(mod_fapar_bet)
##pp_check(mod_fapar_bet)
##pp_check(mod_fapar_bet, type = "ecdf_overlay")

## hypothesis test if correlation >0 Betula nana
post<-posterior_samples(mod_fapar_bet)
hypothesis(post, "rescor__scaleFAPAR__scalethickness > 0") #0.64  

hypothesis(post, "rescor__scaleFAPAR__scaleSLA > 0") #0.66

hypothesis(post, "rescor__scaleFAPAR__scaleLMA > 0") #0.27

hypothesis(post, "rescor__scaleFAPAR__scaleLA > 0") #0.48

hypothesis(post, "rescor__scaleFAPAR__scaleLDMC > 0") #0.38     

hypothesis(post, "rescor__scaleFAPAR__scaleLNC > 0") #0.99

hypothesis(post, "rescor__scaleFAPAR__scaleLPC > 0") #0.97

hypothesis(post, "rescor__scaleFAPAR__scaleCNratio > 0") #0.02 

hypothesis(post, "rescor__scaleFAPAR__scaleNPratio > 0") #0.03

###get rescor summary
#all
library(data.table)
ar_mod_fapar_bet<-as.array(summary(mod_fapar_bet))
resc_mod_fapar_bet<-ar_mod_fapar_bet$rescor_pars
resc_mod_fapar_bet_F<- resc_mod_fapar_bet[row.names(resc_mod_fapar_bet) %like% "FAPAR", ]
resc_mod_fapar_bet_F$y <- c("Thickness","SLA","LMA","LA","LDMC" ,"LNC","LPC","CN ratio","NP ratio")
resc_mod_fapar_bet_F$y
resc_mod_fapar_bet_F$lower<-resc_mod_fapar_bet_F[,3]
resc_mod_fapar_bet_F$upper<-resc_mod_fapar_bet_F[,4]
resc_mod_fapar_bet_F

resc_mod_fapar_bet_F$facet<-"B. nana"
resc_mod_fapar_bet_F$y <- factor(resc_mod_fapar_bet_F$y,levels=c("NP ratio","CN ratio","LPC","LNC","LDMC" ,"LA","LMA","SLA","Thickness"))
plot_corbet<-ggplot(data=resc_mod_fapar_bet_F, aes(x = Estimate, y = y))+
  labs(y="")+
  labs(x="")+
  xlim(c(-0.9,0.9))+
  facet_grid(~facet)+
  geom_vline(xintercept = 0,col="black") +
  geom_linerange(data=resc_mod_fapar_bet_F, aes(xmin = lower, xmax=upper,y = y,col=Estimate))+
  geom_point(data=resc_mod_fapar_bet_F, aes(x = Estimate, y = y,col=Estimate))+
  scale_color_gradient2(low="blue",mid="grey60",high="#C90076")+
  theme_bw()+
  theme(legend.position="n")+
  guides(colour = guide_legend(override.aes = list(alpha=1)))+
  theme(axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        strip.text = element_text(face = "italic"))+
  theme(plot.margin = unit(c(0.1,0.1,0.1, 0.1), "cm"))
plot_corbet

mod_fapar_sal<- brm(bf(mvbind(scale_FAPAR,scale_thickness,
                              scale_SLA,scale_LMA,scale_LA,scale_LDMC,scale_LNC,scale_LPC,scale_CN_ratio,
                              scale_NP_ratio) ~ 1 + 
                         (1| p |blk)) + 
                      set_rescor(TRUE),
                    family = student(),
                    data   = sal, 
                    iter   = 4000, 
                    warmup = 2000,
                    chains = 4,
                    control=list(adapt_delta=0.99,max_treedepth=15)
)

hypothesis test if correlation >0 Salix pulchra

post<-posterior_samples(mod_fapar_sal)
hypothesis(post, "rescor__scaleFAPAR__scalethickness > 0") #0.88   

hypothesis(post, "rescor__scaleFAPAR__scaleSLA > 0") #0.36

hypothesis(post, "rescor__scaleFAPAR__scaleLMA > 0") #0.54 

hypothesis(post, "rescor__scaleFAPAR__scaleLA > 0") # 0.93

hypothesis(post, "rescor__scaleFAPAR__scaleLDMC > 0") #0.3  

hypothesis(post, "rescor__scaleFAPAR__scaleLNC > 0") #0.99

hypothesis(post, "rescor__scaleFAPAR__scaleLPC > 0") #0.54

hypothesis(post, "rescor__scaleFAPAR__scaleCNratio > 0") #0.02

hypothesis(post, "rescor__scaleFAPAR__scaleNPratio > 0") # 0.87

###get rescor summary
#all
library(data.table)
ar_mod_fapar_sal<-as.array(summary(mod_fapar_sal))
resc_mod_fapar_sal<-ar_mod_fapar_sal$rescor_pars
resc_mod_fapar_sal_F<- resc_mod_fapar_sal[row.names(resc_mod_fapar_sal) %like% "FAPAR", ]
resc_mod_fapar_sal_F$y <- c("Thickness","SLA","LMA","LA","LDMC" ,"LNC","LPC","CN ratio","NP ratio")
resc_mod_fapar_sal_F$y
resc_mod_fapar_sal_F$lower<-resc_mod_fapar_sal_F[,3]
resc_mod_fapar_sal_F$upper<-resc_mod_fapar_sal_F[,4]
resc_mod_fapar_sal_F

resc_mod_fapar_sal_F$facet<-"S. pulchra"
resc_mod_fapar_sal_F$y <- factor(resc_mod_fapar_sal_F$y,levels=c("NP ratio","CN ratio","LPC","LNC","LDMC" ,"LA","LMA","SLA","Thickness"))
resc_mod_fapar_sal_F$y
plot_corsal<-ggplot(data=resc_mod_fapar_sal_F, aes(x = Estimate, y = y))+
  labs(y="")+
  labs(x="")+
  xlim(c(-0.9,0.9))+
  facet_grid(~facet)+
  geom_vline(xintercept = 0,col="black") +
  geom_linerange(data=resc_mod_fapar_sal_F, aes(xmin = lower, xmax=upper,y = y,col=Estimate))+
  geom_point(data=resc_mod_fapar_sal_F, aes(x = Estimate, y = y,col=Estimate))+
  scale_color_gradient2(low="blue",mid="grey60",high="#C90076")+
  theme_bw()+
  theme(legend.position="n")+
  guides(colour = guide_legend(override.aes = list(alpha=1)))+
  theme(axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        strip.text = element_text(face = "italic"))+
  theme(plot.margin = unit(c(0.1,0.1,0.1, 0.1), "cm"))
plot_corsal

Preparation for Plots

###plots
img_betula <- readPNG("betula.png")
g_betula <- rasterGrob(img_betula, interpolate=TRUE)
img_salix <- readPNG("salix.png")
g_salix <- rasterGrob(img_salix, interpolate=TRUE)
img_fert <- readPNG("fert.png")
g_fert <- rasterGrob(img_fert, interpolate=TRUE)
img_nofert <- readPNG("nofert.png")
g_nofert <- rasterGrob(img_nofert, interpolate=TRUE)
img_heat <- readPNG("heat.png")
g_heat <- rasterGrob(img_heat, interpolate=TRUE)
img_noheat <- readPNG("noheat.png")
g_noheat <- rasterGrob(img_noheat, interpolate=TRUE)
img_control <- readPNG("control.png")
g_control <- rasterGrob(img_control, interpolate=TRUE)
annotation_custom2 <- 
  function (grob, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, data){ layer(data = data, stat = StatIdentity, position = PositionIdentity, 
                                                                                 geom = ggplot2:::GeomCustomAnn,
                                                                                 inherit.aes = TRUE, params = list(grob = grob, 
                                                                                                                   xmin = xmin, xmax = xmax, 
                                                                                                                   ymin = ymin, ymax = ymax))}
#betula
a1 = annotation_custom2(rasterGrob(img_betula, interpolate=TRUE), xmin=-0.9, xmax=-0.5, ymin=8.5, ymax=9.5, data=resc_mod_fapar_bet_F[1,])
plot_cor2_fin<-plot_corbet + a1
plot_cor2_fin

#plot3
a2 = annotation_custom2(rasterGrob(img_salix, interpolate=TRUE), xmin=-0.9, xmax=-0.5, ymin=8.5, ymax=9.5, data=resc_mod_fapar_sal_F[1,])
plot_cor3_fin<-plot_corsal + a2
plot_cor3_fin

#plot4
plot_cor4

a3 = annotation_custom2(rasterGrob(img_betula, interpolate=TRUE), xmin=-0.9, xmax=-0.5, ymin=8.5, ymax=9.5, data=resc_mod_fapar_cor4_F[1,])
a4 = annotation_custom2(rasterGrob(img_salix, interpolate=TRUE), xmin=-0.55, xmax=-0.15, ymin=8.2, ymax=9.2, data=resc_mod_fapar_cor4_F[1,])
plot_cor4_fin<-plot_cor4 + a4 + a3
plot_cor4_fin

combine plots

library(cowplot)
xlab <- ggdraw() + 
  draw_label("Correlation",  hjust = -0.04,  vjust = -0.5,size=10)

all1<-plot_grid(plot_cor4_fin,plot_cor2_fin,plot_cor3_fin,labels=c("A","B", "C"),ncol=3)

all1

png(file="plot_cor_final.png", width=18, height=8, units="cm", pointsize=12,res=600)
all1
dev.off()