# Electronic Supplement to the article 
# "Context-dependency of trait repeatability and its relevance for management 
#  and conservation of fish populations."
# by Killen SS, Adriaenssens B, Marras S, Claireaux G & Cooke SJ 
# Journal: Conservation Physiology 

# In this supplement we provide the complete R code to
# (1) simulate data used in figure 1
# (2) plot figure 1
# and (3) calculate cross-context correlations and context specific repeatabilities

# comments are blocked with '#' making the entire document fully executable once pasted in R / Rstudio.

# the code relies on the following packages to be installed on your machine
# use the command install.packages("package name") to install packages beforehand

require(lme4) # for simulations and mixed effects models
require(ggplot2) # for graphing
require(ggthemes) # graphical themes
require(Rmisc) # plots the composed figure

# the following package versions were used at the making of this code file
# ggplot2_2.0.0
# lme4_1.1-10
# Rmisc_1.5
# ggthemes_3.0.0

# Start simulating the data:
# set the model the data should follow. Here we set a 
# standard random intercept-random slope model with individual as nesting factor,
# and slopes for a continous contextual variable (context, could be time, temperature, oxygen levels) 
# differing between individuals

form<-as.formula(c("~context+(context|ID)"))

# context is fitted as a fixed effect to account for mean plasticity at
# a group level as well as random effect. The way to interpret this is that 
# all individuals may increase activity with temperature (mean pop effect, fixed effect)
# but some individuals do so more than others (random slope effect)

# In this model there are four (co)variances fitted
# 1) V.ind0 = variance in reaction norm intercepts between individuals
# 2) V.inde = variance in reaction norm slopes between individuals
# 3) COVind0-ind1e = covariance between reaction norm slopes and intercepts
# 4) V.e0 = residual variance

# set the predictor variables:

# number of individuals sampled
N.ind<-40

# number of repeat observations per individual
N.obs<-15

# set.seed command makes the code reproducible and result in the same random draw every time
# the code is ran.
set.seed(25)
simdat <-data.frame(ID=factor(rep(1:N.ind,each=N.obs)),context=runif(N.ind*N.obs,0,1))


# for simplicity we start with simulation of figures b-d and f-h
# after which figures a and e are produced by simplifying the underlying variance structure



##########################################################################
##########################################################################
##########################################################################
# scenario 1: positive covariance between individual slopes and intercepts
# produces figures 1b and 1f

# the predictors
simdat1<-simdat

# the fixed effect parameters
beta<-c(7,2)
names(beta)<-c("(Intercept)","context")
# fixed effects:
# Beta intercept=7
# Population level slope for context=2

# the random effect structure
V.ind0 <- 3
V.inde <- 5
V.err0 <- 3

# scenario 1:
COVind0.ind1e <- 0.7*(sqrt(V.ind0*V.inde))

# variance covariance matrix scenario 1
vcov1<-matrix(c(V.ind0,COVind0.ind1e,COVind0.ind1e,V.inde), 2, 2)

# calculate theta
theta1<-c((chol(vcov1)/sqrt(V.err0))[1,1],
         (chol(vcov1)/sqrt(V.err0))[1,2],
         (chol(vcov1)/sqrt(V.err0))[2,2])
names(theta1)<-c("ID.(Intercept)","ID.context.(Intercept)","ID.context")

# simulate data following above parameters
set.seed(25)
response<-simulate(form,newdata=simdat1,family=gaussian,weights=NULL,
         newparams=list(theta = theta1,beta = beta, sigma = sqrt(V.err0)))

simdat1$resp<-as.vector(response[,1])

# check the simulated data.
model1<-lmer(resp~context+(context|ID),data=simdat1)

# do simulated betas correspond to set beta parameters?
summary(model1)
beta

# do simulated variances correspond to set variance parameters?
as.data.frame(VarCorr(model1))
V.ind0 
V.inde 
COVind0.ind1e
V.err0 

getME(model1,"theta")
theta1

# extract random effects for plotting
fittednorms1<-ranef(model1)$ID[,1:2]
fittednorms1[,1]<-fittednorms1[,1]+getME(model1,"beta")[1]
fittednorms1[,2]<-fittednorms1[,2]+getME(model1,"beta")[2]
colnames(fittednorms1)<-c("intercept","slope") 


########### visualizing data #######################################

# a plot showing slopes and intercepts for all individuals
# 10 random individuals in black, others in grey

set.seed(25)
rand10<-sample(1:N.ind, 10, replace = FALSE)

plot1.1<-ggplot(simdat1,aes(context, resp))+
  geom_rangeframe(data=data.frame(x=c(0, 1), y=c(0, 20)), aes(x, y))+ 
  theme_tufte(base_size = 25)+
  geom_point(colour=NA)+
  geom_segment(aes(x = 0, y = intercept,xend=1, yend=intercept+slope ), data=fittednorms1,colour = "gray95", size = 1)+
  geom_segment(aes(x = 0, y = intercept,xend=1, yend=intercept+slope ), data=fittednorms1[rand10,],colour = "black", size = 1)+
  labs(x = "Context", y="Phenotypic trait",title="slope-intercept correlation=.7")+
  theme(axis.title.x=element_text(vjust=-0.5),axis.title.y=element_text(vjust=1))+
  scale_x_continuous(breaks=seq(0,1,0.2))+
  scale_y_continuous(breaks=seq(0,20,4))
plot1.1


### scenario 1: individual level cross-environmental / 
### contextual correlation between context 0 and 1
### procedure follows Brommer Behav Ecol Sociobiol (2013) 67:1709-1718

context1<-0
context2<-1

k.matrix<-matrix(c(as.data.frame(VarCorr(model1))[1,4],
                as.data.frame(VarCorr(model1))[3,4],
                as.data.frame(VarCorr(model1))[3,4],
                as.data.frame(VarCorr(model1))[2,4]), 2, 2)

phi.matrix<-matrix(c(1,
               1,
               context1,
               context2), 2, 2)

P <- phi.matrix%*%k.matrix%*%t(phi.matrix)

r_between1<-P[1,2]/sqrt(P[1,1]*P[2,2])
r_between1


### scenario 1: context-specific repeatabilities 
### along gradient between context 0 and 1
### procedure follows Martin et al. 2011 Methods in Ecology and Evolution 2, p372

contexts<-seq(context1,context2,0.05)

rept.context<-NULL

for (i in 1:21 ) {
  numerator<-c(as.data.frame(VarCorr(model1))[1,4]+ 
                 2*as.data.frame(VarCorr(model1))[3,4]*contexts[i]+
                 as.data.frame(VarCorr(model1))[2,4]*contexts[i]^2)
    rept.context[i] = numerator/(numerator+as.data.frame(VarCorr(model1))[4,4])
}

repts<-data.frame(context=contexts,rept=rept.context)


### scenario 1: R random intercept model
### repeatability calculated from random intercept model
### ~ neglecting variation in slope 

model1ri<-lmer(resp~context+(1|ID),data=simdat1)
summary(model1ri)
Rri1<-as.data.frame(VarCorr(model1ri))[1,4]/(as.data.frame(VarCorr(model1ri))[1,4]+as.data.frame(VarCorr(model1ri))[2,4])
Rri1


# show graphically
plot1.2<-ggplot(repts,aes(context, rept))+
  geom_rangeframe(data=data.frame(x=c(0, 1), y=c(0, 1)), aes(x, y))+ 
  theme_tufte(base_size = 25)+
  geom_point(colour="red", size=3)+
  geom_hline(aes(yintercept=Rri1),linetype = 2)+
  labs(y="Repeatability",x="Context", title="slope-intercept correlation=.7")+
  theme(axis.title.y=element_text(vjust=1))+
  scale_y_continuous(breaks=seq(0,1,0.2))+
  theme(axis.title.x=element_text(vjust=-0.5),axis.title.y=element_text(vjust=1))
plot1.2



#############################################
#############################################
#############################################
# scenario 2: zero slope-intercept covariance
# produces figures 1c and 1g

# use same scenario as above but set COVind0.ind1e = 0
vcov2<-matrix(c(V.ind0,0,0,V.inde), 2, 2)

theta2<-c((chol(vcov2)/sqrt(V.err0))[1,1],
          (chol(vcov2)/sqrt(V.err0))[1,2],
          (chol(vcov2)/sqrt(V.err0))[2,2])
names(theta2)<-c("ID.(Intercept)","ID.context.(Intercept)","ID.context")


# simulate the data
simdat2<-simdat
set.seed(25)
response2<-simulate(form,newdata=simdat2,family=gaussian,weights=NULL,
                   newparams=list(theta = theta2,beta = beta, sigma = sqrt(V.err0)))

simdat2$resp<-as.vector(response2[,1])

# check the simulated data
model2<-lmer(resp~context+(context|ID),data=simdat2)

# do simulated betas correspond to set beta parameters?
summary(model2)
beta

# simulated variances correspond to set beta parameters?
as.data.frame(VarCorr(model2))
V.ind0 
V.inde 
0
V.err0 

# do simulated theta correspond to set beta parameters?
getME(model2,"theta")
theta2

# extract the random effects for plotting 
fittednorms2<-ranef(model2)$ID[,1:2]
fittednorms2[,1]<-fittednorms2[,1]+getME(model2,"beta")[1]
fittednorms2[,2]<-fittednorms2[,2]+getME(model2,"beta")[2]
colnames(fittednorms2)<-c("intercept","slope")          


########### visualizing data #######################################

# a plot showing slopes and intercepts for all individuals
# 10 random individuals in black, others in grey

set.seed(25)
rand10<-sample(1:N.ind, 10, replace = FALSE)

plot2.1<-ggplot(simdat2,aes(context, resp))+
  geom_rangeframe(data=data.frame(x=c(0, 1), y=c(0, 20)), aes(x, y))+ 
  theme_tufte(base_size = 25)+
  geom_point(colour=NA)+
  geom_segment(aes(x = 0, y = intercept,xend=1, yend=intercept+slope ), data=fittednorms2,colour = "gray95", size = 1)+
  geom_segment(aes(x = 0, y = intercept,xend=1, yend=intercept+slope ), data=fittednorms2[rand10,],colour = "black", size = 1)+
  labs(x = "Context", y="Phenotypic trait",title="slope-intercept correlation=0")+
  theme(axis.title.x=element_text(vjust=-0.5),axis.title.y=element_text(vjust=1))+
  scale_x_continuous(breaks=seq(0,1,0.2))+
  scale_y_continuous(breaks=seq(0,20,4))
plot2.1


### scenario 2: individual level cross-environmental / 
### contextual correlation between context 0 and 1

k.matrix<-matrix(c(as.data.frame(VarCorr(model2))[1,4],
                   as.data.frame(VarCorr(model2))[3,4],
                   as.data.frame(VarCorr(model2))[3,4],
                   as.data.frame(VarCorr(model2))[2,4]), 2, 2)

P <- phi.matrix%*%k.matrix%*%t(phi.matrix)

r_between2<-P[1,2]/sqrt(P[1,1]*P[2,2])
r_between2


### scenario 2: context-specific repeatabilities 
### along gradient between context 0 and 1

contexts<-seq(context1,context2,0.05)

rept.context<-NULL

for (i in 1:21 ) {
  numerator<-c(as.data.frame(VarCorr(model2))[1,4]+ 
                 2*as.data.frame(VarCorr(model2))[3,4]*contexts[i]+
                 as.data.frame(VarCorr(model2))[2,4]*contexts[i]^2)
  rept.context[i] = numerator/(numerator+as.data.frame(VarCorr(model2))[4,4])
}

repts2<-data.frame(context=contexts,rept=rept.context)


#### scenario 2: R random intercept model
# repeatability calculated from random intercept model
# ~ neglecting variation in slope 


model2ri<-lmer(resp~context+(1|ID),data=simdat2)
summary(model2ri)
Rri2<-as.data.frame(VarCorr(model2ri))[1,4]/(as.data.frame(VarCorr(model2ri))[1,4]+as.data.frame(VarCorr(model2ri))[2,4])
Rri2


# show graphically
plot2.2<-ggplot(repts2,aes(context, rept))+
  geom_rangeframe(data=data.frame(x=c(0, 1), y=c(0, 1)), aes(x, y))+ 
  theme_tufte(base_size = 25)+
  geom_point(colour="red", size=3)+
  geom_hline(aes(yintercept=Rri2),linetype = 2)+
  labs(y="Repeatability", x="Context",title="slope-intercept correlation=0")+
  theme(axis.title.y=element_text(vjust=1))+
  scale_y_continuous(breaks=seq(0,1,0.2))+
  theme(axis.title.x=element_text(vjust=-0.5),axis.title.y=element_text(vjust=1))
plot2.2



#############################################
#############################################
#############################################
# scenario 3: zero slope-intercept covariance
# produces figures 1d and 1h

# use same scenario as in one but negative COVind0.ind1e 
vcov3<-matrix(c(V.ind0,-COVind0.ind1e,-COVind0.ind1e,V.inde), 2, 2)

theta3<-c((chol(vcov3)/sqrt(V.err0))[1,1],
          (chol(vcov3)/sqrt(V.err0))[1,2],
          (chol(vcov3)/sqrt(V.err0))[2,2])
names(theta3)<-c("ID.(Intercept)","ID.context.(Intercept)","ID.context")

# simulate the data
simdat3<-simdat
set.seed(25)
response3<-simulate(form,newdata=simdat3,family=gaussian,weights=NULL,
                    newparams=list(theta = theta3,beta = beta, sigma = sqrt(V.err0)))

simdat3$resp<-as.vector(response3[,1])

# check the simulated data
model3<-lmer(resp~context+(context|ID),data=simdat3)

# do simulated betas correspond to set beta parameters?
summary(model3)
beta

# simulated variances correspond to set beta parameters?
as.data.frame(VarCorr(model3))
V.ind0 
V.inde 
-COVind0.ind1e
V.err0 

# simulated theta correspond to set beta parameters?
getME(model3,"theta")
theta3

# extract the random effects for plotting 
fittednorms3<-ranef(model3)$ID[,1:2]
fittednorms3[,1]<-fittednorms3[,1]+getME(model3,"beta")[1]
fittednorms3[,2]<-fittednorms3[,2]+getME(model3,"beta")[2]
colnames(fittednorms3)<-c("intercept","slope")          


########### visualizing data #######################################

# a plot showing slopes and intercepts for all individuals
# 10 random individuals in black, others in grey
set.seed(25)
rand10<-sample(1:N.ind, 10, replace = FALSE)

plot3.1<-ggplot(simdat3,aes(context, resp))+
  geom_rangeframe(data=data.frame(x=c(0, 1), y=c(0, 20)), aes(x, y))+ 
  theme_tufte(base_size = 25)+
  geom_point(colour=NA)+
  geom_segment(aes(x = 0, y = intercept,xend=1, yend=intercept+slope ), data=fittednorms3,colour = "gray95", size = 1)+
  geom_segment(aes(x = 0, y = intercept,xend=1, yend=intercept+slope ), data=fittednorms3[rand10,],colour = "black", size = 1)+
  labs(x = "Context", y="Phenotypic trait",title="slope-intercept correlation=-.7")+
  theme(axis.title.x=element_text(vjust=-0.5),axis.title.y=element_text(vjust=1))+
  scale_x_continuous(breaks=seq(0,1,0.2))+
  scale_y_continuous(breaks=seq(0,20,4))
plot3.1


### scenario 3: individual level cross-environmental / 
### contextual correlation between context 0 and 1

k.matrix<-matrix(c(as.data.frame(VarCorr(model3))[1,4],
                   as.data.frame(VarCorr(model3))[3,4],
                   as.data.frame(VarCorr(model3))[3,4],
                   as.data.frame(VarCorr(model3))[2,4]), 2, 2)

P <- phi.matrix%*%k.matrix%*%t(phi.matrix)

r_between3<-P[1,2]/sqrt(P[1,1]*P[2,2])
r_between3


### scenario 3: context-specific repeatabilities 
### along gradient between context 0 and 1

contexts<-seq(context1,context2,0.05)

rept.context<-NULL

for (i in 1:21 ) {
  numerator<-c(as.data.frame(VarCorr(model3))[1,4]+ 
                 2*as.data.frame(VarCorr(model3))[3,4]*contexts[i]+
                 as.data.frame(VarCorr(model3))[2,4]*contexts[i]^2)
  rept.context[i] = numerator/(numerator+as.data.frame(VarCorr(model3))[4,4])
}

repts3<-data.frame(context=contexts,rept=rept.context)


#### scenario 3: R random intercept model
# repeatability calculated from random intercept model
# ~ neglecting variation in slope 


model3ri<-lmer(resp~context+(1|ID),data=simdat3)
summary(model3ri)
Rri3<-as.data.frame(VarCorr(model3ri))[1,4]/(as.data.frame(VarCorr(model3ri))[1,4]+as.data.frame(VarCorr(model3ri))[2,4])
Rri3

# show graphically
plot3.2<-ggplot(repts3,aes(context, rept))+
  geom_rangeframe(data=data.frame(x=c(0, 1), y=c(0, 1)), aes(x, y))+ 
  theme_tufte(base_size = 25)+
  geom_point(colour="red", size=3)+
  geom_hline(aes(yintercept=Rri3),linetype = 2)+
  labs(y="Repeatability", x="Context",title="slope-intercept correlation=-.7")+
  theme(axis.title.y=element_text(vjust=1))+
  scale_y_continuous(breaks=seq(0,1,0.2))+
  theme(axis.title.x=element_text(vjust=-0.5),axis.title.y=element_text(vjust=1))
plot3.2


#################################
#################################
#################################
# scenario 0: no slope variation
# produces figures 1a and 1e

# the predictors
simdat0<-simdat

# change random effect structure
V.inde <- 0.00000001
COVind0.ind1e <- 0.00000001*(sqrt(V.ind0*V.inde))

# variance covariance matrix scenario 0
vcov0<-matrix(c(V.ind0,COVind0.ind1e,COVind0.ind1e,V.inde), 2, 2)

# calculate theta

theta0<-c((chol(vcov0)/sqrt(V.err0))[1,1],
          (chol(vcov0)/sqrt(V.err0))[1,2],
          (chol(vcov0)/sqrt(V.err0))[2,2])
names(theta0)<-c("ID.(Intercept)","ID.context.(Intercept)","ID.context")

# simulate data following above parameters
set.seed(25)
response<-simulate(form,newdata=simdat0,family=gaussian,weights=NULL,
                   newparams=list(theta = theta0,beta = beta, sigma = sqrt(V.err0)))

simdat0$resp<-as.vector(response[,1])

# check the simulated data.
model0<-lmer(resp~context+(context|ID),data=simdat0)

# simulated betas correspond to set beta parameters?
summary(model0)
beta

# simulated variances correspond to set variance parameters?
as.data.frame(VarCorr(model0))
V.ind0 
V.inde 
COVind0.ind1e
V.err0 

getME(model0,"theta")
theta0

# extract random effects for plotting
fittednorms0<-ranef(model0)$ID[,1:2]
fittednorms0[,1]<-fittednorms0[,1]+getME(model0,"beta")[1]
fittednorms0[,2]<-fittednorms0[,2]+getME(model0,"beta")[2]
colnames(fittednorms0)<-c("intercept","slope") 


########### visualizing data #######################################

# a plot showing slopes and intercepts for all individuals
# 10 random individuals in black, others in grey
set.seed(25)
rand10<-sample(1:N.ind, 10, replace = FALSE)

plot0.1<-ggplot(simdat0,aes(context, resp))+
  geom_rangeframe(data=data.frame(x=c(0, 1), y=c(0, 20)), aes(x, y))+ 
  theme_tufte(base_size = 25)+
  geom_point(colour=NA)+
  geom_segment(aes(x = 0, y = intercept,xend=1, yend=intercept+slope ), data=fittednorms0,colour = "gray90", size = 1)+
  geom_segment(aes(x = 0, y = intercept,xend=1, yend=intercept+slope ), data=fittednorms0[rand10,],colour = "black", size = 1)+
  labs(x = "Context", y="Phenotypic trait",title="no slope variation")+
  theme(axis.title.x=element_text(vjust=-0.5),axis.title.y=element_text(vjust=1))+
  scale_x_continuous(breaks=seq(0,1,0.2))+
  scale_y_continuous(breaks=seq(0,20,4))
plot0.1


### scenario 0: individual level cross-environmental / 
### contextual correlation between context 0 and 1

context1<-0
context2<-1

k.matrix<-matrix(c(as.data.frame(VarCorr(model0))[1,4],
                   as.data.frame(VarCorr(model0))[3,4],
                   as.data.frame(VarCorr(model0))[3,4],
                   as.data.frame(VarCorr(model0))[2,4]), 2, 2)

phi.matrix<-matrix(c(1,
                     1,
                     context1,
                     context2), 2, 2)

P <- phi.matrix%*%k.matrix%*%t(phi.matrix)

r_between0<-P[1,2]/sqrt(P[1,1]*P[2,2])
r_between0

### scenario 0: context-specific repeatabilities 
### along gradient between context 0 and 1

contexts<-seq(context1,context2,0.05)

rept.context<-NULL

for (i in 1:21 ) {
  numerator<-c(as.data.frame(VarCorr(model0))[1,4]+ 
                 2*as.data.frame(VarCorr(model0))[3,4]*contexts[i]+
                 as.data.frame(VarCorr(model0))[2,4]*contexts[i]^2)
  rept.context[i] = numerator/(numerator+as.data.frame(VarCorr(model0))[4,4])
}

repts<-data.frame(context=contexts,rept=rept.context)


#### scenario 0: R random intercept model
# repeatability calculated from random intercept model
# ~ neglecting variation in slope / adjusting for systematic changes

model0ria<-lmer(resp~context+(1|ID),data=simdat0)
summary(model0ria)
Rri0a<-as.data.frame(VarCorr(model0ria))[1,4]/(as.data.frame(VarCorr(model0ria))[1,4]+as.data.frame(VarCorr(model0ria))[2,4])
Rri0a

# repeatability calculated from random intercept model
# ~ neglecting variation in slope / neglecting systematic changes

model0ri<-lmer(resp~1+(1|ID),data=simdat0)
summary(model0ri)
Rri0<-as.data.frame(VarCorr(model0ri))[1,4]/(as.data.frame(VarCorr(model0ri))[1,4]+as.data.frame(VarCorr(model0ri))[2,4])
Rri0

plot0.2<-ggplot(repts,aes(context, rept))+
  geom_rangeframe(data=data.frame(x=c(0, 1), y=c(0, 1)), aes(x, y))+ 
  theme_tufte(base_size = 25)+
  geom_point(colour="red", size=3)+
  geom_hline(aes(yintercept=Rri0a),linetype = 2)+
  geom_hline(aes(yintercept=Rri0),linetype = 1)+
  labs(y="Repeatability",x="Context", title="no slope variation")+
  theme(axis.title.y=element_text(vjust=1))+
  scale_y_continuous(breaks=seq(0,1,0.2))+
  theme(axis.title.x=element_text(vjust=-0.5),axis.title.y=element_text(vjust=1))
plot0.2


# cross-environmental repeatabilities for each of 4 scenario's:
r_between0
r_between1
r_between2
r_between3

# R's from random intercept model
Rri0
Rri1
Rri2
Rri3

# Figures used in paper:
plot1.1
plot1.2
plot2.1
plot2.2
plot3.1
plot3.2
plot0.1
plot0.2


# composed figure

plota<-plot0.1+
  annotate("text", x = 0.1, y = 18, label = "(a)", size=8)+
  labs(title="",x = "")
plote<-plot0.2+
  annotate("text", x = 0.1, y = 0.9, label = "(e)", size=8)+
  labs(x = "",title="")
plotb<-plot1.1+
  annotate("text", x = 0.1, y = 18, label = "(b)", size=8)+
  labs(x = "", y="",title="")
plotf<-plot1.2+
  annotate("text", x = 0.1, y = 0.9, label = "(f)", size=8)+
  labs(y="",title="")+
  theme(axis.title.x=element_text(hjust=1.4))
plotc<-plot2.1+
  annotate("text", x = 0.1, y = 18, label = "(c)", size=8)+
  labs(x = "", y="",title="")
plotg<-plot2.2+
  annotate("text", x = 0.1, y = 0.9, label = "(g)", size=8)+
  labs(x = "",y="",title="")
plotd<-plot3.1+
  annotate("text", x = 0.1, y = 18, label = "(d)", size=8)+
  labs(y="", title="",x = "")
ploth<-plot3.2+
  annotate("text", x = 0.1, y = 0.9, label = "(h)", size=8)+
  labs(x = "",y="", title="")

multiplot(plota, plote, plotb, plotf, plotc, plotg, plotd, ploth,cols=4)