(Italics represent corrected text)
You are making a 'mistake' in your model specification given what you say you want.
Random effects:
Groups Name Variance Std.Dev. Corr
Item (Intercept) 273.508 16.5381
Subject Gramgram 0.000 0.0000
Gramungram 3717.213 60.9689 NaN
Number1 59.361 7.7046 NaN -1.000
You see the numbers there under Corr for Subject? That shows that you are estimating correlations between the random slope of Gramgungram and Gramgram, Number1 and Gramgram, and Number1 and Gramungram by subject. If Gram was a numeric , you could eliminate the random correlation between Gram and Number1 with a model specified by:
m <- lmer(RT ~ Gram*Number + (1|Subject) + (0+Gram|Subject) (0+Number|Subject) + (1|Item),data= data)
You'll notice that any random effect specified in the same set of parentheses yields a random effect correlation. At least that is true for models without the / symbol, I'm not really familiar with the notion for lmer when that is in the mix.
However, given what we saw from a model where you estimated this parameter, I'd advise caution. Moreover, you'll probably note that my code above doesn't work for you.
EDIT
For those of you just now joining our program... for these examples I'll refer to primingHeid
as OP did in the comments, this dataset can be found in languageR
.
library(languageR)
library(lme4)
data(primingHeid)
Why doesn't my code work? It doesn't work because Gram is a factor. Think about it for a little bit... and look at your fixed effects. If a factor has two levels how many parameters must you estimate to explain its effects? Two. Of course, one of the parameters you estimate is the intercept. The interpretation of the intercept will depend on how your factors are coded. In treatment coding (the default in R), the intercept represents the value for a case when all variables are at level 1 (cf. a regression textbook for details of other contrasts). Regardless of your contrasts two parameters are estimated for two levels of a factor. I think what is happening is that when you fail to specify an intercept R is protecting you from yourself and going ahead and estimating two parameters anyway. Try summary(lm(RT ~ 0 + Condition,data=primingHeid))
and you'll see that it went right on ahead and estimated two parameters. So, back to the context of lmer... if you have a factor with two levels, R will gladly estimate two parameters and then correlate them all under the hood. Back again to your comments... estimate lmer(RT ~ Condition +(0+Condition|Subject),data=primingHeid)
and look and the ranef
of that model and you'll see yet again that this is exactly what R has done.
If you wanted to force R to stop that, you'd have to do the factor coding manually by turning Condition into a numeric. The assumptions you'd have to make about the mean value of RT when Condition was at the level you coded as 0 are likely untenable (i.e. that RT is really 0). I won't exclude the possibility that with some careful thought, transformation of the DV (centering on the mean of condition you are setting equal to 0?), and good model specification you might work your way somewhere that made some sense... but, that would be an entirely different question and I can't speak to it at the moment.
\EDIT
I think you probably should step back and think about your model structure a little more (which really is one of the great messages in Barr et al., 2013). Are items crossed with gram and number? How many items occur within a unique arrangement of gram and number per subject?
More general issues now...
I have a huge amount of respect for Barr (no surprises there). However, he is not entirely mainstream on issues related to fitting this type of model. That isn't a bad thing ... but time will tell whether his approach for these models will become the next big thing. I have little doubt that 'keeping it maximal' is great if your data will tolerate it. But sometimes it won't. The backwards selection procedure he published involving the use of non-converged models is a bit unexpected. However, I have to admit now that I've seen his appendices I'm a little less sour on the idea than I was when I first read it. All the same, I'd like to see it be vetted a bit more.
You'll note Barr specifically does not use pvals.fnc() for models that have random correlations. So, only having skimmed the published version of his paper, I'd guess you can only use it under his approach if you can backwards step to a point where you don't have any.
Going now to my training with other stats gurus I feel compelled to say that almost all of this worry is an exercise in p-value fetishism that may be entirely misplaced - especially if you consider that this level of nested decision making yields a test that has a definition that is difficult to define.
This is an old question, but I ran into a similar situation when trying to run a gamm
model using bam
.
I believe what's going on is that s(prev, speaker, bs='re')
in the model above is estimating a variance term for the ~prev:speaker-1
interaction (random slopes), whereas s(time, by=speaker, bs="fs", m=1)
is estimating a separate slope term for each speaker. These slopes do not share a variance term, meaning that they aren't really a random effect (and aren't necessarily centred on zero).
Here's a demonstration of the differences between the two methods, using gamm
to compare them. These simulated data are from a normal distribution, but this should work for whatever distribution you're using:
#Idea: y = [b0 + boi] + x[b1 + b1i] + s(z) + error,
#
# where b0i and b1i are random intercepts and slopes, and s(z) is some curvy spline
library(mgcv)
#Generate some fake data
set.seed(123)
ndat <- 900 #900 data points
ngroup <- 30 #30 groups
gInd <- sample(rep(1:30,length.out=900)) #Group index
xVal <- runif(ndat,-1,1) #Range of x values
zVal <- runif(ndat,-pi,pi) #Range of z values
bo <- -1 #Intercept
b1 <- 2 #Slope of x
b0i <- rnorm(ndat/ngroup,0,0.5) #Random intercepts
b1i <- rnorm(ndat/ngroup,0,0.5) #Random slopes
s_z <- sin(zVal) #Sine curve to estimate by spline
err <- rnorm(ndat,0,0.2) #Noise
yVal <- (bo + b0i[gInd]) + (b1 + b1i[gInd])*xVal + s_z + err #Generate values
plot(xVal,yVal) #Looks OK
plot(zVal,yVal)
dat <- data.frame(y=yVal,x=xVal,z=zVal,group=factor(gInd)) #Assemble into dataframe
str(dat)
#Help file for s(x,bs='re')
# ?smooth.construct.re.smooth.spec
#GAMM MODEL
mod1 <- gamm(y~1+x+s(z,bs='cr'),random=list(group=~1+x),data=dat,
control=list(msVerbose=TRUE))
summary(mod1$lme) #Estimates variance terms well, slope terms well
summary(mod1$gam)
plot(mod1$gam); curve(sin(x),-pi,pi,col='red',add=T) #Pretty good estimation of sine curve
#BAM MODEL
#Using "by" notation within s()
mod2a <- bam(y~1+x+s(z,bs='cr')+s(group,bs='re')+s(x,by=group,bs='re'),data=dat)
gam.vcomp(mod2a) #Note the number of sd terms estimated -- 1 for each group
summary(mod2a)
#Using paired notation
mod2b <- bam(y~1+x+s(z,bs='cr')+s(group,bs='re')+s(x,group,bs='re'),data=dat)
gam.vcomp(mod2b) #Fewer sd terms estimated
summary(mod2b)
#Compare estimates of intercepts and slopes to actual values:
ranInt <- data.frame(actual=b0i,gamm=ranef(mod1$lme)[[2]][,1],
bamA=coef(mod2a)[grepl('s\\(group\\).',names(coef(mod2a)))],
bamB=coef(mod2b)[grepl('s\\(group\\).',names(coef(mod2b)))])
ranSlope <- data.frame(actual=b1i,gamm=ranef(mod1$lme)[[2]][,2],
bamA=coef(mod2a)[grepl('s\\(x\\).',names(coef(mod2a)))],
bamB=coef(mod2b)[grepl('s\\(x,group\\).',names(coef(mod2b)))])
#Good estimates of random intercept terms
pairs(ranInt,upper.panel=function(x,y) {points(x,y); abline(0,1,col='red')},
lower.panel=function(x,y) text(mean(x),mean(y),round(cor(x,y),3),cex=0.1+cor(x,y)*2),
main='Intercepts')
#bamA slope estimates are offset by a constant value when using "by"
pairs(ranSlope,upper.panel=function(x,y) {points(x,y); abline(0,1,col='red')},
lower.panel=function(x,y) text(mean(x),mean(y),round(cor(x,y),3),cex=0.1+cor(x,y)*2),
main='Slopes')
Best Answer
I eventually figured this one out with much help from Doing Bayesian Data Analysis (Kruschke) and Data Analysis Using Regression and Multilevel/Hierarchical Models (Gelman). This model gives varying intercepts and slopes and the correlation between them.
y = dependent variable
sLvl = participant id at each data point
aLvlx = between subjects factor for each id
NaLvl = Number of levels for the between subject's factor
Ntotal = total length of data in long form