3
$\begingroup$

Let's say I have a data set (with a binary DV) which I used a generalized linear mixed-effects model on. I am currently trying to run a simulation-based power analysis to determine the sample size I should try to collect for a replication attempt.

I know about the simulate function for fitted lme4 model objects, which will create simulated responses that follow the size and structure of the actual dataset.

My question is, how can I use this fitted model object to run an a-priori power analysis for my replication?. In other words, is there an easy way to extract the necessary information from the simulate function, vary the sample size, then to simulate responses?

EDIT:

Here is a simple example with ten subjects that should give me a good starting place. The data here are totally fake so the effects are not close to what I have in the pilot data.

The two variables of interest are age in months (between subjects) and whether the stimulus is "similar" or "dissimilar" (within subjects). Each individual is exposed to three items within each similarity group, for a total of six binary responses for each subject. As seen in the model, I am interested in a possible interaction between age and similarity and include a random slope of similarity to do this.

Let me know if I need to modify anything/provide more information to make things easier.

set.seed(28)

fakeData <- data.frame(ID = rep(1:10, each = 6),
                       months = rep(round(runif(10, min = 60, max = 80)), each = 18),
                       similarity = rep(c("Sim", "Dissim"), each = 3),
                       response = sample(c(0, 1), 60, replace = TRUE))

fakeModel <- glmer(response ~ months * similarity + (similarity | ID), family = binomial, data = fakeData)
summary(fakeModel)
$\endgroup$
  • 1
    $\begingroup$ yes: use the simulate function from the development version of lme4. If you give a small reproducible example I'll try to respond. $\endgroup$ – Ben Bolker Jan 17 '14 at 21:31
  • $\begingroup$ I included an example. Let me know if you need anything else. Thanks @BenBolker! $\endgroup$ – dmartin Jan 18 '14 at 1:21
3
$\begingroup$

Generate data:

set.seed(28)
fakeData <- data.frame(ID = rep(1:10, each = 6),
          months = rep(round(runif(10, min = 60, max = 80)), each = 18),
          similarity = rep(c("Sim", "Dissim"), each = 3),
          response = sample(c(0, 1), 60, replace = TRUE))

Fit initial model:

library(lme4)
fakeModel <- glmer(response ~ months * similarity + (similarity | ID), 
                   family = binomial, data = fakeData)
summary(fakeModel)

Now I'm going to write a generic function that replicates your data-generation code above. (I'm a little bit nervous about the way you generated your data by implicit replication of the vectors -- I would prefer to explicitly write out the expected amount of replication for each column. I'm not sure that the code below will work properly if you change the numbers ...)

datfun <- function(nindiv=10,nperindiv=6,
                   nsim=3, monthmin=60, monthmax=80) {
      data.frame(ID=factor(rep(1:nindiv, each=nperindiv)),
      months = rep(round(runif(nindiv, min=monthmin, max=monthmax)),
                    each=nsim*nperindiv),
      similarity = rep(c("Sim","Dissim"), each=nsim),
      response = sample( 0:1, nindiv*nperindiv, replace=TRUE))
}

Once you've done this, however, it's easy (using the development version of lme4) to simulate new responses with the same parameters as the initial model (or to change the parameters as you choose):

simdat <- simulate(formula(fakeModel), newdat=datfun(),
                   newparam=list(theta=getME(fakeModel,"theta"),
                                 beta=getME(fakeModel,"beta")),
                   family=binomial)

The result is a data frame with 1 column. If you want to do this multiple times per data set size you may want to set nsim larger ...

$\endgroup$
  • $\begingroup$ After downloading and calling the development version of lme4 I get: Error in UseMethod("simulate"): no applicable method for 'simulate' applied to an object of class "formula." It runs when I just include the model object (i.e., fakeModel) rather than the formula, but this still simulates the number of responses from the original dataset. So if I have nindiv = 20, simdat still gives me a dataframe the same size as fakeData. $\endgroup$ – dmartin Jan 19 '14 at 19:13
  • $\begingroup$ Additionally, I'm having some problems getting models to converge with the development version (model failing to converge, computed covariance matrix is not non positive-definite). Any help would be appreciated, thanks! $\endgroup$ – dmartin Jan 19 '14 at 19:34
  • $\begingroup$ hmm. For comment #1, how did you install? Was lme4 loaded? I would recommend devtools::install_github("lme4","lme4") if you can manage it. For #2: I'm having a lot of convergence problems too (perhaps slightly mitigated by using the very latest version). I suspect that this fake data set may simply be too small/unstable for reliable fitting? $\endgroup$ – Ben Bolker Jan 19 '14 at 21:10
  • $\begingroup$ Problem #1 was because I was using the whole model, not just the right side, so that's solved. Problem #2 is still there. Both this fake example and my real data fit fine with the CRAN version of lme4. However, they both yield non positive definite matrices with the development version. I also tried the example in the help page of simulate.merMod (with cbpp data), where I simulated a larger N, and it's still giving me convergence problems. Was the optimizer/default iterations number changed? I can start a new question (on Overflow) with my modified example if that's easier. $\endgroup$ – dmartin Jan 20 '14 at 3:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.