[R] fixed effects constant in mcmcsamp

[Daniel Farewell]

I'm fitting a GLMM to some questionnaire data. The structure is J individuals,
nested within I areas, all of whom answer the same K (ordinal) questions. The
model I'm using is based on so-called continuation ratios, so that it can be
fitted using the lme4 package.

The lmer function fits the model just fine, but using mcmcsamp to judge the
variability of the parameter estimates produces some strange results. The
posterior sample is constant for the fixed effects, and the estimates of the
variance components are way out in the tails of their posterior samples.

The model I'm using says (for l = 1, ..., L - 1)

logit P(X[ijk] = l | X[ijk] >= l, U[i], V[j], W[k]) = U[i] + V[j] + W[k] + a[l]

where X[ijk] is the ordinal response to question k for individual j in area i.
The U, V, and W are random effects and the a's are fixed effects. Here's a
function to simulate data which mimics this setup (with a sequence of binary
responses Y[ijkl] = 1 iff X[ijk] = l):

sim <- function(n = c(10, 10, 10), sd = c(0.5, 2, 0.5), a = seq(-5, 1, 2)) {

 u <- rnorm(n[1], 0, sd[1])
 v <- rnorm(prod(n[1:2]), 0, sd[2])
 w <- rnorm(n[3], 0, sd[3])

 i <- factor(rep(1:n[1], each = prod(n[2:3])))
 j <- factor(rep(1:prod(n[1:2]), each = n[3]))
 k <- factor(rep(1:n[3], prod(n[1:2])))

 df <- NULL

 for(l in 1:length(a)) {
    
  y <- rbinom(length(i), 1, plogis(u[i] + v[j] + w[k] + a[l]))
  df <- rbind(df, data.frame(i, j, k, l, y))

  i <- i[!y]
  j <- j[!y]
  k <- k[!y]

 }

 df$l <- factor(df$l)
 return(df)
 
}

And here's a function which shows the difficulties I've been having:

test <- function(seed = 10, ...) {

 require(lme4)
 set.seed(seed)

 df <- sim(...)
 df.lmer <- lmer(y ~ 0 + l + (1 | i) + (1 | j) + (1 | k), family = binomial,
data = df) 
 df.mcmc <- mcmcsamp(df.lmer, 1000, trans = FALSE)

 print(summary(df.lmer))
 print(summary(df.mcmc))
 densityplot(~ as.numeric(df.mcmc) | factor(rep(colnames(df.mcmc), each =
1000)), scales = list(relation = "free"))

}

Running

test()

gives the following:

Generalized linear mixed model fit using PQL
Formula: y ~ 0 + l + (1 | i) + (1 | j) + (1 | k)
Data: df
 Family: binomial(logit link)
      AIC      BIC    logLik deviance
 2316.133 2359.166 -1151.066 2302.133
Random effects:
 Groups Name        Variance Std.Dev.
 j      (Intercept) 4.15914  2.03940
 k      (Intercept) 0.25587  0.50584
 i      (Intercept) 0.56962  0.75473
number of obs: 3455, groups: j, 100; k, 10; i, 10

Estimated scale (compare to 1)  0.8747598

Fixed effects:
   Estimate Std. Error  z value  Pr(>|z|)
l1 -4.50234    0.40697 -11.0632 < 2.2e-16 ***
l2 -3.27643    0.38177  -8.5821 < 2.2e-16 ***
l3 -1.05277    0.36566  -2.8791  0.003988 **
l4  0.76538    0.36832   2.0780  0.037706 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
   l1    l2    l3
l2 0.843
l3 0.841 0.900
l4 0.805 0.865 0.921
       l1               l2               l3               l4
 Min.   :-4.502   Min.   :-3.276   Min.   :-1.053   Min.   :0.7654
 1st Qu.:-4.502   1st Qu.:-3.276   1st Qu.:-1.053   1st Qu.:0.7654
 Median :-4.502   Median :-3.276   Median :-1.053   Median :0.7654
 Mean   :-4.502   Mean   :-3.276   Mean   :-1.053   Mean   :0.7654
 3rd Qu.:-4.502   3rd Qu.:-3.276   3rd Qu.:-1.053   3rd Qu.:0.7654
 Max.   :-4.502   Max.   :-3.276   Max.   :-1.053   Max.   :0.7654
     j.(In)          k.(In)           i.(In)           deviance
 Min.   :1.911   Min.   :0.0509   Min.   :0.06223   Min.   :2011
 1st Qu.:2.549   1st Qu.:0.1310   1st Qu.:0.19550   1st Qu.:2044
 Median :2.789   Median :0.1756   Median :0.25581   Median :2054
 Mean   :2.824   Mean   :0.2085   Mean   :0.29948   Mean   :2054
 3rd Qu.:3.070   3rd Qu.:0.2463   3rd Qu.:0.34640   3rd Qu.:2064
 Max.   :4.615   Max.   :0.8804   Max.   :3.62168   Max.   :2107

As you can see, the posterior samples from the fixed effects are constant (at
the inital estimates) and the estimates of the variance components aren't within
the IQ ranges of their posterior samples.

I understand from various posts that mcmcsamp is still a work in progress, and
may not work on every model. Is this one of those cases? I'm using R 2.3.1 and
lme4 0.995-2 on Windows XP.

Daniel Farewell
Cardiff University