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

[Paul Johnson]

Hi, 

I've moved sim.glmm from Dropbox to github: 

https://github.com/pcdjohnson/sim.glmm.git

The dropbox copy will no longer be updated and may disappear at some point.

Paul



On 30 Sep 2013, at 13:11, Paul Johnson <Paul.Johnson at glasgow.ac.uk> wrote:

> Hi,
> 
> Thanks to those who emailed with feedback on the sim.glmm function. New version here:
> https://www.dropbox.com/s/41xeopjsxt3eerc/functions_for_glmm_power_simulations_v4.2_30September2013.R
> The main change is that I've fixed an error that occurred when supplying a single random effect variance as a vector. 
> 
> Paul
> 
> ________________________________________
> From: r-sig-mixed-models-bounces at r-project.org [r-sig-mixed-models-bounces at r-project.org] On Behalf Of Paul Johnson [paul.johnson at glasgow.ac.uk]
> Sent: 12 September 2013 16:44
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] function for simulating from GLMM
> 
> 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
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models