[R-sig-ME] conditional repeatability from MCMCglmm with random slope

Mon Sep 11 08:41:03 CEST 2017

Hi John, 

Repeatability for models with a random slope are a bit tricky to understand since it is calculated with respect to a particular value of a covariate. In general, it can be calculated as:

R = Vg / (Vg + Vr) 
   =  V_intercept + V_slope * X_i^2 + 2 * Cov( V_intercept + V_slope ) * X_i / (numerator + Vr)

where V_intercept = variance of intercepts, V_slope = variance of slopes, X_i is the covariate at a particular value i, Cov represents the covariance, and Vr = the residual variance. I'm not sure it can be calculated with rptR. 

For more details, see Martin et al. 2011. http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2010.00084.x/abstract (particularly pages 371 - 372)

Best regards
Conor Goold
PhD Student
Phone:        +47 67 23 27 24

Norwegian University of Life Sciences
Campus Ås. www.nmbu.no

________________________________________
From: R-sig-mixed-models <r-sig-mixed-models-bounces at r-project.org> on behalf of John Morrongiello <john.morrongiello at unimelb.edu.au>
Sent: Monday, September 11, 2017 5:26 AM
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] conditional repeatability from MCMCglmm with random slope

Dear list
I would like to estimate conditional repeatability of a behavioural trait from a model including random slopes fit with MCMCglmm. Would someone have some tips for how this can be done? The rptR package offers this option using bootstrapping, based on Johnson's 2014 paper<https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4368045/> (equation 11). Here, the average repeatability is estimated across the distribution of a covariate in question. I can readily estimate raw (intercept only) and a conditional (common slope) repeatability from a MCMCglmm model, but I'm not sure how to get the random slopes repeatability. Any help/ advice is much appreciated.

(Johnson, P.C.D. (2014). Extension of Nakagawa & Schielzeth's R2GLMMRGLMM2 to random slopes models. Methods in Ecology and Evolution 5: 944-946).

library(rptR)
library(MCMCglmm)

#######some data#######
behaviour <-
structure(list(fid = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L,
3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 8L, 8L,
8L, 9L, 9L, 9L, 10L, 10L, 10L, 11L, 11L, 11L, 12L, 12L, 12L,
13L, 13L, 13L, 14L, 14L, 14L, 15L, 15L, 15L, 16L, 16L, 16L, 17L,
17L, 17L, 18L, 18L, 18L), .Label = c("2", "3", "5", "6", "7",
"8", "9", "10", "11", "12", "13", "15", "16", "17", "18", "19",
"20", "21"), class = "factor"), week = c(1L, 2L, 3L, 1L, 2L,
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L,
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L,
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L,
3L), trait1 = c(0.225903614, 0.368017058, 0.805194805, 0.691799861,
0.654561365, 0.815674805, 0.477368175, 0.222362125, 0.09480315,
0.268691229, 0.500809061, 0.262361128, 0.365236524, 0.090701498,
0.344081069, 0.851842461, 0.477558348, 0.324364209, 0.425583266,
0.085414302, 0.311696658, 0.083493708, 0.164316797, 0.104364194,
0.156442474, 0.145193036, 0.237337192, 0.445346658, 0.377818182,
0.256871307, 0.287724075, 0.334800469, 0.324204293, 0.285696594,
0.199901004, 0.227835821, 0.709611452, 0.202676634, 0.084690774,
0.288277218, 0.282670647, 0.408824789, 0.241985203, 0.159210526,
0.0565915, 0.105827529, 0.039312489, 0.049929224, 0.292300806,
0.098316799, 0.15465606, 0.175269499, 0.271719515, 0.230084728
)), .Names = c("fid", "week", "trait1"), row.names = c(NA, -54L
), class = "data.frame")

############# Repeatability from rptR package #####

###intercept only model
rep1 <- rpt(trait1 ~ (1 | fid), grname = "fid", data = behaviour,  datatype = "Gaussian", nboot = 1000, npermut = 0)
###common week slope
rep2 <- rpt(trait1 ~ week+ (1 | fid), grname = "fid", data = behaviour, datatype = "Gaussian", nboot = 1000, npermut = 0)
###individual differences in week slope
rep3 <- rpt(trait1 ~ week+ (week | fid), grname = "fid", data = behaviour,  datatype = "Gaussian", nboot = 1000, npermut = 0)

##raw repeatability: ignore temporal change
print(rep1)
##conditional repeatability 1: shared change over time
print(rep2)
##conditional repeatability 2: change over time that differs among individuals (this is what I want to get from the equivalent MCMCglmm model)
print(rep3)

########### mcmcGLMM code ####
prior1.1 <- list(G = list(G1 = list(V = 1, nu = 0.002)), R = list(V = 1, nu = 0.002))

###intercept only model
MCMC_raw<-MCMCglmm(trait1 ~1,
             random=~fid,family=("gaussian"), prior=prior1.1,
             nitt=110000,thin=25,burnin=10000, data=behaviour, verbose=TRUE)

###common week slope
MCMC_cond1<-MCMCglmm(trait1 ~week,
             random=~fid,family=("gaussian"), prior=prior1.1,
             nitt=110000,thin=25,burnin=10000, data=behaviour, verbose=TRUE)

###common week slope
MCMC_cond2<-MCMCglmm(trait1 ~week,
             random=~us(week):fid,family=("gaussian"), prior=prior1.1,
             nitt=110000,thin=25,burnin=10000, data=behaviour, verbose=TRUE)

####repeatability formula
R= Vg/(Vg + Vr), where R= repeatability, Vg is group level variance and Vr is data level (residual) variance

###MCMCglmm derived raw repeatability
rep4<-(MCMC_raw$VCV[,"fid"]/ (MCMC_raw$VCV[,"fid"]+MCMC_raw$VCV[,"units"]))
posterior.mode(rep4)
HPDinterval(rep4)

###MCMCglmm derived conditional repeatability 1: shared change over time
rep5<-(MCMC_cond$VCV[,"fid"]/ (MCMC_ cond $VCV[,"fid"]+MCMC_ cond $VCV[,"units"]))
posterior.mode(rep5)
HPDinterval(rep5)

####How would I estimate conditional repeatability 2: change over time that differs among individuals using MCMC_cond2?

Cheers
john
--
Dr. John R. Morrongiello
School of BioSciences
University of Melbourne
Victoria 3010, Australia
T: +61 3 8344 8929
M: +61 403 338 554
E: john.morrongiello at unimelb.edu.au<mailto:%20jmorrongiell at unimelb.edu.au>
W: morrongiellolab.com<http://morrongiellolab.com/>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models