Brandon LeBeau

3 minute read

This is a quick note looking for any further feedback on the simglm package prior to CRAN submission later this week. The goal is to submit Thursday or Friday this week. The last few documentation finishing touches are happening now working toward a version 0.5.0 release on CRAN.

For those who have not seen this package yet, the aim is to simulate regression models (single level and multilevel models) as well as employ empirical power analyses based on Monte Carlo simulation. The package is relatively flexible in user control of inputs to generate data.

To install the package and also build the vignettes:

devtools::install_github("lebebr01/simglm", build_vignettes = TRUE)

Then to generate a simple single level data set.

library(simglm)

fixed <- ~ 1 + act + diff + numCourse + act:numCourse
fixed_param <- c(2, 4, 1, 3.5, 2)
cov_param <- list(dist_fun = c('rnorm', 'rnorm', 'rnorm'), 
                  var_type = c("single", "single", "single"),
                  cov_param = list(list(mean = 0, sd = 4),
                                   list(mean = 0, sd = 3),
                                   list(mean = 0, sd = 3)))
n <- 150
error_var <- 3
with_err_gen = 'rnorm'
temp_single <- sim_reg(fixed = fixed, fixed_param = fixed_param, 
                       cov_param = cov_param,
                       n = n, error_var = error_var, 
                       with_err_gen = with_err_gen, 
                       data_str = "single")
head(temp_single)
## # A tibble: 6 x 9
##   X.Intercept.         act       diff   numCourse act.numCourse      Fbeta
##          <dbl>       <dbl>      <dbl>       <dbl>         <dbl>      <dbl>
## 1            1  0.04045215 -0.7738108  0.02183119  0.0008831184  1.4661731
## 2            1  0.09351232  1.3360941  1.34312887  0.1255991021  8.6622926
## 3            1  0.06066363 -1.2230152  1.26423076  0.0766928310  5.5978327
## 4            1  0.01352531  0.3356532 -0.12900022 -0.0017447676  1.9347642
## 5            1 -0.91711660 -0.3569118  0.85737040 -0.7863086260 -0.5971991
## 6            1  1.40470811  1.8956529 -0.83768332 -1.1767005519  4.2291926
## # ... with 3 more variables: err <dbl>, sim_data <dbl>, ID <int>

Then one can extend this to conduct of power analysis. The benefit of this approach is that you are able to generate data that hopefully more closely resembles data that is to be collected and can also evaluate assumption violations, sample size differences, and other conditions directly into the power analysis similar to a Monte Carlo simulation.

fixed <- ~ 1 + act + diff + numCourse + act:numCourse
fixed_param <- c(0.5, 1.1, 0.6, 0.9, 1.1)
cov_param <- list(dist_fun = c('rnorm', 'rnorm', 'rnorm'), 
                  var_type = c("single", "single", "single"),
                  opts = list(list(mean = 0, sd = 2),
                              list(mean = 0, sd = 2),
                              list(mean = 0, sd = 1)))
n <- NULL
error_var <- NULL
with_err_gen <- 'rnorm'
pow_param <- c('(Intercept)', 'act', 'diff', 'numCourse')
alpha <- .01
pow_dist <- "t"
pow_tail <- 2
replicates <- 10
terms_vary <- list(n = c(20, 40, 60, 80, 100), error_var = c(5, 10, 20),
                   fixed_param = list(c(0.5, 1.1, 0.6, 0.9, 1.1), 
                                      c(0.6, 1.1, 0.6, 0.9, 1.1)),
                cov_param = list(list(dist_fun = c('rnorm', 'rnorm', 'rnorm'),
                                       mean = c(0, 0, 0), sd = c(2, 2, 1), 
                                  var_type = c("single", "single", "single")),
                                  list(dist_fun = c('rnorm', 'rnorm', 'rnorm'),
                                       mean = c(0.5, 0, 0), sd = c(2, 2, 1), 
                                  var_type = c("single", "single", "single"))
                                  )
                   )
power_out <- sim_pow(fixed = fixed, fixed_param = fixed_param, 
                     cov_param = cov_param,
                     n = n, error_var = error_var, with_err_gen = with_err_gen, 
                     data_str = "single", pow_param = pow_param, alpha = alpha,
                     pow_dist = pow_dist, pow_tail = pow_tail, 
                     replicates = replicates, terms_vary = terms_vary)
head(power_out)
## # A tibble: 6 x 11
## # Groups:   var, n, error_var, fixed_param [3]
##           var     n error_var         fixed_param
##        <fctr> <dbl>     <dbl>              <fctr>
## 1 (Intercept)    20         5 0.5,1.1,0.6,0.9,1.1
## 2 (Intercept)    20         5 0.5,1.1,0.6,0.9,1.1
## 3 (Intercept)    20         5 0.6,1.1,0.6,0.9,1.1
## 4 (Intercept)    20         5 0.6,1.1,0.6,0.9,1.1
## 5 (Intercept)    20        10 0.5,1.1,0.6,0.9,1.1
## 6 (Intercept)    20        10 0.5,1.1,0.6,0.9,1.1
## # ... with 7 more variables: cov_param <fctr>, avg_test_stat <dbl>,
## #   sd_test_stat <dbl>, power <dbl>, num_reject <dbl>, num_repl <dbl>,
## #   data <list>

Questions and feedback are welcomed by filing an issue on GitHub here: https://github.com/lebebr01/simglm/issues.

comments powered by Disqus