Dear mixed-effects modellers,

I've been testing the importance of a multilevel factor (species, with 6
levels) in a binary mixed-effects model. The (simplified) hypothesis is
that survival varies between species. To test it, I fit a model with
species in it, and then without it, compared the log likelihoods as a test
statistic.

Because log-likelihood ratios are very likely not to be Chi-squared
distributed, I computed an approximate p-value using a simulation approach.
Briefly, I refit the model without species, simulated a response from this
(using simulate) and then refit both models to the simulated data and then
computed the log-likelihood ratios of the models. I repeated the
simulations several times (c. 99) and by comparing the observed and
simulated chisq values got a rough p-value. The results look consistent
with those of an anova(mod1, mod0), although a little more conservative as
I would expect.

However, I was somewhat surprised that many of the simulations resulted in
higher likelihoods for the null model (i.e. the simulated "chi-sq' is < 0).
This got me wondering about whether I'm doing something wrong, but my
searches couldn't find anything on this on the list (but maybe I'm not
looking in the right places).

I've tried to create a repeatable example below. The two functions I wrote
for this,  ChiSim and simAnova are below.

Many thanks in advance for any help.

Best wishes
Robi


## try functions on some test data (functions are below)
library(lme4)
library(mlmRev)
data(Contraception)

fm1 <- glmer(use ~ urban+age+livch+(1|district), Contraception, binomial)
fm0 <- update(fm1, ~.-livch)
## note this next step takes some time - hence the small number of
iterations
test <- simAnova(mod1=fm1, mod0=fm0, nsim=9, use.u=FALSE)

summary(test$simChi) ## often contains negative values



## function to simulate Chisq
ChiSim <- function(mod1, mod0, mod0sim, use.u) {
  converge <- FALSE         ## error trap to deal with non-convergence
issues
  while(converge==FALSE)
    {
      newY <- unlist(simulate(mod0sim, 1, use.u=use.u)) ## simulate new Y
      newmod1 <- try(refit(mod1, newY), silent=TRUE) ## refit H1 model to
null data
      newmod0 <- try(refit(mod0, newY), silent=TRUE) ## refit H0 model to
null data
      ## rerun if either model fit exited abnormally
      if((class(newmod1) == 'try-error') |
         (class(newmod0) == 'try-error'))
        converge==FALSE
      ## accept iteration if model fits converged
      else
        {
          if(length(newmod1@optinfo$warnings) == 0 &
             length(newmod0@optinfo$warnings) ==0)
            converge <- TRUE
        }
    }
  ## calculate the logLik ratio for simulated data
  simChi <- 2*(logLik(newmod1)[1] - logLik(newmod0)[1])
  ##return(list(simChi=simChi, nmod0=newmod0, nmod1=newmod1))
  return(simChi)
}


simAnova <- function(mod1, mod0, nsim, REMLsim=TRUE, use.u=FALSE){
  ## Should the REML fit be used for simulating "new data"
  if(REMLsim)
    mod0sim <- mod0

  ## refit models as ML to allow comparison of fixed effects
  if(isREML(mod1)){
    mod1 <- update(mod1, REML=FALSE)
    mod0 <- update(mod0, REML=FALSE)
  }
  ## if we want to use the ML fit for simulation, then use the updated ML
null model
  ## doesn't change anything if the original model wasn't REML
  if(!REMLsim)
    mod0sim <-  mod0
  ##calculate the observed logLik ratio
  obsChi <- 2*(logLik(mod1)[1] - logLik(mod0)[1])
  ## simulate distribution of Chisq
  sims <- replicate(nsim,
                    ChiSim(mod1=mod1, mod0=mod0, mod0sim=mod0sim,
use.u=use.u),
                    simplify=TRUE)

  ## compare observed Chi to simulated null distribution
  p <-  (sum(sims > obsChi)+1)/(nsim +1)
  ## extract dfs of test

  df <- as.numeric(getME(mod1, 'devcomp')$dims['p'] -
                   getME(mod0, 'devcomp')$dims['p'])

  return(list(ChiSq=obsChi, df=df,  p=p, simChi=sims))}


## Session info in case this is helpful
R version 3.0.2 (2013-09-25)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
 [1] parallel  grid      splines   stats     graphics  grDevices utils
 [8] datasets  methods   base

other attached packages:
 [1] mlmRev_1.0-5      verification_1.35 CircStats_0.2-4   boot_1.3-9
 [5] fields_6.8        maps_2.3-6        spam_0.40-0       waveslim_1.7.1
 [9] survival_2.37-4   sandwich_2.3-0    lmtest_0.9-32     zoo_1.7-10
[13] car_2.0-19        arm_1.6-10        MASS_7.3-29       reshape_0.8.4
[17] plyr_1.8          ggplot2_0.9.3.1   lme4_1.0-5        Matrix_1.1-0
[21] lattice_0.20-24

loaded via a namespace (and not attached):
 [1] abind_1.4-0        AER_1.2-1          bitops_1.0-6       caTools_1.16

 [5] cluster_1.14.4     coda_0.16-1        colorspace_1.2-4
compiler_3.0.2
 [9] dichromat_2.0-0    digest_0.6.3       Formula_1.1-1      gdata_2.13.2

[13] gplots_2.12.1      gtable_0.1.2       gtools_3.1.1       Hmisc_3.13-0

[17] KernSmooth_2.23-10 labeling_0.2       lmerTest_2.0-3     minqa_1.2.1

[21] munsell_0.4.2      nlme_3.1-113       nnet_7.3-7
numDeriv_2012.9-1
[25] pbkrtest_0.3-8     proto_0.3-10       RColorBrewer_1.0-5
reshape2_1.2.2
[29] scales_0.2.3       snow_0.3-13        stringr_0.6.2      tools_3.0.2


-- 
Dr Robert Bagchi
Ecosystem Management Group
Institute for Terrestrial Ecosystems
Universitätstrasse 16
8092 Zurich
Switzerland

tel: +41 (0)44 63 26805

email: bagchi.r@gmail.com (preferred)
robert.bagchi@usys.ethz.ch

web: http://www.ecology.ethz.ch/people/obass/rbagchi

	[[alternative HTML version deleted]]

