首页 > 解决方案 > 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.

标签: r

解决方案


推荐阅读