<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/bootstrap-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/29/bootstrap-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>Now that we’ve seen <a rel="nofollow" href="https://wp.me/p27liZ-TD" target="_blank"><strong>MCHT</strong> basics</a>, <a rel="nofollow" href="https://wp.me/p27liZ-TG" target="_blank">how to make <code>MCHTest()</code> objects self-contained</a>, and <a rel="nofollow" href="https://wp.me/p27liZ-TQ" target="_blank">maximized Monte Carlo (MMC) testing with <strong>MCHT</strong></a>, let’s now talk about bootstrap testing. Not much is different when we’re doing bootstrap testing; the main difference is that the replicates used to generate test statistics depend on the data we feed to the test, and thus are not completely independent of it. You can read more about bootstrap testing in [1].</p>
<p><span></span></p>
<h2>Bootstrap Hypothesis Testing</h2>
<p>Let <img alt="S" title="S" src="https://s0.wp.com/latex.php?latex=S&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> represent our test statistic. For bootstrap hypothesis testing, we will construct <img alt="N" title="N" src="https://s0.wp.com/latex.php?latex=N&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> test statistics from data generated using our sample. Call these test statistics <img alt="S^{*}_1, ldots, S^{*}_N" title="S^{*}_1, ldots, S^{*}_N" src="https://s0.wp.com/latex.php?latex=S%5E%7B%2A%7D_1%2C+%5Cldots%2C+S%5E%7B%2A%7D_N&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. These statistics are generated in such a way that we know that the null hypothesis holds for them. Suppose for the sake of demonstration that large values of <img alt="S" title="S" src="https://s0.wp.com/latex.php?latex=S&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> constitute evidence against the null hypothesis. Then 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 for the bootstrap hypothesis test is</p>
<p><img alt="hat{p} = frac{1}{N} sum_{j = 1}^{N} I_{{S^{*}_j &gt; S}}" title="hat{p} = frac{1}{N} sum_{j = 1}^{N} I_{{S^{*}_j &gt; S}}" src="https://s0.wp.com/latex.php?latex=%5Chat%7Bp%7D+%3D+%5Cfrac%7B1%7D%7BN%7D+%5Csum_%7Bj+%3D+1%7D%5E%7BN%7D+I_%7B%5C%7BS%5E%7B%2A%7D_j+%3E+S%5C%7D%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>Here, <img alt="I" title="I" src="https://s0.wp.com/latex.php?latex=I&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Indicator_function" target="_blank">indicator function</a>.</p>
<p>There are many ways to generate the data used to compute <img alt="S^{*}_j" title="S^{*}_j" src="https://s0.wp.com/latex.php?latex=S%5E%7B%2A%7D_j&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. There’s the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Bootstrapping_(statistics)#Parametric_bootstrap" target="_blank">parametric bootstrap</a>, where the data is used to estimate the parameters of a distribution, then those parameters are plugged into that distribution and then the distribution is used to generate new samples. There’s also the nonparametric bootstrap that doesn’t make such strong assumptions about the data, perhaps sampling from the data itself to generate new samples. Either of these methods can be used in bootstrap testing, and <code>MCHTest()</code> supports both.</p>
<p>Unlike Monte Carlo tests, bootstrap tests cannot claim to be exact tests for any sample size; they’re better for larger sample sizes. That said, they often work well even in small sample sizes and thus are still a good alternative to inference based on asymptotic results. They also could serve as an alternative approach to the nuisance parameter problem, as MMC often has weak power.</p>
<h2>Bootstrap Testing in MCHT</h2>
<p>In <strong>MCHT</strong>, there is little difference between bootstrap testing and Monte Carlo testing. Bootstrap tests need the original dataset to generate replicates; Monte Carlo tests do not. So the difference here is that the function passed to <code>rand_gen</code> needs to accept a parameter <code>x</code> rather than <code>n</code>, with <code>x</code> representing the original dataset, like that passed to <code>test_stat</code>.</p>
<p>That’s the only difference. All else is the same.</p>
<h2>Nonparametric Bootstrap Example</h2>
<p>Suppose we wish to test for the location of the mean. Our nonparametric bootstrap procedure is as follows:</p>
<ol>
<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"> samples of data from the demeaned dataset.</li>
<li>Suppose our mean under the null hypothesis is <img alt="mu_0" title="mu_0" src="https://s0.wp.com/latex.php?latex=%5Cmu_0&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. Add this mean to each generated dataset and compute the <img alt="t" title="t" src="https://s0.wp.com/latex.php?latex=t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> statistic for each of those datasets; these will be the simulated test statistics <img alt="S^{*}_j" title="S^{*}_j" src="https://s0.wp.com/latex.php?latex=S%5E%7B%2A%7D_j&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">.</li>
<li>Compute the test statistic on the main data and 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 a <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.</li>
</ol>
<p>The code below implements this procedure.</p>
<pre>
library(MCHT)
library(doParallel)

registerDoParallel(detectCores())

ts &lt;- function(x, mu = 0) {
  sqrt(length(x)) * (mean(x) - mu)/sd(x)
}

rg &lt;- function(x) {
  x_demeaned &lt;- x - mean(x)
  sample(x_demeaned, replace = TRUE, size = length(x))
}

sg &lt;- function(x, mu = 0) {
  x &lt;- x + mu
  test_stat(x, mu = mu)  # Will be localizing
}

b.t.test &lt;- MCHTest(ts, sg, rg, seed = 123, N = 1000, lock_alternative = FALSE,
                    test_params = &quot;mu&quot;, localize_functions = TRUE)

dat &lt;- c(2.3, 1.1, 8.1, -0.2, -0.8, 4.7, -1.9)

b.t.test(dat, alternative = &quot;two.sided&quot;, mu = 1)
</pre>
<pre>
## 
## 	Monte Carlo Test
## 
## data:  dat
## S = 0.68164, p-value = 0.432
## alternative hypothesis: true mu is not equal to 1
</pre>
<pre>
b.t.test(dat, alternative = &quot;less&quot;, mu = 7)
</pre>
<pre>
## 
## 	Monte Carlo Test
## 
## data:  dat
## S = -3.8626, p-value = 0.025
## alternative hypothesis: true mu is less than 7
</pre>
<h2>Parametric Bootstrap Test</h2>
<p>The parametric bootstrap test assumes that the observed data was generated using a specific distribution, such as the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Normal_distribution" target="_blank">Gaussian distribution</a>. All that’s missing, in essence, is the parameters of that distribution. The procedure thust starts by estimating all nuisance parameters of the assumed distribution using the data. Then the first step of the process mentioned above (which admittedly was specific to a test for the mean but still strongly resembles the general process) is replaced with simulating data from the assumed distribution using any parameters assumed under the null hypothesis and the estimated values of any nuisance parameters. The other two steps of the above process are unchanged.</p>
<p>We can use the parametric bootstrap to test for goodness of fit with the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test#Kolmogorov_distribution" target="_blank">Kolmogorov-Smirnov test</a>. Without going into much detail, suppose that <img alt="F" title="F" src="https://s0.wp.com/latex.php?latex=F&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> represents a distribution that is known except maybe for the values of its parameters. Assume that <img alt="X_1, ldots, X_n" title="X_1, ldots, X_n" src="https://s0.wp.com/latex.php?latex=X_1%2C+%5Cldots%2C+X_n&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is an independently and identically distributed dataset, and we have observed values <img alt="x_1, ldots, x_n" title="x_1, ldots, x_n" src="https://s0.wp.com/latex.php?latex=x_1%2C+%5Cldots%2C+x_n&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. We wish to use the dataset to decide between the hypotheses</p>
<p><img alt="H_0: X_1 sim F" title="H_0: X_1 sim F" src="https://s0.wp.com/latex.php?latex=H_0%3A+X_1+%5Csim+F&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p><img alt="H_A: H_0 text{ is false}" title="H_A: H_0 text{ is false}" src="https://s0.wp.com/latex.php?latex=H_A%3A+H_0+%5Ctext%7B+is+false%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>That is, we want to test whether our data was drawn from the distribution <img alt="F" title="F" src="https://s0.wp.com/latex.php?latex=F&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> or whether it was drawn from a different distribution. This is what the Kolmogorov-Smirnov test checks.</p>
<p>R implements this test in <code>ks.test()</code> but that function does not allow for any nuisance parameters. It will only check for an exact match between distributions. This is often not what we want; we want to check whether out data was drawn from any member of the family of distributions <img alt="F" title="F" src="https://s0.wp.com/latex.php?latex=F&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, not a particular member with a particular combination of parameters. It’s tempting to plug in the estimated values of these parameters, but then 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 needs to be computed differently, not in the way that is prescribed by <code>ks.test()</code>. Thus we will need to approach the test differently.</p>
<p>Since the distribution of the data is known under the null hypothesis, this is a good situation to use a bootstrap test. We’ll use <a rel="nofollow" href="https://en.wikipedia.org/wiki/Maximum_likelihood_estimation" target="_blank">maximum likelihood estimation</a> to estimate the values of the missing parameters, as implemented by <strong>fitdistrplus</strong> (see [2]). Then we generate samples from this distribution using the estimated parameter values and use those samples to generate simulated test statistic values that follow the distribution prescribed by the null hypothesis.</p>
<p>Suppose we wished to test whether the data was drawn from a Weibull distribution. The result is the following test.</p>
<pre>
library(fitdistrplus)

ts &lt;- function(x) {
  param &lt;- coef(fitdist(x, &quot;weibull&quot;))
  shape &lt;- param[['shape']]; scale &lt;- param[['scale']]
  ks.test(x, pweibull, shape = shape, scale = scale,
          alternative = &quot;two.sided&quot;)$statistic[[1]]
}

rg &lt;- function(x) {
  n &lt;- length(x)
  param &lt;- coef(fitdist(x, &quot;weibull&quot;))
  shape &lt;- param[['shape']]; scale &lt;- param[['scale']]
  rweibull(n, shape = shape, scale = scale)
}

b.ks.test &lt;- MCHTest(test_stat = ts, stat_gen = ts, rand_gen = rg,
                     seed = 123, N = 1000)

b.ks.test(rweibull(1000, 2, 2))
</pre>
<pre>
## 
## 	Monte Carlo Test
## 
## data:  rweibull(1000, 2, 2)
## S = 0.021907, p-value = 0.275
</pre>
<pre>
b.ks.test(rbeta(1000, 2, 2))
</pre>
<pre>
## 
## 	Monte Carlo Test
## 
## data:  rbeta(1000, 2, 2)
## S = 0.047165, p-value &lt; 2.2e-16
</pre>
<h2>Conclusion</h2>
<p>Given the choice between a MMC test and a bootstrap test, which should you prefer? If you’re concerned about speed and power, go with the bootstrap test. If you’re concerned about precision and getting an “exact” test that’s at least conservative, then go with a MMC test. I think most of the time, though, the bootstrap test will be good enough, even with small samples, but that’s mostly a hunch.</p>
<p>Next week we will see how we can go beyond one-sample or univariate tests to multi-sample or multivariate tests. <a rel="nofollow" href="https://wp.me/p27liZ-U0" target="_blank">See the next blog post.</a></p>
<h2>References</h2>
<ol>
<li>J. G. MacKinnon, <em>Bootstrap hypothesis testing</em> in <em>Handbook of computational econometrics</em> (2009) pp. 183-213</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>
</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/29/bootstrap-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/4LlPFg1WoWw" height="1"></a></div>