pitakakariki / simr

Power Analysis of Generalised Linear Mixed Models by Simulation
69 stars 19 forks source link

Can we pass 'allow.new.levels' to powerSim? #178

Open JCruk opened 4 years ago

JCruk commented 4 years ago

I am trying to conduct a power analysis on a glmer model where the levels of the grouping variable (random intercept) are not fully crossed.

This results in an error when running powerSim: NAs are not allowed in prediction data for grouping variables unless allow.new.levels is TRUE

Is there a way to pass allow.new.levels = TRUE along to predict.glmer from the powerSim call?

pitakakariki commented 4 years ago

I can have a look at this but I'll need more information. Can you make a reproducible example?

JCruk commented 4 years ago

Great! Thank you!

Here is an example that reproduces the error:

# List of needed packages
Pkgs <- c("tidyverse", "lme4", 
          "lmerTest", "simr")

# Load packages
lapply(Pkgs, require, c = T)

b <- c(-1.18641157630157, 0.45567541002807, -1.60121022845497, -2.61708155295718, 
       -0.382435480880911, 0.161296354182718, 0.485326753074071, 0.781761834854902, 
       1.01764589807192)

V1 <- 0.107061360865083

s <- 0.828580261209386

Data <-
  structure(
    list(
      Condition = structure(
        c(
          1L,
          2L,
          2L,
          1L,
          2L,
          1L,
          1L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          1L,
          1L,
          1L,
          2L,
          1L,
          2L,
          1L,
          1L,
          1L,
          1L,
          1L,
          2L,
          1L,
          2L,
          2L,
          1L,
          2L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          2L,
          1L,
          2L,
          2L,
          2L,
          2L
        ),
        .Label = c("Control", "Intervention"),
        class = "factor"
      ),
      Group = structure(
        c(
          5L,
          1L,
          2L,
          7L,
          8L,
          3L,
          11L,
          9L,
          8L,
          9L,
          1L,
          2L,
          2L,
          9L,
          5L,
          11L,
          3L,
          8L,
          7L,
          1L,
          3L,
          7L,
          5L,
          11L,
          15L,
          16L,
          4L,
          NA,
          10L,
          14L,
          13L,
          6L,
          15L,
          4L,
          6L,
          14L,
          14L,
          10L,
          NA,
          10L,
          13L,
          12L,
          16L
        ),
        .Label = c(
          "1",
          "3",
          "4",
          "6",
          "8",
          "10",
          "12",
          "13",
          "17",
          "19",
          "20",
          "29",
          "33",
          "36",
          "40",
          "41"
        ),
        class = "factor"
      ),
      Session = structure(
        c(
          3L,
          3L,
          2L,
          3L,
          3L,
          2L,
          2L,
          2L,
          1L,
          1L,
          1L,
          1L,
          3L,
          3L,
          2L,
          3L,
          3L,
          2L,
          2L,
          2L,
          1L,
          1L,
          1L,
          1L,
          2L,
          3L,
          3L,
          NA,
          3L,
          2L,
          2L,
          3L,
          1L,
          1L,
          1L,
          1L,
          3L,
          2L,
          NA,
          1L,
          1L,
          1L,
          1L
        ),
        .Label = c("Instruction",
                   "Initial", "Follow-up"),
        class = "factor"
      ),
      Facilitator = structure(
        c(
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          1L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L,
          2L
        ),
        .Label = c("Facilitator 1",
                   "Facilitator 2"),
        class = "factor"
      )
    ),
    class = "data.frame",
    row.names = c(NA,-43L)
  )

Power_Mod <- makeGlmer(y ~ Condition + 
                             Session + 
                             Facilitator + 
                             Condition:Session +
                             Facilitator:Session +
                             (1|Group),
                           family = binomial(link = "logit" ),
                           fixef = b,
                           VarCorr = V1,
                           data = Data)
Error in mkNewReTrms(object, newdata, compReForm, na.action = na.action,  : 
  NAs are not allowed in prediction data for grouping variables unless allow.new.levels is TRUE
pitakakariki commented 4 years ago

Why not just remove the values where Group is NA?

glmer is going to ignore them anyway.