pitakakariki / simr

Power Analysis of Generalised Linear Mixed Models by Simulation
68 stars 19 forks source link

Errors when testing a factorial model in powerSim #267

Open FrVitale opened 9 months ago

FrVitale commented 9 months ago

Hi, I cannot achieve runnig a correct power analysis. When I was doing Power analysis, the result of Power was always 0%. I don't know where the problem is? Here is my code and results My model: m_GO<-lmer(MEP ~ Stim*ToV*Pola +(1|Sujeto)+0, data = GO);

and this is the power analysis: pa<-powerSim(m_GO, test=fixed("Stim"), nsim=100).

And this is the output I get:

"Power for predictor 'Stim', (95% confidence interval):
       0.00% ( 0.00,  3.62)

Test: Likelihood ratio

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

Time elapsed: 0 h 0 m 8 s

nb: result might be an observed power calculation"

I tried to change the fixed effect of "Stim" with this code: fixef(m_GO)["Stim"] <- 1.2

but I get such result:

Error in `fixef<-`(`*tmp*`, value = c(StimICF = 1.8531818181818, StimSICI = 0.497272727272712,  : 
  Stim is not the name of a fixed effect.

Then i try this code, to see which are my fixed effect: fixef(m_GO)

and here are the results:

 StimICF                 StimSICI                   StimSP 
              1.85318182               0.49727273               1.54500000 
                 ToVATTE                  PolaPOS         StimSICI:ToVATTE 
              0.04954545               0.02000000              -0.03727273 
          StimSP:ToVATTE         StimSICI:PolaPOS           StimSP:PolaPOS 
              0.04136364              -0.02818182               0.09954545 
         ToVATTE:PolaPOS StimSICI:ToVATTE:PolaPOS   StimSP:ToVATTE:PolaPOS 
             -0.05454545               0.05636364              -0.08136364 

How can I solve this problem?

Thank you!

pitakakariki commented 9 months ago

It looks like Stim is a factor, with three levels? If you try fixef(m_GO)["StimICF"] <- 1.2 instead that should work, although it will only change the coefficient for that level of the factor - I'd recommend thinking carefully about all 12 coefficients.

The errors in powerSim probably come from the test specification. You can use doTest to find one that works before running the simulations again. In this case you might want something like: doTest(m_GO, fcompare(~ ToV * Pola, "kr")).

If that doesn't work head(pa$errors) should give more information.

FrVitale commented 8 months ago

Thank you for your response. I did run again the simulation with the specified test: pa<-powerSim(m_GO, fixed("Stim","lr"), nsim=100)

And I receive the same error.

Running head(pa$errors), I receive this output, which I actually don't understand: stage index message 1 Testing 1 Test did not return a p-value 2 Testing 2 Test did not return a p-value 3 Testing 3 Test did not return a p-value 4 Testing 4 Test did not return a p-value 5 Testing 5 Test did not return a p-value 6 Testing 6 Test did not return a p-value

pitakakariki commented 8 months ago

You could try one of these tests which should work.

doTest(m_GO, fcompare(~ ToV * Pola, "kr"))
doTest(m_GO, fixed("Stim", "anova"))

Be careful that you're using the test that's appropriate for you research question though - testing main effects when there are interactions is tricky.