MartÃ Casals writes:

The final data employed for the analysis were composed of 2187 cases (27 players * 81 games), having followed a filtering process. I dealt with a balanced study design with repeated measurements given that each player was observed the same number of games, and therefore the player was considered as a random effect. We carried out two mixed models, a Linear Mixed Model (LMM) through the nlme library (Pinheiro) and a Poisson mixed model with lme4 package (Bates).

```
library(lme4.0) ## development version ...
```

```
## Loading required package: Matrix
```

```
## Loading required package: lattice
```

```
library(plyr)
library(ggplot2)
theme_set(theme_bw())
library(reshape2)
set.seed(101)
```

```
simdata <- function(nplayers = 27, ngames = 81, sd.game = 1, sd.player = 1,
sd.obs = 0.5, beta0 = 0.5) {
d <- expand.grid(player = factor(seq(nplayers)), game = factor(seq(ngames)))
u.game <- rnorm(ngames, sd = sd.game)
u.player <- rnorm(nplayers, sd = sd.player)
u.obs <- rnorm(ngames * nplayers, sd = sd.obs)
d <- transform(d, obs = factor(seq(nrow(d))), eta = beta0 + u.game[game] +
u.player[player] + u.obs)
d$y <- rpois(nrow(d), exp(d$eta))
d
}
```

```
d1 <- simdata()
fit1 <- glmer(y ~ (1 | player) + (1 | game) + (1 | obs), data = d1, family = poisson)
```

```
## Number of levels of a grouping factor for the random effects is *equal* to
## n, the number of observations
```

```
unlist(VarCorr(fit1))
```

```
## obs game player
## 0.2345 0.9652 1.0449
```

```
simfun <- function(m = fit1, ...) {
sqrt(unlist(VarCorr(refit(m, simdata(...)$y))))
}
s0 <- simfun()
```

```
nsim <- 500
simres <- raply(nsim, simfun())
```

```
summary(simres)
```

```
## obs game player
## Min. :0.415 Min. :0.757 Min. :0.601
## 1st Qu.:0.482 1st Qu.:0.942 1st Qu.:0.881
## Median :0.495 Median :0.991 Median :0.966
## Mean :0.494 Mean :0.991 Mean :0.976
## 3rd Qu.:0.506 3rd Qu.:1.047 3rd Qu.:1.062
## Max. :0.572 Max. :1.214 Max. :1.353
```

```
library(proto)
source("geom-hboxplot.r")
m <- transform(melt(simres), Var2 = factor(Var2, levels = c("game", "player",
"obs")))
truevals <- data.frame(Var2 = c("game", "obs", "player"), value = c(1, 0.5,
1))
ggplot(m, aes(x = Var2, y = value)) + geom_boxplot(notch = TRUE) + facet_wrap(~Var2,
scale = "free") + geom_hline(aes(yintercept = value), data = truevals, lwd = 2,
alpha = 0.5, colour = "red") + geom_violin(alpha = 0.2, colour = NA, fill = "blue")
```

This looks reasonably good on the whole, although the estimates are *slightly* biased (i.e. notches don't include true values for `player`

and `obs`

); this could
be a finite-size/REML-y issue?

Would be nice to turn the whole thing sideways, but: `ggplot2 does not currently support free scales with a non-cartesian coord or coord_flip.`