Skip to content

Random effects in clmm() of package ordinal

3 messages · Malcolm Fairbrother, Christian Brauner, Jarrod Hadfield

#
Hi Christian,

The ordinal package (which is otherwise very handy) does not allow for
random slopes (only random intercepts), so I don't think you can have
tested what you think you tested using that package.

You could try the MCMCglmm package instead, which allows for ordinal models
*and* random slopes.

Regarding random slopes more generally, Barr et al.'s (2013) paper "Random
effects structure for confirmatory hypothesis testing: Keep it maximal"
shows pretty definitively that not allowing for random slopes can often be
anticonservative. So if there's a gap between the p values you get with and
without random slopes, I'd be more inclined to trust the value *with*
random slopes.

Hope that's useful.

Cheers,
Malcolm

  
  
#
Hi Malcolm,

Interesting. I just read the reference manual for the ordinal package (v.
July 2, 2014) and indeed at the beginning it states "random slopes (random
coefficient models) are not yet implemented" (p. 3). If this i indeed true
then I'm puzzled by some things:
(a) later on he explains the usage of the extractor function VarCorr() for
class "clmm" and states "[...] a model can contain a random intercept (one
column) or a random intercept and a random slope (two columns) for the
same grouping factor" (p. 52).
(b) at first glance this seems inconsistent with the output of clmm():

library(ordinal)
fm2 <- clmm(rating ~ temp + contact + (temp | judge),
            data = wine,
            Hess = TRUE)

Now look at

summary(fm2)

Random effects:
 Groups Name        Variance Std.Dev. Corr
 judge  (Intercept) 1.15608  1.0752
        tempwarm    0.02801  0.1674   0.649
Number of groups:  judge 9

and

ranef(fm2)
$judge
  (Intercept)    tempwarm
1  1.61166859  0.20664909
2 -0.49269459 -0.06317359
3  0.93254466  0.11957142
4 -0.08571746 -0.01099074
5  0.13893821  0.01781474
6  0.44710565  0.05732815
7 -1.77974875 -0.22820043
8 -0.26589478 -0.03409318
9 -0.51044493 -0.06544955

that looks like a random intercept and random slope model.

What I accidently did in my models (illustrating here with his example) is
not to replace "temp" in the fixed part but in the ranom part:
fm2 <- clmm(rating ~ temp + contact + (temp | judge),
            data = wine,
            Hess = TRUE)
fm2a <- clmm(rating ~ temp + contact + (1 | judge),
             data = wine,
             Hess = TRUE)
anova(fm2a, fm2)

The anova output is fairly similar to mine; the p is extremely high as
well.

All the standard stuff works without a warning, e.g.:
fm2b <- clmm(rating ~ temp + contact + (1 | judge) + (temp + 0 | judge),
             data = wine,
             Hess = TRUE)

So what is happening here? Anybody got an idea?

Best,
Christian
On Fri, Aug 29, 2014 at 09:41:26AM -0700, Malcolm Fairbrother wrote:
#
Hi Christian,

I think clmm is dealing with the random 'slopes' correctly (or at  
least it gives the same estimates as MCMCglmm) when temp is  
categorical (see below). Without knowing much about clmm I can imagine  
it might fail when temp is truly continuous because the variance of  
each latent variable varies as a function of the covariate and so (in  
the probit case) the cdf of the normal (i.e. the probit link) would  
need to be evaluated using different standard deviations.

Cheers,

Jarrod

fm2 <- clmm(rating ~ temp + contact + (temp-1 | judge),
? ? ? ? ? ? data = wine, link="probit",
? ? ? ? ? ? Hess = TRUE)

# note that I have changed temp|judge to temp-1|judge because I think  
it makes more sense with categorical variables.  Also I have used  
probit link in order to compare it to MCMCglmm's "threshold".

prior2=list(R=list(V=1, fix=1), G=list(G1=list(V=diag(2), nu=2,  
alpha.V=diag(2)*100, alpha.mu=c(0,0))))

fm2.mcmc <- MCMCglmm(rating ~ temp + contact, random=~us(temp):judge,  
data = wine,family="threshold", prior=prior2, nitt=13000*10,  
thin=10*5, burnin=3000*10, verbose=FALSE)


par(mfrow=c(2,2))
hist(fm2.mcmc$Sol[,2], breaks=50)
abline(v=coef(fm2)[5], col="red")
hist(fm2.mcmc$Sol[,3], breaks=50)
abline(v=coef(fm2)[6], col="red")
hist(fm2.mcmc$VCV[,1], breaks=50)
abline(v=VarCorr(fm2)$judge[1], col="red")
hist(fm2.mcmc$VCV[,4], breaks=50)
abline(v=VarCorr(fm2)$judge[4], col="red")




Quoting Christian Brauner <christianvanbrauner at gmail.com> on Fri, 29  
Aug 2014 19:20:42 +0200: