visreg: plotting random effects from an lme model

[pbreheny]

Communicated via e-mail by Denis Infanger:

Recently, I wanted to visualize the fitted values for a linear mixed effects model that I fitted using nlme. I fitted factor-specific random effects variances (random = list(ID = pdDiag(~g))). Then, I wanted to plot the random effects for each ID using visreg(mod, "xvar", by = "ID", level = 1). Unfortunatley, this gave an error. The reason (I suspect) is that visreg can't extract the data using from this specification of the random effects.

I managed to code a workaround by manually changing the call to something arbitrary, so that the data extraction with get_all_vars() works (fully reproducible example below).

Would this functionality (i.e. plotting from lme with complex random effects) be something you could include in a future release of visreg? Or is this too complex to include?

#-----------------------------------------------------------------------------
# Load packages
#-----------------------------------------------------------------------------

library(lme4)
library(visreg)
library(nlme)

#-----------------------------------------------------------------------------
# Simulate data (taken from Ben Bolker on: http://rpubs.com/bbolker/factorvar)
#-----------------------------------------------------------------------------

set.seed(101)
dsd <- c(1,2,3)
nblock <- 25
ntot <- 750
nfac <- 3
d <- expand.grid(f=letters[1:nfac],g=1:nblock,rep=seq(ntot/(nblock*nfac)))
d$fg <- with(d,interaction(g,f))
delta <- 4:6
u <- rnorm(nblock*nfac,mean=0,sd=rep(dsd,each=nblock))
d$y <- rnorm(nrow(d),mean=delta[d$f],sd=0.2)+u[d$fg]
d$g <- factor(d$g)

#-----------------------------------------------------------------------------
# Fit mixed model using "random slopes" for each group
#-----------------------------------------------------------------------------

mod <- lme(y~f-1, random=list(g = pdDiag(~f)), data=d)

#-----------------------------------------------------------------------------
# Try to call "visreg": Doesn't work
#-----------------------------------------------------------------------------

visreg(mod, "f", by = "g", level = 1)

#-----------------------------------------------------------------------------
# Manually set "random" formula to arbitrary value, as class "call"
#-----------------------------------------------------------------------------

mod$call$random <- as(~1|g, "call")
class(mod$call$random)

#-----------------------------------------------------------------------------
# Works now!
#-----------------------------------------------------------------------------

v <- visreg(mod, "f", by = "g", level = 1)
sub_g <- sample(levels(d$g), 10, replace = FALSE)
vv <- subset(v, g %in% sub_g)
plot(vv, ylab="g", layout=c(5,2))