CODEWITHSUNDEEP
r
llewmills
llewmills

Reputation: 3568

SIMR Package - powerCurve function throws error when equivalent powerSim function does not

I am trying to do a sample size calculation for a grant application. We want to run a simple mixed effects model for repeated measures regression examining change in a continuous variable over 7 days in two groups (treatment vs placebo). In lme4 the model, run on pilot data of 7 participants per group, return the following output

modFull <- lmer(formula = score ~ day*group + (1|id),
                data = longFull)

modFull

# relevant output
Fixed Effects:
     (Intercept)               day      groupplacebo  day:groupplacebo  
         16.0296           -0.7541            1.6035           -0.5561

I want to power the study on the interaction coefficient day:groupplacebo, based on a difference in per-day rate of change between groups of 0.63. I did this via the fixef(foo) <- function.

fixef(modFull)["day:groupplacebo"] <- 0.63

Now when we look at the fixed effects

fixef(modFull)

(Intercept)              day     groupplacebo day:groupplacebo 
  16.0295999       -0.7540535        1.6035141        0.6300000

the original interaction effect has been replaced by the desired interaction effect

Following advice in the simr vignettes and here, in order to get a power estimate the effect of the interaction I used the fcompare function like so (compares full model to model with only main effects, thereby testing the interaction term)

powerSim(modFull, fcompare(~ day + group))

And got the following output

Power for model comparison, (95% confidence interval):==============================================|
      27.80% (25.04, 30.69)

Test: Likelihood ratio
      Comparison to ~day + group + [re]

Based on 1000 simulations, (0 warnings, 0 errors)
alpha = 0.05, nrow = 98

Obviously I need more power than the 27.8% given to me by 14 participants so I use the extend() function like so, this time for the same model with the same effect size but a hypothetical sample size of 40 per group (80 total)

modFullExt <- extend(modFull, along="id", n=80) 
powerSim(modFullExt, fcompare(~ day + group))

# output
    Power for model comparison, (95% confidence interval):==============================================|
      87.20% (84.97, 89.21)

Test: Likelihood ratio
      Comparison to ~day + group + [re]

Based on 1000 simulations, (0 warnings, 0 errors)
alpha = 0.05, nrow = 560

So, all good so far. My problem is that I would like to get a power curve using the powerCurve() function. But when I tried

powerCurve(modFullExt, test = fcompare(~ day + group))

I got the following error

Error in getDefaultXname(fit) : 
  Couldn't automatically determine a default fixed effect for this model.

Now why does the syntax for powerSim() not translate to powerCurve()?

I should note that the following powerSim() syntax

powerSim(modFull, fixed("day:groupplacebo", method = "z"))

returned the same result as powerSim(modFullExt, fcompare(~ day + group)) but when I tried

powerCurve(modFullExt, test = fixed("day:groupplacebo", "z"))

It returned the same error.

Any help greatly appreciated.

Upvotes: 3

Views: 891

Answers (1)

opmorgan
opmorgan

Reputation: 26

This error message may be caused by powerCurve trying to use the default argument for along, which is getDefaultXname(fit). Your code should run as expected if you include an explicit argument for along (the variable whose number of levels will vary in the power analysis, such as number of subjects):

powerCurve(modFullExt, test = fcompare(~ day + group), along = "id")

Upvotes: 1

Related Questions