[R-sig-ME] function for simulating from GLMM

[Paul Johnson]

Dear R mixed modellers,

I've written a function, sim.glmm, to simulate from some GLMMs (Gaussian, binomial, Poisson and negative binomial). I've been using this function for a while, and thought some of you might be interested (I think there was a request for such a function recently). I've use for power analysis and to check bias, CI coverage, type I error, residuals distributions, etc.

The code is here:
https://www.dropbox.com/s/b0h1cptgq3ummjq/functions_for_glmm_power_simulations_v4_16August2013_draft.R

I'm not aware of other R functions that do this, although I guess they must be out there (although I've recently learned that MLwiN can link with R to simulate from GLMMs - I haven't tried this). I'm aware that simulate-mer will simulate from fitted GLMMs, but not, afaik, from a set of parameter values.

Below is an example of how it works. Any comments or suggestions would be greatly appreciated.

Paul Johnson
Glasgow

  # simulate counts of tick burden on grouse chicks in a single year from a Poisson-lognormal GLMM,
  # loosely based on: 
  #   Elston et al. (2001). 
  #   Analysis of aggregation, a worked example: numbers of ticks on red grouse chicks. Parasitology, 122, 563–9.
  #   http://abdn.ac.uk/lambin-group/Papers/Paper%2041%202001%20Elston%20Tick%20aggregation%20Parasitology.pdf
  # chicks are nested within broods, and broods within locations

    # simulate data set that defines sampling of chicks within broods within locations, 
    # assuming 3 chicks per brood and 2 broods per location. N locations = 20.
      
      tickdata <- expand.grid(chick=1:3, brood=1:2, location=1:20)

    # make brood and chick ids unique (otherwise random effects will be simulated as crossed, not nested)

      tickdata$location <- factor(paste("loc", tickdata$location, sep=""))
      tickdata$brood <- factor(paste(tickdata$location, "brd", tickdata$brood, sep=""))
      tickdata$chick <- factor(paste(tickdata$brood, "chk", tickdata$chick, sep=""))
            
    # simulate tick counts with an average burden of 5 ticks per chick 
    # random effect variances are 2, 1 and 0.3 for location, brood and chick respectively
    
      tickdata<-
        sim.glmm(design.data = tickdata, 
          fixed.eff = list(intercept = log(5)),
          rand.V = c(location = 2, brood = 1, chick = 0.3),
          distribution = 'poisson')
      
    # plot counts and fit GLMM

      plot(response~jitter(as.numeric(location),factor=0.5),pch=21,bg=as.numeric(brood),data=tickdata)
      library(lme4)
      glmer(response~(1|location)+(1|brood)+(1|chick),family='poisson',data=tickdata)

    # adding categorical or continuous fixed effects is straightforward, e.g 
    # if tickdata has altiude (in km and centred on zero) and sex: 
    # fixed.eff = list(intercept = log(5), altitude.km = log(2), sex=log(c(M = 1, F = 0.8)))
    # random effects can have non-zero covariances if rand.V is supplied as a variance-covariance matrix