Очень высокие требования к памяти для моделирования мощности с использованием функции R simr :: powerSim - PullRequest
0 голосов
/ 13 апреля 2020

Я пытаюсь провести анализ мощности для проекта, где мы ожидаем 40 000 наблюдений (по крайней мере) на 40-60 "единицах".

Когда я делаю этот анализ мощности, используя гораздо меньшее число наблюдения в единицах (скажем, 100 на единицу), моделирование работает нормально. Однако требования к памяти, кажется, взрываются, когда я увеличиваю количество наблюдений на единицу до 300 или 500. Мне не удалось запустить симуляцию даже на компьютере с пределом памяти 70 ГБ.

Что происходит на? Почему требования к памяти так высоки и как я могу сделать их более управляемыми?

Ниже приведен общий скрипт имитации мощности c для справки.

#psacr001_power.R: Run a job using specifications determined using a shell script

#Get the arguments passed from the shell script
args <- commandArgs(trailingOnly=TRUE)
job <- as.numeric(args[[1]])
sample_size <- as.numeric(args[[2]])
number_of_labs <- as.numeric(args[[3]])
icc_size <- as.numeric(args[[4]])
slope_size <- as.numeric(args[[5]])
effect <- as.numeric(args[[6]])

#Required R libraries
library(lme4)
library(simr)
library(pbkrtest)
library(tibble)

#create dataset
simulated.df <- tibble(
                       x = sample(rep(0:1, sample_size * number_of_labs / 2)), #group assignment
                       g = rep(1:number_of_labs, sample_size), #number of clusters
                       y = rnorm(sample_size * number_of_labs) #outcome
                      )

#Create simple model
model.of.interest <- lmer(y ~ x + (x|g), data=simulated.df)

#Create a small effect of interest
fixef(model.of.interest)['x'] <- effect

#create various ICC intercept
VarCorr(model.of.interest)["g"][["g"]][1] <- icc_size

#create varioous slopes
VarCorr(model.of.interest)["g"][["g"]][4] <- slope_size

#try not to be singular
attr(VarCorr(model.of.interest)[["g"]], "correlation")[2:3] <- .3

power_summary <- tibble(
                        job = numeric(),
                        successes = numeric(),
                        trials = numeric(),
                        mean = numeric(),
                        lower = numeric(),
                        upper = numeric()
                       )

#simulate those models
temp <- powerSim(fit = model.of.interest,
                 nsim = 200,
                 progress=FALSE
                )

power_summary[1 , 1] <- job
power_summary[1 , 2:6] <- summary(temp)

write.csv(power_summary, paste(c("res_", job, ".csv"), collapse=""), row.names=FALSE)
...