Power with simglm - legacy code

2019-05-30

Power Analysis with simglm

The simglm package allows the ability to conduct a power analysis through simulation. This will be particularly helpful with multilevel models and generalized linear models. To show the process, we will start with basic regression models.

Single Level Power Analysis

Let’s look at a simple single level regression example to get started:

Much of the output here is the same from the sim_reg function. The additional arguments, pow_param represents the terms to conduct a power analysis for and must be a subset of the fixed argument, alpha represents the per term level of significance, pow_dist represents the sampling distribution to refer to, either ‘z’ or ‘t’, pow_tail represents whether a one or two tailed hypothesis is being tested, and replicates represents the number of simulations to conduct. Note, to do a power analysis for the intercept, ‘(Intercept)’ must be used. By default, if pow_param is not specified power is conducted for all terms.

Finally, looking at the output from the above call:

## # A tibble: 4 x 9
##   term  avg_test_stat sd_test_stat avg_std_err sd_std_err power num_reject
##   <chr>         <dbl>        <dbl>       <dbl>      <dbl> <dbl>      <dbl>
## 1 (Int…         0.497        0.398       0.368     0.0215     0          0
## 2 act           1.12         0.179       0.187     0.0148     0          0
## 3 diff          0.634        0.192       0.184     0.0134     0          0
## 4 numC…         0.972        0.347       0.375     0.0286     0          0
## # … with 2 more variables: num_repl <dbl>, data <list>

The output contains the variable name, the average test statistic, the standard deviation of the test statistic, the power rate, the number of null hypotheses rejects, and the total number of replications. Increasing the number of replications would increase the precision of the power analysis, however may significantly increase the computational time.

Standardized Coefficients

By default, the simglm package uses unstandardized regression coefficients when doing the simulation. A way to use standardized coefficients however, would be to generate standardized variables. For example:

## # A tibble: 4 x 9
##   term  avg_test_stat sd_test_stat avg_std_err sd_std_err power num_reject
##   <chr>         <dbl>        <dbl>       <dbl>      <dbl> <dbl>      <dbl>
## 1 (Int…         0.198       0.0970      0.0829    0.00539     0          0
## 2 act           0.400       0.0867      0.0837    0.00772     0          0
## 3 diff          0.262       0.0758      0.0834    0.00704     0          0
## 4 numC…         0.707       0.0852      0.0837    0.00716     0          0
## # … with 2 more variables: num_repl <dbl>, data <list>

Model Misspecification

It is also possible to specify a model different than the generating model for evaluation of power. This may be useful if it is thought another variable is important in explaining variation, however due to design issues is not possible to collect this variable. As a result, there would likely be additional variation due to this variable that will not be able to be explained. Building this into the generating model can help provide additional information about the impact on power.

An example using single level power analysis is shown below.

## # A tibble: 3 x 9
##   term  avg_test_stat sd_test_stat avg_std_err sd_std_err power num_reject
##   <chr>         <dbl>        <dbl>       <dbl>      <dbl> <dbl>      <dbl>
## 1 (Int…         0.470        0.472       0.416     0.0295     0          0
## 2 act           1.11         0.263       0.209     0.0168     0          0
## 3 diff          0.609        0.192       0.207     0.0193     0          0
## # … with 2 more variables: num_repl <dbl>, data <list>

You’ll notice, compared to above, the power is somewhat reduced for the other three fixed effects when not modeling the numCourse variable. You can also build in the lm_fit_mod argument into the design via the terms_vary argument.

## # A tibble: 6 x 13
## # Groups:   term, n, error_var, fixed_param [3]
##   term      n error_var fixed_param lm_fit_mod avg_test_stat sd_test_stat
##   <chr> <dbl>     <dbl> <fct>       <fct>              <dbl>        <dbl>
## 1 (Int…    20         5 0.5,1.1,0.… ~,sim_dat…         0.433        0.717
## 2 (Int…    20         5 0.5,1.1,0.… ~,sim_dat…         0.624        0.798
## 3 (Int…    20         5 0.6,1.1,0.… ~,sim_dat…         0.534        0.774
## 4 (Int…    20         5 0.6,1.1,0.… ~,sim_dat…         0.720        0.814
## 5 (Int…    20        10 0.5,1.1,0.… ~,sim_dat…         0.471        0.886
## 6 (Int…    20        10 0.5,1.1,0.… ~,sim_dat…         0.531        0.892
## # … with 6 more variables: avg_std_err <dbl>, sd_std_err <dbl>,
## #   power <dbl>, num_reject <dbl>, num_repl <dbl>, data <list>

Nested Data

Extending the power analysis to nested data is fairly straightforward. Below is a three level example

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00962439
## (tol = 0.002, component 1)
## # A tibble: 4 x 9
##   term  avg_test_stat sd_test_stat avg_std_err sd_std_err power num_reject
##   <chr>         <dbl>        <dbl>       <dbl>      <dbl> <dbl>      <dbl>
## 1 act            2.32       0.0150      0.0352   0.00237  0              0
## 2 actC…          7.14       0.255       0.300    0.0485   1              3
## 3 diff           5.99       0.0351      0.0251   0.000674 1              3
## 4 time           2.08       0.686       0.344    0.0745   0.333          1
## # … with 2 more variables: num_repl <dbl>, data <list>

Generalized Power Analysis

## # A tibble: 3 x 9
##   term  avg_test_stat sd_test_stat avg_std_err sd_std_err power num_reject
##   <chr>         <dbl>        <dbl>       <dbl>      <dbl> <dbl>      <dbl>
## 1 (Int…        0.0883        0.566       0.387     0.111      0          0
## 2 act          0.673         0.538       0.260     0.121      0          0
## 3 diff         0.386         0.243       0.134     0.0727     0          0
## # … with 2 more variables: num_repl <dbl>, data <list>

Vary Arguments

## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## # A tibble: 6 x 11
## # Groups:   term, n [3]
##   term      n fixed_param avg_test_stat sd_test_stat avg_std_err sd_std_err
##   <chr> <dbl> <fct>               <dbl>        <dbl>       <dbl>      <dbl>
## 1 (Int…    20 0.5,0.1,0.2         0.456        0.952       0.702     0.188 
## 2 (Int…    20 0.6,0.1,0.2         5.23        16.4      7809.    20316.    
## 3 (Int…    40 0.5,0.1,0.2         0.808        0.362       0.488     0.0681
## 4 (Int…    40 0.6,0.1,0.2         0.804        0.502       0.498     0.104 
## 5 (Int…    60 0.5,0.1,0.2         0.501        0.434       0.366     0.0423
## 6 (Int…    60 0.6,0.1,0.2         0.690        0.500       0.367     0.0504
## # … with 4 more variables: power <dbl>, num_reject <dbl>, num_repl <dbl>,
## #   data <list>