# [R] fixed effects constant in mcmcsamp

From: Daniel Farewell <farewelld_at_Cardiff.ac.uk>
Date: Tue 08 Aug 2006 - 19:58:12 EST

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, 0, sd)
v <- rnorm(prod(n[1:2]), 0, sd)
w <- rnorm(n, 0, sd)

i <- factor(rep(1:n, each = prod(n[2:3])))
```
j <- factor(rep(1:prod(n[1:2]), each = n))  k <- factor(rep(1:n, 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

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

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help