<!DOCTYPE html>
<!-- Generated by pkgdown: do not edit by hand --><html lang="en">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>hmclearn:  Logistic Regression Example • hmclearn</title>
<!-- jquery --><script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/3.4.1/jquery.min.js" integrity="sha256-CSXorXvZcTkaix6Yvo6HppcZGetbYMGWSFlBw8HfCJo=" crossorigin="anonymous"></script><!-- Bootstrap --><link href="https://cdnjs.cloudflare.com/ajax/libs/bootswatch/3.4.0/sandstone/bootstrap.min.css" rel="stylesheet" crossorigin="anonymous">
<script src="https://cdnjs.cloudflare.com/ajax/libs/twitter-bootstrap/3.4.1/js/bootstrap.min.js" integrity="sha256-nuL8/2cJ5NDSSwnKD8VqreErSWHtnEP9E7AySL+1ev4=" crossorigin="anonymous"></script><!-- bootstrap-toc --><link rel="stylesheet" href="../bootstrap-toc.css">
<script src="../bootstrap-toc.js"></script><!-- Font Awesome icons --><link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/all.min.css" integrity="sha256-mmgLkCYLUQbXn0B1SRqzHar6dCnv9oZFPEC1g1cwlkk=" crossorigin="anonymous">
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/v4-shims.min.css" integrity="sha256-wZjR52fzng1pJHwx4aV2AO3yyTOXrcDW7jBpJtTwVxw=" crossorigin="anonymous">
<!-- clipboard.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/2.0.6/clipboard.min.js" integrity="sha256-inc5kl9MA1hkeYUt+EC3BhlIgyp/2jDIyBLS6k3UxPI=" crossorigin="anonymous"></script><!-- headroom.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/headroom.min.js" integrity="sha256-AsUX4SJE1+yuDu5+mAVzJbuYNPHj/WroHuZ8Ir/CkE0=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/jQuery.headroom.min.js" integrity="sha256-ZX/yNShbjqsohH1k95liqY9Gd8uOiE1S4vZc+9KQ1K4=" crossorigin="anonymous"></script><!-- pkgdown --><link href="../pkgdown.css" rel="stylesheet">
<script src="../pkgdown.js"></script><meta property="og:title" content="hmclearn:  Logistic Regression Example">
<meta property="og:description" content="hmclearn">
<!-- mathjax --><script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js" integrity="sha256-nvJJv9wWKEm88qvoQl9ekL2J+k/RWIsaSScxxlsrv8k=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/config/TeX-AMS-MML_HTMLorMML.js" integrity="sha256-84DKXVJXs0/F8OTMzX4UR909+jtl4G7SPypPavF+GfA=" crossorigin="anonymous"></script><!--[if lt IE 9]>
<script src="https://oss.maxcdn.com/html5shiv/3.7.3/html5shiv.min.js"></script>
<script src="https://oss.maxcdn.com/respond/1.4.2/respond.min.js"></script>
<![endif]-->
</head>
<body data-spy="scroll" data-target="#toc">
    <div class="container template-article">
      <header><div class="navbar navbar-default navbar-fixed-top" role="navigation">
  <div class="container">
    <div class="navbar-header">
      <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar" aria-expanded="false">
        <span class="sr-only">Toggle navigation</span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
      </button>
      <span class="navbar-brand">
        <a class="navbar-link" href="../index.html">hmclearn</a>
        <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.0.2</span>
      </span>
    </div>

    <div id="navbar" class="navbar-collapse collapse">
      <ul class="nav navbar-nav">
<li>
  <a href="../index.html">
    <span class="fas fa fas fa-home fa-lg"></span>
     
  </a>
</li>
<li>
  <a href="../reference/index.html">Reference</a>
</li>
<li class="dropdown">
  <a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
    Articles
     
    <span class="caret"></span>
  </a>
  <ul class="dropdown-menu" role="menu">
<li>
      <a href="../articles/linear_mixed_effects_hmclearn.html">hmclearn package:  Linear Mixed Effects Regression Example</a>
    </li>
    <li>
      <a href="../articles/linear_regression_hmclearn.html">hmclearn package:  Linear Regression Example</a>
    </li>
    <li>
      <a href="../articles/logistic_mixed_effects_hmclearn.html">hmclearn package:  Logistic Mixed Effects Regression Example</a>
    </li>
    <li>
      <a href="../articles/logistic_regression_hmclearn.html">hmclearn:  Logistic Regression Example</a>
    </li>
    <li>
      <a href="../articles/poisson_regression_hmclearn.html">hmclearn:  Poisson Regression Example</a>
    </li>
  </ul>
</li>
<li>
  <a href="../news/index.html">Changelog</a>
</li>
      </ul>
<ul class="nav navbar-nav navbar-right"></ul>
</div>
<!--/.nav-collapse -->
  </div>
<!--/.container -->
</div>
<!--/.navbar -->

      

      </header><div class="row">
  <div class="col-md-9 contents">
    <div class="page-header toc-ignore">
      <h1 data-toc-skip>hmclearn: Logistic Regression Example</h1>
                        <h4 class="author">Samuel Thomas</h4>
            
            <h4 class="date">2020-06-01</h4>
      
      
      <div class="hidden name"><code>logistic_regression_hmclearn.Rmd</code></div>

    </div>

    
    
<div id="introduction" class="section level2">
<h2 class="hasAnchor">
<a href="#introduction" class="anchor"></a>Introduction</h2>
<p>This vignette demonstrates fitting a logistic regression model via Hamiltonian Monte Carlo (HMC) using the <strong>hmclearn</strong> package.</p>
<div class="sourceCode" id="cb1"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/r/base/library.html">library</a></span>(<span class="no">hmclearn</span>)
<span class="fu"><a href="https://rdrr.io/r/base/library.html">library</a></span>(<span class="no">MCMCpack</span>)</pre></body></html></div>
<p>For binary response, we let</p>
<p><span class="math display">\[
p = Pr(\mathbf{y} = 1 | \mathbf{X}) = [1 + e^{-\mathbf{X}\boldsymbol\beta}]^{-1}.
\]</span></p>
<p>The vector of responses is <span class="math inline">\(\mathbf{y} = (y_1, ..., y_n)^T\)</span>. The covariate values for the th subject are <span class="math inline">\(\mathbf{x}_i^T = (x_{i0}, ..., x_{iq})\)</span> for <span class="math inline">\(q\)</span> covariates plus an intercept. We write the full design matrix as <span class="math inline">\(\mathbf{X} = (\mathbf{x}_1^T, ..., \mathbf{x}_n^T)^T \in \mathbb{R}^{n\times(q+1)}\)</span> for <span class="math inline">\(n\)</span> observations. The regression coefficients are a vector of length <span class="math inline">\(q + 1\)</span>, <span class="math inline">\(\boldsymbol\beta = (\beta_0, ..., \beta_q)^T\)</span>.</p>
</div>
<div id="derive-log-posterior-and-gradient-for-hmc" class="section level2">
<h2 class="hasAnchor">
<a href="#derive-log-posterior-and-gradient-for-hmc" class="anchor"></a>Derive log posterior and gradient for HMC</h2>
<p>We develop the likelihood and log likelihood,</p>
<p><span class="math display">\[
\begin{aligned}
f(\mathbf{y} | \mathbf{X},\boldsymbol\beta) &amp;= \prod_{i=1}^n p^{y_i} (1-p)^{1-y_i}, \\
&amp;= \prod_{i=1}^{n} \left(\frac{1}{1+e^{-\mathbf{x}_i^T\boldsymbol\beta}}\right)^{y_i} 
\left(\frac{e^{-\mathbf{x}_i^T\boldsymbol\beta}}{1+e^{-\mathbf{x}_i^T\boldsymbol\beta}}\right)^{1-y_i},   \\
\log f(\mathbf{y} | \mathbf{X}, \boldsymbol\beta) &amp;= \sum_{i=1}^n -y_i\log(1+e^{-\mathbf{x}_i^T\boldsymbol\beta}) + (1-y_i)(-\mathbf{x}_i\boldsymbol\beta - \log(1+e^{-\mathbf{x}_i^T\boldsymbol\beta})),  \\
&amp;= \sum_{i=1}^n -\log(1+e^{-\mathbf{x}_i\boldsymbol\beta}) - \mathbf{x}_i^T\boldsymbol\beta(1 - y_i),  \\
&amp;= \sum_{i=1}^n \mathbf{x}_i^T\boldsymbol\beta(y_i - 1) - \log(1 + e^{-\mathbf{x}_i^T\boldsymbol\beta}),  \\
&amp;= \boldsymbol\beta^T\mathbf{X}^T(\mathbf{y} - \mathbf{1}_n) - \mathbf{1}_n^T [ \log( 1 + e^{-\mathbf{X}\boldsymbol\beta})]. 
\end{aligned}
\]</span></p>
<p>We set a multivariate normal prior for <span class="math inline">\(\boldsymbol\beta\)</span></p>
<p><span class="math display">\[
\begin{aligned}
\boldsymbol\beta &amp;\sim N(0, \sigma_\beta^2 \mathbf{I}), 
\end{aligned}
\]</span></p>
<p>with pdf, omitting constants</p>
<p><span class="math display">\[
\begin{aligned}
\pi(\boldsymbol\beta | \sigma_\beta^2) &amp;= \frac{1}{\sqrt{\lvert 2\pi \sigma_\beta^2 \rvert }}e^{-\frac{1}{2}\boldsymbol\beta^T \boldsymbol\beta / \sigma_\beta^2},  \\
\log \pi(\boldsymbol\beta | \sigma_\beta^2) &amp;= -\frac{1}{2}\log(2\pi \sigma_\beta^2) - \frac{1}{2}\boldsymbol\beta^T \boldsymbol\beta / \sigma_\beta^2, \\
&amp;\propto -\frac{1}{2}\log \sigma_\beta^2 - \frac{\boldsymbol\beta^T\boldsymbol\beta}{2\sigma_\beta^2}.
\end{aligned}
\]</span></p>
<p>Next, we derive the log posterior, omitting constants,</p>
<p><span class="math display">\[
\begin{aligned}
f(\boldsymbol\beta | \mathbf{X}, \mathbf{y}, \sigma_\beta^2) &amp;\propto f(\mathbf{y} | \mathbf{X}, \boldsymbol\beta)  \pi(\boldsymbol\beta | \sigma_\beta^2), \\
\log f(\boldsymbol\beta | \mathbf{X}, \mathbf{y}, \sigma_\beta^2) &amp; \propto \log f(\mathbf{y} | \mathbf{X}, \boldsymbol\beta) + \log \pi(\boldsymbol\beta|\sigma_\beta^2), \\
&amp;\propto \sum_{i=1}^n \mathbf{x}_i^T\boldsymbol\beta(y_i - 1) - \log(1 + e^{-\mathbf{x}_i^T\boldsymbol\beta}) - \frac{1}{2}\boldsymbol\beta^T\boldsymbol\beta / \sigma_\beta^2,  \\
&amp;\propto \boldsymbol\beta^T\mathbf{X}^T(\mathbf{y} - \mathbf{1}_n) - \mathbf{1}_n^T[\log( 1 + e^{-\mathbf{X}\boldsymbol\beta})] - \frac{\boldsymbol\beta^T\boldsymbol\beta}{2\sigma_\beta^2}.
\end{aligned}
\]</span></p>
<p>Next, we need to derive the gradient of the log posterior for the leapfrog function</p>
<p><span class="math display">\[
\begin{aligned}
\nabla_{\boldsymbol\beta} \log f(\boldsymbol\beta, \mathbf{X}, \mathbf{y}, \sigma_\beta^2)  &amp;\propto \mathbf{X}^T \left ( \mathbf{y} - \mathbf{1}_n+ \frac{e^{-\mathbf{X}\boldsymbol\beta}}{1 + e^{-\mathbf{X}\boldsymbol\beta}}\right) - \boldsymbol\beta / \sigma_\beta^2
\end{aligned}
\]</span></p>
</div>
<div id="logistic-regression-example-data" class="section level2">
<h2 class="hasAnchor">
<a href="#logistic-regression-example-data" class="anchor"></a>Logistic Regression Example Data</h2>
<p>The user must define provide the design matrix directly for use in <strong>hmclearn</strong>. Our first step is to load the data and store the design matrix <span class="math inline">\(X\)</span> and dependent variable vector <span class="math inline">\(y\)</span>.</p>
<p>First, we load the Endometrial cancer data set (Heinze and Schember 2002) and create <span class="math inline">\(X\)</span> and <span class="math inline">\(y\)</span>. This example also appears in Agresti (2015), and we compare results to his.</p>
<div class="sourceCode" id="cb2"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/r/utils/data.html">data</a></span>(<span class="no">Endometrial</span>)

<span class="co"># data prep</span>
<span class="no">Endometrial</span>$<span class="no">PI2</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/with.html">with</a></span>(<span class="no">Endometrial</span>, (<span class="no">PI</span> - <span class="fu"><a href="https://rdrr.io/r/base/mean.html">mean</a></span>(<span class="no">PI</span>)) / <span class="fu"><a href="https://rdrr.io/r/stats/sd.html">sd</a></span>(<span class="no">PI</span>))
<span class="no">Endometrial</span>$<span class="no">EH2</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/with.html">with</a></span>(<span class="no">Endometrial</span>, (<span class="no">EH</span> - <span class="fu"><a href="https://rdrr.io/r/base/mean.html">mean</a></span>(<span class="no">EH</span>)) / <span class="fu"><a href="https://rdrr.io/r/stats/sd.html">sd</a></span>(<span class="no">EH</span>))
<span class="no">Endometrial</span>$<span class="no">NV2</span> <span class="kw">&lt;-</span> <span class="no">Endometrial</span>$<span class="no">NV</span> - <span class="fl">0.5</span>

<span class="no">X</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/cbind.html">cbind</a></span>(<span class="fl">1</span>, <span class="fu"><a href="https://rdrr.io/r/base/matrix.html">as.matrix</a></span>(<span class="no">Endometrial</span>[, <span class="fu"><a href="https://rdrr.io/r/base/which.html">which</a></span>(<span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">Endometrial</span>)
            <span class="kw">%in%</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="st">"PI2"</span>, <span class="st">"EH2"</span>, <span class="st">"NV2"</span>))]))
<span class="no">y</span> <span class="kw">&lt;-</span> <span class="no">Endometrial</span>$<span class="no">HG</span>

<span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">X</span>) <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="st">"(Intercept)"</span>, <span class="st">"PI2"</span>, <span class="st">"EH2"</span>, <span class="st">"NV2"</span>)</pre></body></html></div>
</div>
<div id="comparison-models---frequentist-and-metropolis-hastings" class="section level2">
<h2 class="hasAnchor">
<a href="#comparison-models---frequentist-and-metropolis-hastings" class="anchor"></a>Comparison models - Frequentist and Metropolis-Hastings</h2>
<p>To compare results, we first fit a standard linear model using the frequentist function <em>glm</em>. Note the high standard error estimates for the Intercept and NV2. This occurs because all 13 subjects with neovasculation present also have <span class="math inline">\(y_i = 1\)</span>. See section 11.1.8 of Agresti (2015).</p>
<div class="sourceCode" id="cb3"><html><body><pre class="r"><span class="no">f</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/stats/glm.html">glm</a></span>(<span class="no">y</span> ~ <span class="no">X</span>-<span class="fl">1</span>, <span class="kw">family</span> <span class="kw">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/family.html">binomial</a></span>())
<span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span>(<span class="no">f</span>)
<span class="co">#&gt; </span>
<span class="co">#&gt; Call:</span>
<span class="co">#&gt; glm(formula = y ~ X - 1, family = binomial())</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; Deviance Residuals: </span>
<span class="co">#&gt;      Min        1Q    Median        3Q       Max  </span>
<span class="co">#&gt; -1.50137  -0.64108  -0.29432   0.00016   2.72777  </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; Coefficients:</span>
<span class="co">#&gt;               Estimate Std. Error z value Pr(&gt;|z|)    </span>
<span class="co">#&gt; X(Intercept)    7.8411   857.8755   0.009 0.992707    </span>
<span class="co">#&gt; XPI2           -0.4217     0.4432  -0.952 0.341333    </span>
<span class="co">#&gt; XEH2           -1.9219     0.5599  -3.433 0.000597 ***</span>
<span class="co">#&gt; XNV2           18.1856  1715.7509   0.011 0.991543    </span>
<span class="co">#&gt; ---</span>
<span class="co">#&gt; Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; (Dispersion parameter for binomial family taken to be 1)</span>
<span class="co">#&gt; </span>
<span class="co">#&gt;     Null deviance: 109.517  on 79  degrees of freedom</span>
<span class="co">#&gt; Residual deviance:  55.393  on 75  degrees of freedom</span>
<span class="co">#&gt; AIC: 63.393</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; Number of Fisher Scoring iterations: 17</span></pre></body></html></div>
<p>As a result of challenges with the frequentist fit, Agresti (2015) in section 10.3.2 also fits this logistic regression model using a random walk Metropolis algorithm from the  package (Martin, Quinn, Park 2011). We fit the model with the same number of MCMC iterations, prior distributions, and hyperparameters as in the text.</p>
<p>This model also assigns a normal prior for <span class="math inline">\(\boldsymbol\beta\)</span> where <span class="math inline">\(\sigma_\beta^2 = 100\)</span>.  uses a precision parameterization, where <span class="math inline">\(B0 = 1/\sigma_\beta^2\)</span>. The mean hyperparameter is 0 as specified in <span class="math inline">\(b0\)</span>.</p>
<p>Note the high number of Metropolis iterations specified in  and the time required to fit, even with optimized compiled code. Our abbreviation MH refers to the general Metropolis-Hastings algorithm.</p>
<div class="sourceCode" id="cb4"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/r/base/Random.html">set.seed</a></span>(<span class="fl">921</span>)
<span class="no">t1.mh</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/Sys.time.html">Sys.time</a></span>()
<span class="no">fitMH</span> <span class="kw">&lt;-</span> <span class="kw pkg">MCMCpack</span><span class="kw ns">::</span><span class="fu"><a href="https://rdrr.io/pkg/MCMCpack/man/MCMClogit.html">MCMClogit</a></span>(<span class="no">HG</span> ~ <span class="no">NV2</span> + <span class="no">PI2</span> + <span class="no">EH2</span>, <span class="kw">data</span><span class="kw">=</span><span class="no">Endometrial</span>,
                   <span class="kw">mcmc</span><span class="kw">=</span><span class="fl">10000000</span>, <span class="kw">b0</span><span class="kw">=</span><span class="fl">0</span>, <span class="kw">B0</span><span class="kw">=</span><span class="fl">0.01</span>)
<span class="no">t2.mh</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/Sys.time.html">Sys.time</a></span>()</pre></body></html></div>
<p>The time to fit the logistic regression using MH is</p>
<div class="sourceCode" id="cb5"><html><body><pre class="r"><span class="no">t2.mh</span> - <span class="no">t1.mh</span>
<span class="co">#&gt; Time difference of 2.569034 mins</span></pre></body></html></div>
<p>We store the Metropolis parameter estimates in .</p>
<div class="sourceCode" id="cb6"><html><body><pre class="r">(<span class="no">summMH</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span>(<span class="no">fitMH</span>))</pre></body></html></div>
<pre><code>#&gt; 
#&gt; Iterations = 1001:10001000
#&gt; Thinning interval = 1 
#&gt; Number of chains = 1 
#&gt; Sample size per chain = 1e+07 
#&gt; 
#&gt; 1. Empirical mean and standard deviation for each variable,
#&gt;    plus standard error of the mean:
#&gt; 
#&gt;                Mean     SD  Naive SE Time-series SE
#&gt; (Intercept)  3.2114 2.5585 0.0008091      0.0026811
#&gt; NV2          9.1134 5.0921 0.0016103      0.0053685
#&gt; PI2         -0.4723 0.4540 0.0001436      0.0006411
#&gt; EH2         -2.1381 0.5935 0.0001877      0.0008742
#&gt; 
#&gt; 2. Quantiles for each variable:
#&gt; 
#&gt;                2.5%     25%     50%     75%   97.5%
#&gt; (Intercept) -0.3453  1.2675  2.7207  4.6828  9.3351
#&gt; NV2          2.1065  5.2285  8.1258 12.0409 21.3190
#&gt; PI2         -1.4132 -0.7664 -0.4545 -0.1584  0.3664
#&gt; EH2         -3.4015 -2.5152 -2.1004 -1.7222 -1.0815</code></pre>
</div>
<div id="fit-model-using-hmc" class="section level2">
<h2 class="hasAnchor">
<a href="#fit-model-using-hmc" class="anchor"></a>Fit model using <em>hmc</em>
</h2>
<p>Next, we fit the logistic regression model using HMC with the same hyperprior for <span class="math inline">\(\sigma_\beta^2\)</span> as in Agresti’s example.</p>
<div class="sourceCode" id="cb8"><html><body><pre class="r"><span class="no">N</span> <span class="kw">&lt;-</span> <span class="fl">1e4</span>

<span class="fu"><a href="https://rdrr.io/r/base/Random.html">set.seed</a></span>(<span class="fl">412</span>)
<span class="no">t1.hmc</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/Sys.time.html">Sys.time</a></span>()
 <span class="no">f_hmc</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/hmc.html">hmc</a></span>(<span class="kw">N</span> <span class="kw">=</span> <span class="no">N</span>, <span class="kw">theta.init</span> <span class="kw">=</span> <span class="fu"><a href="https://rdrr.io/r/base/rep.html">rep</a></span>(<span class="fl">1</span>, <span class="fl">4</span>),
            <span class="kw">epsilon</span> <span class="kw">=</span> <span class="fl">0.1</span>, <span class="kw">L</span> <span class="kw">=</span> <span class="fl">20</span>,
            <span class="kw">logPOSTERIOR</span> <span class="kw">=</span> <span class="no">logistic_posterior</span>,
            <span class="kw">glogPOSTERIOR</span> <span class="kw">=</span> <span class="no">g_logistic_posterior</span>,
            <span class="kw">varnames</span> <span class="kw">=</span> <span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">X</span>),
            <span class="kw">param</span><span class="kw">=</span><span class="fu"><a href="https://rdrr.io/r/base/list.html">list</a></span>(<span class="kw">y</span> <span class="kw">=</span> <span class="no">y</span>, <span class="kw">X</span><span class="kw">=</span><span class="no">X</span>, <span class="kw">sig2beta</span><span class="kw">=</span><span class="fl">100</span>),
            <span class="kw">parallel</span><span class="kw">=</span><span class="fl">TRUE</span>, <span class="kw">chains</span><span class="kw">=</span><span class="fl">2</span>)
<span class="no">t2.hmc</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/Sys.time.html">Sys.time</a></span>()</pre></body></html></div>
<p>HMC requires fewer iterations and lower computation time than MH in this example.</p>
<div class="sourceCode" id="cb9"><html><body><pre class="r"><span class="no">t2.hmc</span> - <span class="no">t1.hmc</span>
<span class="co">#&gt; Time difference of 20.98491 secs</span></pre></body></html></div>
</div>
<div id="mcmc-summary-and-diagnostics" class="section level2">
<h2 class="hasAnchor">
<a href="#mcmc-summary-and-diagnostics" class="anchor"></a>MCMC summary and diagnostics</h2>
<p>The acceptance ratio for each of the HMC chains is sufficiently high for an efficient simulation.</p>
<div class="sourceCode" id="cb10"><html><body><pre class="r"><span class="no">f_hmc</span>$<span class="no">accept</span>/<span class="no">N</span>
<span class="co">#&gt; [1] 0.9175 0.9186</span></pre></body></html></div>
<p>The posterior quantiles are summarized after removing an initial <em>burnin</em> period. The <span class="math inline">\(\hat{R}\)</span> statistics provide an indication of convergence. Values close to one indicate that the multiple MCMC chains coverged to the same distribution, while values above 1.1 indicate possible convergence problems. All <span class="math inline">\(\hat{R}\)</span> values in our example are close to one.</p>
<div class="sourceCode" id="cb11"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span>(<span class="no">f_hmc</span>, <span class="kw">burnin</span><span class="kw">=</span><span class="fl">3000</span>)
<span class="co">#&gt; Summary of MCMC simulation</span>
<span class="co">#&gt;                   2.5%          5%       25%        50%        75%        95%</span>
<span class="co">#&gt; (Intercept) -0.3074032 -0.01396594  1.253784  2.6564096  4.5602601  7.8440969</span>
<span class="co">#&gt; PI2         -1.4050342 -1.24701482 -0.771735 -0.4684764 -0.1716551  0.2268731</span>
<span class="co">#&gt; EH2         -3.3370587 -3.12960806 -2.488422 -2.0969937 -1.7318884 -1.2628641</span>
<span class="co">#&gt; NV2          2.0706893  2.68949683  5.217250  8.0376322 11.8207796 18.3763529</span>
<span class="co">#&gt;                  97.5%      rhat</span>
<span class="co">#&gt; (Intercept)  9.1070020 1.0018737</span>
<span class="co">#&gt; PI2          0.3561943 0.9999416</span>
<span class="co">#&gt; EH2         -1.1126949 0.9999392</span>
<span class="co">#&gt; NV2         20.8586412 1.0018367</span></pre></body></html></div>
<p>Trace plots provide a visual indication of stationarity. These plots indicate that the MCMC chains are reasonably stationary.</p>
<div class="sourceCode" id="cb12"><html><body><pre class="r"><span class="fu"><a href="../reference/hmclearn-plots.html">mcmc_trace</a></span>(<span class="no">f_hmc</span>, <span class="kw">burnin</span><span class="kw">=</span><span class="fl">3000</span>)</pre></body></html></div>
<p><img src="logistic_regression_hmclearn_files/figure-html/unnamed-chunk-13-1.png" width="576"></p>
<p>Histograms of the posterior distribution show that HMC parameter estimates align with Metropolis estimates from <strong>MCMCpack</strong>.</p>
<div class="sourceCode" id="cb13"><html><body><pre class="r"><span class="no">beta.mh</span> <span class="kw">&lt;-</span> <span class="no">summMH</span>$<span class="no">quantiles</span>[, <span class="st">"50%"</span>][<span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">X</span>)]
<span class="fu"><a href="../reference/diagplots.html">diagplots</a></span>(<span class="no">f_hmc</span>, <span class="kw">burnin</span><span class="kw">=</span><span class="fl">3000</span>, <span class="kw">comparison.theta</span><span class="kw">=</span><span class="no">beta.mh</span>)
<span class="co">#&gt; $histogram</span></pre></body></html></div>
<p><img src="logistic_regression_hmclearn_files/figure-html/unnamed-chunk-14-1.png" width="576"></p>
</div>
<div id="source" class="section level2">
<h2 class="hasAnchor">
<a href="#source" class="anchor"></a>Source</h2>
<p>Heinze, G., &amp; Schemper, M. (2002). <em>A solution to the problem of separation in logistic regression</em>. Statistics in medicine, 21(16), 2409-2419.</p>
</div>
<div id="reference" class="section level2">
<h2 class="hasAnchor">
<a href="#reference" class="anchor"></a>Reference</h2>
<p>Agresti, A. (2015). <em>Foundations of linear and generalized linear models</em>. John Wiley &amp; Sons. ISBN: 978-1-118-73003-4</p>
<p>Martin AD, Quinn KM, Park JH (2011). “MCMCpack: Markov Chain Monte Carlo in R.” <em>Journal of Statistical</em> Software, 42(9), 22. <a href="http://www.jstatsoft.org/v42/i09/" class="uri">http://www.jstatsoft.org/v42/i09/</a>.</p>
</div>
  </div>

  <div class="col-md-3 hidden-xs hidden-sm" id="pkgdown-sidebar">

        <nav id="toc" data-toggle="toc"><h2 data-toc-skip>Contents</h2>
    </nav>
</div>

</div>



      <footer><div class="copyright">
  <p>Developed by Samuel Thomas.</p>
</div>

<div class="pkgdown">
  <p>Site built with <a href="https://pkgdown.r-lib.org/">pkgdown</a> 1.5.1.</p>
</div>

      </footer>
</div>

  


  </body>
</html>
