[R] Adding Nested Random Effects to MCMCglmm

Marian L. Firke mfirke1 at swarthmore.edu
Sat Apr 12 21:29:37 CEST 2014


Dear R Experts, 

Does anyone have advice for how to nest random effects in the MCMCglmm package? 

Right now, I am running my models with Individual as a random effect; however, the 25 individuals in my study are birds taken from 13 different nests, so I feel that a more accurate model would nest the variable Individual (called RNR in my code) within the variable Nest ID. 

I am trying to run both univariate and bivariate mixed models using the MCMCglmm package. (My interest is in calculating repeatability estimates using the univariate models, and in partitioning out the within-individual and between-individual correlations using the bivariate models.) I am also trying to include a fixed effect (SMI) but am not sure that I am doing this correctly for either case. 



Example Univariate MM: 
Cort0.mcmc<-MCMCglmm(CORT0_LOG10_Z~SMI_CEN, random=~indiv, family="gaussian", data=cort0.data,verbose=FALSE)

CORT0_LOG10_Z is the variable of interest (which has been measured repeatedly for each given individual) 
RNR is individual ID 
SMI_CEN is a fixed effect that varies both within and between individuals

Example piece of data frame: 

  RNR HPA CORT0_LOG10_Z  SMI_CEN  NEST_ID
1  AR04273   1   -3.29844317 1.236170087  3
2  AR04273   2   -0.51328002 0.122035902  3
3  AR04281   1    0.83631319 1.212914325  8


Example Bivariate MM: 
cort0_cort15.bivarC.mcmc<-MCMCglmm(cbind(CORT0_std,CORT15_std)~(SMI_cen), random=~us(trait):RNR, rcov=~us(trait):units, family=c("gaussian","gaussian"), prior=prior.bivar, nitt=1300000,thin=1000,burnin=300000, data=cort0_cort15,verbose=FALSE)

CORT0_std and CORT15_std are the 2 variables of interest
RNR is individual ID 
SMI_cen is a fixed effect that varies both within and between individuals

Example piece of data frame: 
RNR      SMI  SMI_avg      SMI_dev  CORT0 CORT15      SMI_cen SMI_avg_cen   CORT0_std   CORT15_std  NEST_ID
1  AR04273 18.03640 17.48987  0.546522102  0.460  4.590  1.236170087  0.73332960 -3.80922656 -3.170123623  3
2  AR04273 16.94335 17.48987 -0.546522102  3.493 13.604  0.122035902  0.73332960 -0.04862913 -0.480133012  3
3  AR04281 15.63367 15.86610 -0.232427804  4.875 13.507 -1.212914325 -1.05393597  0.56974795 -0.497849822  8

I am basing my code off of the code supplied in the 2013 paper by Dingemanse and Dochtermann about applying linear mixed modelling to behavioral ecology. 

Does anyone have advice for how to nest RNR within Nest ID, as well as insight into whether I have coded the fixed effect SMI correctly? 


Many thanks, 
Marian




More information about the R-help mailing list