0

I would like to do a power simulation in R to plan the sample size for my master thesis study in Psychology. I'm using a linear mixed model with two fixed effects, conflict and openness, and random effects for subjects (persons) and item (posts).

The simulation produces different error messages and I don't know how to fix it:

library(lme4)

### effect sizes
## random effects
sigma <- 0.59 # Standard deviation
ssubj <- 0.16  # person
sscen <- 0.05  # item
theta <- c(sscen, ssubj)

# intercept, main effect conflict, main effect openness, interaction
beta <- c(4.305,-0.27,-0.11,-0.3)

par <- list(beta = beta, theta = theta, sigma = sigma)




n <- 200 # persons
n_post = 21 # number of posts

newdat <- data.frame(subjectID = rep(1:n, each = n_post),
                     postID = 1:n_post,
             conflict = runif(n,0,8),  
             openness = runif(n,1,5)
)



library(optimx)  
pval <- replicate(200, {
  sim <- simulate(~ conflict * openness + (1 | postID) + (1 | subjectID), newparams = par, newdata = newdat, REML=FALSE, control = lmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')))$sim_1
  m0  <- lmer(sim ~ conflict + openness + (1 | postID) + (1 | subjectID), newdat, REML=FALSE, control = lmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')))
  m1  <- lmer(sim ~ conflict * openness + (1 | postID) + (1 | subjectID), newdat, REML=FALSE, control = lmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')))
  anova(m0, m1)["m1", "Pr(>Chisq)"]
} 
)

#theta parameter vector not named: assuming same order as internal vector
#beta parameter vector not named: assuming same order as internal vector
#boundary (singular) fit: see ?isSingular
#There were 50 or more warnings (use warnings() to see the first 50)

summary(m1)

hist(pval)

# Simulated power at a = 5%
# -> goal 0.8
mean(pval < .05)
# 1 ?!

I expect to find out how big the sample size has to be to have a power of 0.8, at the moment I always get a power of 1, which can't be correct.

I tried different numbers for n (persons) and the number of replications in the simulation. I also tried different convergence controls:

control = lmerControl(optimizer = "nloptwrap")
control = lmerControl(optimizer ="Nelder_Mead")
control = lmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')))
control = lmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))

I get different warning messages, e.g.:

Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
Model failed to converge with max|grad| = 0.0350272 (tol = 0.002, component 1)
Anne B
  • 1
  • 1
  • This is as much (or more) of a statistical question as a programming/coding one. So maybe consider taking this to [CrossValidated](stats.stackexchange.com). Regarding the convergence issues, have a look at [this helpful page](https://rstudio-pubs-static.s3.amazonaws.com/33653_57fc7b8e5d484c909b615d8633c01d51.html) in case you don't know it yet. Also, one simple reason for convergence trouble can be big differences between the scales of predictors, e.g. if one predictor ranges 0-2 and others 2000-10'000. – benimwolfspelz Dec 05 '22 at 13:09

0 Answers0