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

Mon Sep 11 11:06:42 CEST 2017

Hi Conor
Thanks for getting back to me. I'll have a close read of the Martin etal paper (I'm familiar with the Goldstein etal 2002 paper they cite in this section). In regards to rptR being able to calculate a conditional repeatability involving random slopes, they provide an example of this towards the end of their vignette (https://cran.r-project.org/web/packages/rptR/vignettes/rptR.html). The trick (which is a little beyond my coding ability) is to properly average repeatability estimates across the range of the covariate as indicated in your formula below from MCMCglmm. I've had a look at the source code for rpt and I can't see where this function is estimating each covariate specific random effect variance to then calculate the mean random effect variance.

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
W: morrongiellolab.com

-----Original Message-----
From: Conor Michael Goold [mailto:conor.goold at nmbu.no] 
Sent: Monday, 11 September 2017 4:49 PM
To: John Morrongiello <john.morrongiello at unimelb.edu.au>; r-sig-mixed-models at r-project.org
Subject: Re: conditional repeatability from MCMCglmm with random slope

Hi again John, 

In the repeatability formula I worte, it should be Cov( V_intercept, V_slope ) with an addition sign!

Conor
________________________________________
From: Conor Michael Goold
Sent: Monday, September 11, 2017 8:41 AM
To: John Morrongiello; r-sig-mixed-models at r-project.org
Subject: Re: conditional repeatability from MCMCglmm with random slope

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. http://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: http://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