# Power/assessment of Poisson GLMMs

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.`