No applicable method for 'test' applied to an object of class "c('glmerMod', 'merMod')"

[JCruk]

In fitting some gaussian log models, I am getting an error for each simulation in powerSim:
no applicable method for 'test' applied to an object of class "c('glmerMod', 'merMod')"

The model is relatively straightforward:

Mod <- glmer(Y ~ A + B + C + A:B + B:C + A:C + (1|G),
data = data,
family = gaussian(link = "log"),
mustart = pmax(data$Y),
glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000))

I then call powerSIm with

powerSim(Mod, test = fixed("A"), nsim = 100)

I get an estimated power of 0.00 with 100 errors (as above)

Has there been a change in how to call powerSim for glmer models?

[pitakakariki]

That looks like maybe simr's heuristics are choosing the wrong test.

Does it work if you manually specify the test? i.e. fixed("A", "???").

Is A numeric or categorical here?

[JCruk]

A is categorical here. I tried powerSim(Mod, test = fixed("A", "z"), nsim = 25) with the same no applicable method... result

[pitakakariki]

How about if you use simr::fixed in place of fixed?

[JCruk]

That's what I was thinking too - and it was a function masking issue. I usually load simr last, but didn't here. Sorry for the bother!

I'm now getting "Test did not return a p-value" with all interaction terms. I suspect that's due to model specification rather than anything in simr. Is that a correct assumption?

[pitakakariki]

That's a catch-all in the test wrapper, so could be a number of things, but most likely an NA p-value.

You should be able to bypass the wrapper with simr::lrtest(Mod, "A") to see what's happening.

[JCruk]

Thanks!
I get an error trying to bypass the wrapper:
Error: 'lrtest' is not an exported object from 'namespace:simr

[pitakakariki]

Sorry, triple-double dots: simr:::lrtest(Mod, "A").

[JCruk]

So there is an NA

simr:::lrtest(Mod "A")
[1] NA

Is there anything I can do about that?

[pitakakariki]

Try drop1(Mod, ~A, test="Chisq"), which is the test simr's calling from lme4. If that doesn't work, it's probably a problem with your model, but you might get a more helpful error messsge.

[JCruk]

That does work. Is there a way to combine the single term deletion with simr so I can plot out a powerCurve?

[pitakakariki]

fixed("A", "lr") or even just fixed("A") should be using that drop1 result. If drop1 is working though, maybe the wrapper's having trouble extracting the p-value.

What does the drop1 output look like?

[JCruk]

In the model with my actual data (where I get the no p-value errors):
powerSim(Valence_glmer, test = fixed("Category", "lr"), nsim = 10),

drop1(Valence_glmer, ~Category, test="Chisq") returns:

Model:
valence ~ condition + filter_type + Category + condition:filter_type + 
    condition:Category + Category:filter_type + (1 | initials)
         Df   AIC        LRT Pr(Chi)
<none>      27427                   
Category  0 27427 1.4767e-08 

[pitakakariki]

Aha, drop1 isn't returning a p-value. I suspect that the problem is with the zero degrees of freedom, although I have no idea what could cause that.

[JCruk]

It's likely due to something in specification. When I run powerSim on a "typical" glmer specification I get cannot find valid starting values: please specify some for every simulation - even though the model runs just fine. So I added in mustart = pmax(data$y) to the model call, but then I get the no p-value errors.

Have you encountered the cannot find valid starting values: please specify some error previously?

[JCruk]

I've tried a few more things and while my model runs fine in glmer with the specified values for mustart, it seems something is happening when I pass the model over to powerCurve.

When I try plot(powerCurve(Mod, test = fixed("A), along = "G", nsim = 10)),

I get a flat plot at zero and lastResult()$errors shows:
variable lengths differ (found for '(mustart)') for each simulation

[JCruk]

In case it's helpful to be able to reproduce my issue:

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

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

## Build a simulated data set
Sub_Group <- data.frame(Subject = LETTERS[1:20], 
                        Group = c(rep("Y", 10), rep("Z", 10)))

Cat_Item <- data.frame(expand.grid(Category = c(rep("U", 25), rep("V", 25)), 
                                        Item = factor(rep(c(1:25),2))))

Cond <- data.frame(Condition = letters[1:5])

Cond_Cat_Item <- merge(Cond, Cat_Item, by = NULL)

Data <- merge(Sub_Group, Cond_Cat_Item, by = NULL) %>%
        mutate(Y = round(rtnorm(nrow(.), 4.5, 2, 1, 9), digits = 0))

## Run initial glmer without any specified starting values

summary(Mod1 <- glmer(Y ~ Category + Condition + Group +
                             (1 + Group|Subject),
                     data = Data,
        family = gaussian(link = "log")))

plot(powerCurve(extend(Mod1, along = "Subject", n = 50), along = "Subject", nsim = 10))

# Errors suggesting to provide starting values
lastResult()$errors
# "cannot find valid starting values: please specify some"

## Second glmer with starting values
summary(Mod2 <- glmer(Y ~ Category + Condition + Group +
                             (1 + Group|Subject),
                     data = Data,
                     mustart = Data$Y,
                     family = gaussian(link = "log")))

plot(powerCurve(extend(Mod2, along = "Subject", n = 50), along = "Subject", nsim = 10))

## Errors regarding variable lengths
lastResult()$errors
# "variable lengths differ (found for '(mustart)')"

[chri4354]

I am encountering the same problem. I am fairly new to simr, so pardon me.

I have the following model:
m <- glmer('error ~ type+ (1 | p_id:group) + (1|Q)', data = df, family=gaussian(link="log"))

where type is a categorical variable with 4 levels, p_id is the participant id, group is the group id.

When I run
powerSim(m, test = simr::fixed("type"), nsim=20 )

I get the error "cannot find valid starting values: please specify some" in all 20 simulations, even if the glmer function runs smoothly. Running drop1(m, ~type, test="Chisq") also does not give any problem and returns a p-value.

Any hint to how to fix it?

[pitakakariki]

It looks like the problem is happening at the simulation stage (it works for simulate.merMod but simr uses simulate.formula).

Workaround for now would be a custom simulation function, e.g.:

p <- list(beta=fixef(m), theta=getME(m, "theta"), sigma=sigma(m))
s <- function() simulate(m, newparams=p)[[1]]

powerSim(m, sim=s, etc)