[R-sig-ME] variance explained by fixed effects in MCMCglmm

[rafter sass ferguson]

For any who comes looking:

I asked this same question on researchgate here:
https://www.researchgate.net/post/How_can_I_calculate_R2_for_an_Bayesian_MCMC_multilevel_model

Happily, it was answered by Shinichi Nakagawa himself (author of the
original article), who provided updated
code based on the same example models but using MCMCglmm.
I'm pasting in his response below:

Hi, Rafter
>
> I have added MCMCglmm to our supplementary R code. Note that I have only
> done it for the Gaussian model (see below) - it comes with credible
> intervals for R2 as well. The next week or so I will do this for binary and
> Poisson models (binary model one is a bit complicated). Once finished, I
> will upload it to my website (http://www.i-deel.org/
> <https://www.researchgate.net/go.Deref.html?url=http%3A%2F%2Fwww.i-deel.org%2F>)
> linking it with the paper.
>
> Thanks
>
> Shinichi
>
> ####################################################
>
> # A. Preparation
>
> ####################################################
>
> # Note that data generation appears below the analysis section.
>
> # You can use the simulated data table from the supplementary files to
> reproduce exactly the same results as presented in the paper.
> # Set the work directy that is used for rading/saving data tables
>
> # setwd("/Users/R2")
>
>
> # load R required packages
>
> # If this is done for the first time, it might need to first download and
> install the package
>
> # install.packages("arm")
>
> library(arm)
>
> # install.packages("lme4")
>
> # the verson 1.0-5
>
> library(lme4)
>
> # install.packages("MCMCglmm")
>
> library(MCMCglmm)
> ####################################################
>
> # B. Analysis
>
> ####################################################
> # 1. Analysis of body size (Gaussian mixed models)
>
> #---------------------------------------------------
> # Clear memory
>
> rm(list = ls())
> # Read body length data (Gaussian, available for both sexes)
>
> Data <- read.csv("BeetlesBody.csv")
> # Fit null model without fixed effects (but including all random effects)
>
> m0 <- lmer(BodyL ~ 1 + (1 | Population) + (1 | Container), data = Data)
> # MCMCglmm model
> mm0<-MCMCglmm(BodyL ~ 1 , random = ~ Population + Container, data=Data)
> # Fit alternative model including fixed and all random effects
>
> mF <- lmer(BodyL ~ Sex + Treatment + Habitat + (1 | Population) + (1 |
> Container), data = Data)
> # MCMCglmm model
> mmF<-MCMCglmm(BodyL ~ Sex + Treatment + Habitat , random = ~ Population +
> Container, data=Data)
> # View model fits for both models
>
> summary(m0)
>
> summary(mm0)
> summary(mF)
>
> summary(mmF)
> # Extraction of fitted value for the alternative model
>
> # fixef() extracts coefficents for fixed effects
>
> # mF at pp$X returns fixed effect design matrix
>
> Fixed <- fixef(mF)[2] * mF at pp$X[, 2] + fixef(mF)[3] * mF at pp$X[, 3] +
> fixef(mF)[4] * mF at pp$X[, 4]
>
>
> # MCMCglmm (it is probably better to get a posterior distribuiton of R2
> rather than getting each varaince component - we do this below as an
> alternative)
>
> mFixed <- mean(mmF$Sol[,2]) * mmF$X[, 2] + mean(mmF$Sol[, 3]) * mmF$X[, 3]
> + mean(mmF$Sol[ ,4]) * mmF$X[, 4]
>
>
> # Calculation of the variance in fitted values
>
> VarF <- var(Fixed)
>
> mVarF<- var(mFixed)
>
>
> # An alternative way for getting the same result
>
> VarF <- var(as.vector(fixef(mF) %*% t(mF at pp$X)))
>
> mVarF <- var(as.vector(apply(mmF$Sol,2,mean) %*% t(mmF$X)))
> # R2GLMM(m) - marginal R2GLMM
>
> # Equ. 26, 29 and 30
>
> # VarCorr() extracts variance components
>
> # attr(VarCorr(lmer.model),'sc')^2 extracts the residual variance
>
> VarF/(VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] +
> attr(VarCorr(mF), "sc")^2)
>
>
> # MCMCglmm - marginal
> mVarF/(mVarF+sum(apply(mmF$VCV,2,mean)))
>
>
> # alternative with crebile intervals
> vmVarF<-numeric(1000)
>
> for(i in 1:1000){
>
> Var<-var(as.vector(mmF$Sol[i,] %*% t(mmF$X)))
>
> vmVarF[i]<-Var}
>
>
> R2m<-vmVarF/(vmVarF+mmF$VCV[,1]+mmF$VCV[,2]+mmF$VCV[,3])
>
> mean(R2m)
>
> posterior.mode(R2m)
>
> HPDinterval(R2m)
> # R2GLMM(c) - conditional R2GLMM for full model
>
> # Equ. 30
>
> (VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1])/(VarF +
> VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + (attr(VarCorr(mF),
> "sc")^2))
>
>
> # MCMCglmm - conditional
>
> (mVarF+sum(apply(mmF$VCV,2,mean)[-3]))/(mVarF+sum(apply(mmF$VCV,2,mean)))
> # alternative with crebile intervals
>
>
> R2c<-(vmVarF+mmF$VCV[,1]+mmF$VCV[,2])/(vmVarF+mmF$VCV[,1]+mmF$VCV[,2]+mmF$VCV[,3])
>
> mean(R2c)
>
> posterior.mode(R2c)
>
> HPDinterval(R2c)
>




Rafter Sass Ferguson, MS
PhD Candidate | Crop Sciences Department
University of Illinois in Urbana-Champaign
liberationecology.org
518 567 7407

On Thu, Jul 16, 2015 at 1:05 AM, rafter sass ferguson <
liberationecology at gmail.com> wrote:

> I've been searching for ways to calculate some R^2-like statistics for a
> multi-level multi-response model fit with MCMCglmm (with 3 Gaussian
> responses).
>
> I can see from the Course Notes and elsewhere that it's very
> straightforward to calculate the variance explained by the random effects,
> but after much searching I haven't found a discussion for fixed effects. It
> looks like one option would be refitting with all my fixed effects as
> random - but I'm concerned that might be bonkers, or that there might be an
> easier way.
>
> For previous multilevel modeling with lmer, I've used the approach
> following Nakagawa et al. 2013 (A general and simple method for obtaining
> R2 from generalized linear mixed-effects models. Methods in Ecology and
> Evolution, 4(2), 133–142. http://doi.org/10.1111/j.2041-210x.2012.00261.x)
> to calculate conditional and marginal R^2...
> but I'm not sure how to adapt it for the (for me) brave new world of
> MCMCglmm.
>
> I would be grateful for any advice!
>
> Thanks so much for your time.
>
> Warmly,
> Rafter
>
>
> Rafter Sass Ferguson, MS
> PhD Candidate | Crop Sciences Department
> University of Illinois in Urbana-Champaign
> liberationecology.org
> 518 567 7407
>

	[[alternative HTML version deleted]]