pitakakariki / simr

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

Couldn't automatically determine a default fixed effect for this model #96

Closed rsiugzda closed 4 years ago

rsiugzda commented 6 years ago

library(simr) fit <- lmer(oxy ~ 1 + (1|subject), data = df) power <- powerSim(fit,nsim = 200)

we got an error: Error in getDefaultXname(fit) : Couldn't automatically determine a default fixed effect for this model.​

any ideas what did go wrong?

rsiugzda commented 6 years ago

I saw previous comments, i need to specify test= variable.

pitakakariki commented 6 years ago

Looks like I probably need to write a better error message :)

rsiugzda commented 6 years ago

hi, we tried in a different way: fit <- lmer(oxy ~ (task|subject) + hemi, data = df)

and it gives the same error... help!!

pitakakariki commented 6 years ago

Have you specified a test= argument?

rsiugzda commented 6 years ago

yes, and now we got another error:

power <- powerSim(lm_ABC, test = "task",nsim=200) Warning message: In observedPowerWarning(sim) : This appears to be an "observed power" calculation​

pitakakariki commented 6 years ago

That's a warning rather than an error. It will go away if you specify the relevant parameter (e.g. with fixef<-)

rsiugzda commented 6 years ago

sorry, could you be more specific, pls? if I had this: power <- powerSim(lm_ABC, test = "task", nsim=200)

how it can be changed? if in lm_ABC i have 3 fixed effects, but I want to test only one (task).

pitakakariki commented 6 years ago

You want something like: power <- powerSim(lm_ABC, test=fixed("task"), nsim=200)

You'll also want something like fixef(lm_ABC) <- ??? based on the size of the effect you'd like to be able to detect.

rsiugzda commented 6 years ago

thanks we manage to make it work ;)

Porites commented 6 years ago

I also received this error message and below is my code for both the model and the powerCurve function. I realize that this is a more complex model than usual so perhaps powerCurve doesn't work at this level.

fit1<-lmer(Coral ~ WYear*Frames*Points + (1+WYear|Trans2) +(1|Year), data=Benthic)
library(simr)
fixef(fit1)["WYear"] <- -0.25  #Changes the size of the fixed effect for this variable
pc1 <- powerCurve(fit1, test=fixed("WYear"), nsim=50)

Oddly enough, powerSim seems to work with the same syntax, but I do receive a large number of errors. BTW, I am expecting very low power among factor levels due to the small differences among years.

## Power for predictor 'WYear', (95% confidence
## interval):===========================================================|
##       0.00% ( 0.00,  7.11)
## Test: Kenward Roger (package pbkrtest)
##    Effect size for WYear is -0.25
## Based on 50 simulations, (50 warnings, 50 errors)
## alpha = 0.05, nrow = 4320
## Time elapsed: 0 h 1 m 49 s 

Any suggestions on the syntax or what might be causing the error message? Thanks

pitakakariki commented 6 years ago

What does head(pc1$errors) say?

I suspect the problem is testing the main effect WYear when there are interactions. Perhaps you want test=fcompare(~ Frames*Points)?

Porites commented 6 years ago

Peter,

Thanks for the quick and personal reply. R responses are below in red with comments in blue. In case you are color blind I added a dash (-) for R responses and an asterisk (*) for my comments.

What does head(pc1$errors) say?

I suspect the problem is testing the main effect WYear when there are interactions. Perhaps you want test=fcompare(~ Frames*Points)?

Hope this info helps and I appreciate the feedback.

Eric

On Tue, Jul 31, 2018 at 12:47 PM, Peter Green notifications@github.com wrote:

What does head(pc1$errors) say?

I suspect the problem is testing the main effect WYear when there are interactions. Perhaps you want test=fcompare(~ Frames*Points)?

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

-- Eric Brown, Ph.D. Marine Ecologist Kalaupapa NHP 7 Puahi Street P.O. Box 2222 Kalaupapa, HI 96742 (808)-567-6802, x1502 FAX (808)-567-6682

pitakakariki commented 6 years ago

NB: it's better to respond on the gihub website, as things like colour get stripped if you try to reply to emails.

Looks like the problem's fairly simple - you haven't specified an along argument, and simr's very simple heuristics can't guess one when the first term has interactions.

Does it work if you specify along="WYear" (or whichever variable is appropriate in your case)?

Porites commented 6 years ago

The next iteration of the code is pc1 <- powerCurve(fit1, extend(along="WYear"), nsim=25) but this yielded the same error message.

I also tried pc1 <- extend(fit1, along="WYear", n=20) pc2 <- powerCurve(pc1, fixed("WYear"), nsim=25) thinking that I have to sequence the steps. pc1 was generated, but the powerCurve function returned the same error message. Other suggestions?

Porites commented 6 years ago

Peter,

I also responded in gitHub so people could follow the thread.

The next iteration of the code is pc1 <- powerCurve(fit1, extend(along="WYear"), nsim=25) but this yielded the same error message.

I also tried pc1 <- extend(fit1, along="WYear", n=20) pc2 <- powerCurve(pc1, fixed("WYear"), nsim=25) thinking that I have to sequence the steps. pc1 was generated, but the powerCurve function returned the same error message. Other suggestions?

Eric

On Tue, Jul 31, 2018 at 3:27 PM, Peter Green notifications@github.com wrote:

NB: it's better to respond on the gihub website, as things like colour get stripped if you try to reply to emails.

Looks like the problem's fairly simple - you haven't specified an along argument, and simr's very simple heuristics can't guess one when the first term has interactions.

Does it work if you specify along="WYear" (or whichever variable is appropriate in your case)?

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

-- Eric Brown, Ph.D. Marine Ecologist Kalaupapa NHP 7 Puahi Street P.O. Box 2222 Kalaupapa, HI 96742 (808)-567-6802, x1502 FAX (808)-567-6682

pitakakariki commented 6 years ago
pc1 <- powerCurve(fit1, along="WYear", test=..., nsim=25)
Porites commented 6 years ago

This code worked

pc1 <- powerCurve(fit1, along="WYear", test=fcompare(~ Frames*Points), nsim=25)

Thanks

On Tue, Jul 31, 2018 at 4:22 PM, Peter Green notifications@github.com wrote:

pc1 <- powerCurve(fit1, along="WYear", test=..., nsim=25)

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

-- Eric Brown, Ph.D. Marine Ecologist Kalaupapa NHP 7 Puahi Street P.O. Box 2222 Kalaupapa, HI 96742 (808)-567-6802, x1502 FAX (808)-567-6682

AbsoluteBeginners commented 6 years ago

Peter thanks for all this info, it is very helpful. I have a problem when running a model with 4 factors and interactions

fit <- lmer(model ~ A*B*C*D+ (1|Subject), data = Example)

When I am running the power analysis for any main effect or interaction, for example:

power <- powerSim(fit, test=fixed("A*B"), nsim = 200)

I receive this:

Power for predictor 'A*B', (95% confidence interval):
       0.00% ( 0.00,  1.83)
Test: Kenward Roger (package pbkrtest)
Based on 200 simulations, (0 warnings, 200 errors)
alpha = 0.05, nrow = 984
Time elapsed: 0 h 0 m 19 s
nb: result might be an observed power calculation

Note that the A*B interaction is significant (t=4,29)

What am I doing wrong?

AB

pitakakariki commented 6 years ago

Have a look at the names in summary(fit)$coef to see which strings are valid for fixed.

Note that you'll probably need to specify a method too, e.g.:

fixed("Asomething:Bsomething", "z")
jrode11 commented 5 years ago

Hello,

I'm also getting the error: "Couldn't automatically determine a default fixed effect for this model."

I ran this mixed model from the lmerTest package: mod2 <- lmer(anxiety ~ output*time*session + (time|id), data=long2.t, REML=TRUE) ,

where "output" "time" and "session" are all factor variables, and there is a random slope and intercept for time.

The powerSim function works: works<-powerSim(mod2, test=fixed("outputspss:timepost:session2", method="t"), nsim=20)

But the powerCurve function doesn't work: doesnt<-powerCurve(mod2, test=fixed("outputspss:timepost:session2", method="t"), nsim=20)

Sorry if I'm missing something here, but looking for some clarification.

Thank you!! Jacob

pitakakariki commented 5 years ago

I think you just need to specify an along argument for powerCurve.

jrode11 commented 5 years ago

That worked, thank you!

MattK99 commented 5 years ago

G'day

I have a similar issue to the one above by AbsoluteBeginners.

Model<-lmer(log(Catch2)~B_A+C_I+B_A:C_I+(1 | id)+(1 | Period),data=AllData)

pow<-powerSim(fixef(Model), test=fixed("B_A:C_I", "anova"), nsim=100, alpha=0.1) ##https://www.rdocumentation.org/packages/simr/versions/1.0.4/topics/powerSim Warning message:====| In observedPowerWarning(sim) : This appears to be an "observed power" calculation pow Power for predictor 'B_A:C_I', (95% confidence interval): 0.00% ( 0.00, 3.62)

Test: Type-I F-test

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

Time elapsed: 0 h 0 m 0 s

nb: result might be an observed power calculation

I tried summary(Model) Which gave me the names B_ABefore:C_IImpacted

Using those made no difference. Any help would be appreciated.

Thanks Matt

pitakakariki commented 5 years ago

Try:

pow <- powerSim(Model, ... 

not pow<-powerSim(fixef(Model), ....

MattK99 commented 5 years ago

Thanks for that.

I haven't seen it noted anywhere, but are you aware of a conflict with the package stringr? It results in the same problem I described above.

Thanks Matt

pitakakariki commented 5 years ago

Thanks, yes I'm aware of the clash.

benjaminwnelson commented 5 years ago

Hello, Thanks for all your helpful responses. I got the following error

Power for predictor 'group*Task', (95% confidence interval): 0.00% ( 0.00, 1.83) Test: unknown test Based on 200 simulations, (0 warnings, 200 errors) alpha = 0.05, nrow = NA Time elapsed: 0 h 0 m 12 s nb: result might be an observed power calculation

Here is the model hr_unadjusted <- lmer(HeartRate ~ quiet_hr_bl_centered + group*Task + (1 | subject_id), REML = "FALSE", contrasts = "TRUE", data = heart_rate_analyses)

Here is the power code power <- powerSim(hr_unadjusted, test = fixed("group*Task"), nsim = 200) power

Any idea how to fix this? Thanks!

pitakakariki commented 5 years ago

Looks like the problem is your test specification. Maybe e.g. test=fixed("group:Task", "kr")?

benjaminwnelson commented 5 years ago

This gives the following error message

Warning message:===========================================================================================| In observedPowerWarning(sim) : This appears to be an "observed power" calculation

Power for predictor 'group:Task', (95% confidence interval): 0.00% ( 0.00, 1.83)

Test: Kenward Roger (package pbkrtest)

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

Any ideas?

pitakakariki commented 5 years ago

Try some of the examples here: https://cran.r-project.org/web/packages/simr/vignettes/examples.html

MariaHolcekova commented 4 years ago

Hi,

I've been getting the same error when trying powerSim on my model:

model <- glmer(outcome ~ poly(predictor,2) + control + (1 | study), data = data, family = "binomial")

I've read the previous examples and tried the following: powerSim(model, test = fixed("predictor"), nsim = 200) powerSim(model, test ="predictor", nsim = 200) powerSim(model, nsim = 200)

I've also tried it with a model without a polynomial term but keep getting power either 0% or 100% depending on the specifications.

Any suggestions on how to make it work?

Thanks.

pitakakariki commented 4 years ago

Have a look at the Cake example in the link above. You'll probably need fcompare.

DanielUnn commented 4 years ago

Hi @pitakakariki

my model is

m4 <- glmer(formula='Revenue ~ VarA + VarB + VarC + (1 | VarD)', data=df_sample, family=gaussian(link="log"), nAGQ=0, control=glmerControl(optimizer="nloptwrap"), mustart=df_sample$Revenue)

I wanna powerSim(m4, nsim=10) to detect 5% or 10% treatment effect on Revenue, but I keep getting the same result as below.

Power for predictor 'VarA', (95% confidence interval):=========================================================================|
      100.0% (69.15, 100.0)

Test: Likelihood ratio

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

Time elapsed: 0 h 0 m 35 s

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

How should I do? Please guide me here, thanks.

I have manually added 5% treatment to 50% of the data randomly at column Revenue. 50% for treatment and 50% is not.