#Solution 2 attempt
set.seed(101)
dd <- data.frame(x1= rnorm(1920), x2=rnorm(1920), x3=rnorm(1920), x4=rnorm(1920),
                 treatment = factor(rep(1:2, each=3)),
                 replicate = factor(rep(1:3, each=1)),
                 stage = factor(rep(1:5, each=384)),
                 country = factor(rep(1:4, each=96)),
                 plot = factor(rep(1:10, each=24)),
                 chamber = factor(rep(1:6, each=1)),
                 n = 1920)

library(lme4)
dd$y <- simulate(~ x1 + x2 + x3 + (1|plot),
                 family = binomial,
                 weights = dd$n,
                 newdata = dd,
                 newparams = list(beta = c(1,1,1,1),
                                theta = 1))[[1]]

# my real response variable 'y' has a poisson distribution, but I had difficulty figuring
# out how to simulate a poisson distribution so I left the bionomial.

m0 <- lmer(y~ x1 + x2 + x3 + x4 + treatment*replicate*stage + (1|chamber) + (1|country/plot),
            data=dd,
            na.action = "na.fail",
           REML = F,
            lmercontrol = glmerControl(optimizer="bobyqa"))

nullmodel <- MuMIn:::.nullFitRE(m0)
dredge(m0, m.lim = c(0,5), rank = "AIC", extra =list(R2 = function(x)     {
  r.squaredGLMM(x, null = nullmodel)["delta", ]}))