---
title: |
  | R analysis
  | How do developmental and parental exposures to predation affect personality and immediate behavioural plasticity in the snail <i>Physa acuta?</i>
  | J. Tariel, S. Plénet and É. Luquet
  | Proceedings of the Royal Society B
  | doi:10.1098/rspb.2020.1761
author: "Juliette Tariel"
date: "29/10/2020"
output: 
  pdf_document:
    toc: true
    toc_depth: 2
    fig_width: 4
    fig_height: 4
editor_options: 
  chunk_output_type: console
---

# Load packages and data
## Load packages
```{r, message=FALSE, warning=FALSE, eval = F}
# data manipulation
library(dplyr)
library(tidyr)
# graphics
library(ggplot2)
library(cowplot)
theme_set(theme_classic())
# statistics
library(lme4)
library(lmerTest)
library(emmeans)
library(MCMCglmm)
se <- function(x) sd(x, na.rm = TRUE)/sqrt(sum(!is.na(x)))
# table output in word format
library(flextable) 
library(officer) 
library(measurements) # for the function 'conv_unit'
```

## Load data
```{r, eval = F}
dat <- read.table("data.csv", header = T, sep=";", dec=".",stringsAsFactors =T)
str(dat)
nrow(dat) # total number of trials
dat %>% group_by(Developmental, Parental) %>% summarise(n = n_distinct(ID)) # number of individuals

# log10 transformation of time values to achieve normality
# and time values multiplied by -1
dat$logTime <- -log10(dat$Time)

# Mass standardization
dat$sMass <- scale(dat$Mass)

# Labels of graphics
xlab = "Immediate environment"

# Recode immediate, developmental and parental environments as numeric to have the intercept 
#in a theoretical mean environment, i.e. a center environmental gradient 
# see Dingemanse, N. J., A. J. N. Kazem, D. Réale, and J. Wright. 2010. 
#Behavioural reaction norms: animal personality meets individual plasticity. 
#Trends Ecol. Evol. 25:81–89. https://doi.org/10.1016/j.tree.2009.07.013. 
dat$Immediate2    <- dat$Immediate  
dat$Developmental2<- dat$Developmental
dat$Parental2     <- dat$Parental
dat$Immediate     <- as.numeric(dat$Immediate) - 1.5
dat$Developmental <- as.numeric(dat$Developmental) - 1.5
dat$Parental      <- as.numeric(dat$Parental) - 1.5
```

# Group diversity in personality and immediate plasticity
## Model equations
In all LMMs, $\mathrm{Y}_{jkl}$ is a single measurement of escape behaviour (time to crawl-out the water in seconds) of:

+ the individual j (j spreads from 1 to 80),
+ in the immediate environment k (k = 1 for immediate environment control C and k = 2 for immediate environment predator-cue P),
+ at the trial number l (l = 1 for the first trial, l = 2 for the second trial in the immediate environment k).

$\mathrm{E}_{jkl}$ is the residual value of the measurement $\mathrm{Y}_{jkl}$. $\sigma_{res}^2$ is the residual variance. 

In our paper, fixed effects parameters were estimated with the LMM3. They are then only detailed in this model. 

### Linear model 0 (LMM0, null model)

$-log10(\mathrm{Y}_{jkl}) = \beta_{0} + \mathrm{fixed} \ \mathrm{effects} + \mathrm{E}_{jkl}$

$[\mathrm{E}_{jkl}] \stackrel{i.i.d.}{\sim}\mathcal{N}(0,\sigma^{2}_{res})$

$\beta_{0}$ is the parameter for average intercept. 

### Linear mixed model 1 (LMM1, random intercept model)

$-log10(\mathrm{Y}_{jkl}) = \beta_{0} + \mathrm{fixed} \ \mathrm{effects} + \mathrm{U_{0}}_{j} + \mathrm{E}_{jkl}$

$[\mathrm{E}_{jkl}] \stackrel{i.i.d.}{\sim}\mathcal{N}(0,\sigma^{2}_{res})$

$[\mathrm{U_{0}}_{j}] \stackrel{i.i.d.}{\sim}\mathcal{N}(0,\sigma^{2}_{i})$

$\mathrm{U_{0}}_{j}$ is the random intercept value of the individual j. $\sigma_i^2$ is the variance in intercept, representing group diversity in personality. 

### Linear mixed model 2 (LMM2, random slope model)

$-log10(\mathrm{Y}_{jkl}) = \beta_{0} + \mathrm{fixed} \ \mathrm{effects} + \mathrm{U_{0}}_{j} + \mathrm{U_{1}}_{j}\times x_{jkl}^{(I)}  + \mathrm{E}_{jk}$

$[\mathrm{E}_{jkl}] \stackrel{i.i.d.}{\sim}\mathcal{N}(0,\sigma^{2}_{res})$

$\begin{aligned}
\left[\begin{array}{l} 
\mathrm{U_{0}}_{j} \\
\mathrm{U_{1}}_{j}
\end{array}\right]
\stackrel{i.i.d.}{\sim}
\mathcal{N}
(
\left[\begin{array}{l} 
0 \\
0
\end{array}\right]
,
\left[\begin{array}{ll} 
\sigma^{2}_{i} & cov_{is} \\
cov_{is} & \sigma^{2}_{s}
\end{array}\right]
)
\end{aligned}$ 

$x_{jk}^{(I)}$ is the immediate environment for the individual j with:

+ $x_{jk}^{(I)}$ = -0.5 for k = 1 (immediate environment control C)
+ $x_{jk}^{(I)}$ =  0.5 for k = 2 (immediate environment predator-cue P)

$\mathrm{U_{1}}_{j}$ is the random slope value of the individual j. $\sigma^{2}_{s}$ is the variance in slope, representing group diversity in immediate plasticity. 

$cov_{is}$ is the covariance between intercept and slopeand it uses for the calculation of the correlation between intercept and slope:  $cor_{is} = \frac{cov_{is}}{\sqrt{\sigma^{2}_{i} \times \sigma^{2}_{s}}}$. 

## R code
### Load models and extract variance components
```{r, eval = F}
## LMM0 (null model)
LMM0 <- lm(logTime ~ sMass + Developmental*Parental*Immediate, data = dat)

## LMM1 (random intercept model)
LMM1 <- lmer(logTime ~ sMass + Developmental*Parental*Immediate+ (1|ID), data = dat)
### Verification
plot(LMM1)
qqnorm(residuals(LMM1))
qqline(residuals(LMM1))
qqnorm(ranef(LMM1)$ID[,1])
qqline(ranef(LMM1)$ID[,1])
### Extract variance
(Var.LMM1 <- data.frame(VarCorr(LMM1)))
### Extract confidence intervals
#(confint.LMM1 <- confint(LMM1, method="boot", nsim =2000,oldNames = FALSE))
#save(confint.LMM1, file="confint.LMM1")
load("confint.LMM1")

## LMM2 (random slope model)
LMM2 <- lmer(logTime ~ sMass + Developmental*Parental*Immediate + (1+ Immediate|ID), data = dat)
### Verification
plot(LMM2)
qqnorm(residuals(LMM2))
qqline(residuals(LMM2))
qqnorm(ranef(LMM2)$ID[,1])
qqline(ranef(LMM2)$ID[,1])
qqnorm(ranef(LMM2)$ID[,2])
qqline(ranef(LMM2)$ID[,2])
### Extract variance 
(Var.LMM2 <- data.frame(VarCorr(LMM2)))
### Extract confidence intervals
#(confint.LMM2 <- confint(LMM2, method="boot", nsim =2000,oldNames = FALSE))
#save(confint.LMM2, file="confint.LMM2")
load("confint.LMM2")
```

### Repeatability calculation
```{r, eval = F}
round(Var.LMM1[1,4]/(Var.LMM1[1,4]+Var.LMM1[2,4]),2)
round(confint.LMM1[1,1]^2/(confint.LMM1[1,1]^2+confint.LMM1[2,1]^2),2) ### lower CI
round(confint.LMM1[1,2]^2/(confint.LMM1[1,1]^2+confint.LMM1[2,2]^2),2) ### upper CI
```

### Likelihood ratio tests
```{r, eval = F}
(anova1 <- anova(LMM1,LMM0))
(anova2 <- anova(LMM2,LMM1,refit=FALSE))
```

## Table and Figure
### Table 2 - group diversity in personality and immediate plasticity
This code create a beautiful final table in word format with all variances and their confidence intervals, in addition to the results of the likelihood ratio tests
```{r, eval = F}
# function to calculate confidence interval of the residual variance of the LMM0 (null model)
V0 <- function(sigma, n, k, alpha = 0.05){
  # lm(y ~ x1 + x2)
  # n <- length(y)
  # k = number of parameters
  lower <- (n-(k+1))*sigma^2/qchisq(alpha/2, df = n-(k+1), lower.tail = FALSE)
  upper <- (n-(k+1))*sigma^2/qchisq(1-alpha/2, df = n-(k+1), lower.tail = FALSE)
  (temp <- format(round(c(sigma^2, lower,upper), digits = 3)))
  (temp <- paste0(temp[1], " [",temp[2],", ",temp[3],"]"))
}

# function to concatenate variance with confidence interval in brackets
V1 <- function(variance, confint){
  (temp <- c(variance, confint[1],confint[2]))
  (temp <- format(round(temp, digits = 3)))
  (temp <- paste0(temp[1], " [",temp[2],", ",temp[3],"]"))
  return(temp)
}

(Vresidual = c(V0(sigma(LMM0), length(dat$logTime),8),
               V1(Var.LMM1[2,4], confint.LMM1[2,]^2),
               V1(Var.LMM2[4,4], confint.LMM2[4,]^2)))

(Vintercept = c("",
                V1(Var.LMM1[1,4], confint.LMM1[1,]^2),
                V1(Var.LMM2[1,4], confint.LMM2[1,]^2)))

(Vslope = c("","", V1(Var.LMM2[2,4], confint.LMM2[3,]^2)))

(Coris = c("","", V1(Var.LMM2[3,5], confint.LMM2[2,])))

(`Chisqdf` = c("", format(round(c(anova1$Chisq[2], 
                                  anova2$Chisq[2]), digits = 2))))

(P = c("","<0.0001",format(round( anova2$`Pr(>Chisq)`[2], digits = 3))))


table <- data.frame(Model = c("LM0","LMM1","LMM2"),
           Vresidual = Vresidual,
           Vintercept = Vintercept,
           Vslope = Vslope,
           Coris = Coris,
           `Tested models` = c("","0 vs 1","1 vs 2"),
           `Chisqdf` = `Chisqdf`,
           P = P) %>% 
  flextable(.)%>%
  flextable::align(.,align = "left", part = "all") %>%
  hline(i =3,border = fp_border())%>%
  bg(i=2, bg = grey(0.9))%>%
  font(fontname = "Times", part = "all")%>%
  fontsize(size = 12, part = "all")%>%
  width(j = 1, width = conv_unit(1.6, "cm", "inch"))%>%
  width(j = 2:4, width = conv_unit(3.9, "cm", "inch"))%>%
  width(j = 5, width = conv_unit(4.3, "cm", "inch"))%>%
  width(j = 6:7, width = conv_unit(1.6, "cm", "inch"))%>%
  width(j = 8, width = conv_unit(1.8, "cm", "inch"))%>%
  bold( ~ P == "<0.0001", ~ P)
table

# run if you want to create the output in word format
#doc <- read_docx() %>%
#  body_end_section_continuous() %>% 
#  body_add_flextable(value = table) %>%
#  body_end_section_landscape() 
# print(doc, target = "Table 2.docx") 

rm(Vresidual,Vintercept,Vslope,Coris,Chisqdf,P,table)
```

### Figure 3 - group diversity in immediate plasticity
```{r, eval = F}
fig3 <- dat %>% 
  group_by(ID, Immediate2,Parental2,Developmental2) %>% 
  dplyr::summarise(Value = mean(logTime)) %>%
  ggplot(., aes(x = Immediate2, y = Value, group = ID)) +
  geom_line(alpha = 0.7) +
  xlab(xlab)+
  ylab("Espace behaviour (sec) \nSlow                                                       Fast")+
  scale_y_continuous(breaks = c(- 2.39794,-2.198657,-2,-1.799341), 
                     labels =c(250,158,100,63), 
                     limits = c(-2.48,-1.62))
fig3

# run if you want to create the output
#ggsave(fig3, filename = "Figure 3.pdf", height = 9, width = 9, units="cm")
```

***

# Effects of developmental and parental environments
## Model equation
### Linear mixed model 3 (LMM3, model with random intercepts for each combination of developmental and parental environments)

$-log10(\mathrm{Y}_{jkl}) = \beta_{0} + \beta_{1} \times x_{j}^{(W)} + \beta_{2} \times x_{jk}^{(I)} + \beta_{3} \times x_{j}^{(D)} + \beta_{4} \times x_{j}^{(P)}  + \gamma_{1} \times x_{jk}^{(I)} \times x_{j}^{(D)} + \gamma_{2} \times x_{jk}^{(I)} \times x_{j}^{(P)} + \gamma_{3} \times x_{j}^{(D)} \times x_{j}^{(P)} + \delta_{1} \times x_{jk}^{(I)} \times x_{j}^{(D)} \times x_{j}^{(P)} + \mathrm{U_{0}}_{j}^{\mathrm{CC}} \times \mathrm{I^{CC}} + \mathrm{U_{0}}_{j}^{\mathrm{CP}} \times \mathrm{I^{CP}} + \mathrm{U_{0}}_{j}^{\mathrm{PC}} \times \mathrm{I^{PC}} + \mathrm{U_{0}}_{j}^{\mathrm{PP}} \times \mathrm{I^{PP}} + \mathrm{E}_{jkl}$

$[\mathrm{E}_{jkl}] \stackrel{i.i.d.}{\sim}\mathcal{N}(0,\sigma^{2}_{res})$

$[\mathrm{U_{0}}_{j}^{\mathrm{CC}}] \stackrel{i.i.d.}{\sim}\mathcal{N}(0, \sigma^{2, \mathrm{CC}}_{i})$

$[\mathrm{U_{0}}_{j}^{\mathrm{CP}}] \stackrel{i.i.d.}{\sim}\mathcal{N}(0, \sigma^{2, \mathrm{CP}}_{i})$

$[\mathrm{U_{0}}_{j}^{\mathrm{PC}}] \stackrel{i.i.d.}{\sim}\mathcal{N}(0, \sigma^{2, \mathrm{PC}}_{i})$

$[\mathrm{U_{0}}_{j}^{\mathrm{PP}}] \stackrel{i.i.d.}{\sim}\mathcal{N}(0, \sigma^{2, \mathrm{PP}}_{i})$


\underline{Fixed effects:}

$\beta$ are parameters for main fixed effects, $\gamma$ are parameters for double-interaction terms and $\delta$ is the parameter for the triple-interaction term. 

$x_{j}^{(W)}$ is the snail total weight for the individual j.

$x_{j}^{(D)}$ is the developmental environment for the individual j with:

+ $x_{j}^{(D)}$ = -0.5 for unexposed snails
+ $x_{j}^{(D)}$ =  0.5 for predator-exposed snails

$x_{j}^{(P)}$ is the parental environment for the individual j with:

+ $x_{j}^{(P)}$ = -0.5 for snails from unexposed parents
+ $x_{j}^{(P)}$ =  0.5 for snails from predator-exposed parents

\underline{Random effects:}

For unexposed snails from predator-exposed parents, $\mathrm{U_{0}}_{i}^{\mathrm{PC}}$ is the random intercept value of the individual j.$\sigma^{2, \mathrm{PC}}_{i}$ is the variance in intercept, representing group diversity in personality when snails were not exposed to predator cues and from predator-exposed parents. The same logic applied for CC, CP and PP snails.   

$\mathrm{I^{CC}}$, $\mathrm{I^{CP}}$, $\mathrm{I^{PC}}$, $\mathrm{I^{PP}}$ are dummy variables:

+ $\mathrm{I^{CC}}$ = 1 if individual j was unexposed from unexposed parents, 0 otherwise
+ $\mathrm{I^{CP}}$ = 1 if individual j was predator-exposed from unexposed parents, 0 otherwise
+ $\mathrm{I^{PC}}$ = 1 if individual j was unexposed from predator-exposed parents, 0 otherwise
+ $\mathrm{I^{PP}}$ = 1 if individual j was predator-exposed from predator-exposed parents, 0 otherwise


## R code and Bayesian procedure
### Priors
We used a non-informative prior with an inverse-Wishart distribution (expected variance V = 1; degree of belief $\nu$ = 0.002; corresponding to a Gamma distribution with shape = 1 and scale = 1).

```{r, eval = F}
prior.normal <-list(R=list(V=1, nu=0.002), G=list(G1=list(V=1, nu=0.002),G2=list(V=1, nu=0.002),
                                                  G3=list(V=1, nu=0.002),G4=list(V=1, nu=0.002)))
```

### Models
We used a chain of 700 000 iterations with a thin of 100 iterations (1 out of 100 iterations kept) and a burnin of 200 000 iterations (the first 200 000 iterations were discarded). Posterior distribution of each fixed effects parameters and random variances were then made on the 5 000 remaining iterations after thinning and burning.

```{r, eval = F}
# Create dummy variables
dat$CC=as.numeric(dat$Parental==-0.5&dat$Developmental==-0.5)
dat$CP=as.numeric(dat$Parental==-0.5&dat$Developmental==0.5)
dat$PC=as.numeric(dat$Parental==0.5&dat$Developmental==-0.5)
dat$PP=as.numeric(dat$Parental==0.5&dat$Developmental==0.5)

# Run model

#LMM3 <-MCMCglmm(logTime ~ sMass + Immediate*Developmental*Parental,
#                random=~us(CC):ID + (CP):ID + us(PC):ID + us(PP):ID, 
#                rcov=~units, data=dat, prior=prior.normal, 
#                pr=F, nitt=700000, thin=100, burnin=200000)
# save(file = "LMM3",LMM3)
load("LMM3")
summary(LMM3)$solutions
summary(LMM3)$Gcovariances
```

### Check of MCMC algorithms

Convergence of MCMC algorithms was checked by: 

(1) running separately three other chains and analysing between-chain variance with Gelman & Rubin diagnostic implemented in the package coda (Plummer et al. 2006); 

```{r, eval = F}
#LMM3.chain2 <- MCMCglmm(logTime ~ sMass + Immediate*Developmental*Parental, 
#                        random=~us(CC):ID + us(CP):ID + us(PC):ID + us(PP):ID,
#                        rcov=~units, data=dat, prior=prior.normal, 
#                        pr=F, nitt=700000, thin=100, burnin=200000)
#save(LMM3.chain2 ,file = "LMM3.chain2")
load("LMM3.chain2")
#LMM3.chain3 <- MCMCglmm(logTime ~ sMass + Immediate*Developmental*Parental, 
#                        random=~us(CC):ID + us(CP):ID + us(PC):ID + us(PP):ID, 
#                        rcov=~units, data=dat, prior=prior.normal, 
#                        pr=F, nitt=700000, thin=100, burnin=200000)
#save(LMM3.chain3 ,file = "LMM3.chain3")
load("LMM3.chain3")
#LMM3.chain4 <- MCMCglmm(logTime ~ sMass + Immediate*Developmental*Parental, 
#                        random=~us(CC):ID + us(CP):ID + us(PC):ID + us(PP):ID, 
#                        rcov=~units, data=dat, prior=prior.normal, 
#                        pr=F, nitt=700000, thin=100, burnin=200000)
#save(LMM3.chain4 ,file = "LMM3.chain4")
load("LMM3.chain4")

# Convergence of fixed effects parameters
gelman.diag(list(LMM3$Sol,LMM3.chain2$Sol,LMM3.chain3$Sol,LMM3.chain4$Sol))
gelman.plot(list(LMM3$Sol,LMM3.chain2$Sol,LMM3.chain3$Sol,LMM3.chain4$Sol)) 

# Convergence of random effects parameters
gelman.diag(list(LMM3$VCV,LMM3.chain2$VCV,LMM3.chain3$VCV,LMM3.chain4$VCV))
gelman.plot(list(LMM3$VCV,LMM3.chain2$VCV,LMM3.chain3$VCV,LMM3.chain4$VCV))
```

(2) looking at the trace and shape of posterior distributions;

```{r, eval = F}
plot(LMM3$Sol)
plot(LMM3$VCV)
```

(3) looking at  95 % CI of cumulative values of iterations for each parameter and at each iteration;

```{r, eval = F}
cumuplot(LMM3$Sol)
cumuplot(LMM3$VCV)
```

(4) calculating effective size of our parameters (5 000 independent iterations for variances and at least 4 500 independent iterations for fixed effects parameters). 

```{r, eval = F}
effectiveSize(LMM3$Sol)
effectiveSize(LMM3$VCV)
```

### Sensitivity to prior

Sensitivity to prior was checked by ran two other chains with an improper prior (inverse-Wishart with V = 0 and $\nu$ = -2) and a flat prior (inverse-Wishart with V = 1 and $\nu$ = 0).

```{r, eval = F}
prior.flat <-list(R=list(V=1, nu=0), G=list(G1=list(V=1, nu=0),G2=list(V=1, nu=0),
                                            G3=list(V=1, nu=0),G4=list(V=1, nu=0)))
prior.improper <-list(R=list(V=1e-16, nu=-2), G=list(G1=list(V=1e-16, nu=-2),G2=list(V=1e-16, nu=-2),
                                                     G3=list(V=1e-16, nu=-2),G4=list(V=1e-16, nu=-2)))

#LMM3.flat <-MCMCglmm(logTime ~ sMass + Immediate*Developmental*Parental, 
#                     random=~us(CC):ID + us(CP):ID + us(PC):ID + us(PP):ID,
#                     rcov=~units, data=dat,prior=prior.flat,
#                     pr=F,nitt=700000, thin=100, burnin=200000)
#save(LMM3.flat ,file = "LMM3.flat")
load("LMM3.flat")

#LMM3.improper <-MCMCglmm(logTime ~ sMass + Immediate*Developmental*Parental, 
#                         random=~us(CC):ID + us(CP):ID + us(PC):ID + us(PP):ID,
#                         rcov=~units, data=dat,prior=prior.improper,
#                         pr=F,nitt=700000, thin=100, burnin=200000)
#save(LMM3.improper ,file = "LMM3.improper")
load("LMM3.improper")

summary(LMM3)$solutions
summary(LMM3.flat)$solutions
summary(LMM3.improper)$solutions

posterior.mode(LMM3$VCV)
posterior.mode(LMM3.flat$VCV)
posterior.mode(LMM3.improper$VCV)

gelman.diag(list(LMM3$Sol, LMM3.flat$Sol, LMM3.improper$Sol))
gelman.plot(list(LMM3$VCV, LMM3.flat$VCV, LMM3.improper$VCV))

plot(LMM3.flat$VCV)
plot(LMM3.improper$VCV)
cumuplot(LMM3.flat$VCV)
cumuplot(LMM3.improper$VCV)
effectiveSize(LMM3.flat$VCV)
effectiveSize(LMM3.improper$VCV)
```

## Post-hoc tests
### Supplementary material - contrasts for significance of group mean in immediate plasticity for each treatment
```{r, eval =F}
contrasts <- as.data.frame(emmeans(LMM3, pairwise ~ Immediate | Developmental*Parental, data=dat, 
                                   at = list(Immediate=c(-0.5,0.5),Developmental=c(-0.5,0.5),Parental=c(-0.5,0.5), sMass = 0))$contrasts)%>%
  mutate(Developmental = ifelse(Developmental==-0.5,"C","P"), 
         Parental = ifelse(Parental==-0.5,"C","P"),
         `Treatment` = paste0(Parental,Developmental),
         estimate = format(round(-estimate, digits = 3)),
         lower.HPD = format(round(-lower.HPD, digits = 3)),
         upper.HPD = format(round(-upper.HPD, digits = 3)))%>%
  mutate(Estimate = paste0(estimate, " (",upper.HPD,", ",lower.HPD,")"))%>%
  select(`Treatment`,Estimate) %>% 
  flextable() %>% 
  autofit() %>% 
  bg(i=c(2,4), bg = grey(0.9))
contrasts
```

### Contrasts to test the effects of parental and developmental environments on group mean in personality in each immediate environments

```{r, eval = F}
# Immediate env = C 
emmeans(LMM3,  pairwise~ Parental, data=dat, at = list(Immediate = -0.5, sMass = 0))$contrasts
emmeans(LMM3,  pairwise~ Developmental, data=dat, at = list(Immediate = -0.5, sMass = 0))$contrasts 
# Immediate env = P   
emmeans(LMM3,  pairwise~ Parental, data=dat, at = list(Immediate = 0.5, sMass = 0))$contrasts
emmeans(LMM3,  pairwise~ Developmental, data=dat, at = list(Immediate = 0.5, sMass = 0))$contrasts 
```

### Effects of developmental environment on group diversity in personality
```{r, eval =F}
# more group diversity in personality in the C developmental group than in the P developmental group? 
round(table(LMM3$VCV[,1]+ LMM3$VCV[,3] >= LMM3$VCV[,2] + LMM3$VCV[,4])/dim(LMM3$VCV)[1],3)

# more group diversity in personality in the C parental group than in the P parental group? 
round(table(LMM3$VCV[,1]+ LMM3$VCV[,2] >= LMM3$VCV[,3] + LMM3$VCV[,4])/dim(LMM3$VCV)[1],3)
```

## Table and figures
### Table 1 - effects on developmental and parental environments on group mean in personality and immediate plasticity
```{r, eval =F}
Fixed  <- summary(LMM3)$solutions %>% 
  as.data.frame(.)%>% round(., digits =3) %>% 
  format(., digits =3) %>% 
  mutate(`beta [95% CI]` = paste0(post.mean," [",`l-95% CI`,", ",`u-95% CI`,"]"),
         pMCMC =replace(pMCMC, pMCMC=="0.000","<0.001"),
         `Fixed effects` = c("Intercept","Mass","Immediate (Im.) ","Developmental (Dev.)",
                             "Parental (Par.)","Im. x Dev.","Im. x Par.","Dev. x Par.",
                             "Im. x Dev. x Par.")) %>% 
  select(`Fixed effects`,`beta [95% CI]`,pMCMC)
Fixed 

# Final
table1 <- flextable(Fixed) %>%
  flextable::align(align = "left", part = "all" )%>%
  bold( ~ pMCMC <= 0.05, ~ pMCMC) %>%
  bg(i=c(2,4,6,8), bg = grey(0.9))%>%
  font(fontname = "Times", part = "all")%>%
  fontsize(size = 12, part = "all")%>%
  width(j = 1, width = conv_unit(4.2, "cm", "inch"))%>%
  width(j = 2, width = conv_unit(4.3, "cm", "inch"))%>%
  width(j = 3, width = conv_unit(1.9, "cm", "inch"))
table1

# run if you want to create the output in word format
#doc <- read_docx()
#doc <- body_add_flextable(doc, value = table1)
#print(doc, target = "Table 1.docx")
```

### Figure 2 - effects of developmental and parental environments on group mean in personality and immediate plasticity

```{r, eval =F}
dodgi =0.2 # setting dodge between dots

fig2 <- dat %>% group_by(Immediate2, Developmental2, Parental2) %>% 
  dplyr::summarise(Temps = mean(logTime), SE = se(logTime)) %>%
  ggplot(., aes(x = Immediate2, y =Temps, 
                group = interaction(Parental2,Developmental2), 
                shape  = interaction(Parental2,Developmental2),
                fill = Parental2))+
    geom_line(aes(linetype= Parental2, color = Developmental2),
              size = 0.5, 
              position = position_dodge(width = dodgi))+
    geom_errorbar(aes(ymin=Temps-SE, ymax=Temps +SE), 
                  width= 0.1, size = 0.1, 
                  position = position_dodge(width = dodgi))+
    geom_point(aes(color = Developmental2), size = 4, position = position_dodge(width = dodgi))+
    xlab(xlab)+
    ylab("Espace behaviour (sec) \nSlow                                                       Fast")+
    scale_shape_manual(values =c(21,16,22,15))+
    scale_linetype_manual(values =c(2,1))+
    scale_fill_manual(values =c("white","black"))+
    scale_color_manual(values = c("#3182bd","#de2d26"))+
    scale_y_continuous(breaks = c(-2.158362,-2.079181,-2,-1.920819,-1.841638), 
                       labels =c(144,120,100,83,69), 
                       limits = c(-2.175,-1.825)) +
    guides(linetype = "none",shape = "none",fill = "none", color = "none")
fig2

# run if you want to create the output
# ggsave(fig2, filename = "Figure 2 input.pdf", height = 9, width = 9, units="cm")
```

### Figure 4 - effect of developmental environment on group diversity in personality

```{r, eval =F}
variance <- summary(LMM3)$Gcovariances %>% 
  as.data.frame()%>%
  mutate(history = c("CC","CP","PC","PP"), min = `l-95% CI`, max = `u-95% CI`)
variance

fig4 <- ggplot(variance,aes(x = history, y = post.mean, 
                            shape = history, 
                            fill = history, 
                            color = history))+
  geom_errorbar(aes(ymin = min, ymax = max), width = 0)+
    geom_point(size = 3)+
  scale_shape_manual(values =c(21,22,16,15))+
  scale_fill_manual(values =c("white","white","black","black"))+
  scale_color_manual(values = c("#3182bd","#de2d26","#3182bd","#de2d26"))+
  theme(legend.position = "none")+
  xlab("Combination of parental and\ndevelopmental environments")+
  ylab("Variance in escape behaviour")+
  ylim(0,0.05)
fig4

# run if you want to create the output
# ggsave(fig4, filename = "Figure 4.pdf", height = 9, width = 9, units="cm")
```
