[R-sig-ME] conditional repeatability from MCMCglmm with random slope
Conor Michael Goold
conor.goold at nmbu.no
Mon Sep 11 11:52:56 CEST 2017
Ah, I see. It's my understanding that you can just take the range of values for the covariate, and apply the formula below. I'm not sure what you mean by the 'covariate specific random effect variance'...the variance of the random effects (intercept and slope variance terms) do not change with the covariate, but the value of the ICC will change depending on which value of the covariate it is estimated. If X is some covariate, I'd just do something like:
#########
library(rethinking) # for the HPDI interval function
# get the unique values of your covariate X
covariate_values <- unique( X )
# assuming V_intercept, V_slope, Cov_intercept_slope and Vr are vectors of posterior probability distributions from the fitted model
ICC_vector <- sapply( covariate_values,
function(z) {
numerator <- V_intercept + V_slope * z^2 + 2 * Cov_intercept_slope * z
denominator <- numerator + Vr
numerator / denominator
}
)
# ICC_vector is now a matrix with dimensions N (the number of posterior draws) x P (the length of covariate_values)
# create data frame to hold ICC summary statistics across the values of covariate_values
ICC_df <- data.frame(covariate_values = covariate_values,
ICC_mu = apply(ICC_vector , 2 , mean ) ,
hdi_low = apply(ICC_vector , 2 , function(z) HPDI(z,0.95) )[1,] ,
hdi_high = apply(ICC_vector , 2 , function(z) HPDI(z,0.95) )[2,]
)
##############
>From there, you can take the average of the ICC values in the ICC_df data frame.
Hope this makes sense,
Conor
________________________________________
From: John Morrongiello <john.morrongiello at unimelb.edu.au>
Sent: Monday, September 11, 2017 11:06 AM
To: Conor Michael Goold; r-sig-mixed-models at r-project.org
Subject: RE: conditional repeatability from MCMCglmm with random slope
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
More information about the R-sig-mixed-models
mailing list