pitakakariki / simr

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

Unhelpful error messages for not-quite-a-model objects #285

Open bronson-hui opened 4 days ago

bronson-hui commented 4 days ago

I am struggling to get the power analysis for an interaction term (categorical : categorical) to run. Here is the model summary:

> mod.sum
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: dv ~ stim_version.e * group + (1 | sub_id) + (1 | stim_id)
   Data: d

     AIC      BIC   logLik deviance df.resid 
 -554855  -554781   277435  -554869   319993 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.753 -0.674  0.000  0.675  4.386 

Random effects:
 Groups   Name        Variance Std.Dev.
 stim_id  (Intercept) 0.000205 0.0143  
 sub_id   (Intercept) 0.001799 0.0424  
 Residual             0.010246 0.1012  
Number of obs: 320000, groups:  stim_id, 800; sub_id, 200

Fixed effects:
                            Estimate    Std. Error            df t value Pr(>|t|)    
(Intercept)                 0.297101      0.004286    207.154789   69.32   <2e-16 ***
stim_version.e             -0.100160      0.000506 319000.994191 -197.90   <2e-16 ***
groupNS                     0.001655      0.006019    201.421052    0.27     0.78    
stim_version.e:groupNS     -0.399668      0.000716 319000.994213 -558.39   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) stm_v. gropNS
stim_versn. -0.059              
groupNS     -0.702  0.042       
stm_vrs.:NS  0.042 -0.707 -0.059

Here is the code I use for 'powerSim'

sim1 <- powerSim(mod.sum, 
                 nsim= 100, 
                 test = simr::fixed('stim_version.e:group', 'z'),
                 progress = T)
sim1

The output is always 0% power with all simulations resulted in errors

> sim1
Power for predictor 'stim_version.e:group', (95% confidence interval):
       0.00% ( 0.00,  3.62)

Test: z-test

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

Time elapsed: 0 h 0 m 0 s

nb: result might be an observed power calculation

The errors are as follow:

                                                                                                           message
1   no applicable method for 'simulate' applied to an object of class "c('summary.lmerModLmerTest', 'summary.merMod')"

Any pointers are appreciated.

pitakakariki commented 4 days ago

Looks like you're passing the model summary to powerSim rather than the model itself?

bronson-hui commented 4 days ago

I am so embarrassed. OMG! But, it now gives me a different error

sim1 <- powerSim(mod, 
                 nsim= 10, 
                 test = simr::fixed('stim_version.e:group', 'z'),
                 progress = T)

same issue with all sims resulted in errors:

sim1
Power for predictor 'stim_version.e:group', (95% confidence interval):
       0.00% ( 0.00, 30.85)

Test: z-test

Based on 10 simulations, (0 warnings, 10 errors)
alpha = 0.05, nrow = 320000

Time elapsed: 0 h 1 m 30 s

nb: result might be an observed power calculation

The error is different though:

> sim1$errors
     stage index                 message
1  Testing     1 subscript out of bounds
pitakakariki commented 4 days ago

Looks like your interaction term involves multiple fixed effects parameters - in which case a z-test won't work. Try 'lr' for the test instead of 'z'?

bronson-hui commented 3 days ago

it works.. thanks for the support!