Zero power in complex models?

[chengafni]

Hi, I tried using the powerSim function on a linear mixed effects model with interaction terms, but I always get 0% power, even if the tested effects are highly significant, and even if I examine only the power of main effects. When I re-fitted the model to include only the single most significant main effect and no interactions, I got 100% power. Does powerSim have problems with complex models? If not, how can I use it to test the power of various effects? Thanks, Chen

[pitakakariki]

0% power could mean that there are errors during the simulations, so you should check lastResult()$errors.

Most likely problem with interactions is you haven't specified a test. Examples here: https://cran.r-project.org/web/packages/simr/vignettes/examples.html

[chengafni]

Thank you. Indeed, there were many errors of the type "Test did not return a p-value", but I'm not sure why. I tried running doTest on the model with anova test and f test. Neither of them worked (f test took a really long time, too). My model is of the form y~x1*x2*x3+(1|g1)+(1|g2)+(1|g3)+(1|g1:g2)+(1|g1:g3)+(1|g4). y is continuous, and all factors are categorical except for g4 which is continuous. I would like to test the power of the 3-way interaction of the fixed effects, and also of 2-way interactions, if possible. How should I phrase the command?

[pitakakariki]

I would probably start with test=fixed("???", "z") where ??? is the name generated for the 3-way interaction.

The random effect structure of your model doesn't look right to me. Would you mind posting the output from summary?

[chengafni]

It's funny, I tried again doTest with anova, and it worked, and the powerSim gave me 100% power for the 3-way interaction, though I'm a bit skeptical about it... Also, I got a non-significant result for a two-way interaction that came out significant when running the anova function on the model. What's the difference between the two methods? When I tried running doTest with a z test I got an error: "Error in a[xname, "t value"] : subscript out of bounds". So, I'm a bit confused about the process. I would appreciate any pointers... See below the output of summary. It's quite long though... Thank you!###

Linear mixed model fit by REML t-tests use Satterthwaite approximations to
  degrees of freedom [lmerMod]
Formula: adjRT ~ Transposition * Rotation * Length + (1 | Subject) + (1 |  
    Order) + (1 | Word) + (1 | Subject:Word) + (1 | Subject:Order) +  
    (1 | MBF_con_freq)
   Data: nonwordsRTData

REML criterion at convergence: 900.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.1734 -0.6030 -0.0413  0.5776  3.6731 

Random effects:
 Groups        Name        Variance  Std.Dev.
 Subject:Word  (Intercept) 0.0033186 0.05761 
 Subject:Order (Intercept) 0.0081528 0.09029 
 MBF_con_freq  (Intercept) 0.0038482 0.06203 
 Subject       (Intercept) 0.0551161 0.23477 
 Word          (Intercept) 0.0001632 0.01278 
 Order         (Intercept) 0.0020329 0.04509 
 Residual                  0.0511840 0.22624 
Number of obs: 8200, groups:  
Subject:Word, 2126; Subject:Order, 1060; MBF_con_freq, 120; Subject, 76; Word, 40; Order, 14

Fixed effects:
                                           Estimate Std. Error         df
(Intercept)                              -1.537e+00  3.798e-02  1.960e+02
TranspositionInitial                      2.303e-02  3.307e-02  1.630e+02
TranspositionMiddle                       1.255e-01  3.445e-02  1.890e+02
Rotation-35                               3.887e-02  2.289e-02  6.292e+03
Rotation0                                -5.891e-04  2.257e-02  6.163e+03
Rotation35                                5.723e-02  2.296e-02  6.286e+03
Rotation70                               -1.154e-03  2.235e-02  6.105e+03
Length5                                   6.130e-02  3.208e-02  2.510e+02
Length6                                  -1.611e-03  3.841e-02  2.480e+02
TranspositionInitial:Rotation-35         -1.166e-01  3.200e-02  6.292e+03
TranspositionMiddle:Rotation-35          -1.005e-01  3.439e-02  6.431e+03
TranspositionInitial:Rotation0           -8.651e-02  3.161e-02  6.107e+03
TranspositionMiddle:Rotation0            -8.070e-02  3.386e-02  6.271e+03
TranspositionInitial:Rotation35          -1.492e-01  3.174e-02  6.231e+03
TranspositionMiddle:Rotation35           -1.144e-01  3.433e-02  6.374e+03
TranspositionInitial:Rotation70          -3.815e-02  3.201e-02  6.207e+03
TranspositionMiddle:Rotation70            5.600e-04  3.434e-02  6.319e+03
TranspositionInitial:Length5             -4.809e-02  4.455e-02  1.630e+02
TranspositionMiddle:Length5              -3.673e-02  4.716e-02  2.040e+02
TranspositionInitial:Length6             -3.761e-03  5.346e-02  1.640e+02
TranspositionMiddle:Length6              -4.734e-02  6.121e-02  2.620e+02
Rotation-35:Length5                      -8.122e-02  3.116e-02  6.290e+03
Rotation0:Length5                        -8.159e-02  3.078e-02  6.251e+03
Rotation35:Length5                       -1.210e-01  3.093e-02  6.254e+03
Rotation70:Length5                        3.231e-02  3.090e-02  6.232e+03
Rotation-35:Length6                      -1.782e-02  3.630e-02  6.149e+03
Rotation0:Length6                        -5.546e-02  3.663e-02  6.270e+03
Rotation35:Length6                       -1.040e-01  3.748e-02  6.400e+03
Rotation70:Length6                        3.286e-02  3.656e-02  6.082e+03
TranspositionInitial:Rotation-35:Length5  1.305e-01  4.355e-02  6.317e+03
TranspositionMiddle:Rotation-35:Length5   1.444e-01  4.801e-02  6.429e+03
TranspositionInitial:Rotation0:Length5    9.556e-02  4.272e-02  6.065e+03
TranspositionMiddle:Rotation0:Length5     1.384e-01  4.758e-02  6.413e+03
TranspositionInitial:Rotation35:Length5   1.661e-01  4.322e-02  6.257e+03
TranspositionMiddle:Rotation35:Length5    1.221e-01  4.782e-02  6.414e+03
TranspositionInitial:Rotation70:Length5   1.825e-02  4.344e-02  6.190e+03
TranspositionMiddle:Rotation70:Length5   -3.400e-02  4.829e-02  6.468e+03
TranspositionInitial:Rotation-35:Length6  1.367e-03  5.036e-02  6.106e+03
TranspositionMiddle:Rotation-35:Length6   1.449e-01  6.463e-02  6.744e+03
TranspositionInitial:Rotation0:Length6    9.397e-02  5.138e-02  6.284e+03
TranspositionMiddle:Rotation0:Length6     1.265e-01  6.641e-02  6.924e+03
TranspositionInitial:Rotation35:Length6   1.650e-01  5.129e-02  6.271e+03
TranspositionMiddle:Rotation35:Length6    2.002e-01  6.724e-02  7.098e+03
TranspositionInitial:Rotation70:Length6  -2.042e-02  5.152e-02  6.186e+03
TranspositionMiddle:Rotation70:Length6    1.512e-02  6.396e-02  6.669e+03
                                         t value Pr(>|t|)    
(Intercept)                              -40.477  < 2e-16 ***
TranspositionInitial                       0.696 0.487194    
TranspositionMiddle                        3.645 0.000346 ***
Rotation-35                                1.698 0.089501 .  
Rotation0                                 -0.026 0.979177    
Rotation35                                 2.493 0.012691 *  
Rotation70                                -0.052 0.958837    
Length5                                    1.911 0.057183 .  
Length6                                   -0.042 0.966567    
TranspositionInitial:Rotation-35          -3.644 0.000271 ***
TranspositionMiddle:Rotation-35           -2.922 0.003491 ** 
TranspositionInitial:Rotation0            -2.737 0.006222 ** 
TranspositionMiddle:Rotation0             -2.383 0.017196 *  
TranspositionInitial:Rotation35           -4.700 2.65e-06 ***
TranspositionMiddle:Rotation35            -3.333 0.000864 ***
TranspositionInitial:Rotation70           -1.192 0.233375    
TranspositionMiddle:Rotation70             0.016 0.986988    
TranspositionInitial:Length5              -1.079 0.282046    
TranspositionMiddle:Length5               -0.779 0.437039    
TranspositionInitial:Length6              -0.070 0.944003    
TranspositionMiddle:Length6               -0.773 0.440061    
Rotation-35:Length5                       -2.607 0.009160 ** 
Rotation0:Length5                         -2.650 0.008064 ** 
Rotation35:Length5                        -3.913 9.21e-05 ***
Rotation70:Length5                         1.046 0.295781    
Rotation-35:Length6                       -0.491 0.623611    
Rotation0:Length6                         -1.514 0.130078    
Rotation35:Length6                        -2.776 0.005527 ** 
Rotation70:Length6                         0.899 0.368809    
TranspositionInitial:Rotation-35:Length5   2.996 0.002748 ** 
TranspositionMiddle:Rotation-35:Length5    3.009 0.002634 ** 
TranspositionInitial:Rotation0:Length5     2.237 0.025314 *  
TranspositionMiddle:Rotation0:Length5      2.910 0.003630 ** 
TranspositionInitial:Rotation35:Length5    3.844 0.000122 ***
TranspositionMiddle:Rotation35:Length5     2.553 0.010691 *  
TranspositionInitial:Rotation70:Length5    0.420 0.674498    
TranspositionMiddle:Rotation70:Length5    -0.704 0.481476    
TranspositionInitial:Rotation-35:Length6   0.027 0.978348    
TranspositionMiddle:Rotation-35:Length6    2.242 0.024987 *  
TranspositionInitial:Rotation0:Length6     1.829 0.067449 .  
TranspositionMiddle:Rotation0:Length6      1.904 0.056905 .  
TranspositionInitial:Rotation35:Length6    3.218 0.001299 ** 
TranspositionMiddle:Rotation35:Length6     2.978 0.002914 ** 
TranspositionInitial:Rotation70:Length6   -0.396 0.691931    
TranspositionMiddle:Rotation70:Length6     0.236 0.813166    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 45 > 12.
Use print(x, correlation=TRUE)  or
         vcov(x)         if you need it

[pitakakariki]

High power is not surprising with this amount of data, unless the effect is very small.

The tests used by fixed look at one 3-way interation at a time. It probably makes more sense to test them all simultaneously, for example:

test=fcompare(~ Transposition*Rotation*Length - Transposition:Rotation:Length)

[chengafni]

Thank you! Some follow-up questions, if I may:

1) I use the anova function from the lmerTest package to test effects in my model. I would expect the doTest function with the "anova" option to give similar results, but it doesn't. Any idea why?

2) I performed powerSim on my model. I tested the highest interaction term ('Transposition:Rotation:Length') and one two-way interaction term ('Transposition:Rotation'). The 3-way interaction got 99.5% power and the 2-way interaction received only 74% power. I would expect it to be the other way around. Any thoughts?

3) During powerSim the following message appeared: "anova from lme4 is returned some computational error has occurred in lmerTest", and then the simulation process seemed to start from the beginning. What would be the implications of the switch from lmerTest to lme4?

[pitakakariki]

1. At a guess this may be to do with the type of test? I think lmerTest defaults to Type III but simr defaults to Type II to match Anova from the car package.

2. I would also expect this to be the other way around, but it's not implausible depending on the relative effect sizes.

3. This looks like lmerTest using message when warning might be more appropriate. I suspect the progress bar is being interrupted by the message, rather than the computations restarting?

[chengafni]

Thank you for the clarifications.

1. Is it possible to make simr use Type III test? I didn't quite get it from the manual.

2. the vignettes page states that "Main effects should not be tested when they appear in an interaction term." Does this apply to powerSim as well? My model has a 3-way interaction term and I wanted to test the power of 2-way terms. It seems to work fine, I just wanted to make sure the results are reliable for lower order terms.

3. on the same note, I tried powerCurve on the 2-way term to get the required number of subjects, as follows: powerCurve(fit, test=fixed("Transposition:Rotation"), along = "Subject" , nsim = 50, seed=1). The simulation returned 0% power on all numbers of levels of Subject. powerSim returned 74% for the 2-way term. Any idea what could be the problem? I basically want to determine how many subjects I need for 80% power. Thanks again!

[pitakakariki]

1. I think you can use lmerTestType=3 as an option in powerSim and powerCurve.

2. This is an interesting question, normally I'd expect this to be fine but you seem to have unusually large 3-way interactions. You might want to try running powerSim on a model without 3-way interactions as well.

3. 0% usually means a lot of errors. Did anything show up in lastResult()$errors?

[chengafni]

1. yes, that works.

2. it's a weird situation. I am actually interested in a 2-way interaction which becomes statistically significant only when there is also a 3-way term in the model (I have no idea why). When I remove the 3-way term I can run simPower on the 2-way term with Type III test with no errors. When I use Type III test on the 2-way term in a model that also includes a 3-way term, I get the following message several times:

   anova from lme4 is returned 
   some computational error has occurred in lmerTest"

   The failure results from: Test did not return a p-value, which is not very informative. Nevertheless, the summary of results indicates that Type III test was performed after all...

3. powerCurve reports: Models have either equal fixed mean stucture or are not nested (note the typo in stucture). I guess this is why I get 0%. How can I fix this issue?

Thanks!

[pitakakariki]

That error message comes from KRmodcomp in pbkrtest, so it's probably something to do with your fcompare specifications.

What are the formulas for the model you're testing and for the argument to fcompare?

[chengafni]

I actually don't use fcompare (and don't know what it does). These are the commands I use:

Model: fit<-lmer("adjRT ~ Transposition*Rotation*Length+(1|Subject)+(1|Order)+(1|Word)+(1|Subject:Word)+(1|Subject:Order)+(1|MBF_con_freq)",data=nonwordsRTData)

PowerCurve: pc<- powerCurve(fit, test=fixed("Transposition:Rotation"), along = "Subject" , nsim = 50, seed=1,lmerTestType=3)

[pitakakariki]

I think you want test=fixed("Transposition:Rotation", method="anova") here.

[chengafni]

Yes!! I still get many of those Test did not return a p-value errors, but the simulation seems to work. Many many thanks!!!

[chengafni]

Hi again, I have another 0% power issue. This time it seems to be associated with "random effects specified in re.form that were not present in original model". I found a couple of discussions on this error, but I couldn't find a solution that works for me. This is what my data looks like:

a
# A tibble: 3,039 x 3
   Subject Condition  adjRT
   <fct>   <fct>      <dbl>
 1 MIE1_10 noroot    -1.11 
 2 MIE1_10 root      -0.935
 3 MIE1_10 root      -0.845
 4 MIE1_10 root      -0.783
 5 MIE1_10 noroot    -0.927
 6 MIE1_10 root      -0.767
 7 MIE1_10 noroot    -1.07 
 8 MIE1_10 noroot    -0.955
 9 MIE1_10 noroot    -0.914
10 MIE1_10 noroot    -0.863
# ... with 3,029 more rows

I fitted a model with the following formula: fit<-lmer("adjRT~Condition+(Condition|Subject)",data=a)

and then tried to calculate the power: power <- powerSim(fit, test=fixed("Condition", method="anova"),lmerTestType=3), which, as I mentioned, yielded 0%. Each iteration of the simulation returned the aforementioned error. Also, the results of powerSim included nrow = NA. I'm guessing it's part of the problem.

Any idea how to solve this? Chen

[pitakakariki]

I can't replicate the error. What versions of simr and lmerTest?

[chengafni]

lmerTest: 2.0.36 simr: 1.0.3 R: 3.4.3 (64-bit)

[pitakakariki]

Try upgrading both, see if that fixes it?

simr: 1.0.4 lmerTest: 3.0-1

[chengafni]

I updated R to version 3.5.0 and also all the packages. I still get 0% power with the same data but this time with a different error: object of type 'symbol' is not subsettable. Also, the nrow parameter of powerSim is calculated correctly now. Any thoughts?

[pitakakariki]

I don't get any errors when I run something similar. Are you running any other commands first?

[chengafni]

I load a larger data set from an Excel file. Then I perform some manipulations on the data and use the subset function to remove some rows and columns I don't need.

[pitakakariki]

Probably the only way I can solve this is if you send me the data and script. Email at https://cran.r-project.org/web/packages/simr/

[pitakakariki]

Looks like some sort of compatibility problem between simr and lmerTest.

If you comment out the 'library(lmerTest)` line it should work.

[chengafni]

Wow, I never would have thought trying that! Thank you so much!

[chengafni]

Sorry to bother again, but now I have a problem with powerSim of a glmer object: model not of class 'lmerMod': cannot coerce to class 'lmerModLmerTest

Any idea what that might be about?

[pitakakariki]

I don't think the lmerTest package supports GLMMs, so I think you'll need to try a different test.

[chengafni]

Oh, didn't know that. Thanks!

[jrode11]

Hello,

I'm also having trouble detecting power for an interaction. I have the following model:

mod3.s1<-glmer(true.other ~ c.ideology*pol.state.contr*c.truth2+(pol.state.contr*c.truth2|id)+(c.ideology|number), data = facts.long, family = binomial, control=glmerControl(optimizer="bobyqa"), nAGQ = 1

It’s a three-way interaction (continuous/between-subjects, categorical/within-subjects, continuous/within-subjects) with by-subjects random slopes and intercepts for the two within-subjects variables, and by-items random slopes and intercepts for the between-subjects variable. I have around 1300 participants with 24 observations each.

I’m interested in the power of one of the interactions: c.ideology:pol.state.contr

I set the effect size of the interaction to be small (.05) and tried to calculate power:

model.small<-mod3.s1 fixef(model.small)['c.ideology:pol.state.contr'] <- .05 model.small.test<-extend(model.small, along="id", n=2000)

try<-powerCurve(model.small.test, along="id", test=fixed("c.ideology:pol.state.contr", "z"), nsim=20)

But I get the following errors when I check lastResult()$errors:

1 Fitting 1 object 'true.other' not found 2 Fitting 2 object 'true.other' not found 3 Fitting 3 object 'true.other' not found 4 Fitting 4 object 'true.other' not found 5 Fitting 5 object 'true.other' not found

…all the way until the last simulation number. 'true.other' is my dependent variable

Is there a reason that it can’t find my dependent variable? The original model I ran (mod3.s1) works fine.

[pitakakariki]

What sort of binomial variable is true.other? i.e. is it TRUE / FALSE, 1 / 0, or a successes / failures matrix?

[jrode11]

It's a 0/1 numeric variable in the data frame

[pitakakariki]

No idea what might be causing this. Make sure you have the latest version from CRAN. If that doesn't help, if you email me the data I can have a closer look.

[jrode11]

Just figured out it has to do with missing data, sorry I didn't realize this earlier! I made sure none of the predictors nor the DV had missing data before running the glmer model and it worked. Thanks!

[schwerdt-aut]

Hi, I want to chime in here, because I also experience problems using simr. The lmer-model looks as follows: lmer(dv~a+b+c+d+e+f*g*h+(1|id). g is a factorized variable (0,1,2,3,), h a dichotomous and f a continuous variable. I am particularly interested in the power of the three-way interaction with a particular factor level of g (say g2). Running simr I get the following output: powerSim(data, test=fixed("f:g2:h", "z"), nsim=1000, seed=1)

> 
_Power for predictor 'f:g2:h', (95% confidence interval):==============================================|
       0.00% ( 0.00,  0.37)
Test: z-test
      Effect size for f:g2:h is 0.070
Based on 1000 simulations, (0 warnings, 1000 errors)
alpha = 0.05, nrow = NA
Time elapsed: 0 h 0 m 0 s

nb: result might be an observed power calculation
Warning message:
In observedPowerWarning(sim) :
This appears to be an "observed power" calculation_
> 

Entering lastResult()$errors I get the following error: argument is of length zero

The same error happens when using a much simpler model (say, dv~a+(1|id)). The dataset does not contain missing values. I am using R 3.6.1, lme4 _1.1-21, and simr_1.0.5. Any help is much appreciated.

Best regards, Andreas

[pitakakariki]

This might be a silly question but are you giving powerSim your data instead of your model as its first argument?

[schwerdt-aut]

Sorry for the confusion. The first argument refers to the lmer model.

[pitakakariki]

I'd need more info to guess what's happening. What stage is it telling you has an error? What test are you using for the simpler model and what error are you getting?

[schwerdt-aut]

Thank you Peter for your help. Meanwhile, I figured out what went wrong. First, I was analyzing data imported from SPSS (read.spss). While lmer worked perfectly, simr seems to need data assigned to a data frame (data.frame). Thereafter, simr worked well. I still got some error messages at first, but these were due to one predictor variable still having missing data. Finally, the three-way interaction was nicely simulated using "z" as the test.

Best regards, Andreas

[mrahimit]

Hi All, I have a lme model and I wanna run power analysis. Linear mixed-effects model fit by REML Data: NULL Log-restricted-likelihood: -212.0827 Fixed: br ~ type * imp (Intercept) type2 imp2 imp3 type2:imp2 type2:imp3 0.3346635 1.5904707 0.1059542 0.1010908 -1.3444451 -1.3160416

Random effects: Formula: ~1 | subj (Intercept) Residual StdDev: 0.09222212 0.1949634

Variance function: Structure: Different standard deviations per stratum Formula: ~1 | type imp Parameter estimates: 02 22 03 23 01 2*1 1.000000 1.695103 1.022635 1.905845 1.006305 13.732514 Number of Observations: 2906 Number of Groups: 20 When I use powerSim(modl1), I keep getting the error: could not find function "powerSim" I updated everything. Do you know what is going on?

[pitakakariki]

simr is based on lme4, and does not support nlme models.

[amanol]

Hi,

I am getting 0.00% power and errors because of the levels/contrasts in the categorical factor I have in my model. How can I resolve this?

Thanks, Athina

[io-inspace]

I'm posting this for anyone else troubleshooting power analysis in simr for interaction effects: if you're getting power for test="t" but not "chisq", "anova", or "f", or vice versa, check how you're specifying the interaction term in the test function. For the "t" test the interaction effect should appear as it looks in summary(yourmodel). For the other tests, it should appear as it looks in Anova(yourmodel). This tripped me up quite a bit, hope it helps someone else!

[pitakakariki]

Thanks, even I tend to forget which is which. Maybe 2022 is the year I finally update the documentation!

[obada-alzoubi]

For those who have the same issue. Use lme4 instead of lmerTest package for model fitting. All factors have to be numeric factors i.e. as.numeric(as.factor(df$variable)). This how it was solved for me. Otherwise, I was getting zero power.