This RPub is an RMarkdown version of the supplementary code provided with the below paper, along with some additional in-line comments.

“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 Conservation Physiology 4 (1): cow007 open access link

In this publication we discuss how changing contexts (age, temperature, lab-field,…) can cause individual differences in trait values to become less/more pronounced and individual rankings for that trait to be reshuffled. As a result, the fish with the highest metabolism,fastest swim speed or most bold behaviour in one context not necessarily beats the others in another context. Apart from discussing the implications of this for conservation and management of fish populations we use simulated data to show how the extent of individual phenotypic plastiticy and the interplay between individual reaction norm variation and individual mean trait values determine how rankings will shuffle and repeatabilities will change with context. We also provide code for two metrics that become of interest when wishing to study variable traits across different contexts: context-specific repeatabilities and cross-context correlations.

Some explanation: the different panels in the figure illustrate four different scenario’s:

a & e show a scenario in which individuals differ in mean trait values but all share the same response to changing contexts. In b & f individuals differ in mean trait values and trait plasticity to changing contexts but their mean trait values and plastic responses are positively correlated (e.g. bold individuals also increase their boldness most with context). In c & g individuals again differ in mean trait values and trait plasticity to changing contexts but their mean trait values and plastic responses are not correlated. Finally, in d & g individuals differ in mean trait values and trait plasticity to changing contexts but their mean trait values and plastic responses are negatively correlated (e.g. shy individuals increase most in boldness with context).

The upper panels always show individual specific slopes and intercepts and the lower panel shows three kinds of repeatabilities:

R, full black line in panel e only: repeatability / intraclass correlation coefficient ignoring variation in slope and systematic changes for scenario 0

Radj, dashed black lines: ‘adjusted’ repeatabilities ignoring variation in slope, but accounting for systematic changes for scenario 0

Rcontext, red dotted lines: context-specific repeatabilities calculating repeatability at any givenpoint along a gradient

The code relies on the following packages to be installed on your machine. You can use the command install.packages(“package name”) to install packages beforehand.

require(lme4) # for simulations and mixed effects models
## Loading required package: lme4
## Loading required package: Matrix
require(ggplot2) # for graphing
## Loading required package: ggplot2
require(ggthemes) # graphical themes
## Loading required package: ggthemes
require(Rmisc) # plots the composed figure
## Loading required package: Rmisc
## Loading required package: lattice
## Loading required package: plyr

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

First we start simulating the data by setting the model the data should follow. I here use the same approach to simulate data following a random slope model as explained here.

In this example I 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

Next set the predictor variables for N.ind individuals…

N.ind<-40

…and N.obs repeat observations per individual

N.obs<-15

The 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 I start with simulation of panels b & d, c & g and f & h after which panels a & e are produced by simplifying the underlying variance structure.

SCENARIO 1

panels b & f: positive covariance between individual slopes and intercepts

The predictors:

simdat1<-simdat

The fixed effect parameters:

beta<-c(7,2)
names(beta)<-c("(Intercept)","context")

Generates as fixed effects a Beta intercept=7 and population level slope for context=2.

The random effect structure:

V.ind0 <- 3
V.inde <- 5
V.err0 <- 3

Sets variances for random effects.

COVind0.ind1e <- 0.7*(sqrt(V.ind0*V.inde))

Sets covariance 0.7 between slope and intercept.

Generate the variance covariance matrix for 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)
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ context + (context | ID)
##    Data: simdat1
## 
## REML criterion at convergence: 2507.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.86745 -0.64037 -0.01198  0.67059  2.70785 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  ID       (Intercept) 3.432    1.853        
##           context     4.623    2.150    0.76
##  Residual             2.933    1.713        
## Number of obs: 600, groups:  ID, 40
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.7063     0.3261  20.564
## context       1.8034     0.4216   4.277
## 
## Correlation of Fixed Effects:
##         (Intr)
## context 0.328
beta
## (Intercept)     context 
##           7           2

Do simulated variances correspond to set variance parameters?

as.data.frame(VarCorr(model1))
##        grp        var1    var2     vcov     sdcor
## 1       ID (Intercept)    <NA> 3.432159 1.8526086
## 2       ID     context    <NA> 4.623206 2.1501643
## 3       ID (Intercept) context 3.046502 0.7647969
## 4 Residual        <NA>    <NA> 2.932810 1.7125448
V.ind0 
## [1] 3
V.inde 
## [1] 5
COVind0.ind1e
## [1] 2.711088
V.err0 
## [1] 3
getME(model1,"theta")
##         ID.(Intercept) ID.context.(Intercept)             ID.context 
##              1.0817869              0.9602312              0.8089070
theta1
##         ID.(Intercept) ID.context.(Intercept)             ID.context 
##              1.0000000              0.9036961              0.9219544

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")

Visualize the data with a plot showing slopes and intercepts for all individuals. 10 random individuals are sketched in black and the remainder 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

Calculate the individual level cross-contextual correlation between context 0 and 1 This follows procedures from 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
## [1] 0.9297118

Calculate Rcontext or context-specific repeatabilities along gradient between context 0 and 1 (red dots in lower figure panels). This 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)

Now calculate the Radj (dashed line on lower figure panels) i.e. repeatability neglecting individual variation in slope by fitting arandom intercept model.

model1ri<-lmer(resp~context+(1|ID),data=simdat1)
summary(model1ri)
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ context + (1 | ID)
##    Data: simdat1
## 
## REML criterion at convergence: 2561.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.76587 -0.63735  0.00699  0.66116  2.81793 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 7.488    2.736   
##  Residual             3.311    1.819   
## Number of obs: 600, groups:  ID, 40
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.7051     0.4586  14.621
## context       1.8201     0.2641   6.892
## 
## Correlation of Fixed Effects:
##         (Intr)
## context -0.289
Radj1<-as.data.frame(VarCorr(model1ri))[1,4]/(as.data.frame(VarCorr(model1ri))[1,4]+as.data.frame(VarCorr(model1ri))[2,4])
Radj1
## [1] 0.6934343

Plot

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=Radj1),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

panels c & g: zero slope-intercept covariance

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)
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ context + (context | ID)
##    Data: simdat2
## 
## REML criterion at convergence: 2508.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.87680 -0.60104 -0.01203  0.64908  2.78129 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  ID       (Intercept) 3.306    1.818        
##           context     4.610    2.147    0.00
##  Residual             2.935    1.713        
## Number of obs: 600, groups:  ID, 40
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.6890     0.3222  20.761
## context       1.5940     0.4225   3.773
## 
## Correlation of Fixed Effects:
##         (Intr)
## context -0.233
beta
## (Intercept)     context 
##           7           2

Do simulated variances correspond to set beta parameters?

as.data.frame(VarCorr(model2))
##        grp        var1    var2       vcov       sdcor
## 1       ID (Intercept)    <NA> 3.30611934 1.818273726
## 2       ID     context    <NA> 4.61011758 2.147118437
## 3       ID (Intercept) context 0.00656179 0.001680765
## 4 Residual        <NA>    <NA> 2.93530348 1.713272741
V.ind0 
## [1] 3
V.inde 
## [1] 5
0
## [1] 0
V.err0 
## [1] 3

Do simulated theta correspond to set beta parameters?

getME(model2,"theta")
##         ID.(Intercept) ID.context.(Intercept)             ID.context 
##            1.061286789            0.002106379            1.253224517
theta2
##         ID.(Intercept) ID.context.(Intercept)             ID.context 
##               1.000000               0.000000               1.290994

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 with 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

Individual level cross-environmental / contextual correlation between context 0 and 1 for scenario 2.

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
## [1] 0.6469955

Calculate Rcontext along gradient between context 0 and 1 for scenario 2.

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)

Now calculate the Radj i.e. repeatability neglecting individual variation in slope by fitting a random intercept model.

model2ri<-lmer(resp~context+(1|ID),data=simdat2)
summary(model2ri)
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ context + (1 | ID)
##    Data: simdat2
## 
## REML criterion at convergence: 2544.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.91831 -0.65181  0.00259  0.64995  2.50221 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 4.495    2.120   
##  Residual             3.320    1.822   
## Number of obs: 600, groups:  ID, 40
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.7221     0.3681  18.259
## context       1.5286     0.2643   5.785
## 
## Correlation of Fixed Effects:
##         (Intr)
## context -0.361
Radj2<-as.data.frame(VarCorr(model2ri))[1,4]/(as.data.frame(VarCorr(model2ri))[1,4]+as.data.frame(VarCorr(model2ri))[2,4])
Radj2
## [1] 0.5751858

Plot

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=Radj2),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

panels d & h: negative slope-intercept covariance

Use same scenario as in scenario 1 but use 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)
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ context + (context | ID)
##    Data: simdat3
## 
## REML criterion at convergence: 2481.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.88631 -0.61744 -0.03651  0.67280  2.85803 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID       (Intercept) 3.265    1.807         
##           context     5.067    2.251    -0.71
##  Residual             2.938    1.714         
## Number of obs: 600, groups:  ID, 40
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.6868     0.3209   20.84
## context       1.8876     0.4359    4.33
## 
## Correlation of Fixed Effects:
##         (Intr)
## context -0.747
beta
## (Intercept)     context 
##           7           2

Do simulated variances correspond to set beta parameters?

as.data.frame(VarCorr(model3))
##        grp        var1    var2      vcov      sdcor
## 1       ID (Intercept)    <NA>  3.264972  1.8069234
## 2       ID     context    <NA>  5.067132  2.2510291
## 3       ID (Intercept) context -2.897221 -0.7122965
## 4 Residual        <NA>    <NA>  2.937754  1.7139878
V.ind0 
## [1] 3
V.inde 
## [1] 5
-COVind0.ind1e
## [1] -2.711088
V.err0 
## [1] 3

Do simulated theta correspond to set beta parameters?

getME(model3,"theta")
##         ID.(Intercept) ID.context.(Intercept)             ID.context 
##              1.0542218             -0.9354793              0.9217973
theta3
##         ID.(Intercept) ID.context.(Intercept)             ID.context 
##              1.0000000             -0.9036961              0.9219544

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 with 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

Individual level cross-environmental / contextual correlation between context 0 and 1 for scenario 3.

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
## [1] 0.1277607

Calculate Rcontext along gradient between context 0 and 1 for scenario 3.

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)

Now calculate the Radj i.e. repeatability neglecting individual variation in slope by fitting a random intercept model.

model3ri<-lmer(resp~context+(1|ID),data=simdat3)
summary(model3ri)
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ context + (1 | ID)
##    Data: simdat3
## 
## REML criterion at convergence: 2518.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.91578 -0.66264 -0.01383  0.64903  2.59906 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 1.727    1.314   
##  Residual             3.375    1.837   
## Number of obs: 600, groups:  ID, 40
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.7382     0.2581  26.108
## context       1.7846     0.2656   6.719
## 
## Correlation of Fixed Effects:
##         (Intr)
## context -0.517
Radj3<-as.data.frame(VarCorr(model3ri))[1,4]/(as.data.frame(VarCorr(model3ri))[1,4]+as.data.frame(VarCorr(model3ri))[2,4])
Radj3
## [1] 0.338551

Plot

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=Radj3),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

panels a & e: no individual slope variation

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)

Do simulated betas correspond to set beta parameters?

summary(model0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ context + (context | ID)
##    Data: simdat0
## 
## REML criterion at convergence: 2461.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.95089 -0.65687 -0.02382  0.65663  2.83893 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  ID       (Intercept) 3.5610982 1.88709       
##           context     0.0007204 0.02684  -1.00
##  Residual             2.9122465 1.70653       
## Number of obs: 600, groups:  ID, 40
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.7141     0.3307    20.3
## context       2.4991     0.2475    10.1
## 
## Correlation of Fixed Effects:
##         (Intr)
## context -0.391
beta
## (Intercept)     context 
##           7           2

Do simulated variances correspond to set variance parameters?

as.data.frame(VarCorr(model0))
##        grp        var1    var2          vcov       sdcor
## 1       ID (Intercept)    <NA>  3.5610981834  1.88708722
## 2       ID     context    <NA>  0.0007203892  0.02684007
## 3       ID (Intercept) context -0.0506495469 -1.00000000
## 4 Residual        <NA>    <NA>  2.9122464855  1.70653054
V.ind0 
## [1] 3
V.inde 
## [1] 1e-08
COVind0.ind1e
## [1] 1.732051e-12
V.err0 
## [1] 3
getME(model0,"theta")
##         ID.(Intercept) ID.context.(Intercept)             ID.context 
##           1.105803e+00          -1.572786e-02           7.777867e-08
theta0
##         ID.(Intercept) ID.context.(Intercept)             ID.context 
##           1.000000e+00           5.773503e-13           5.773503e-05

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 with a plot showing slopes and intercepts for all individuals for 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

Individual level cross-environmental / contextual correlation between context 0 and 1 for scenario 0.

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
## [1] 1

Calculate Rcontext along gradient between context 0 and 1 for scenario 0.

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)

Now calculate the Radj i.e. repeatability ignoring individual variation in slope by fitting random intercept model.

model0ri<-lmer(resp~context+(1|ID),data=simdat0)
summary(model0ri)
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ context + (1 | ID)
##    Data: simdat0
## 
## REML criterion at convergence: 2461.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.94936 -0.65897 -0.02333  0.66121  2.83932 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 3.513    1.874   
##  Residual             2.912    1.707   
## Number of obs: 600, groups:  ID, 40
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.7145     0.3288   20.42
## context       2.4982     0.2474   10.10
## 
## Correlation of Fixed Effects:
##         (Intr)
## context -0.378
Radj0<-as.data.frame(VarCorr(model0ri))[1,4]/(as.data.frame(VarCorr(model0ri))[1,4]+as.data.frame(VarCorr(model0ri))[2,4])
Radj0
## [1] 0.5467254

Calculate R or repeatability calculated from random intercept model ignoring variation in slope and ignoring systematic changes (shown as full black line in figure panel e).

model0ria<-lmer(resp~1+(1|ID),data=simdat0)
summary(model0ria)
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ 1 + (1 | ID)
##    Data: simdat0
## 
## REML criterion at convergence: 2554.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.83970 -0.70448 -0.01216  0.62766  2.51756 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 3.557    1.886   
##  Residual             3.433    1.853   
## Number of obs: 600, groups:  ID, 40
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   7.9694     0.3077    25.9
Radj0a<-as.data.frame(VarCorr(model0ria))[1,4]/(as.data.frame(VarCorr(model0ria))[1,4]+as.data.frame(VarCorr(model0ria))[2,4])
Radj0a
## [1] 0.508898
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=Radj0),linetype = 2)+
  geom_hline(aes(yintercept=Radj0a),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

Wrapping up

To summarize, the cross-environmental repeatabilities for each of 4 scenario’s:

r_between0
## [1] 1
r_between1
## [1] 0.9297118
r_between2
## [1] 0.6469955
r_between3
## [1] 0.1277607

Radj i.e. repeatability ignoring individual variation in slope for each of 4 scenario’s:

Radj0
## [1] 0.5467254
Radj1
## [1] 0.6934343
Radj2
## [1] 0.5751858
Radj3
## [1] 0.338551

And the 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)