Issue with estimating power for an interaction term.

[cv0104]

Hello,

I'm attempting to estimate power for an interaction term but seem to be running into some issues.

My model: diss_model2<-glmer(Log_rt~1+v*a+(1|Subject), family = gaussian(link="log"), data=data.situation) summary(diss_model2)

V represents valence, and is coded 0/1. A represents Agreeableness, and is a continuous predictor. Log_rt is log-transformed reaction time. The outcome was a result of a repeated measures design, so Subject serves as the random effect parameter. I've double checked to be sure that all variables are coded correctly (i.e., V is read as a factor, A as numeric, etc.).

My code for the power analysis: powerSim(diss_model2, fixed("v1:a", "z"), seed=123, nsim=20)

The error call states that for each simulation, the function "cannot find valid starting values: please specify some". lastResult()$errors

Is it possible to specify start values? Additionally, is there something wrong in my code that is leading to the error messages?

Thanks in advance for your time!

[pitakakariki]

My first guess is that you really want the model:

diss_model2 <- lmer(Log_rt ~ v * a + (1|Subject), data=data.situation) 

If that's not the case, could you please let me know what summary(diss_model2) looks like?

[cv0104]

Thanks for the quick reply!

I've run both lmer() and glmer() models and the results are below:

GLMER

> diss_model2<-glmer(Log_rt~1+v*a+(1|Subject), family = gaussian(link="log"), data=data.situation)
> summary(diss_model2)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: gaussian  ( log )
Formula: Log_rt ~ 1 + v * a + (1 | Subject)
   Data: data.situation

     AIC      BIC   logLik deviance df.resid 
120168.5 120209.9 -60078.3 120156.5     7314 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5012 -0.6504 -0.1285  0.6235  3.9153 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 143916   379.4   
 Residual             581477   762.5   
Number of obs: 7320, groups:  Subject, 122

Fixed effects:
            Estimate Std. Error t value Pr(>|z|)    
(Intercept)  7.20635    0.04540  158.75  < 2e-16 ***
v1           0.08552    0.01059    8.08 6.63e-16 ***
a            0.01355    0.04541    0.30    0.765    
v1:a        -0.04647    0.01038   -4.48 7.63e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr) v1     a     
v1   -0.125              
a    -0.002  0.011       
v1:a  0.011 -0.046 -0.126

LMER

> diss_model2_lm<-lmer(Log_rt~1+v*a+(1|Subject), data=data.situation)
> summary(diss_model2_lm)
Linear mixed model fit by REML ['lmerMod']
Formula: Log_rt ~ 1 + v * a + (1 | Subject)
   Data: data.situation

REML criterion at convergence: 118467

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.4407 -0.6542 -0.1261  0.6298  3.7591 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 426100   652.8   
 Residual             590261   768.3   
Number of obs: 7320, groups:  Subject, 122

Fixed effects:
            Estimate Std. Error t value
(Intercept)  1489.38      60.45  24.639
v1            155.50      17.96   8.658
a              55.62      60.45   0.920
v1:a          -81.76      17.96  -4.552

Correlation of Fixed Effects:
     (Intr) v1     a     
v1   -0.149              
a     0.000  0.000       
v1:a  0.000  0.000 -0.149

Here are the results of the power analyses for the respective models:

GLMER MODEL

Power for predictor 'v1:a', (95% confidence interval):
       0.00% ( 0.00,  3.62)

Test: z-test
      Effect size for v1:a is -0.046

Based on 100 simulations, (0 warnings, 100 errors)
alpha = 0.05, nrow = 7320

Time elapsed: 0 h 0 m 10 s

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

LMER MODEL

powerSim(diss_model2_lm, fixed("v1:a", "z"), seed=123, nsim=100)

Power for predictor 'v1:a', (95% confidence interval):
      100.0% (96.38, 100.0)

Test: z-test
      Effect size for v1:a is -82.

Based on 100 simulations, (0 warnings, 0 errors)
alpha = 0.05, nrow = 7320

Time elapsed: 0 h 0 m 12 s

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

There appears to be something off with the power calculations for the lmer model as well, given that the results indicate 100% power to detect the interaction term.

[pitakakariki]

You have a large number of observations, and a very large estimated effect size for the interaction. So the power estimate is probably correct for this scenario. You could try a smaller effect size, e.g.

fixef(diss_model2)["v1:a"] <- -10 

I'm not quite sure what you're trying to do with the glmer model. It looks to me like you're logging rt twice? It's a pretty unusual model.

What's happening in the power simulations is that some of the simulated responses end up being negative, which makes it very hard for lme4 to fit. If you have a very good reason for preferring that model I can investigate further, but I suspect that the lmer model is more appropriate?

[cv0104]

Ahh, I see, thanks for pointing out the issues with my glmer model. I'll go ahead and tinker around with smaller effect sizes for the lmer model. Thanks so much for the help and developing the package!

[mengzhu-nz]

My linear mixed effects model had an interaction and random intercepts:

fit<-lmer(rttrans~Scondition*TargetType+(1|ID)+(1|probe),data=data1)

I am trying to estimate power for the above model. And I used the following code and got warnings and errors:

powerSim(fit, test=fixed("SconditionScleftO:TargetTypecontrol"),nsim = 30) [ I have deleted the NA datapoints.] Power for predictor 'SconditionScleftO:TargetTypecontrol', (95% confidence interval):============================| 0.00% ( 0.00, 11.57)

Test: unknown test Effect size for SconditionScleftO:TargetTypecontrol is 0.13

Based on 30 simulations, (0 warnings, 30 errors) alpha = 0.05, nrow = 5535

Time elapsed: 0 h 0 m 16 s

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

lastResult()$errors stage index message 1 Testing 1 scope is not a subset of term labels 2 Testing 2 scope is not a subset of term labels ...

powerSim(fit, test=fixed("TargetType"),nsim = 30) Power for predictor 'TargetType', (95% confidence interval):=====================================================| 0.00% ( 0.00, 11.57)

Test: Likelihood ratio

Based on 30 simulations, (30 warnings, 30 errors) alpha = 0.05, nrow = 5535

Time elapsed: 0 h 0 m 26 s

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

lastResult()$errors stage index message 1 Testing 1 Test did not return a p-value 2 Testing 2 Test did not return a p-value ...

But when I only had the simple effects in the model (without their interaction), the power analysis seems to return a resealable power but still with the same warning.

fit<-lmer(rttrans~Scondition+TargetType+(1|ID)+(1|probe),data=pilot) powerSim(fit, test=fixed("TargetType"),nsim = 30) Power for predictor 'TargetType', (95% confidence interval):=====================================================| 100.0% (88.43, 100.0)

Test: Likelihood ratio

Based on 30 simulations, (0 warnings, 0 errors) alpha = 0.05, nrow = 5535

Time elapsed: 0 h 0 m 24 s

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

I would really appreciate if you could help me with solving this problem.

[mengzhu-nz]

More info to my model: fit<-lmer(rttrans~Scondition*TargetType+(1|ID)+(1|probe),data=data1)

both Scondition and TargetType are categorical variables.

I further tried the following, but I got a different error. Please advise if you have any ideas. Thank you very much in advance. @pitakakariki @palday @cv0104

`> powerSim(fit, test=fixed("Scondition:TargetType", "anova")) Power for predictor 'Scondition:TargetType', (95% confidence interval):==========================================| 0.00% ( 0.00, 0.37)

Test: Type-2 F-test with Satterthwaite degrees of freedom (package lmerTest)

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

Time elapsed: 0 h 2 m 56 s

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

lastResult()$errors stage index message 1 Simulating 1 non-conformable arrays 2 Simulating 2 non-conformable arrays 3 Simulating 3 non-conformable arrays 4 Simulating 4 non-conformable arrays`

[pitakakariki]

Try specifying a particular test, e.g. fixed("SconditionScleftO:TargetTypecontrol", "z"),

[mengzhu-nz]

Thanks for your reply. I just tried the z test, but I got an error too:

powerSim(fit, test=fixed("SconditionScleftO:TargetTypecontrol","z")) Power for predictor 'SconditionScleftO:TargetTypecontrol', (95% confidence interval):====================================| 0.00% ( 0.00, 0.37)

Test: z-test Effect size for SconditionScleftO:TargetTypecontrol is 0.13

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

And I make the interaction a single term using "df$AB = interaction(df$A, df$B)", and I got the followings with warnings but not errors. Do they think this is a correct way to do? `> powerSim(fit, test=fixed("Q"))

Power for predictor 'Q', (95% confidence interval):======================================================================|

  100.0% (99.63, 100.0)

Test: Likelihood ratio

Based on 1000 simulations, (26 warnings, 0 errors)

alpha = 0.05, nrow = 5577

Time elapsed: 0 h 9 m 1 s

nb: result might be an observed power calculation

Warning message:

In observedPowerWarning(sim) :

This appears to be an "observed power" calculation

powerCurve(fit, test=fixed("Q"))

Calculating power at 10 sample sizes along Q

Power for predictor 'Q', (95% confidence interval),======================================================================|

by number of levels in Q:

  3: 51.70% (48.55, 54.84) - 1401 rows

  4: 67.50% (64.50, 70.40) - 1855 rows

  5: 75.70% (72.92, 78.33) - 2316 rows

  6: 73.90% (71.06, 76.60) - 2779 rows

  7: 99.80% (99.28, 99.98) - 3246 rows

  8: 99.80% (99.28, 99.98) - 3718 rows

  9: 100.0% (99.63, 100.0) - 4186 rows

 10: 100.0% (99.63, 100.0) - 4651 rows

 11: 100.0% (99.63, 100.0) - 5113 rows

 12: 100.0% (99.63, 100.0) - 5577 rows

Time elapsed: 0 h 49 m 58 s`

[pitakakariki]

The warning is because you're taking the effect size from your pilot data. Think about what sort of effect size you need to detect, and then set that with fixef<- before running the simulation.

More info here: http://core.ecu.edu/psyc/wuenschk/StatHelp/Power-Retrospective.htm

[mengzhu-nz]

Thank you very much. I really appreciate your help!

I set a medium effect size with fixef following the code in your SimR paper. and I got the following. I am now a litttle bit confused after reading the link that you sent to me and specifying the effect size. In your paper, page 494, you put it:"Retrospective ‘observed power’ calculations, where the target effect size comes from the data, give misleading results (Hoenig&Heisey 2001)."

My aim was to detect whether my model has a power of at least 80% given the data that I collected for my experiment. Why I cannot take the effect size from my data? As I have 3 and 4 levels for each factor, so I would need to specify the effect size for each of them, which becomes a bit confused...

powerSim(fit, test=fixed("SconditionScleftO:TargetTypecontrol","z")) Power for predictor 'SconditionScleftO:TargetTypecontrol', (95% confidence interval):====================================| 44.60% (41.49, 47.74)

Test: z-test Effect size for SconditionScleftO:TargetTypecontrol is 0.50

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

Time elapsed: 0 h 4 m 45 s

On Fri, Apr 5, 2019 at 5:31 PM Peter Green notifications@github.com wrote:

The warning is because you're taking the effect size from your pilot data. Think about what sort of effect size you need to detect, and then set that with fixef<- before running the simulation.

More info here: http://core.ecu.edu/psyc/wuenschk/StatHelp/Power-Retrospective.htm

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/pitakakariki/simr/issues/91#issuecomment-480144801, or mute the thread https://github.com/notifications/unsubscribe-auth/AuU0vgLeTPbqaUTnubZ4Beyrwf-ouJihks5vdtGugaJpZM4Qz7G9 .

[pitakakariki]

You can take the effect size from your data, it's just that it's usually a bad idea, which is why you get a warning.

The reason it's usually a bad idea is explained in Hoenig & Heisey.

[mengzhu-nz]

Thanks for the feedback. Am I right that when we would like to test the power of an interaction, we have to test each level of the interaction one by one?

I used the following code to get the power for the interaction Q, but I got the warning message of 'observed power calculation'.

pilot$Q = interaction(pilot$Scondition, pilot$TargetType) pilot <- droplevels(pilot) fit<-lmer(rttrans~Q+(1|ID)+(1|probe),data=pilot) Anova(fit)

powerSim(fit, test=fixed("Q"),nsim = 1000) Power for predictor 'Q', (95% confidence interval):=======================================================================================| 100.0% (99.63, 100.0)

Test: Likelihood ratio

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

Time elapsed: 0 h 8 m 44 s

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

And I tried to set an effect size with fixef, I got 16%, 75.40% and 100% for the different effect sizes that I had. You mentioned in your paper that the bigger the effect size is, the higher the power is. However, is there a standard effect size that we should expect? [if I need to detect a big effect size like 1, the power would be 100%].

In the following case, the actual estimate/effect size from the model was 0.18.

fixef(fit)["QScleftS.Contrastive"]<-0.18 powerSim(fit, test=fixed("QScleftS.Contrastive","z"),nsim = 1000) Power for predictor 'QScleftS.Contrastive', (95% confidence interval):====================================================================|

• 16.00% (13.78, 18.42)*

Test: z-test Effect size for QScleftS.Contrastive is 0.18

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

Time elapsed: 0 h 4 m 30 s

fixef(fit)["QScleftS.Contrastive"]<-0.5 powerSim(fit, test=fixed("QScleftS.Contrastive","z"),nsim = 1000) Power for predictor 'QScleftS.Contrastive', (95% confidence interval):====================================================================|

• 75.40% (72.61, 78.04)*

Test: z-test Effect size for QScleftS.Contrastive is 0.50

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

Time elapsed: 0 h 4 m 44 s

fixef(fit)["QScleftS.Contrastive"]<-1 powerSim(fit, test=fixed("QScleftS.Contrastive","z"),nsim = 1000) Power for predictor 'QScleftS.Contrastive', (95% confidence interval):====================================================================|

• 100.0% (99.63, 100.0)*

Test: z-test Effect size for QScleftS.Contrastive is 1.0

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

Time elapsed: 0 h 4 m 38 s

On Tue, Apr 9, 2019 at 2:56 PM Peter Green notifications@github.com wrote:

You can take the effect size from your data, it's just that it's usually a bad idea, which is why you get a warning.

The reason it's usually a bad idea is explained in Hoenig & Heisey.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/pitakakariki/simr/issues/91#issuecomment-481082203, or mute the thread https://github.com/notifications/unsubscribe-auth/AuU0vsf_bwBtcyLiP7SJdzOCwbY7_Tueks5vfAF3gaJpZM4Qz7G9 .

[pitakakariki]

I would normally test the entire interaction, but it really depends on what you're trying to do.

You could try a range of effect sizes until you get about 80% power - this would give you an idea of your Minimal Detectable Effect.

[maureenbowers]

Hello, I am trying to estimate sample size needed from pilot data (n=8). There are two conditions: Reward and Cue, each with two levels.

Below is my model and the summary:

random_intercept_ACC_mlm_lmer = lmer(Acc_FE ~ (1|id_num) + Reward.EC * Cue.EC , data = df_flank.Reward, na.action = "na.exclude") summary(random_intercept_ACC_mlm_lmer) REML criterion at convergence: -19.3

Scaled residuals: Min 1Q Median 3Q Max -2.03963 -0.65864 0.03047 0.63006 1.68580

Random effects: Groups Name Variance Std.Dev. id_num (Intercept) 0.02639 0.1625
Residual 0.01145 0.1070
Number of obs: 36, groups: id_num, 9

Fixed effects: Estimate Std. Error df t value Pr(>|t|)
(Intercept) -0.280093 0.057013 7.999995 -4.913 0.00117 * Reward.EC 0.009259 0.017837 24.000002 0.519 0.60844
Cue.EC 0.016204 0.017837 24.000002 0.908 0.37268
Reward.EC:Cue.EC -0.037037 0.017837 24.000002 -2.076 0.04873

I then ran the power analysis below and it seems like the power is around 53%

powerSim(random_intercept_ACC_mlm_lmer, test=fixed("Reward.EC:Cue.EC", method="z")) Power for predictor 'Reward.EC:Cue.EC', (95% confidence interval): 53.30% (50.15, 56.43)

Test: z-test Effect size for Reward.EC:Cue.EC is -0.037

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

Time elapsed: 0 h 1 m 38 s

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

When I try to extend the model to see the power I would get with 80 participants using the below code, I get an error:

model2<-extend(random_intercept_ACC_mlm_lmer, along="id_num", n=80) powerSim(model2) Error in getDefaultXname(fit) : Couldn't automatically determine a default fixed effect for this model.

When I try to specify the test, I still get an error

model2<-extend(random_intercept_ACC_mlm_lmer, along="id_num",test=fixed("Reward.EC:Cue.EC", method="z"), n=80) Error in extend(random_intercept_ACC_mlm_lmer, along = "id_num", test = fixed("Reward.EC:Cue.EC", : unused argument (test = fixed("Reward.EC:Cue.EC", method = "z"))

Any ideas? Thanks in advance!

[pitakakariki]

You need to specify the test in powerSim, not in extend.

[maureenbowers]

Great! Thank you so much. That worked.

I know have the following code:

Running the Power Analysis

powerSim(random_intercept_ACC_mlm_lmer, test=fixed("Reward.EC:Cue.EC", method="z"))

Increasing the sample size

model2<-extend(random_intercept_ACC_mlm_lmer, along="id_num", n=80) powerSim(model2, test=fixed("Reward.EC:Cue.EC", method="z"))

Then, I was trying to produce a powercurve. When I do this: pc2<-powerCurve(model2, breaks=c(10,20,30,40,50,60,70,80))

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

Do I need to specify the test in powerCurve too? Even when I try: pc2<-powerCurve(model2,test=fixed("Reward.EC:Cue.EC", method="z"), breaks=c(10,20,30,40,50,60,70,80)) OR pc2<-powerCurve(model2,test=fixed("Reward.EC:Cue.EC"), breaks=c(10,20,30,40,50,60,70,80))

I still get the above error about not automatically determining a fixed effect.

This is my first time using your package (with somewhat limited R knowledge), so I appreciate your help!

[pitakakariki]

I think you need to manually specify along=id_num for this model.

[MichelMcD]

I have a model with two factors: a (4 levels) and b (3 levels). Each participant receives two problems (id as a random effect). I want to estimate sample size based on pilot data to give me sufficient power to detect an interaction effect.

The model would be as follows:

mod1 <- lmer(outcome ~ a*b+ (1 | id), data).

If I run a powerSim on the interaction given the current data (without changing effect sizes), the power is ~ 66% using the following:

powerSim(mod1, fixed(“a:b”, "lr"), nsim=100)

And is similar if I were to test the interaction using fcompare (to the main effects):

powerSim(mod1, fcompare(~ a+b), nsim=100)

Now I know that estimating power based on one’s own data is not good, so I want to change the effect size estimates. However, I am unsure whether I need to change each fixed effect (e.g., a1:b1, a2:b1…) separately, or if there is a way to test the whole interaction with one test (e.g., like the fcompare test above)? If I do need to change the fixed effect estimate for each of the interaction terms, do I do so individually, estimate the power/sample size, and repeat for each fixed effect or do all changes in the same model (e.g., change all fixed effects)?

e.g.,:

fixef(mod1)["a1:b1"] <- 2 fixef(mod1)["a1:b2"] <- 2.1 ...

The next step would then be to extend the model and increase the sample within each group, for example:

mod2 <- extend(mod1, within=“a+b”, n=100)

I have seen the interaction term tested when one of the predictors is continuous, but not when both are factors. Any help would be much appreciated.

[palday]

@MichelMcD For a 2x3 design, there are two interaction terms in the model. Since the fixef(model) bit works by modifying terms in a model, you have to modify those individually.

Depending on your contrast coding, you may also want to change the size of your main effects as well so that your actually manipulating the part(s) of the interaction you care about

[palday]

@pitakakariki This thread is slowly but surely turning into a general help forum, I think it might be best to close it and suggest an alternative venue for general help vs. issues in (verses in using) the package.

[MichelMcD]

@pitakakariki Thank you. That was my hunch but I was unsure whether there was an alternative for the whole interaction term (versus the individual fixed effects for the interactions). I will take this approach.

[adrose]

I would normally test the entire interaction, but it really depends on what you're trying to do.

You could try a range of effect sizes until you get about 80% power - this would give you an idea of your Minimal Detectable Effect.

can you expand on what you mean by this?