<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/maximized-monte-carlo-testing-with-mcht/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://ntguardian.wordpress.com/2018/10/22/maximized-monte-carlo-testing-mcht/"> R – Curtis Miller's Personal Website</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Introduction</h2>
<p>I <a rel="nofollow" href="https://wp.me/p27liZ-TD" target="_blank">introduced <strong>MCHT</strong></a> two weeks ago and presented it as a package for Monte Carlo and boostrap hypothesis testing. Last week, I delved into important technical details and showed how to <a rel="nofollow" href="https://wp.me/p27liZ-TG" target="_blank">make self-contained <code>MCHTest</code> objects</a> that don’t suffer side effects from changes in the global namespace. In this article I show how to perform maximized Monte Carlo hypothesis testing using <strong>MCHT</strong>, as described in [1].</p>
<p><span></span></p>
<h2>Maximized Monte Carlo Hypothesis Testing</h2>
<p>The usual procedure for Monte Carlo hypothesis testing is:</p>
<ol>
<li>Compute a test statistic for the data on which you wish to test a hypothesis</li>
<li>Generate <img alt="N" title="N" src="https://s0.wp.com/latex.php?latex=N&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> random datasets like the one of interest but with the data generating process (DGP) being the one prescribed by the null hypothesis, and compute the test statistic on each of these datasets</li>
<li>Use the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Empirical_distribution_function" target="_blank">empirical distribution function</a> of the simulated test statistics to compute the <img alt="p" title="p" src="https://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-value of the test</li>
</ol>
<p>Monte Carlo tests often make strong distributional assumptions, such as what distribution generated the dataset being tested, but when those assumptions hold, they are exact tests (see [2]). They are not as powerful as if we had the exact distribution of the test statistic under those assumptions but the power increases with <img alt="N" title="N" src="https://s0.wp.com/latex.php?latex=N&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> (see [3]) and given the power of modern computers getting a large <img alt="N" title="N" src="https://s0.wp.com/latex.php?latex=N&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is usually not a problem. Thus Monte Carlo tests are attractive in small sample situations where we do not want to rely on an asymptotic distribution for inference.</p>
<p>However, the procedure outlined above does not allow for nuisance parameters (parameters that are not the subject of interest but whose values are needed in order to conduct inference). In the introductory blog post, while one may view the population standard deviation <img alt="sigma" title="sigma" src="https://s0.wp.com/latex.php?latex=%5Csigma&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> as a nuisance parameter, the test statistic <img alt="frac{bar{X} - mu_0}{frac{S}{sqrt{n}}}" title="frac{bar{X} - mu_0}{frac{S}{sqrt{n}}}" src="https://s0.wp.com/latex.php?latex=%5Cfrac%7B%5Cbar%7BX%7D+-+%5Cmu_0%7D%7B%5Cfrac%7BS%7D%7B%5Csqrt%7Bn%7D%7D%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> does not depend on it when the data follows a <a rel="nofollow" href="https://en.wikipedia.org/wiki/Normal_distribution" target="_blank">Gaussian distribution</a> so there was no need to worry about it. In the case when we switched to data following the exponential distribution, it was still not a problem since its value was specified under the null hypothesis (<img alt="sigma = mu_0" title="sigma = mu_0" src="https://s0.wp.com/latex.php?latex=%5Csigma+%3D+%5Cmu_0&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">). Thus it was no longer a nuisance parameter.</p>
<p>That said, nuisance parameters can still appear when we need to perform inference. Suppose, for example, that our data follows a <a rel="nofollow" href="https://en.wikipedia.org/wiki/Weibull_distribution" target="_blank">Weibull distribution</a>, denoted by <img alt="X_i sim text{WEI}(k, lambda)" title="X_i sim text{WEI}(k, lambda)" src="https://s0.wp.com/latex.php?latex=X_i+%5Csim+%5Ctext%7BWEI%7D%28k%2C+%5Clambda%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, with <img alt="k" title="k" src="https://s0.wp.com/latex.php?latex=k&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> being the shape parameter and <img alt="lambda" title="lambda" src="https://s0.wp.com/latex.php?latex=%5Clambda&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> the scale parameter. We want to test the set of hypotheses:</p>
<p><img alt="H_0: lambda = lambda_0" title="H_0: lambda = lambda_0" src="https://s0.wp.com/latex.php?latex=H_0%3A+%5Clambda+%3D+%5Clambda_0&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p><img alt="H_A: lambda neq lambda_0" title="H_A: lambda neq lambda_0" src="https://s0.wp.com/latex.php?latex=H_A%3A+%5Clambda+%5Cneq+%5Clambda_0&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>We can use the generalized likelihood ratio test to form a test statistic (which I won’t repeat hear but does appear in the code below). While <a rel="nofollow" href="https://en.wikipedia.org/wiki/Wilks%27_theorem" target="_blank">Wilks’ theorem</a> tells us about the asymptotic distribution of the test statistic, it says nothing about the exact distribution of the test statistic at a particular sample size, and it’s not given that the test statistic is <a rel="nofollow" href="https://en.wikipedia.org/wiki/Pivotal_quantity" target="_blank">pivotal</a> and thus independent of the value of nuisance parameters under the null hypothesis (the nuisance parameter, in this case, being the shape parameter <img alt="k" title="k" src="https://s0.wp.com/latex.php?latex=k&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">).</p>
<p>What then can we do? A bootstrap test would estimate the value of the nuisance parameter under the null hypothesis and use that estimate as the actual value when generating new, simulated test statistics. Bootstrap tests, however, are not exact tests (see [2]) and we’ve decided that we want a test with stronger guarantees.</p>
<p>[1] introduced the maximized Monte Carlo (MMC) test, which proceeds as described below:</p>
<ol>
<li>Compute the test statistic from the data.</li>
<li>Generate <img alt="N" title="N" src="https://s0.wp.com/latex.php?latex=N&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> collections of random numbers, such as uniformly distributed random numbers, and use those random numbers for generating random copies of test statistics that depend on the values of nuisance parameters (notice that the random numbers are effectively <em>not</em> discarded)</li>
<li>Use an optimization procedure ([1] suggested simulated annealing) to pick values for the nuisance parameters such that the <img alt="p" title="p" src="https://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-value is maximized; the maximally chosen <img alt="p" title="p" src="https://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-value is the <img alt="p" title="p" src="https://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-value of the test</li>
</ol>
<p>[1] showed that this procedure yields <img alt="p" title="p" src="https://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-values that, while not as precise as if we knew the values of the nuisance parameters that produced the data, are at least <em>conservative</em>, in the sense that they’re no smaller than they should be (thus biasing our conclusions in favor of the null hypothesis). This is the best we can hope for in this context.</p>
<p>MMC is intuitive and compelling, and the theoretical guarantee gives us confidence in our conclusions, but it’s not a panacea. First, the optimization procedure is costly in work and time. Second (and, in my opinion, the biggest problem), the procedure may be <em>too</em> conservative. There’s a strong chance that the procedure will find <em>some</em> values for nuisance parameters that yield a large <img alt="p" title="p" src="https://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-value, perhaps a combination not at all resembling the actual values of the nuisance parameters that produced the data. In short, MMC can be severely lacking in <a rel="nofollow" href="https://en.wikipedia.org/wiki/Power_(statistics)" target="_blank">power</a>. When it does reject the null hypothesis, it’s compelling, but otherwise it’s not convincing that the alternative hypothesis is false.</p>
<h2>MMC in MCHT</h2>
<p>Creating an implementation of MMC in R was my original goal in developing <strong>MCHT</strong>, and all that needs to be done to perform MMC is pass a value the <code>nuisance_params</code> parameter and an appropriate list to <code>optim_control</code>.</p>
<p>Let’s take the hypothesis test mentioned above and prepare to implement it using <strong>MCHT</strong>. I will be using <strong>fitdistrplus</strong> for maximum likelihood estimation, as required by the test statistic (see <a rel="nofollow" href="https://www.jstatsoft.org/article/view/v064i04" target="_blank">[4]</a>).</p>
<pre>
library(MCHT)
</pre>
<pre>
## .------..------..------..------.
## |M.--. ||C.--. ||H.--. ||T.--. |
## | (\/) || :/\: || :/\: || :/\: |
## | :\/: || :\/: || (__) || (__) |
## | '--'M|| '--'C|| '--'H|| '--'T|
## `------'`------'`------'`------' v. 0.1.0
## Type citation(&quot;MCHT&quot;) for citing this R package in publications
</pre>
<pre>
library(fitdistrplus)
</pre>
<pre>
registerDoParallel(detectCores())

# To be passed to test_stat
ts &lt;- function(x, scale = 1) {
  fit_null &lt;- coef(fitdist(x, &quot;weibull&quot;, fix.arg = list(&quot;scale&quot; = scale)))
  kt &lt;- fit_null[[&quot;shape&quot;]]
  l0 &lt;- scale
  fit_all &lt;- coef(fitdist(x, &quot;weibull&quot;))
  kh &lt;- fit_all[[&quot;shape&quot;]]
  lh &lt;- fit_all[[&quot;scale&quot;]]
  n &lt;- length(x)

  # Test statistic, based on the negative-log-likelihood ratio
  suppressWarnings(n * ((kt - 1) * log(l0) - (kh - 1) * log(lh) -
      log(kt/kh) - log(lh/l0)) - (kt - kh) * sum(log(x)) + l0^(-kt) *
      sum(x^kt) - lh^(-kh) * sum(x^kh))
}

# To be passed to stat_gen; localize_functions will be TRUE
sg &lt;- function(x, scale = 1, shape = 1) {
  x &lt;- qweibull(x, shape = shape, scale = scale)
  test_stat(x, scale = scale)
}
</pre>
<p>The <code>MCHTest()</code> parameter <code>nuisance_params</code> accepts a character vector giving the names of nuisance parameters the distribution of the test statistic may depend upon, and those names must be among the arguments of the function passed to <code>stat_gen</code>; that function is expected to know how to handle those parameters. In this case, <code>rand_gen</code> will not be specified since by default it gives uniformly distributed random variables. It’s a well-known fact in probability that the inverse of the CDF of a random variable (which are the <code>q</code> functions in R) applied to a uniformly distributed random variable yields a random variable that follows the distribution prescribed by the CDF. Hence the use of <code>qweibull()</code> above, which is being applied to datasets of uniformly distributed random variables that are effectively fixed when <code>stat_gen</code> will be called. Then the test statistic computed will be computed from data that follows the scale parameter <img alt="lambda" title="lambda" src="https://s0.wp.com/latex.php?latex=%5Clambda&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> prescribed by the null hypothesis but for some set value of <img alt="k" title="k" src="https://s0.wp.com/latex.php?latex=k&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, the shape parameter.</p>
<p>The <code>MCHTest</code> object will then perform simulated annealing to choose the value of the nuisance parameter that maximizes the <img alt="p" title="p" src="https://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-value under the null hypothesis for the given dataset. Simulated annealing is implemented in the <code>GenSA()</code> function provided by the <strong>GenSA</strong> package (see <a rel="nofollow" href="https://journal.r-project.org/archive/2013/RJ-2013-002/index.html" target="_blank">[5]</a>). <code>GenSA()</code> needs a description of what set of parameter values over which to optimize and there is no general method for choosing this, so <code>MCHTest()</code> will require that a list be passed to <code>optim_control</code> that effectively contains the parameters that will be passed to <code>GenSA()</code>. At minimum, this list must contain an <code>upper</code> and a <code>lower</code> element, each of which are numeric vectors with names exactly like the character vector passed to <code>nuisance_params</code>; these vectors specify the space <code>GenSA()</code> will search to find the optima. Other elements can be passed to control <code>GenSA()</code>, and I highly recommend reading the function’s documentation for more details.</p>
<p>There’s an additional parameter of <code>MCHTest()</code>, <code>threshold_pval</code>, that matters to the optimization. <code>GenSA()</code> will take many steps to make sure it reaches a good optimal value, perhaps taking too long. The authors of <strong>GenSA</strong> recommend specifying an additional terminating condition to speed up the process. <code>threshold_pval</code> will alter the <code>threshold.stop</code> parameter in the <code>control</code> list of <code>optim_control</code> so that the algorithm terminates when the estimated <img alt="p" title="p" src="https://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-value crosses <code>threshold_pval</code>‘s value. Effectively, this means that we know that whatever the true <img alt="p" title="p" src="https://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-value of the test is, it’s larger than that threshold, and if the threshold is chosen appropriately, we know that we should not reject the null hypothesis based on the results of this test.</p>
<p>While giving <code>threshold_pval</code> a number less than 1 can help terminate the algorithm in the case of not rejecting the null hypothesis, the algorithm can still run for a long time if the test will eventually return a statistically significant result. For this reason, I recommend that <code>optim_list</code> contain a list called <code>control</code> and that the list should have a <code>max.time</code> element telling the algorithm the maximum running time (in seconds) it should have.</p>
<p>With all this in mind, we create the <code>MCHTest</code> object below:</p>
<pre>
mc.wei.scale.test &lt;- MCHTest(ts, sg, N = 1000, seed = 123,
                             test_params = &quot;scale&quot;, nuisance_params = &quot;shape&quot;,
                             optim_control = list(&quot;lower&quot; = c(&quot;shape&quot; = 0),
                                                  &quot;upper&quot; = c(&quot;shape&quot; = 100),
                                                  &quot;control&quot; = list(
                                                    &quot;max.time&quot; = 10
                                                  )), threshold_pval = .2,
                             localize_functions = TRUE)

mc.wei.scale.test(rweibull(10, 2, 2), scale = 4)
</pre>
<pre>
## 
## 	Monte Carlo Test
## 
## data:  rweibull(10, 2, 2)
## S = 7.2983, p-value &lt; 2.2e-16
</pre>
<h2>Conclusion</h2>
<p>The MMC procedure is intersting and I don’t think any package implements it to the level I have in <strong>MCHT</strong>. The power of the procedure itself concerns me, though. Fortunately, the package also supports bootstrap testing, which I will discuss <a rel="nofollow" href="https://wp.me/p27liZ-TU" target="_blank">next week</a>.</p>
<h2>References</h2>
<ol>
<li>J-M Dufour, <em>Monte Carlo tests with nuisance parameters: A general approach to finite-sample inference and nonstandard asymptotics</em>, Journal of Econometrics, vol. 133 no. 2 (2006) pp. 443-477</li>
<li>J. G. MacKinnon, <em>Bootstrap hypothesis testing</em> in <em>Handbook of computational econometrics</em> (2009) pp. 183-213</li>
<li>A. C. A. Hope, <em>A simplified Monte Carlo test procedure</em>, JRSSB, vol. 30 (1968) pp. 582-598</li>
<li>M. L. Delignette-Muller and C. Dulag, <em>fitdistrplus: an R package for fitting distributions</em>, J. Stat. Soft., vol. 64 no. 4 (2015)</li>
<li>Y. Xiang et. al., <em>Generalized simulated annealing for global optimization: the GenSA package</em>, R Journal, vol. 5 no. 1 (2013) pp. 13-28</li>
</ol>
<hr>
<p><a rel="nofollow" href="https://www.packtpub.com/" target="_blank">Packt Publishing</a> published a book for me entitled <a rel="nofollow" href="https://ntguardian.wordpress.com/books/hands-on-data-analysis-with-numpy-and-pandas/" target="_blank"><em>Hands-On Data Analysis with NumPy and Pandas</em></a>, a book based on my video course <a rel="nofollow" href="https://ntguardian.wordpress.com/video-courses/unpacking-numpy-and-pandas/" target="_blank"><em>Unpacking NumPy and Pandas</em></a>. This book covers the basics of setting up a Python environment for data analysis with Anaconda, using Jupyter notebooks, and using NumPy and pandas. If you are starting out using Python for data analysis or know someone who is, please consider <a rel="nofollow" href="https://www.packtpub.com/big-data-and-business-intelligence/hands-data-analysis-numpy-and-pandas" target="_blank">buying my book</a> or at least spreading the word about it. You can buy the book directly or purchase a <a rel="nofollow" href="https://mapt.packtpub.com/?gclid=Cj0KCQjwuuHdBRCvARIsAELQRQGrXDoR4ox216bUlSD6mTQYibN24CA5bYF7QikiQQ6j3cGm63HRXsEaAjJjEALw_wcB" target="_blank">subscription to Mapt</a> and read it there.</p>
<p>If you like my blog and would like to support it, spread the word (if not get a copy yourself)!</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://ntguardian.wordpress.com/2018/10/22/maximized-monte-carlo-testing-mcht/"> R – Curtis Miller's Personal Website</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 rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/eImQD2q5fHU" height="1"></a></div>