pitakakariki / simr

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

powerSim freezing #150

Open charlijk opened 5 years ago

charlijk commented 5 years ago

I have been using powerSim on my data and it is working perfectly when I am looking at the power of categorical variables with two or three levels to them, however continuous variables and a categorical variable with 14 levels cause powerSim to freeze and will result in an error message from R. I was wondering if anyone had a solution for this? I have a large data set with around 14,000 results for each variable.

Thank you.

mengzhu-nz commented 5 years ago

Sorry this is not an answer. I am shamelessly using the discussion here to ask a question that I asked about a week ago with no feedback yet: https://github.com/pitakakariki/simr/issues/91.

As powerSim is working perfectly for you with categorical variables, could I ask a question that relates to this? I have a linear mixed effects model with two categorical variables (3 and 4 levels) and their interaction in the fixed effects. I am trying to estimate the power for this two-way interaction term, but I always have an error. Do you have an idea how to fix this? Thanks in advance!

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

charlijk commented 5 years ago

Hello, I have removed my interactions from my model as I found it resulted in a similar issue, even when they weren't the variables I was testing the power for, so I would also be happy to have an answer to your question as it's another issue I've faced!

mengzhu-nz commented 5 years ago

Hi Thank you so much for your prompt reply! This is reallly helpful. Yes, when I removed the interaction and tested the power for the single effects one by one, powerSim returned a reasonable power but with a warning about "observed power" calculation. Did you have a similar warning as well?

This is my code and the output: 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

mengzhu-nz commented 5 years ago

Hi if you make the interaction a single term using "df$AB = interaction(df$A, df$B)", the problem should be solved!

pitakakariki commented 5 years ago

Are you sure it's freezing? Try e.g. nsim=3 to see if it's just taking a long time.

What is the error message you are getting?

camilo-r-r commented 5 years ago

Hi everyone,

I guess I'll also use this thread to ask a related question. running powerSim is freezing. Even when i try to run 2 simulations it doenst work. I have a model that looks like this:

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: logratio ~ status.numeric * verb.numeric + (verb.numeric + status.numeric ||  
    subj) + (status.numeric || item)
   Data: data3

     AIC      BIC   logLik deviance df.resid 
 14017.8  14091.9  -6998.9  13997.8    12182 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.98777 -0.59118 -0.01071  0.77250  2.38581 

Random effects:
 Groups   Name           Variance  Std.Dev.
 item     status.numeric 0.0070597 0.08402 
 item.1   (Intercept)    0.0055717 0.07464 
 subj     status.numeric 0.0022048 0.04696 
 subj.1   verb.numeric   0.0009728 0.03119 
 subj.2   (Intercept)    0.0120178 0.10963 
 Residual                0.1804459 0.42479 
Number of obs: 12192, groups:  item, 36; subj, 12

Fixed effects:
                                Estimate   Std. Error           df t value            Pr(>|t|)    
(Intercept)                     0.067156     0.034229    15.648670   1.962              0.0678 .  
status.numeric                  0.048668     0.019879    32.816828   2.448              0.0199 *  
verb.numeric                   -0.020590     0.009818    11.248360  -2.097              0.0593 .  
status.numeric:verb.numeric    -0.064259     0.003916 12155.395970 -16.408 <0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) stts.n vrb.nm
status.nmrc -0.001              
verb.numerc  0.002 -0.002       
stts.nmrc:. -0.002  0.011 -0.021

I wanna run a simulation on the interaction term:

power_1 <- powerSim(m.onset1.sim, fixed("status.numeric:verb.numeric"),nsim = 2)

but the progress bar doesnt move at all and my computer starts purring :) Any ideas on what I'm doing wrong?

pitakakariki commented 5 years ago

I'd start by running through a single simulation to see which step your computer is struggling with:

y <- doSim(m.onset1.sim)
z <- doFit(y, m.onset1.sim)
doTest(z, fixed("status.numeric:verb.numeric"))

If the model fitting in lme4 is working fine, then I suspect it will be the testing step that's the problem. In that case you might want to specify a less computationally intensive test e.g.:

fixed("status.numeric:verb.numeric", "z")
jhr87 commented 3 years ago

I have been using powerSim on my data and it is working perfectly when I am looking at the power of categorical variables with two or three levels to them, however continuous variables and a categorical variable with 14 levels cause powerSim to freeze and will result in an error message from R. I was wondering if anyone had a solution for this? I have a large data set with around 14,000 results for each variable.

Thank you.

Hello,

I had the same problem and I fixed it by removing the NA values:

df <- df[!is.na(df$var1),]
model <- lmer(dv ~ var1 + var2 + (1|subject) + (1|item), data = df)
powerSim(model, test=fixed("var1"))