[prev in list] [next in list] [prev in thread] [next in thread] 

List:       r-sig-mixed-models
Subject:    Re: [R-sig-ME] function for simulating from GLMM
From:       Paul Johnson <paul.johnson () glasgow ! ac ! uk>
Date:       2013-09-30 12:11:34
Message-ID: DAB55EFEBDAD1F4888EDB2E8F5180BC95E065070EF () CMS03 ! campus ! gla ! ac ! uk
[Download RAW message or body]

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@r-project.org \
[r-sig-mixed-models-bounces@r-project.org] On Behalf Of Paul Johnson \
                [paul.johnson@glasgow.ac.uk]
Sent: 12 September 2013 16:44
To: r-sig-mixed-models@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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
_______________________________________________
R-sig-mixed-models@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


[prev in list] [next in list] [prev in thread] [next in thread]