pitakakariki / simr

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

Problems with building power sim from scratch #102

Closed psyargh closed 2 years ago

psyargh commented 6 years ago

Hello,

I've had problems building my power simulation from scratch.

This is my model:

model1 <- makeLmer(y ~ (x1+x2)*x3*x4 + (1|subject) + (1|stimset) + (1|stimuli), fixef=b, VarCorr=V2, sigma=s, data=mydata)

These are the parameters I estimated:

b <- c(1.04, 1.73, -0.17, 0.3, 0.4, 0.2, 0.3, 0.1, -0.2, -0.1, 0.3, 1.2) # fixed intercept and slope
V2 <- matrix(c(0.12, 0.06, 0.06, 0.06, 0.12, 0.06, 0.06, 0.06, 0.12), 3) # random intercept and slope variance-covariance matrix
s <- 0.28 # residual variance

This is the error message I received:

Error in setParams(object, newparams) : length mismatch in theta (6!=3)

Could you please help me identify and solve the problem?

Thank you!

pitakakariki commented 6 years ago

You have three independent random intercepts in your formula, so V2 should be a list of three random effect variances. e.g.:

V2 <- list(0.12, 0.12, 0.12)
psyargh commented 6 years ago

Thanks! However, once I do that, I get the following error:

Error in chol.default(V, pivot = TRUE) : 'a' must be a square matrix

pitakakariki commented 6 years ago

What does your V2 look like now?

psyargh commented 6 years ago

Apologies-- I had retained the matrix command, but I replaced it with "list" and now it works. Thanks!

However, when I do powerSim, it gives me the power of 0% no matter was I specify as my parameters. Here's how it looks:

fixef(model1) ["(x1+x2)*x3*x4"] <- 0.05
powerSim(model1, fixed("(x1+x2)*x3*x4"), nsim=10)
pitakakariki commented 6 years ago

You need to specify effect names like this:

fixef(model1) ["x1"] <- ...
fixef(model1) ["x3:x4"] <- ...

And the same applies to fixed. Did you not get an error from your fixef<- call?

psyargh commented 6 years ago

Thanks! I did get an error from the fixef <- call.

Do I need to specify effect names for all main effects and 2-way interactions separately, even though I only want my power analysis to be based on the potential to capture the effect of the 3-way interaction?

pitakakariki commented 6 years ago

You should be okay just specifying the effect sizes you're interested in. Note that you have two 3-way interactions though, "x1:x3:x4" and "x2:x3:x4".