Сколько времени должно занять моделирование анализа мощности по методу Монте-Карло в R? Это потенциально часы? (1000 повторений, 1000 загрузок)

#r #simulation #montecarlo #statistics-bootstrap #power-analysis

#r #Симуляция #монтекарло #статистика-bootstrap #анализ мощности

Вопрос:

Я использую моделирование по методу Монте-Карло для выполнения анализа мощности для модели продольного посредничества. Я использую функцию power.boot из пакета bmem (lavaan).

Я проверил код только с помощью 5 повторений / 5 загрузок, чтобы убедиться, что он работает, и это сработало.

Затем я запустил код с 1000 повторений, 1000 загрузок, как рекомендует документация пакета.

Прошло уже больше часа, и он все еще работает — это нормально? Сколько времени слишком долго?

 powermodel1 <-'

x2 ~ start(.6)*x1   x*x1 
x3 ~ start(.6)*x2   x*x2

m2 ~ start(.15)*x1   a*x1   start(.3)*m1   m*m1 
m3 ~ start(.15)*x2   a*x2   start(.3)*m2   m*m2

y2 ~ start(.5)*m1   b*m2   start(.3)*y1   y*y1
y3 ~ start(.5)*m2   b*m2   start(.3)*y2   y*y2   start(0.05)*x1   c*x1

x1 ~~ start(.15)*m1

x1 ~~ start(.15)*y1


y1 ~~ start(.5)*m1

'
indirect <- 'ab:=a*b'

N<-200

system.time(bootstrap<-power.boot(powermodel1, indirect, N, nrep=1000, nboot=1000, parallel = 'multicore'))

            summary(bootstrap)

1. Каков «большой O» вашего алгоритма? heatonresearch.com/2015/03/21/r_big_o.html

2. Я не знаю, как это найти — я использовал код из опубликованной статьи, в которой приведены примеры того, как запускать эти типы анализа мощности, но они ничего не упоминали об этом

Ответ №1:

К сожалению, похоже, что это займет некоторое время; ~ 8 часов в моей системе:

 library(bmem)

powermodel1 <-'

x2 ~ start(.6)*x1   x*x1 
x3 ~ start(.6)*x2   x*x2

m2 ~ start(.15)*x1   a*x1   start(.3)*m1   m*m1 
m3 ~ start(.15)*x2   a*x2   start(.3)*m2   m*m2

y2 ~ start(.5)*m1   b*m2   start(.3)*y1   y*y1
y3 ~ start(.5)*m2   b*m2   start(.3)*y2   y*y2   start(0.05)*x1   c*x1

x1 ~~ start(.15)*m1

x1 ~~ start(.15)*y1


y1 ~~ start(.5)*m1

'
indirect <- 'ab:=a*b'

N<-200

system.time(bootstrap<-bmem::power.boot(powermodel1, indirect, N, nrep = 10, nboot = 10, parallel = 'multicore'))
system.time(bootstrap<-bmem::power.boot(powermodel1, indirect, N, nrep = 30, nboot = 30, parallel = 'multicore'))
system.time(bootstrap<-bmem::power.boot(powermodel1, indirect, N, nrep = 60, nboot = 60, parallel = 'multicore'))
system.time(bootstrap<-bmem::power.boot(powermodel1, indirect, N, nrep = 100, nboot = 100, parallel = 'multicore'))

library(tidyverse)
# Load the times from above into a dataframe
benchmark <- tibble(bootstraps = c(10, 30, 60, 100),
                    times = c(4.021, 30.122, 121.103, 311.236)) 

# Plot the points and fit a curve
ggplot(benchmark, aes(x = bootstraps, y = times))  
  geom_point()  
  geom_smooth(se = FALSE, span = 5)

 # Fit a model
fit <- lm(data = benchmark, times~poly(bootstraps,
                                       2, raw=TRUE))
newtimes <- data.frame(bootstraps = seq(100, 1000, length = 4))

# Predict the time it will take for larger bootstrap/rep values
predict(fit, newdata = newtimes)
>        1          2          3          4 
>  311.6829  4568.3812 13789.6754 27975.5655 

# Convert from seconds to hours
print(27975.5655/60/60)
>[1] 7.77099