Dear list
I have a mixed model with crossed random effects and a random slope ((x|ID) + (1|year)) for which I'd like to produce plots of predicted responses with CIs by continuous fixed effects (akin to plots returned by plotLMER.fnc in the languageR library). Unfortunately, pvals.fnc does not seem to support a random slope within group with correlated intercept (x|ID) that my model has so plotLMER.fnc won't work. I saw on the DRAFT r-sig-mixed-models FAQ <http://glmm.wikidot.com/faq>  webpage a method for generating similar plots using the orthodont data, but I couldn't get this to work for crossed random effects or for the partial effects of multiple continuous fixed effects (orthodont example has continuous and categorical).  Would someone be able to help me out with this?

Cheers
John

I've provided an example below:
#####example dataset
nID<-20
nwithin<-10
prosw<-nID*nwithin
growth<-vector(length=prosw)
age<-seq(2,11, length=nwithin)
a<-rnorm(prosw,mean=2)
b<-120;c<-10;d<-5;e<- -2;rho<--0.9
for (i in 1:nID)
{
sf<-0.1*rnorm(1);sc<-2*rnorm(1);ss<-rnorm(1,mean=2);err<-0.1*rnorm(1)
for (j in 1:nwithin)
{
growth[nwithin*(i-1)+j]<-sf +  b + (c+sc)*a[nwithin*(i-1)+j] - (d+ss)*age[j] +sf*a[nwithin*(i-1)+j]+
                   err
err<-err*rho + 0.1*rnorm(1)
}
}
a<-seq((1980+sample(0:10,1)),length.out=10,by=1)
b<-seq((1980+sample(0:10,1)),length.out=10,by=1)
c<-seq((1980+sample(0:10,1)),length.out=10,by=1)
d<-seq((1980+sample(0:10,1)),length.out=10,by=1)
e<-seq((1980+sample(0:10,1)),length.out=10,by=1)
f<-seq((1980+sample(0:10,1)),length.out=10,by=1)
g<-seq((1980+sample(0:10,1)),length.out=10,by=1)
h<-seq((1980+sample(0:10,1)),length.out=10,by=1)
i<-seq((1980+sample(0:10,1)),length.out=10,by=1)
j<-seq((1980+sample(0:10,1)),length.out=10,by=1)
k<-seq((1980+sample(0:10,1)),length.out=10,by=1)
l<-seq((1980+sample(0:10,1)),length.out=10,by=1)
m<-seq((1980+sample(0:10,1)),length.out=10,by=1)
n<-seq((1980+sample(0:10,1)),length.out=10,by=1)
o<-seq((1980+sample(0:10,1)),length.out=10,by=1)
p<-seq((1980+sample(0:10,1)),length.out=10,by=1)
q<-seq((1980+sample(0:10,1)),length.out=10,by=1)
r<-seq((1980+sample(0:10,1)),length.out=10,by=1)
s<-seq((1980+sample(0:10,1)),length.out=10,by=1)
t<-seq((1980+sample(0:10,1)),length.out=10,by=1)
year<-c(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t)##sorry couldn't get for loop to work for this
age<-rep(age, nID)
ID<-rep(seq(1,nID), each=nwithin)
growthdata<-data.frame(ID,age,year,growth)
growthdata$ID<-as.factor(growthdata$ID)
temp<-data.frame(temp=c(12,12.3,13.2,12,14,15.2,13,14,15,17,15.2,14.7,16,16.2,15,17.7,18.2,16,18),year=seq(1980,1998,by=1))
growthdata<-merge(growthdata,temp)

###I illustrate where I've got to using three models
M1<-lmer(growth~age+temp+(age|ID)+(1|as.factor(year)),growthdata)###model I'd like to produce plots for
M2<-lmer(growth~age+temp+(0+age|ID)+(1|ID)+(1|as.factor(year)),growthdata)
M3<-lmer(growth~age+temp+(age|ID),growthdata)

##error when use pvals.fnc with M1
modelmcmc<- pvals.fnc(M1, nsim=1000, withMCMC=TRUE)
##returns:
Error in pvals.fnc(model.lmer, nsim = 1000, withMCMC = TRUE) :
MCMC sampling is not yet implemented in lme4_0.999375
for models with random correlation parameters

###using M2 with uncorrelated random effects, plotLMER.fnc produces similar graphs to what I want
modelmcmc<- pvals.fnc(M2, nsim=1000, withMCMC=TRUE)
par(mfrow=c(2,1), mar=c(5,5,1,1))
plotLMER.fnc(M2, mcmcMat=modelmcmc$mcmc, fun=exp, ylabel = "growth",pred='age')
plotLMER.fnc(M2, mcmcMat=modelmcmc$mcmc, fun=exp, ylabel = "growth",pred='temp')

###modification of orthodont example on DRAFT r-sig-mixed-models FAQ <http://glmm.wikidot.com/faq> webpage using M3, couldn't work out how to do this with crossed random effects from M1
newdat <- expand.grid(
    age=seq(2,11,by=1)
    , temp=seq(12,18.2,by=.1)
    , growth = 0
)
mm = model.matrix(terms(M3),newdat)
newdat$growth = mm %*% fixef(M3)
pvar1 <- diag(mm %*% tcrossprod(vcov(M3),mm))
tvar1 <- pvar1+VarCorr(M3)$ID[1]
newdat <- data.frame(
    newdat
    , plo = newdat$growth-2*sqrt(pvar1)
    , phi = newdat$growth+2*sqrt(pvar1)
    )
#plot with CIs (ideally, I would like a plot with just one curve for each partial main effect, these treated as continuous as in plotLMER.fnc)
ggplot(newdat)+
layer(
    geom = 'point'
    , mapping = aes(
        x = temp
        , y = growth
        , colour = age
    )
)+
layer(
    geom = 'errorbar'
    , mapping = aes(
        x = temp
        , ymin = plo
        , ymax = phi
        , colour = age
    )
)





	[[alternative HTML version deleted]]


