<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/make-power-fun-again/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://educate-r.org//2017/02/24/csp2017.html"> Educate-R - R</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<section>
<h1>Make Power Fun (Again?)</h1>
<h2>Brandon LeBeau</h2>
<h3>University of Iowa</h3>
</section>
<section>
<h1>Overview</h1>
<ol>
<li>(G)LMMs</li>
<li>Power</li>
<li><code>simglm</code> package</li>
<li>Shiny Demo – Broken!<br>

</li>
</ol>
<section>
<h1>Linear Mixed Model (LMM)</h1>
<p><img alt="" src="https://i2.wp.com/educate-r.org/figs/equations.png?w=456" data-recalc-dims="1"><br>
</section>

<section>
<h1>Power</h1>
<ul>
<li>Power is the ability to statistically detect a true effect (i.e. non-zero population effect).</li>
<li>For simple models (e.g. t-tests, regression) there are closed form equations for generating power.
<ul>
<li>R has routines for these: <code>power.t.test, power.anova.test</code></li>
<li>Gpower3<br>

</li>
</ul>
</li>
</ul>
<section>
<h1>Power Example</h1>
<pre><code>n &lt;- seq(4, 1000, 2)
power &lt;- sapply(seq_along(n), function(i) 
  power.t.test(n = n[i], delta = .15, sd = 1, type = 'two.sample')$power)
</code></pre>
<p><img alt="" src="https://i1.wp.com/educate-r.org/figs/power_plot-1.png?w=456" data-recalc-dims="1"><br>
</section>

<section>
<h1>Power for (G)LMM</h1>
<ul>
<li>Power for more complex models is not as straightforward;
<ul>
<li>particularly with messy real world data.</li>
</ul>
</li>
<li>There is software for GLMM models to generate power:
<ul>
<li>Optimal Design: <a target="_blank" rel="nofollow" href="http://hlmsoft.net/od/">http://hlmsoft.net/od/</a></li>
<li>MLPowSim: <a target="_blank" rel="nofollow" href="http://www.bristol.ac.uk/cmm/software/mlpowsim/">http://www.bristol.ac.uk/cmm/software/mlpowsim/</a></li>
<li>Snijders, <em>Power and Sample Size in Multilevel Linear Models</em>.<br>

</li>
</ul>
</li>
</ul>
<section>
<h1>Power is hard</h1>
<ul>
<li>In practice, power is hard.</li>
<li>Need to make many assumptions on data that has not been collected.
<ul>
<li>Therefore, data assumptions made for power computations will likely differ from collected sample.</li>
</ul>
</li>
<li>A power analysis needs to be flexible, exploratory, and well thought out.<br>

</li>
</ul>
<section>
<h1>Power is Fun?</h1>
<ul>
<li>Three common reasons to do power analysis:
<ol>
<li>Power evidence for grant/planning</li>
<li>Post Hoc to explore insignificant results</li>
<li>Monte Carlo studies<br>

</li>
</ol>
</li>
</ul>
<section>
<h1><code>simglm</code> Overview</h1>
<ul>
<li><code>simglm</code> aims to simulate (G)LMMs with up to three levels of nesting (aim to add more later).</li>
<li>Flexible data generation allows:
<ul>
<li>any number of covariates and discrete covariates</li>
<li>change distribution of continuous covariates</li>
<li>change random distribution</li>
<li>unbalanced data</li>
<li>missing data</li>
<li>serial correlation<br>

</li>
</ul>
</li>
</ul>
<section>
<h1>Power with <code>simglm</code></h1>
<ul>
<li>Power with <code>simglm</code> takes on a Monte Carlo approach
<ul>
<li>This can provide a more thorough analysis/understanding of power.</li>
</ul>
</li>
<li>Always outputs a data frame
<ul>
<li>Useful for plotting</li>
<li>Data manipulation</li>
<li>etc.</li>
</ul>
</li>
<li>Serves as a wrapper around data generation process.<br>

</li>
</ul>
<section>
<h1>Power Analysis with <code>simglm</code></h1>
<ul>
<li>
<p>Factorial Design:</p>
<ol>
<li>Idenfity factors that influences power</li>
<li>Determine number of replications</li>
<li>Explore results</li>
</ol>
</li>
<li>
<p>Future Development</p>
<ol>
<li>Add ability for data generation and power model to differ<br>

</li>
</ol>
</li>
</ul>
<section>
<h1>Simple Example</h1>
<ul>
<li>Suppose we wished to generate data for a simple logistic regression.</li>
</ul>
<pre><code>library(simglm)

fixed &lt;- ~ 1 + act + diff
fixed_param &lt;- c(0.1, 0.5, 0.3)
cov_param &lt;- list(dist_fun = c('rnorm', 'rnorm'),
                  var_type = c(&quot;single&quot;, &quot;single&quot;),
                  opts = list(list(mean = 0, sd = 2),
                              list(mean = 0, sd = 4)))
n &lt;- 50
temp_single &lt;- sim_glm(fixed = fixed, fixed_param = fixed_param, 
                      cov_param = cov_param, 
                      n = n, data_str = &quot;single&quot;)
</code></pre>
</section>
<section>
<h1>Output</h1>
<pre><code>head(temp_single)
</code></pre>
<pre><code>##   X.Intercept.         act       diff       Fbeta  logistic sim_data ID
## 1            1 -0.02913722 -0.4430546 -0.04748497 0.4881310        1  1
## 2            1  0.66199364  2.1443743  1.07430910 0.7454155        1  2
## 3            1  1.44621026 -1.1909231  0.46582819 0.6143959        0  3
## 4            1 -0.26011629  3.4395304  1.00180096 0.7314125        0  4
## 5            1 -0.09984213  0.8485436  0.30464201 0.5755769        1  5
## 6            1 -2.72704127  3.3246515 -0.26612517 0.4338586        0  6
</code></pre>
</section>
<section>
<h1>Simple Power Analysis</h1>
<ul>
<li>Suppose we wish to use the same generating model for a power analysis</li>
</ul>
<pre><code>pow_param &lt;- c('(Intercept)', 'act', 'diff')
alpha &lt;- .01
pow_dist &lt;- &quot;z&quot;
pow_tail &lt;- 2
replicates &lt;- 100

power_out &lt;- sim_pow_glm(fixed = fixed, fixed_param = fixed_param, 
                         cov_param = cov_param, 
                         n = n, data_str = &quot;single&quot;, 
                         pow_param = pow_param, alpha = alpha,
                         pow_dist = pow_dist, pow_tail = pow_tail, 
                         replicates = replicates)
</code></pre>
</section>
<section>
<h1>Output</h1>
<pre><code>power_out
</code></pre>
<pre><code>## # A tibble: 3 × 6
##           var avg_test_stat sd_test_stat power num_reject num_repl
##        &lt;fctr>         &lt;dbl>        &lt;dbl> &lt;dbl>      &lt;dbl>    &lt;dbl>
## 1 (Intercept)      0.878713    0.6709319  0.01          1      100
## 2         act      2.342617    0.5777646  0.34         34      100
## 3        diff      2.609432    0.5506204  0.56         56      100
</code></pre>
</section>
<section>
<h1>Varying Arguments</h1>
<ul>
<li>Now suppose we wish to vary the following arguments:
<ul>
<li>Vary n – 50 vs 150</li>
<li>vary effect size on diff – .3 vs .45</li>
</ul>
</li>
</ul>
<pre><code>terms_vary &lt;- list(n = c(50, 150),
                   fixed_param = list(c(0.1, 0.5, 0.3), 
                                      c(0.1, 0.5, 0.45)))

power_out &lt;- sim_pow_glm(fixed = fixed, fixed_param = fixed_param, 
                         cov_param = cov_param, 
                         n = n, data_str = &quot;single&quot;, 
                         pow_param = pow_param, alpha = alpha,
                         pow_dist = pow_dist, pow_tail = pow_tail, 
                         replicates = replicates, 
                         terms_vary = terms_vary)
</code></pre>
</section>
<section>
<h1>Output</h1>
<pre><code>power_out
</code></pre>
<pre><code>## Source: local data frame [12 x 8]
## Groups: var, n [?]
## 
##            var     n  fixed_param avg_test_stat sd_test_stat power
##         &lt;fctr> &lt;dbl>       &lt;fctr>         &lt;dbl>        &lt;dbl> &lt;dbl>
## 1  (Intercept)    50  0.1,0.5,0.3     0.7778328    0.5863240  0.00
## 2  (Intercept)    50 0.1,0.5,0.45     0.8364212    0.6377631  0.01
## 3  (Intercept)   150  0.1,0.5,0.3     0.8629973    0.5814426  0.00
## 4  (Intercept)   150 0.1,0.5,0.45     0.9183353    0.6879182  0.01
## 5          act    50  0.1,0.5,0.3     2.4246997    0.6222346  0.44
## 6          act    50 0.1,0.5,0.45     2.2247451    0.6688308  0.34
## 7          act   150  0.1,0.5,0.3     4.3196568    0.6233962  0.99
## 8          act   150 0.1,0.5,0.45     3.9515646    0.6332452  0.97
## 9         diff    50  0.1,0.5,0.3     2.7887204    0.4892985  0.73
## 10        diff    50 0.1,0.5,0.45     3.0747886    0.3988745  0.89
## 11        diff   150  0.1,0.5,0.3     4.7892881    0.5025082  1.00
## 12        diff   150 0.1,0.5,0.45     5.6060130    0.2823105  1.00
## # ... with 2 more variables: num_reject &lt;dbl>, num_repl &lt;dbl>
</code></pre>
</section>
<section>
<h1>Move to Mixed Models</h1>
<ul>
<li>It is simple to move from single level to multilevel or mixed models.</li>
</ul>
<pre><code>fixed &lt;- ~1 + time + diff + act + time:act
random &lt;- ~1 + time
fixed_param &lt;- c(0, 0.2, 0.1, 0.3, 0.05)
random_param &lt;- list(random_var = c(3, 2), rand_gen = &quot;rnorm&quot;)
cov_param &lt;- list(dist_fun = c('rnorm', 'rnorm'),
                  var_type = c(&quot;lvl1&quot;, &quot;lvl2&quot;),
                  opts = list(list(mean = 0, sd = 3),
                              list(mean = 0, sd = 2)))
n &lt;- 50
p &lt;- 6
data_str &lt;- &quot;long&quot;

temp_long &lt;- sim_glm(fixed = fixed, random = random, fixed_param = fixed_param,
                     random_param = random_param, cov_param = cov_param,
                     n = n, p = p, k = NULL, data_str = data_str)
</code></pre>
</section>
<section>
<h1>Output</h1>
<pre><code>head(temp_long)
</code></pre>
<pre><code>##   X.Intercept. time        diff        act   time.act        b0        b1
## 1            1    0 -6.76572749 -0.3932853  0.0000000 -1.947485 -2.295427
## 2            1    1  0.15530420 -0.3932853 -0.3932853 -1.947485 -2.295427
## 3            1    2  0.07605058 -0.3932853 -0.7865707 -1.947485 -2.295427
## 4            1    3 -1.11192544 -0.3932853 -1.1798560 -1.947485 -2.295427
## 5            1    4 -4.17141062 -0.3932853 -1.5731413 -1.947485 -2.295427
## 6            1    5  4.77024867 -0.3932853 -1.9664267 -1.947485 -2.295427
##         Fbeta    randEff   logistic         prob sim_data withinID clustID
## 1 -0.79455835  -1.947485  -2.742044 6.053757e-02        0        1       1
## 2  0.07788055  -4.242913  -4.165032 1.529175e-02        0        2       1
## 3  0.25029093  -6.538340  -6.288049 1.854935e-03        0        3       1
## 4  0.31182906  -8.833767  -8.521938 1.990136e-04        0        4       1
## 5  0.18621627 -11.129195 -10.942978 1.768142e-05        0        5       1
## 6  1.26071793 -13.424622 -12.163904 5.215325e-06        0        6       1
</code></pre>
</section>
<section>
<h1>Doing Power</h1>
<ul>
<li>Power is also easily extended.</li>
</ul>
<pre><code>pow_param &lt;- c('time', 'diff', 'act')
alpha &lt;- .01
pow_dist &lt;- &quot;z&quot;
pow_tail &lt;- 2
replicates &lt;- 20

power_out &lt;- sim_pow_glm(fixed = fixed, random = random, 
                     fixed_param = fixed_param, 
                     random_param = random_param, cov_param = cov_param, 
                     k = NULL, n = n, p = p,
                     data_str = data_str, unbal = FALSE, pow_param = pow_param, 
                     alpha = alpha, pow_dist = pow_dist, pow_tail = pow_tail,
                     replicates = replicates)
</code></pre>
</section>
<section>
<h1>Output</h1>
<pre><code>power_out
</code></pre>
<pre><code>## # A tibble: 3 × 6
##      var avg_test_stat sd_test_stat power num_reject num_repl
##   &lt;fctr>         &lt;dbl>        &lt;dbl> &lt;dbl>      &lt;dbl>    &lt;dbl>
## 1    act      12.06197     46.70227  0.20          4       20
## 2   diff      11.89673     45.13827  0.25          5       20
## 3   time      18.78877     79.36869  0.05          1       20
</code></pre>
</section>
<section>
<h1>Vary Arguments</h1>
<ul>
<li>Perhaps our effect size estimate is conservative.</li>
</ul>
<pre><code>terms_vary &lt;- list(fixed_param = list(c(0, 0.2, 0.1, 0.3, 0.05), 
                                      c(0, 0.2, 0.3, 0.3, 0.05)))

power_out &lt;- sim_pow_glm(fixed = fixed, random = random, 
                     fixed_param = fixed_param, 
                     random_param = random_param, cov_param = cov_param, 
                     k = NULL, n = n, p = p,
                     data_str = data_str, unbal = FALSE, pow_param = pow_param, 
                     alpha = alpha, pow_dist = pow_dist, pow_tail = pow_tail,
                     replicates = replicates, 
                     terms_vary = terms_vary)
</code></pre>
</section>
<section>
<h1>Output</h1>
<pre><code>power_out
</code></pre>
<pre><code>## Source: local data frame [6 x 7]
## Groups: var [?]
## 
##      var        fixed_param avg_test_stat sd_test_stat power num_reject
##   &lt;fctr>             &lt;fctr>         &lt;dbl>        &lt;dbl> &lt;dbl>      &lt;dbl>
## 1    act 0,0.2,0.1,0.3,0.05     1.1914255    0.8114762  0.10          2
## 2    act 0,0.2,0.3,0.3,0.05    22.9059014   96.3531136  0.15          3
## 3   diff 0,0.2,0.1,0.3,0.05     1.3071639    0.8681348  0.05          1
## 4   diff 0,0.2,0.3,0.3,0.05    17.4774138   62.2814403  0.95         19
## 5   time 0,0.2,0.1,0.3,0.05     0.9281452    0.7670600  0.05          1
## 6   time 0,0.2,0.3,0.3,0.05    12.1678311   49.9607401  0.05          1
## # ... with 1 more variables: num_repl &lt;dbl>
</code></pre>
</section>
<section>
<h1>Shiny App</h1>
<ul>
<li>Note: This app currently looks nice, but is utterly broken!</li>
</ul>
<pre><code>shiny::runGitHub('simglm', username = 'lebebr01', subdir = 'inst/shiny_examples/demo')
</code></pre>
<p>or</p>
<pre><code>devtools::install_github('lebebr01/simglm')
library(simglm)
run_shiny()
</code></pre>
<ul>
<li>Must have following packages installed: <code>simglm, shiny, shinydashboard, ggplot2, lme4, DT</code>.<br>

</li>
</ul>
<section>
<h1><code>simglm</code> timeline</h1>
<ul>
<li>Aim to have this package submitted to CRAN by the end of March.</li>
<li>Fix Shiny application.</li>
<li>For now look for the package on GitHub <a target="_blank" rel="nofollow" href="http://github.com/lebebr01/simglm">http://github.com/lebebr01/simglm</a><br>

</li>
</ul>
<section>
<h1>Questions?</h1>
<ul>
<li>Twitter: @blebeau11</li>
<li>Website: <a target="_blank" rel="nofollow" href="http://educate-r.org/">http://educate-r.org</a></li>
<li>Slides: <a target="_blank" rel="nofollow" href="http://educate-r.org/2017/02/24/csp2017.html">http://educate-r.org/2017/02/24/csp2017.html</a></li>
<li>GitHub: <a target="_blank" rel="nofollow" href="http://github.com/lebebr01/simglm">http://github.com/lebebr01/simglm</a><br>

</li>
</ul>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://educate-r.org//2017/02/24/csp2017.html"> Educate-R - R</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/PUKFI5cn2tU"></a></div></section></section></section></section></section></section></section></section></section></section></section>