[R-sig-ME] design matrix for R2

[Ben Bolker]

  try getME(mF,"X") (and more generally see ?getME)

On 16-05-01 03:47 PM, Rob Griffin wrote:
> Dear list,
> I am trying to implement the R2 method of Nakagawa and Schielzeth for a model in lme4 but having problems. Using the documentation provided with the paper (including R script and dummy data: http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210x.2012.00261.x/suppinfo) I run in to the same problem as when I try it on my own data. The full script from supplementary file S4 is pasted at the end of this message for the first example in the paper (I've slightly modified to make self contained dummy data and dropping the "Condition" variable from the model, which is absent from the dummy data file [S1] and not mentioned in the paper). The issue arises when running:
>> Fixed <- fixef(mF)[2] * mF at X[, 2] + fixef(mF)[3] * mF at X[, 3] + fixef(mF)[4] * mF at X[, 4]Error: no slot of name "X" for this object of class "lmerMod"
> My suspicion is that there has been a change in the structure of the models in lmer such that "mF at X" no longer relates to the design matrix, but reading the latest lme4 documentation I can't solve the issue... any suggestions? I'm using R version 3.2.4 (2016-03-16), and lme4 version 1.1-12.
> Many thanks,Robert
> -----
> # Clear memoryrm(list = ls())
> # Read body length data (Gaussian, available for both sexes)Data <- data.frame(rep(1:12, each = 80), rep(1:120, each = 8), rep(c("Female","Male"), each = 8), rep(c("A","B"), each = 4), rep(c("Cont","Exp")), rep(sample(c(9,10,11), replace=T, size = 120), each = 8)+rnorm(960,0,0.5))colnames(Data) = c("Population","Container","Sex","Habitat","Treatment","BodyL")
> # Fit null model without fixed effects (but including all random effects)m0 <- lmer(BodyL ~ 1 + (1 | Population) + (1 | Container), data = Data)
> # Fit alternative model including fixed and all random effectsmF <- lmer(BodyL ~ Sex + Treatment + (1 | Population) + (1 | Container), data = Data)
> # View model fits for both modelssummary(m0)summary(mF)
> # Extraction of fitted value for the alternative model# fixef() extracts coefficents for fixed effects# mF at X returns fixed effect design matrixFixed <- fixef(mF)[2] * mF at X[, 2] + fixef(mF)[3] * mF at X[, 3] + fixef(mF)[4] * mF at X[, 4]
> # Calculation of the variance in fitted valuesVarF <- var(Fixed)
> # An alternative way for getting the same resultVarF <- var(as.vector(fixef(mF) %*% t(mF at X)))
> # R2GLMM(m) - marginal R2GLMM# Equ. 26, 29 and 30# VarCorr() extracts variance components# attr(VarCorr(lmer.model),'sc')^2 extracts the residual varianceVarF/(VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + attr(VarCorr(mF), "sc")^2)
> # R2GLMM(c) - conditional R2GLMM for full model# Equ. XXX, XXX(VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1])/(VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + (attr(VarCorr(mF), "sc")^2))
> # AIC and BIC needs to be calcualted with ML not REML in body size modelsm0ML <- lmer(BodyL ~ 1 + (1 | Population) + (1 | Container), data = Data, REML = FALSE)mFML <- lmer(BodyL ~ Sex + Treatment + (1 | Population) + (1 | Container), data = Data, REML = FALSE)
> # View model fits for both models fitted by MLsummary(m0ML)summary(mFML)
>  		 	   		  
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>