<!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 package:  Logistic Mixed Effects 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 package:  Logistic Mixed Effects 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 package: Logistic Mixed Effects Regression Example</h1>
                        <h4 class="author">Samuel Thomas</h4>
            
            <h4 class="date">2020-06-01</h4>
      
      
      <div class="hidden name"><code>logistic_mixed_effects_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 mixed effects regression model via Hamiltonian Monte Carlo (HMC) using the <strong>hmclearn</strong> package.</p>
<p>For a mixed effects model with binary response, we let</p>
<p><span class="math display">\[
Pr(\mathbf{y} = \mathbf{1} | \mathbf{X}, \mathbf{Z}) = [1 + e^{-\mathbf{X}\boldsymbol\beta-\mathbf{Z} \mathbf{u} } ]^{-1},
\]</span> or</p>
<p><span class="math display">\[
\begin{aligned}
\text{logit}[P(\mathbf{y} = 1 | \mathbf{u})] &amp;= \mathbf{X}\boldsymbol\beta + \mathbf{Z}\mathbf{u}, \\
\mathbf{u} &amp;\sim N(0, \mathbf{G}).
\end{aligned}
\]</span></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>)</pre></body></html></div>
<p>The response for each subject is a vector <span class="math inline">\(\mathbf{y} = (\mathbf{y}_1, ..., \mathbf{y}_n)\)</span> for <span class="math inline">\(n\)</span> subjects <span class="math inline">\(i= 1, ..., n\)</span>. Each subject has <span class="math inline">\(d\)</span> observations <span class="math inline">\(\mathbf{y}_i = (y_{i1}, ..., y_{id})\)</span> and we let <span class="math inline">\(j = 1, ..., d\)</span>. The fixed effect design matrix is composed of matrices for each subject, <span class="math inline">\(\mathbf{X} = (\mathbf{X}_1, ..., \mathbf{X}_n)\)</span>, and <span class="math inline">\(\mathbf{X}_i \in \mathbb{R}^{d\times (q+1)}\)</span> for the fixed effects parameters <span class="math inline">\(\boldsymbol\beta = (\beta_0, ..., \beta_q)\)</span>. The full fixed effects design matrix is therefore <span class="math inline">\(\mathbf{X} \in \mathbb{R}^{nd \times (q+1)}\)</span>.</p>
<p>For random effects, <span class="math inline">\(\mathbf{Z} = \text{diag}(\mathbf{Z}_1, ..., \mathbf{Z}_n)\)</span>, with individual random effects matrices <span class="math inline">\(\mathbf{Z}_i\)</span> for each of the <span class="math inline">\(i\)</span> subjects. A random intercept model specifies <span class="math inline">\(\mathbf{Z}_i\)</span> as a column vector of ones where <span class="math inline">\(\mathbf{Z}_i = \mathbf{z}_i = \mathbf{1}_d\)</span>. The full random effects design matrix <span class="math inline">\(\mathbf{Z} \in \mathbb{R}^{nd\times n}\)</span>. The parameterization for random effects is <span class="math inline">\(\mathbf{u} = (\mathbf{u}_1, ..., \mathbf{u}_n)^T\)</span> with vectors <span class="math inline">\(\mathbf{u}_i\)</span> for each subject. A random intercept model is somewhat simplified where <span class="math inline">\(\mathbf{u}_i = u_i\)</span> denotes a single random intercept parameter for each subject <span class="math inline">\(i\)</span>, and <span class="math inline">\(\mathbf{u} = (u_1, ..., u_n)^T\)</span>.</p>
<p>We set <span class="math inline">\(\mathbf{u}\)</span> as one of our priors, following a multivariate normal distribution, <span class="math inline">\(\mathbf{u} \sim N(0, \mathbf{G})\)</span>. For our random intercept model, the specification of the covariance matrix <span class="math inline">\(\mathbf{G}\)</span> is expanded to facilitate efficient sampling using HMC. We let <span class="math inline">\(\mathbf{u} = \mathbf{G}^{1/2}\boldsymbol\tau\)</span> where <span class="math inline">\(\mathbf{G}^{1/2} = \lambda \mathbf{I}_n\)</span>. An additional parameter <span class="math inline">\(\boldsymbol\tau = (\tau_1, ..., \tau_n)^T\)</span> where each of these parameters is standard normal <span class="math inline">\(\tau_i \sim N(0, 1)\)</span>. The full covariance matrix is then <span class="math inline">\(\mathbf{G} = \lambda^2 \mathbf{I}_n \boldsymbol\tau\)</span>.</p>
<p>The parameterization approach for this model uses a strategy recommended by Betancourt, Girolami (2013) to facilitate more efficient sampling in HMC.</p>
<p>Further, we select a half-t family of distributions appropriate for hierarchical models per Gelman (2006). This parameterization is well-behaved around 0, in contrast to inverse gamma, and provides flexibility for informed priors.</p>
<p>We select a parameterization of <span class="math inline">\(\mathbf{G}\)</span> such that the likelihood and its gradient can be derived for HMC. To this end, we uses LDL decomposition of <span class="math inline">\(\mathbf{G}\)</span> to form a flexibile parameterization that can easily handle restrictions (Chan, Jelizkov 2009).</p>
<p><span class="math display">\[
\begin{aligned}
\mathbf{u} &amp;\sim N(0, \mathbf{G}),  \\
\mathbf{G} &amp;= \mathbf{L} \mathbf{D} \mathbf{L}^T,  \\
&amp;= \mathbf{L} \mathbf{D}^{1/2} \mathbf{D}^{1/2} \mathbf{L}^T. \\
\end{aligned}
\]</span></p>
<p>Let <span class="math inline">\(\boldsymbol{\lambda} = (\lambda_1, ..., \lambda_p)\)</span> denote the diagonal elements of <span class="math inline">\(\mathbf{D}^{1/2}\)</span> where <span class="math inline">\(p\)</span> indicates the number of random effect parameters, specified as <em>nrandom</em> in <strong>hmclearn</strong>. A future release of <strong>hmclearn</strong> will allow prior specification for the off-diagonal elements of <span class="math inline">\(L\)</span>. For the current version, we let <span class="math inline">\(L = I_{p}\)</span>.</p>
<p><span class="math display">\[
\mathbf{D}^{1/2} := 
\begin{pmatrix}
\lambda_1 &amp; 0 &amp; ... &amp; 0 \\
0 &amp; \lambda_2 &amp; 0 ... &amp; 0 \\
... &amp; ... &amp; ... &amp; ... \\
0 &amp; 0 &amp; ... &amp; \lambda_p
\end{pmatrix}, \quad
\mathbf{L} :=
\begin{pmatrix}
1 &amp; 0 &amp; 0 &amp; ... &amp; 0 \\
0 &amp; 1 &amp; 0 &amp; ... &amp; 0 \\
0 &amp; 0 &amp; 1 &amp; ... &amp; ... \\
... &amp; ... &amp; ... &amp; ... &amp; ... \\
0 &amp; 0 &amp; ... &amp; ... &amp; 1 \\
\end{pmatrix}.
\]</span></p>
<p>We set the prior for <span class="math inline">\(\boldsymbol\beta\)</span> as multivariate normal with variance <span class="math inline">\(\sigma_\beta^2\)</span>, a hyperparameter set by the analyst. The priors for <span class="math inline">\(\boldsymbol\lambda\)</span> are half-t per Gelman (2006) on hierarchical models.</p>
<p><span class="math display">\[
\begin{aligned}
\pi(\boldsymbol\beta | \sigma_\beta^2) &amp;\propto N(0, \sigma_\beta^2 \mathbf{I}), \\
\pi(\boldsymbol\lambda) &amp;\sim  \left(1 + \frac{1}{\nu_\lambda}\left(\frac{\boldsymbol\lambda}{A_\lambda} \right)^2 \right)^{-(\nu+1)/2}.
\end{aligned}
\]</span></p>
<p>We want proposals of <span class="math inline">\(\boldsymbol\lambda\)</span> over the real number line. Therfore we derive the distribution of the transformed parameter <span class="math inline">\(\boldsymbol\xi\)</span> based on a change of variable</p>
<p><span class="math display">\[
\begin{aligned}
\boldsymbol\xi &amp;:= \log\boldsymbol\lambda,  \\
\boldsymbol\lambda &amp;:= e^{\boldsymbol\xi}. 
\end{aligned}
\]</span></p>
<p>We need to compute the Jacobian of the transformation</p>
<p><span class="math display">\[
\begin{aligned}
\pi_{\boldsymbol\xi}(\boldsymbol\xi) &amp;= \pi_{\boldsymbol\xi}(g^{-1}(\boldsymbol\xi)) \left\lvert \frac{d\boldsymbol\lambda}{d\boldsymbol\xi}  \right\rvert, \\
&amp;= \pi_{\boldsymbol\lambda} (e^{\boldsymbol\xi})\lvert e^{\boldsymbol\xi}\rvert, \\
&amp;= \left(1 + \frac{1}{\nu_\xi} \frac{e^{2\boldsymbol\xi}}{A_{\xi}^2}\right)^{-\frac{\nu_{\xi}+1}{2}}e^{\boldsymbol\xi}, \\
\log \pi(\boldsymbol\xi) &amp;\propto -\frac{\nu_{\xi}+1}{2}\log\left(1 + \frac{1}{\nu_{\xi}} \frac{e^{2\boldsymbol\xi}}{A_{\xi}^2}\right)+ \boldsymbol\xi.
\end{aligned}
\]</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>First, we derive the likelihood for our logistic mixed effects regression model.</p>
<p><span class="math display">\[
\begin{aligned}
f(\mathbf{y}|\mathbf{X}, \mathbf{Z}, \boldsymbol\beta, \mathbf{u}) &amp;= \prod_{i=1}^n \prod_{j=1}^d p(y_{ij})^{y_{ij}}(1-p(y_{ij}))^{1-y_{ij}}, \\
&amp;= \prod_{i=1}^n \prod_{j=1}^d \left(\frac{1}{1 + e^{-\mathbf{x}_{ij}^T\boldsymbol\beta - \mathbf{z}_{ij}^T\mathbf{u}_i}} \right)^{y_{ij}}  \left( \frac{e^{-\mathbf{x}_{ij}^T\boldsymbol\beta - \mathbf{z}_{ij}^T\mathbf{u}_i}}{1+e^{-\mathbf{x}_{ij}^T\boldsymbol\beta - \mathbf{z}_{ij}^T\mathbf{u}_i}}  \right)^{1-y_{ij}},
 \\
\log f(\mathbf{y}|\mathbf{X}, \mathbf{Z}, \boldsymbol\beta, \mathbf{u}) &amp;= \sum_{i=1}^n \sum_{j=1}^d -y_{ij} \log\left(1+e^{-\mathbf{x}_{ij}^T\boldsymbol\beta - \mathbf{z}_{ij}^T\mathbf{u}_i} \right) + (1 - y_{ij})\log e^{-\mathbf{x}_{ij}^T\boldsymbol\beta - \mathbf{z}_{ij}^T\mathbf{u}_i} - (1 - y_{ij})\log\left(1+e^{-\mathbf{x}_{ij}^T\boldsymbol\beta - \mathbf{z}_{ij}^T\mathbf{u}_i} \right), \\
&amp;= \sum_{i=1}^n\sum_{j=1}^d -y_{ij} \log\left(1+e^{-\mathbf{x}_{ij}^T\boldsymbol\beta - \mathbf{z}_{ij}\mathbf{u}_i} \right) -\mathbf{x}_{ij}\boldsymbol\beta -\mathbf{z}_{ij}\mathbf{u}_i + y_{ij}\mathbf{x}_{ij}^T\boldsymbol\beta + y_{ij}\mathbf{z}_{ij}^T\mathbf{u}_i - \\
&amp;\qquad \log\left(1+e^{-\mathbf{x}_{ij}^T\boldsymbol\beta - \mathbf{z}_{ij}^T\mathbf{u}_i} \right) + y_{ij} \log\left(1+e^{-\mathbf{x}_{ij}^T\boldsymbol\beta - \mathbf{z}_{ij}^T\mathbf{u}_i} \right), \\
&amp;= \sum_{i=1}^n \sum_{j=1}^d -\mathbf{x}_{ij}^T\boldsymbol\beta - \mathbf{z}_{ij}^T\mathbf{u}_i + y_{ij}\mathbf{x}_{ij}^T\boldsymbol\beta + y_{ij}\mathbf{z}_{ij}^T\mathbf{u}_i - \log\left(1+e^{-\mathbf{x}_{ij}^T\boldsymbol\beta - \mathbf{z}_{ij}^T\mathbf{u}_i} \right), \\
&amp;= \sum_{i=1}^n \sum_{j=1}^d (y_{ij}-1 )(\mathbf{x}_{ij}^T\boldsymbol\beta + \mathbf{z}_{ij}^T\mathbf{u}_i) - \log\left(1+e^{-\mathbf{x}_{ij}^T\boldsymbol\beta - \mathbf{z}_{ij}^T\mathbf{u}_i} \right), \\
&amp;= (\mathbf{y} - \mathbf{1}_{nd})^T (\mathbf{X}\boldsymbol\beta + \mathbf{Z}\mathbf{u}) - \mathbf{1}_{nd}\log (1 + e^{-\mathbf{X}\boldsymbol\beta - \mathbf{Z}\mathbf{u}}).
\end{aligned}
\]</span></p>
<p>Additional notation is added to simplify the log likelihood, log posterior, and gradient formulations, <span class="math inline">\(\widetilde{\mathbf{D}}^{1/2} = \mathbf{I}_n \otimes \mathbf{D}^{1/2}\)</span>. Note that <span class="math inline">\(\widetilde{\mathbf{D}}^{1/2} = \left(\widetilde{\mathbf{D}}^{1/2}\right ) ^T\)</span> due to symmetry.</p>
<p>We write the re-parameterized log likelihood,</p>
<p><span class="math display">\[
\begin{aligned}
\log f(\mathbf{y}|\mathbf{X}, \mathbf{Z}, \boldsymbol\beta, \xi, \boldsymbol\tau) &amp;= (\mathbf{y} - \mathbf{1}_{nd})^T (\mathbf{X}\boldsymbol\beta + \mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau) - \mathbf{1}_{nd}\log (1 + e^{-\mathbf{X}\boldsymbol\beta - \mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau}).
\end{aligned}
\]</span></p>
<p>Next, we express the log priors that we use with transformations, omitting constants. Note that the log densities of the priors with log densities include an additive term for the transformed distribution (i.e. <span class="math inline">\(\xi\)</span> and <span class="math inline">\(\gamma\)</span>),</p>
<p><span class="math display">\[
\begin{aligned}
\log \pi(\boldsymbol\beta | \sigma_\beta^2) &amp;\propto -\frac{\boldsymbol\beta^T \boldsymbol\beta}{2\sigma_\beta^2}, \\
\log \pi(\boldsymbol\xi | \nu_\xi, A_\xi) &amp;\propto -\frac{\nu_\xi + 1}{2} \log \left( 1 + \frac{1}{\nu_\xi} \left(\frac{e^\boldsymbol\xi}{A_\xi} \right)^2 \right) + \boldsymbol\xi, \\
\log \pi(\boldsymbol\tau) &amp;\propto -\frac{1}{2}\boldsymbol\tau^T \boldsymbol\tau.
\end{aligned}
\]</span></p>
<p>The full log posterior with transformed variables is the log likelihood plus the log prior. We develop the log posterior omitting constants,</p>
<p><span class="math display">\[
\begin{aligned}
\log f(\boldsymbol\beta, \xi, \boldsymbol\tau|\mathbf{y}, \mathbf{X}, \mathbf{Z}) &amp;\propto \log f(\mathbf{y}|\mathbf{X}, \mathbf{Z}, \boldsymbol\beta, \xi, \boldsymbol\tau) + \log f(\boldsymbol\beta, \xi, \boldsymbol\tau | \sigma_\beta^2, \nu_\xi, A_\xi), \\
&amp;\propto \log f(\mathbf{y}|\mathbf{X}, \mathbf{Z}, \boldsymbol\beta, \xi, \boldsymbol\tau) + \log f(\boldsymbol\beta | \sigma_\beta^2) + \log f(\xi | \nu_\xi, A_\xi) + \log f(\boldsymbol\tau), \\
&amp;\propto  (\mathbf{y} - \mathbf{1}_{nd})^T (\mathbf{X}\boldsymbol\beta + \mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau) - \mathbf{1}_{nd}^T\log (1 + e^{-\mathbf{X}\boldsymbol\beta - \mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau})  - \\ &amp;\qquad \frac{\boldsymbol\beta^T\boldsymbol\beta}{2\sigma_\beta^2}-\frac{\nu_\xi + 1}{2} \log \left( 1 + \frac{1}{\nu_\xi} \left(\frac{e^\xi}{A_\xi} \right)^2 \right) + \xi-\frac{1}{2}\boldsymbol\tau^T \boldsymbol\tau.
\end{aligned}
\]</span></p>
<p>Next we derive the gradient of the log posterior, comprised of partial derivatives for each of our parameters.</p>
<p>We omit <span class="math inline">\(\mathbf{L}\)</span> since this is currently defined as the identity matrix.</p>
<p>We derive the partial derivative of <span class="math inline">\(\boldsymbol\beta\)</span>,</p>
<p><span class="math display">\[
\begin{aligned}
\nabla_\beta \log f(\boldsymbol\beta, \xi, \boldsymbol\tau | \mathbf{y}, \mathbf{X}, \mathbf{Z}, \sigma_\beta^2, \nu_\xi, A_\xi) &amp;\propto \mathbf{X}^T (\mathbf{y} - \mathbf{1}_{nd})+ \mathbf{X}^T \left[ 1 + e^{-\mathbf{X}\boldsymbol\beta - \mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau } \right]^{-1}e^{-\mathbf{X}\boldsymbol\beta-\mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau} - \boldsymbol\beta / \sigma_\beta^2, \\
&amp;\propto \mathbf{X}^T(\mathbf{y} - \mathbf{1}_{nd}) + \mathbf{X}^T \left(\frac{e^{-\mathbf{X}\boldsymbol\beta-\mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau}}{1 + e^{-\mathbf{X}\boldsymbol\beta-\mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau}} \right) - \boldsymbol\beta / \sigma_\beta^2, \\
&amp;\propto  \mathbf{X}^T(\mathbf{y} - \mathbf{1}_{nd}) + \mathbf{X}^T \left(\frac{1}{1 + e^{\mathbf{X}\boldsymbol\beta+\mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau}} \right) - \boldsymbol\beta / \sigma_\beta^2,  \\
&amp;\propto \mathbf{X}^T \left( \mathbf{y} - \mathbf{1}_{nd} +\frac{1}{1 + e^{\mathbf{X}\boldsymbol\beta+\mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau}} \right) - \boldsymbol\beta / \sigma_\beta^2.
\end{aligned}
\]</span></p>
<p>Next we derive the partial derivative of each parameter <span class="math inline">\(\xi_1, ..., \xi_p\)</span> where <span class="math inline">\(jj = 1, ..., p\)</span>,</p>
<p><span class="math display">\[
\begin{aligned}
\nabla_{\xi_{jj} }\log f(\boldsymbol\beta, \xi_{jj}, \boldsymbol\tau| \mathbf{y}, \mathbf{X}, \mathbf{Z}, \sigma_\beta^2, \nu_\xi, A_\xi) 
&amp;\propto e^{\xi_{jj}} \boldsymbol\tau^T \mathbf{Z}^T (\mathbf{y} - \mathbf{1}_{nd})^T + e^{\xi_{jj}}\boldsymbol\tau^T\mathbf{Z}^T\left[1 + e^{-\mathbf{X}\boldsymbol\beta - e^{\xi_{jj}}\mathbf{Z}\boldsymbol\tau} \right ]^{-1} e^{-\mathbf{X}\boldsymbol\beta - e^{\xi_{jj}}\mathbf{Z}\boldsymbol\beta}  
- \frac{\nu_\xi + 1}{1 + \nu_\xi A_\xi^2 e^{-2\xi_{jj}}} + 1, \\
&amp;\propto e^{\xi_{jj}} \boldsymbol\tau^T \mathbf{Z}^T(\mathbf{y} - \mathbf{1}_{nd})^T + e^{\xi_{jj}}\boldsymbol\tau^T\mathbf{Z}^T \left ( \frac{e^{-\mathbf{X}\boldsymbol\beta - e^{\xi_{jj}}\mathbf{Z}\boldsymbol\tau} }{1 + e^{-\mathbf{X}\boldsymbol\beta - e^{\xi_{jj}}\mathbf{Z}\boldsymbol\tau}}  \right)
- \frac{\nu_\xi + 1}{1 + \nu_\xi A_\xi^2 e^{-2\xi_{jj}}} + 1, \\
&amp;\propto e^{\xi_{jj}}\boldsymbol\tau^T \mathbf{Z}^T \left(\mathbf{y} - \mathbf{1}_{nd} +  \frac{1 }{1 + e^{\mathbf{X}\boldsymbol\beta + e^{\xi_{jj}}\mathbf{Z}\boldsymbol\tau}} \right )
- \frac{\nu_\xi + 1}{1 + \nu_\xi A_\xi^2 e^{-2\xi_{jj}}} + 1.
\end{aligned}
\]</span></p>
<p>Next, we derive the partial derivative of <span class="math inline">\(\boldsymbol\tau\)</span>,</p>
<p><span class="math display">\[
\begin{aligned}
\nabla_\tau \log f(\boldsymbol\beta, \xi, \boldsymbol\tau | \mathbf{y}, \mathbf{X}, \mathbf{Z}, \sigma_\beta^2, \nu_\xi, A_\xi)
&amp;\propto \widetilde{\mathbf{D}}^{1/2}\mathbf{Z}^T (\mathbf{y} - \mathbf{1}_{nd}) + \widetilde{\mathbf{D}}^{1/2}\mathbf{Z}^T \left(\frac{e^{-\mathbf{X}\boldsymbol\beta-\mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau}}{1 + e^{-\mathbf{X}\boldsymbol\beta-\mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau}} \right)  - \tau \\ 
&amp;\propto \widetilde{\mathbf{D}}^{1/2}\mathbf{Z}^T \left(\mathbf{y} - \mathbf{1}_{nd} + \frac{1}{1 + e^{\mathbf{X}\boldsymbol\beta + \mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau}} \right) - \tau .
\end{aligned}
\]</span></p>
<p>The gradient of full log posterior can now be specified,</p>
<p><span class="math display">\[
\begin{aligned}
\nabla_\beta \log f(\boldsymbol\beta, \xi, \boldsymbol\tau, \gamma | \mathbf{y}, \mathbf{X}, \mathbf{Z})
&amp;\propto \mathbf{X}^T \left( \mathbf{y} - \mathbf{1}_{nd} +\frac{1}{1 + e^{\mathbf{X}\boldsymbol\beta+\mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau}} \right) - \boldsymbol\beta / \sigma_\beta^2, \\
\nabla_{\xi_{jj}} \log f(\boldsymbol\beta, \xi, \boldsymbol\tau, \gamma | \mathbf{y}, \mathbf{X}, \mathbf{Z}) 
&amp;\propto e^{\xi_{jj}}\mathbf{Z}\boldsymbol\tau \left(\mathbf{y} - \mathbf{1}_{nd} +  \frac{1 }{1 + e^{\mathbf{X}\boldsymbol\beta + e^{\xi_{jj}}\mathbf{Z}\boldsymbol\tau}} \right )
- \frac{\nu_\xi + 1}{1 + \nu_\xi A_\xi^2 e^{-2\xi_{jj}}} + 1,\quad \forall (\xi_1, ..., \xi_p) \in \boldsymbol\xi \\
\nabla_\tau \log f(\boldsymbol\beta, \xi, \boldsymbol\tau, \gamma | \mathbf{y}, \mathbf{X}, \mathbf{Z}) &amp;\propto \widetilde{\mathbf{D}}^{1/2}\mathbf{Z}^T \left(\mathbf{y} - \mathbf{1}_{nd} + \frac{1}{1 + e^{\mathbf{X}\boldsymbol\beta + \mathbf{Z}\widetilde{\mathbf{D}}^{1/2}\boldsymbol\tau}} \right) - \tau.
\end{aligned}
\]</span></p>
<p>Note that a random intercept model only has a single <span class="math inline">\(\xi\)</span> parameter, which simplifies the log posterior and gradient formulations. For a random intercept model, <span class="math inline">\(\widetilde{\mathbf{D}}^{1/2} = e^{\xi}\mathbf{I}_n\)</span>.</p>
</div>
<div id="logistic-mixed-effects-model-example-data" class="section level2">
<h2 class="hasAnchor">
<a href="#logistic-mixed-effects-model-example-data" class="anchor"></a>Logistic mixed effects model 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 fixed effect design matrix <span class="math inline">\(\mathbf{X}\)</span>, random effects design matrix <span class="math inline">\(\mathbf{Z}\)</span>, and dependent variable vector <span class="math inline">\(y\)</span>.</p>
<p>We load drug Contraception data (Bates, et. al. 2014) and create the design matrices <span class="math inline">\(\mathbf{X}\)</span> and <span class="math inline">\(\mathbf{Z}\)</span> and dependent vector <span class="math inline">\(\mathbf{y}\)</span>. For this model, the random effects design matrix <span class="math inline">\(\mathbf{Z}\)</span> is specified for a random intercept model.</p>
<div class="sourceCode" id="cb2"><html><body><pre class="r"><span class="no">Contraception</span> <span class="kw">&lt;-</span> <span class="kw pkg">mlmRev</span><span class="kw ns">::</span><span class="no"><a href="https://rdrr.io/pkg/mlmRev/man/Contraception.html">Contraception</a></span>

<span class="no">Contraception</span>$<span class="no">liv2</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/ifelse.html">ifelse</a></span>(<span class="no">Contraception</span>$<span class="no">livch</span> <span class="kw">==</span> <span class="st">"0"</span>, <span class="fl">0</span>, <span class="fl">1</span>)

<span class="co">##########</span>
<span class="co"># block diagonal</span>
<span class="no">Zi.lst</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/split.html">split</a></span>(<span class="fu"><a href="https://rdrr.io/r/base/rep.html">rep</a></span>(<span class="fl">1</span>, <span class="fu"><a href="https://rdrr.io/r/base/nrow.html">nrow</a></span>(<span class="no">Contraception</span>)), <span class="no">Contraception</span>$<span class="no">district</span>)
<span class="no">Zi.lst</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/lapply.html">lapply</a></span>(<span class="no">Zi.lst</span>, <span class="no">as.matrix</span>)
<span class="no">Z</span> <span class="kw">&lt;-</span> <span class="kw pkg">Matrix</span><span class="kw ns">::</span><span class="fu"><a href="https://rdrr.io/pkg/Matrix/man/bdiag.html">bdiag</a></span>(<span class="no">Zi.lst</span>)
<span class="no">Z</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/matrix.html">as.matrix</a></span>(<span class="no">Z</span>)

<span class="no">urban</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/ifelse.html">ifelse</a></span>(<span class="no">Contraception</span>$<span class="no">urban</span> <span class="kw">==</span> <span class="st">"Y"</span>, <span class="fl">1</span>, <span class="fl">0</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="no">Contraception</span>$<span class="no">age</span>, <span class="no">Contraception</span>$<span class="no">age</span>^<span class="fl">2</span>, <span class="no">urban</span>, <span class="no">Contraception</span>$<span class="no">liv2</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">"int"</span>, <span class="st">"age"</span>, <span class="st">"age_sq"</span>, <span class="st">"urban"</span>, <span class="st">"liv2"</span>)
<span class="no">y</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/ifelse.html">ifelse</a></span>(<span class="no">Contraception</span>$<span class="no">use</span> <span class="kw">==</span> <span class="st">"Y"</span>, <span class="fl">1</span>, <span class="fl">0</span>)</pre></body></html></div>
</div>
<div id="qr-decomposition-of-design-matrix" class="section level2">
<h2 class="hasAnchor">
<a href="#qr-decomposition-of-design-matrix" class="anchor"></a>QR decomposition of design matrix</h2>
<p>To facilitate a more efficient fitting of the model, we apply QR decomposition to the fixed effects design matrix <span class="math inline">\(\mathbf{X}\)</span>.</p>
<p>Let <span class="math inline">\(\widetilde{\boldsymbol\beta} = R^*\boldsymbol\beta\)</span>, where <span class="math inline">\(\widetilde{\boldsymbol\beta}\)</span> is the transformed vector of <span class="math inline">\(\boldsymbol\beta\)</span>. The HMC estimates <span class="math inline">\(\widetilde{\boldsymbol\beta}\)</span>, from which we can use the deterministic formula to determine <span class="math inline">\(\boldsymbol\beta\)</span>.</p>
<p><span class="math display">\[
\begin{aligned}
\mathbf{X} &amp;= \mathbf{Q}^* \mathbf{R}^*, \quad \mathbf{Q}^* = \mathbf{Q} \cdot \sqrt{nd-1}, \quad \mathbf{R}^* = \frac{1}{\sqrt{nd-1}}\mathbf{R}, \\
\mathbf{X}\boldsymbol\beta &amp;= \mathbf{Q}^* \mathbf{R}^* \boldsymbol\beta \\
\boldsymbol\beta &amp;= \mathbf{R}^{*^{-1}}\widetilde{\boldsymbol\beta}, \quad \widetilde{\boldsymbol\beta} = \mathbf{R}^*\boldsymbol\beta
\end{aligned}
\]</span></p>
<p>We use the <em>qr</em> function from base <strong>R</strong> for our parameter transformation,</p>
<div class="sourceCode" id="cb3"><html><body><pre class="r"><span class="no">xqr</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/qr.html">qr</a></span>(<span class="no">X</span>)
<span class="no">Q</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/qraux.html">qr.Q</a></span>(<span class="no">xqr</span>)
<span class="no">R</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/qraux.html">qr.R</a></span>(<span class="no">xqr</span>)

<span class="no">n</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/nrow.html">nrow</a></span>(<span class="no">X</span>)
<span class="no">X2</span> <span class="kw">&lt;-</span> <span class="no">Q</span> * <span class="fu"><a href="https://rdrr.io/r/base/MathFun.html">sqrt</a></span>(<span class="no">n</span>-<span class="fl">1</span>)
<span class="no">Rstar</span> <span class="kw">&lt;-</span> <span class="no">R</span> / <span class="fu"><a href="https://rdrr.io/r/base/MathFun.html">sqrt</a></span>(<span class="no">n</span>-<span class="fl">1</span>)
<span class="no">Rstar_inv</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/solve.html">solve</a></span>(<span class="no">Rstar</span>)
<span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">X2</span>) <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="st">"int"</span>, <span class="st">"age"</span>, <span class="st">"age_sq"</span>, <span class="st">"urban"</span>, <span class="st">"liv2"</span>)</pre></body></html></div>
<p>We could also apply QR decomposition to <span class="math inline">\(\matbf{Z}\)</span>, but this is not necessary for a random intercept model where all elements of <span class="math inline">\(\mathbf{Z}\)</span> are 0 or 1.</p>
</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 linear mixed effects regression model using HMC. A vector of <em>tuning parameter</em> <span class="math inline">\(\epsilon\)</span> values are specified to align with the data. The hyperparameters for <span class="math inline">\(\nu_\xi\)</span> and <span class="math inline">\(A_\xi\)</span> are set to defaults in <strong>hmclearn</strong> (Gelman 2006). The hyperprior <span class="math inline">\(\sigma_\beta^2\)</span> is set lower than the default based on the dependent variable variance for this hierarchical model.</p>
<div class="sourceCode" id="cb4"><html><body><pre class="r"><span class="no">N</span> <span class="kw">&lt;-</span> <span class="fl">2e3</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">initvals</span><span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="fu"><a href="https://rdrr.io/r/base/rep.html">rep</a></span>(<span class="fl">0</span>, <span class="fl">5</span>), <span class="co"># fixed effects</span>
                <span class="fu"><a href="https://rdrr.io/r/base/rep.html">rep</a></span>(<span class="fl">0</span>, <span class="fl">60</span>), <span class="co"># random intercepts</span>
                <span class="fl">0</span>) <span class="co"># variance of random intercepts</span>

<span class="no">vnames</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></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="https://rdrr.io/r/base/paste.html">paste0</a></span>(<span class="st">"tau_int"</span>, <span class="fl">1</span>:<span class="fl">60</span>),
            <span class="st">"xi"</span>)

<span class="no">epsvals</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="fl">5e-2</span>, <span class="fu"><a href="https://rdrr.io/r/base/rep.html">rep</a></span>(<span class="fl">1e-2</span>, <span class="fl">4</span>), <span class="fu"><a href="https://rdrr.io/r/base/rep.html">rep</a></span>(<span class="fl">5e-2</span>, <span class="fl">61</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="no">initvals</span>,
             <span class="kw">epsilon</span> <span class="kw">=</span> <span class="no">epsvals</span>, <span class="kw">L</span> <span class="kw">=</span> <span class="fl">10</span>,
             <span class="kw">logPOSTERIOR</span> <span class="kw">=</span> <span class="no">glmm_bin_posterior</span>,
             <span class="kw">glogPOSTERIOR</span> <span class="kw">=</span> <span class="no">g_glmm_bin_posterior</span>,
             <span class="kw">varnames</span> <span class="kw">=</span> <span class="no">vnames</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="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">X2</span>, <span class="kw">Z</span><span class="kw">=</span><span class="no">Z</span>, <span class="kw">n</span><span class="kw">=</span><span class="fl">60</span>, <span class="kw">nrandom</span><span class="kw">=</span><span class="fl">1</span>, <span class="kw">sig2beta</span><span class="kw">=</span><span class="fl">5</span>,
                        <span class="kw">nuxi</span><span class="kw">=</span><span class="fl">1</span>, <span class="kw">Axi</span><span class="kw">=</span><span class="fl">25</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>()
<span class="no">t2.hmc</span> - <span class="no">t1.hmc</span>
<span class="co">#&gt; Time difference of 7.115439 mins</span></pre></body></html></div>
<p>The acceptance ratio for each of the HMC chains is sufficiently high for an efficient simulation.</p>
<div class="sourceCode" id="cb5"><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.6640 0.6595</span></pre></body></html></div>
<p>Since we used QR decomposition to transform <span class="math inline">\(\beta\)</span> prior to fitting via HMC, we need to reverse the transformation to obtain the original parameter scale.</p>
<div class="sourceCode" id="cb6"><html><body><pre class="r"><span class="co"># restore beta from Rstar in QR decomposition</span>
<span class="no">calc_beta</span> <span class="kw">&lt;-</span> <span class="kw">function</span>(<span class="no">beta_tilde_param</span>, <span class="no">Rstarinv</span>) {
  <span class="fu"><a href="https://rdrr.io/r/base/numeric.html">as.numeric</a></span>(<span class="no">Rstarinv</span> <span class="kw">%*%</span> <span class="no">beta_tilde_param</span>)
}

<span class="co"># reverse qr decomposition</span>
<span class="no">f_hmc2</span> <span class="kw">&lt;-</span> <span class="no">f_hmc</span>
<span class="no">f_hmc2</span>$<span class="no">thetaCombined</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/lapply.html">lapply</a></span>(<span class="no">f_hmc</span>$<span class="no">thetaCombined</span>, <span class="kw">function</span>(<span class="no">xx</span>) {
  <span class="no">xx</span>[, <span class="fl">1</span>:<span class="fl">5</span>] <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/t.html">t</a></span>(<span class="fu"><a href="https://rdrr.io/r/base/apply.html">apply</a></span>(<span class="no">xx</span>[, <span class="fl">1</span>:<span class="fl">5</span>], <span class="fl">1</span>, <span class="no">calc_beta</span>, <span class="kw">Rstarinv</span><span class="kw">=</span><span class="no">Rstar_inv</span>))
  <span class="no">xx</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 are close to one, indicating that both HMC chains converged to the same distribution. The <span class="math inline">\(\hat{R}\)</span> statistics provide an indication of convergence. Values close to one indicate that the multiple MCMC chains converged 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="cb7"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span>(<span class="no">f_hmc2</span>, <span class="kw">burnin</span><span class="kw">=</span><span class="fl">1000</span>)
<span class="co">#&gt; Summary of MCMC simulation</span>
<span class="co">#&gt;                   2.5%           5%          25%          50%          75%</span>
<span class="co">#&gt; int       -1.339765267 -1.292385942 -1.128009491 -1.019416922 -0.907314990</span>
<span class="co">#&gt; age       -0.009673678 -0.007411716  0.000777096  0.006236717  0.011598903</span>
<span class="co">#&gt; age_sq    -0.006298813 -0.005942879 -0.005188464 -0.004660974 -0.004154788</span>
<span class="co">#&gt; urban      0.461242801  0.490947335  0.604555712  0.690351383  0.771730252</span>
<span class="co">#&gt; liv2       0.565051987  0.617225111  0.762088447  0.861900584  0.966908404</span>
<span class="co">#&gt; tau_int1  -2.477301164 -2.295846320 -1.831002643 -1.485545520 -1.164350389</span>
<span class="co">#&gt; tau_int2  -1.251975957 -1.050339233 -0.403339467  0.012510607  0.444257310</span>
<span class="co">#&gt; tau_int3  -1.496266849 -1.125830130 -0.181008471  0.512069690  1.099865165</span>
<span class="co">#&gt; tau_int4  -0.966503731 -0.712255547 -0.068020500  0.374592021  0.822052629</span>
<span class="co">#&gt; tau_int5  -0.827741170 -0.641609874 -0.144274565  0.245582568  0.573794386</span>
<span class="co">#&gt; tau_int6  -1.478551092 -1.351362204 -0.786215231 -0.480988538 -0.144205377</span>
<span class="co">#&gt; tau_int7  -1.737657080 -1.527971313 -0.825051831 -0.329389586  0.168319141</span>
<span class="co">#&gt; tau_int8  -0.997424538 -0.777771580 -0.133135491  0.297942975  0.724042032</span>
<span class="co">#&gt; tau_int9  -1.643903931 -1.464987371 -0.762294182 -0.337359661  0.034850365</span>
<span class="co">#&gt; tau_int10 -2.359678741 -2.105265779 -1.327151722 -0.715741817 -0.161124706</span>
<span class="co">#&gt; tau_int11 -2.975526923 -2.717381167 -2.051897535 -1.536501259 -1.007861538</span>
<span class="co">#&gt; tau_int12 -1.398553510 -1.196599051 -0.616111099 -0.207894822  0.168799803</span>
<span class="co">#&gt; tau_int13 -1.032358207 -0.844696528 -0.228153558  0.188442621  0.631796254</span>
<span class="co">#&gt; tau_int14  0.559458337  0.721712505  1.069203730  1.342730088  1.623509380</span>
<span class="co">#&gt; tau_int15 -1.365090935 -1.145741575 -0.512050262 -0.067173547  0.368870520</span>
<span class="co">#&gt; tau_int16 -0.064201495  0.135029841  0.829913139  1.295959783  1.734034213</span>
<span class="co">#&gt; tau_int17 -1.592079894 -1.400937404 -0.744838461 -0.321090887  0.120132412</span>
<span class="co">#&gt; tau_int18 -1.153543773 -0.974146585 -0.482445878 -0.089331621  0.238488623</span>
<span class="co">#&gt; tau_int19 -1.303914398 -1.054172616 -0.314905282  0.136838466  0.534635479</span>
<span class="co">#&gt; tau_int20 -1.067261225 -0.829041944  0.009314237  0.594828736  1.108766609</span>
<span class="co">#&gt; tau_int21 -1.309736941 -1.120473575 -0.447187052  0.062560124  0.460796795</span>
<span class="co">#&gt; tau_int22 -2.408175260 -2.126023008 -1.308798064 -0.759845544 -0.292215780</span>
<span class="co">#&gt; tau_int23 -1.749373324 -1.446025651 -0.828086758 -0.349566893  0.113698516</span>
<span class="co">#&gt; tau_int24 -2.329122632 -2.141378285 -1.478914611 -1.000907613 -0.490455706</span>
<span class="co">#&gt; tau_int25 -0.540036541 -0.422645376  0.016679088  0.304140922  0.627940994</span>
<span class="co">#&gt; tau_int26 -1.351668511 -1.142046677 -0.440593643  0.061109897  0.545336103</span>
<span class="co">#&gt; tau_int27 -2.145374899 -1.954037586 -1.344605341 -0.988112387 -0.607366834</span>
<span class="co">#&gt; tau_int28 -1.820291883 -1.626725227 -1.015316805 -0.645003049 -0.280632825</span>
<span class="co">#&gt; tau_int29 -1.491347848 -1.216564131 -0.541549508 -0.162957885  0.195636426</span>
<span class="co">#&gt; tau_int30 -0.047293333  0.106622876  0.582074062  0.911978092  1.219722545</span>
<span class="co">#&gt; tau_int31 -0.321135200 -0.109160266  0.366222564  0.776824770  1.185801478</span>
<span class="co">#&gt; tau_int32 -2.254218808 -2.134555625 -1.487004438 -0.974969538 -0.506001715</span>
<span class="co">#&gt; tau_int33 -1.635821042 -1.372336868 -0.615861808 -0.104346197  0.406080821</span>
<span class="co">#&gt; tau_int34  0.254477324  0.394040586  1.008690765  1.411154022  1.775017112</span>
<span class="co">#&gt; tau_int35 -0.490136768 -0.335186646  0.227358260  0.573131204  0.895082541</span>
<span class="co">#&gt; tau_int36 -1.661440235 -1.450058926 -0.771555053 -0.254392894  0.293770981</span>
<span class="co">#&gt; tau_int37 -0.904188684 -0.493811623  0.328799567  0.839476446  1.338323458</span>
<span class="co">#&gt; tau_int38 -1.669368277 -1.351369655 -0.608118957 -0.136361127  0.362715112</span>
<span class="co">#&gt; tau_int39 -0.352821461 -0.073515331  0.489814711  0.966400114  1.413938187</span>
<span class="co">#&gt; tau_int40 -1.004600501 -0.779743166 -0.222373162  0.164659958  0.581337584</span>
<span class="co">#&gt; tau_int41 -0.685159504 -0.529584717  0.084617820  0.537188788  0.967631007</span>
<span class="co">#&gt; tau_int42 -1.087748990 -0.862336968 -0.171009371  0.359401202  0.872360584</span>
<span class="co">#&gt; tau_int43 -0.188416379  0.006353872  0.583026811  0.973056859  1.336939929</span>
<span class="co">#&gt; tau_int44 -1.920893963 -1.775253474 -1.075238768 -0.612325712 -0.204538828</span>
<span class="co">#&gt; tau_int45 -1.375816390 -1.214593860 -0.661396619 -0.247330696  0.106653253</span>
<span class="co">#&gt; tau_int46  0.260955774  0.385831544  0.828625933  1.098098656  1.413419865</span>
<span class="co">#&gt; tau_int47 -1.065313590 -0.846629294 -0.074804657  0.417714248  1.010515010</span>
<span class="co">#&gt; tau_int48 -0.154266238  0.032249004  0.495837288  0.904583872  1.287996119</span>
<span class="co">#&gt; tau_int49 -2.292638576 -2.074924835 -1.120632986 -0.465344095  0.254080874</span>
<span class="co">#&gt; tau_int50 -0.694975585 -0.502580514  0.147183021  0.593669294  1.049450921</span>
<span class="co">#&gt; tau_int51 -0.850862185 -0.706340487 -0.149812867  0.231557500  0.603255298</span>
<span class="co">#&gt; tau_int52 -0.400872580 -0.234890160  0.218027332  0.504078904  0.838352856</span>
<span class="co">#&gt; tau_int53 -1.865292382 -1.490577454 -0.792623129 -0.353175268  0.106754012</span>
<span class="co">#&gt; tau_int54 -2.637911673 -2.331394595 -1.440817339 -0.796543076 -0.119464717</span>
<span class="co">#&gt; tau_int55  0.234680515  0.412128967  0.895952077  1.248927955  1.631155334</span>
<span class="co">#&gt; tau_int56 -2.367275531 -2.087173030 -1.505090725 -1.040664297 -0.649099371</span>
<span class="co">#&gt; tau_int57 -0.642412799 -0.473892204  0.199193278  0.629168374  1.010768769</span>
<span class="co">#&gt; tau_int58 -2.479820044 -2.213332561 -1.457421029 -0.844450385 -0.306002122</span>
<span class="co">#&gt; tau_int59 -2.075535059 -1.837203938 -1.237516378 -0.810022666 -0.360490885</span>
<span class="co">#&gt; tau_int60 -2.132356364 -1.975811873 -1.407047094 -1.022341822 -0.687099338</span>
<span class="co">#&gt; xi        -1.025640064 -0.972884948 -0.806205309 -0.688589885 -0.579083896</span>
<span class="co">#&gt;                    95%        97.5%      rhat</span>
<span class="co">#&gt; int       -0.730945734 -0.694148399 1.0004496</span>
<span class="co">#&gt; age        0.019824477  0.021691001 1.0001230</span>
<span class="co">#&gt; age_sq    -0.003448999 -0.003151154 0.9996103</span>
<span class="co">#&gt; urban      0.892366455  0.925163643 0.9999805</span>
<span class="co">#&gt; liv2       1.100691135  1.151162473 0.9995266</span>
<span class="co">#&gt; tau_int1  -0.799589589 -0.669744207 1.0002275</span>
<span class="co">#&gt; tau_int2   1.181023361  1.348978663 1.0123033</span>
<span class="co">#&gt; tau_int3   1.796161785  2.090884947 1.0018321</span>
<span class="co">#&gt; tau_int4   1.399780513  1.618444831 1.0059114</span>
<span class="co">#&gt; tau_int5   1.094726502  1.252372909 1.0038215</span>
<span class="co">#&gt; tau_int6   0.296673780  0.425040708 1.0186125</span>
<span class="co">#&gt; tau_int7   0.843221047  1.022866044 1.0003223</span>
<span class="co">#&gt; tau_int8   1.307405711  1.456389568 1.0117206</span>
<span class="co">#&gt; tau_int9   0.648364858  0.867039474 1.0249419</span>
<span class="co">#&gt; tau_int10  0.711639294  0.997120565 1.0042454</span>
<span class="co">#&gt; tau_int11 -0.305128022 -0.044896361 1.0054944</span>
<span class="co">#&gt; tau_int12  0.759889717  0.954494483 1.0028436</span>
<span class="co">#&gt; tau_int13  1.280758654  1.460665846 1.0051275</span>
<span class="co">#&gt; tau_int14  2.146233237  2.343043362 1.0003173</span>
<span class="co">#&gt; tau_int15  1.042446268  1.215022601 0.9999793</span>
<span class="co">#&gt; tau_int16  2.375195988  2.607269118 1.0007880</span>
<span class="co">#&gt; tau_int17  0.714932777  0.931174987 1.0045396</span>
<span class="co">#&gt; tau_int18  0.718577986  0.934903492 1.0047520</span>
<span class="co">#&gt; tau_int19  1.191775401  1.393233342 1.0079420</span>
<span class="co">#&gt; tau_int20  1.809275450  1.926635301 0.9998867</span>
<span class="co">#&gt; tau_int21  1.168546901  1.544283141 1.0006834</span>
<span class="co">#&gt; tau_int22  0.421293857  0.640374363 1.0000198</span>
<span class="co">#&gt; tau_int23  1.120706910  1.399107469 1.0145157</span>
<span class="co">#&gt; tau_int24  0.195017088  0.497331539 1.0028844</span>
<span class="co">#&gt; tau_int25  1.158575285  1.285652174 1.0133637</span>
<span class="co">#&gt; tau_int26  1.356310609  1.594988909 1.0060216</span>
<span class="co">#&gt; tau_int27  0.005555420  0.161674130 0.9995001</span>
<span class="co">#&gt; tau_int28  0.229115852  0.446685154 1.0009935</span>
<span class="co">#&gt; tau_int29  0.789414127  0.977283620 0.9998314</span>
<span class="co">#&gt; tau_int30  1.689971972  1.829473728 0.9998212</span>
<span class="co">#&gt; tau_int31  1.738535476  1.963106919 1.0003311</span>
<span class="co">#&gt; tau_int32  0.115847158  0.356481074 0.9997730</span>
<span class="co">#&gt; tau_int33  1.105972373  1.326058437 0.9999655</span>
<span class="co">#&gt; tau_int34  2.366490502  2.518571484 1.0060984</span>
<span class="co">#&gt; tau_int35  1.485976781  1.731665495 0.9997770</span>
<span class="co">#&gt; tau_int36  0.936032250  1.076176241 1.0030959</span>
<span class="co">#&gt; tau_int37  2.008294757  2.343421693 1.0002054</span>
<span class="co">#&gt; tau_int38  1.043362204  1.148039747 1.0104009</span>
<span class="co">#&gt; tau_int39  2.115662668  2.260788790 0.9995159</span>
<span class="co">#&gt; tau_int40  1.129694872  1.369096048 0.9998209</span>
<span class="co">#&gt; tau_int41  1.539274586  1.778110612 1.0015194</span>
<span class="co">#&gt; tau_int42  1.586178401  1.850719607 1.0166520</span>
<span class="co">#&gt; tau_int43  1.962268814  2.230524064 1.0081014</span>
<span class="co">#&gt; tau_int44  0.351470732  0.494781239 0.9997364</span>
<span class="co">#&gt; tau_int45  0.667483591  0.837574465 1.0013793</span>
<span class="co">#&gt; tau_int46  1.810090021  1.984648247 1.0044925</span>
<span class="co">#&gt; tau_int47  1.795064682  2.063204056 0.9995069</span>
<span class="co">#&gt; tau_int48  1.849395364  2.066445196 1.0057786</span>
<span class="co">#&gt; tau_int49  1.040068288  1.252388248 1.0065090</span>
<span class="co">#&gt; tau_int50  1.788166907  2.027244709 0.9997836</span>
<span class="co">#&gt; tau_int51  1.133122335  1.300567668 1.0091042</span>
<span class="co">#&gt; tau_int52  1.301499630  1.480462042 0.9999622</span>
<span class="co">#&gt; tau_int53  0.723585762  0.953061185 0.9995559</span>
<span class="co">#&gt; tau_int54  0.659315950  0.833145130 1.0034914</span>
<span class="co">#&gt; tau_int55  2.209654121  2.400583051 1.0001205</span>
<span class="co">#&gt; tau_int56 -0.032734459  0.152406047 1.0268956</span>
<span class="co">#&gt; tau_int57  1.648366407  1.821616636 1.0042383</span>
<span class="co">#&gt; tau_int58  0.512621248  0.721277596 1.0049889</span>
<span class="co">#&gt; tau_int59  0.209610312  0.419462276 0.9998181</span>
<span class="co">#&gt; tau_int60 -0.100502679  0.050438431 0.9995032</span>
<span class="co">#&gt; xi        -0.410063707 -0.369595618 1.0031427</span></pre></body></html></div>
<p>We create trace plots on the transformed simulation data.</p>
<div class="sourceCode" id="cb8"><html><body><pre class="r"><span class="fu"><a href="../reference/hmclearn-plots.html">mcmc_trace</a></span>(<span class="no">f_hmc2</span>, <span class="kw">burnin</span><span class="kw">=</span><span class="fl">1000</span>, <span class="kw">pars</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>))</pre></body></html></div>
<p><img src="logistic_mixed_effects_hmclearn_files/figure-html/unnamed-chunk-8-1.png" width="576"></p>
</div>
<div id="comparison-model---frequentist" class="section level2">
<h2 class="hasAnchor">
<a href="#comparison-model---frequentist" class="anchor"></a>Comparison model - Frequentist</h2>
<p>To compare results, we first fit a logistic mixed effects model using the frequentist package <strong>lme4</strong> (Bates et. al. 2015).</p>
<div class="sourceCode" id="cb9"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/r/base/library.html">library</a></span>(<span class="no">lme4</span>)
<span class="co">#&gt; Loading required package: Matrix</span>

<span class="no">fm1</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/pkg/lme4/man/glmer.html">glmer</a></span>(<span class="no">use</span> ~ <span class="no">age</span> + <span class="fu"><a href="https://rdrr.io/r/base/AsIs.html">I</a></span>(<span class="no">age</span>^<span class="fl">2</span>) + <span class="no">urban</span> + <span class="no">liv2</span> + (<span class="fl">1</span> <span class="kw">|</span> <span class="no">district</span>),
           <span class="kw">data</span><span class="kw">=</span><span class="no">Contraception</span>, <span class="kw">family</span><span class="kw">=</span><span class="no">binomial</span>,
           <span class="kw">control</span><span class="kw">=</span><span class="fu"><a href="https://rdrr.io/pkg/lme4/man/lmerControl.html">glmerControl</a></span>(<span class="kw">optCtrl</span><span class="kw">=</span><span class="fu"><a href="https://rdrr.io/r/base/list.html">list</a></span>(<span class="kw">maxfun</span><span class="kw">=</span><span class="fl">20000</span>)))
<span class="co">#&gt; Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue</span>
<span class="co">#&gt;  - Rescale variables?</span>
<span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span>(<span class="no">fm1</span>)
<span class="co">#&gt; Generalized linear mixed model fit by maximum likelihood (Laplace</span>
<span class="co">#&gt;   Approximation) [glmerMod]</span>
<span class="co">#&gt;  Family: binomial  ( logit )</span>
<span class="co">#&gt; Formula: use ~ age + I(age^2) + urban + liv2 + (1 | district)</span>
<span class="co">#&gt;    Data: Contraception</span>
<span class="co">#&gt; Control: glmerControl(optCtrl = list(maxfun = 20000))</span>
<span class="co">#&gt; </span>
<span class="co">#&gt;      AIC      BIC   logLik deviance df.resid </span>
<span class="co">#&gt;   2385.2   2418.6  -1186.6   2373.2     1928 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; Scaled residuals: </span>
<span class="co">#&gt;     Min      1Q  Median      3Q     Max </span>
<span class="co">#&gt; -1.8151 -0.7620 -0.4619  0.9518  3.1033 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; Random effects:</span>
<span class="co">#&gt;  Groups   Name        Variance Std.Dev.</span>
<span class="co">#&gt;  district (Intercept) 0.2247   0.474   </span>
<span class="co">#&gt; Number of obs: 1934, groups:  district, 60</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; Fixed effects:</span>
<span class="co">#&gt;               Estimate Std. Error z value Pr(&gt;|z|)    </span>
<span class="co">#&gt; (Intercept) -1.0063737  0.1691107  -5.951 2.67e-09 ***</span>
<span class="co">#&gt; age          0.0062561  0.0078848   0.793    0.428    </span>
<span class="co">#&gt; I(age^2)    -0.0046353  0.0007207  -6.432 1.26e-10 ***</span>
<span class="co">#&gt; urbanY       0.6929220  0.1206566   5.743 9.31e-09 ***</span>
<span class="co">#&gt; liv2         0.8603821  0.1483014   5.802 6.57e-09 ***</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; Correlation of Fixed Effects:</span>
<span class="co">#&gt;          (Intr) age    I(g^2) urbanY</span>
<span class="co">#&gt; age       0.529                     </span>
<span class="co">#&gt; I(age^2) -0.533 -0.506              </span>
<span class="co">#&gt; urbanY   -0.259 -0.032  0.020       </span>
<span class="co">#&gt; liv2     -0.801 -0.565  0.345  0.094</span>
<span class="co">#&gt; convergence code: 0</span>
<span class="co">#&gt; Model is nearly unidentifiable: very large eigenvalue</span>
<span class="co">#&gt;  - Rescale variables?</span></pre></body></html></div>
<div class="sourceCode" id="cb10"><html><body><pre class="r">
<span class="no">freqvals</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="fu"><a href="https://rdrr.io/r/base/numeric.html">as.numeric</a></span>(<span class="fu"><a href="https://rdrr.io/pkg/lme4/man/fixef.html">fixef</a></span>(<span class="no">fm1</span>)),
              <span class="fu"><a href="https://rdrr.io/r/base/numeric.html">as.numeric</a></span>(<span class="fu"><a href="https://rdrr.io/pkg/lme4/man/ranef.html">ranef</a></span>(<span class="no">fm1</span>)$<span class="no">district</span>[, <span class="fl">1</span>]),
              <span class="fu"><a href="https://rdrr.io/r/base/Log.html">log</a></span>(<span class="fu"><a href="https://rdrr.io/r/base/MathFun.html">sqrt</a></span>(<span class="fu"><a href="https://rdrr.io/r/base/numeric.html">as.numeric</a></span>(<span class="fu"><a href="https://rdrr.io/pkg/lme4/man/VarCorr.html">VarCorr</a></span>(<span class="no">fm1</span>)[<span class="fl">1</span>]))))</pre></body></html></div>
<p>Histograms of the posterior distribution show that Bayesian parameter estimates align with frequentist estimates. The <em>cols</em> parameter specifies the parameters to be displayed in <em>diagplots</em>, based on the order provided to the <em>hmc</em> function.</p>
<div class="sourceCode" id="cb11"><html><body><pre class="r"><span class="fu"><a href="../reference/diagplots.html">diagplots</a></span>(<span class="no">f_hmc2</span>, <span class="kw">burnin</span><span class="kw">=</span><span class="fl">1000</span>,
          <span class="kw">comparison.theta</span><span class="kw">=</span><span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="no">freqvals</span>[<span class="fl">1</span>:<span class="fu"><a href="https://rdrr.io/r/base/nrow.html">ncol</a></span>(<span class="no">X</span>)],
                             <span class="no">freqvals</span>[<span class="fu"><a href="https://rdrr.io/r/base/length.html">length</a></span>(<span class="no">freqvals</span>)]),
          <span class="kw">cols</span><span class="kw">=</span><span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="fl">1</span>:<span class="fu"><a href="https://rdrr.io/r/base/nrow.html">ncol</a></span>(<span class="no">X</span>), <span class="fu"><a href="https://rdrr.io/r/base/length.html">length</a></span>(<span class="no">initvals</span>)))
<span class="co">#&gt; $histogram</span></pre></body></html></div>
<p><img src="logistic_mixed_effects_hmclearn_files/figure-html/unnamed-chunk-11-1.png" width="576"></p>
<p>We also compare the random effects parameter estimates with <strong>lme4</strong>. We apply the linear transformtion back to <span class="math inline">\(\mathbf{u}\)</span> for comparison.</p>
<div class="sourceCode" id="cb12"><html><body><pre class="r"><span class="no">u.freq</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/pkg/lme4/man/ranef.html">ranef</a></span>(<span class="no">fm1</span>)$<span class="no">district</span>[, <span class="fl">1</span>]
<span class="no">lambda.freq</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/MathFun.html">sqrt</a></span>(<span class="fu"><a href="https://rdrr.io/pkg/lme4/man/VarCorr.html">VarCorr</a></span>(<span class="no">fm1</span>)$<span class="no">district</span>[<span class="fl">1</span>])

<span class="co"># transform parameters back to original scale</span>
<span class="no">f_hmc</span>$<span class="no">thetaCombined</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/lapply.html">lapply</a></span>(<span class="no">f_hmc</span>$<span class="no">thetaCombined</span>, <span class="kw">function</span>(<span class="no">xx</span>) {
  <span class="no">tau_mx</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/matrix.html">as.matrix</a></span>(<span class="no">xx</span>[, <span class="fu"><a href="https://rdrr.io/r/base/grep.html">grepl</a></span>(<span class="st">"tau"</span>, <span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">xx</span>))])
  <span class="no">u_mx</span> <span class="kw">&lt;-</span> <span class="no">tau_mx</span> * <span class="fu"><a href="https://rdrr.io/r/base/Log.html">exp</a></span>(<span class="no">xx</span>[, <span class="st">"xi"</span>])
  <span class="no">u_df</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/as.data.frame.html">as.data.frame</a></span>(<span class="no">u_mx</span>)
  <span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">u_df</span>) <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/paste.html">paste0</a></span>(<span class="st">"u"</span>, <span class="fl">1</span>:<span class="fu"><a href="https://rdrr.io/r/base/nrow.html">ncol</a></span>(<span class="no">u_df</span>))
  <span class="no">xx</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/cbind.html">cbind</a></span>(<span class="no">xx</span>, <span class="no">u_df</span>, <span class="fu"><a href="https://rdrr.io/r/base/Log.html">exp</a></span>(<span class="no">xx</span>[, <span class="st">"xi"</span>]))
  <span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">xx</span>)[<span class="fu"><a href="https://rdrr.io/r/base/nrow.html">ncol</a></span>(<span class="no">xx</span>)] <span class="kw">&lt;-</span> <span class="st">"lambda"</span>
  <span class="no">xx</span>
})</pre></body></html></div>
<p>Since we have 60 random effect parameters plus <span class="math inline">\(\lambda\)</span>, we split the plots into more manageable chunks. The random effects parameters align with frequentist estimates.</p>
<div class="sourceCode" id="cb13"><html><body><pre class="r"><span class="co"># histograms with lines for frequentist estimates</span>
<span class="co"># ucols &lt;- which(grepl("^[u][0-9]{1,2}", colnames(f_hmc$thetaCombined[[1]])))</span>
<span class="no">ucols</span> <span class="kw">&lt;-</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/grep.html">grepl</a></span>(<span class="st">"^[u](?:[1-9]|0[1-9]|1[0-9]|20)$"</span>,
                     <span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">f_hmc</span>$<span class="no">thetaCombined</span><span class="kw">[[</span><span class="fl">1</span>]])))
<span class="no">lambdacol</span> <span class="kw">&lt;-</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/grep.html">grepl</a></span>(<span class="st">"^lambda"</span>, <span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">f_hmc</span>$<span class="no">thetaCombined</span><span class="kw">[[</span><span class="fl">1</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">1000</span>,
          <span class="kw">comparison.theta</span> <span class="kw">=</span> <span class="no">u.freq</span>[<span class="fl">1</span>:<span class="fl">20</span>], <span class="kw">cols</span> <span class="kw">=</span> <span class="no">ucols</span>)
<span class="co">#&gt; $histogram</span></pre></body></html></div>
<p><img src="logistic_mixed_effects_hmclearn_files/figure-html/unnamed-chunk-13-1.png" width="768"></p>
<div class="sourceCode" id="cb14"><html><body><pre class="r"><span class="co"># histograms with lines for frequentist estimates</span>
<span class="co"># ucols &lt;- which(grepl("^[u][0-9]{1,2}", colnames(f_hmc$thetaCombined[[1]])))</span>
<span class="no">ucols</span> <span class="kw">&lt;-</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/grep.html">grepl</a></span>(<span class="st">"^[u](?:2[1-9]|3[0-9]|40)$"</span>,
                     <span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">f_hmc</span>$<span class="no">thetaCombined</span><span class="kw">[[</span><span class="fl">1</span>]])))
<span class="no">lambdacol</span> <span class="kw">&lt;-</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/grep.html">grepl</a></span>(<span class="st">"^lambda"</span>, <span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">f_hmc</span>$<span class="no">thetaCombined</span><span class="kw">[[</span><span class="fl">1</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">1000</span>,
          <span class="kw">comparison.theta</span> <span class="kw">=</span> <span class="no">u.freq</span>[<span class="fl">21</span>:<span class="fl">40</span>], <span class="kw">cols</span> <span class="kw">=</span> <span class="no">ucols</span>)
<span class="co">#&gt; $histogram</span></pre></body></html></div>
<p><img src="logistic_mixed_effects_hmclearn_files/figure-html/unnamed-chunk-14-1.png" width="768"></p>
<div class="sourceCode" id="cb15"><html><body><pre class="r"><span class="co"># histograms with lines for frequentist estimates</span>
<span class="co"># ucols &lt;- which(grepl("^[u][0-9]{1,2}", colnames(f_hmc$thetaCombined[[1]])))</span>
<span class="no">ucols</span> <span class="kw">&lt;-</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/grep.html">grepl</a></span>(<span class="st">"^[u](?:4[1-9]|5[0-9]|60)$"</span>,
                     <span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">f_hmc</span>$<span class="no">thetaCombined</span><span class="kw">[[</span><span class="fl">1</span>]])))
<span class="no">lambdacol</span> <span class="kw">&lt;-</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/grep.html">grepl</a></span>(<span class="st">"^lambda"</span>, <span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">f_hmc</span>$<span class="no">thetaCombined</span><span class="kw">[[</span><span class="fl">1</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">1000</span>,
          <span class="kw">comparison.theta</span> <span class="kw">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="no">u.freq</span>[<span class="fl">41</span>:<span class="fl">60</span>], <span class="no">lambda.freq</span>), <span class="kw">cols</span> <span class="kw">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="no">ucols</span>, <span class="no">lambdacol</span>))
<span class="co">#&gt; $histogram</span></pre></body></html></div>
<p><img src="logistic_mixed_effects_hmclearn_files/figure-html/unnamed-chunk-15-1.png" width="768"></p>
</div>
<div id="source" class="section level2">
<h2 class="hasAnchor">
<a href="#source" class="anchor"></a>Source</h2>
<p>Steele, F., Diamond, I. And Amin, S. (1996). Immunization uptake in rural Bangladesh: a multilevel analysis. <em>Journal of the Royal Statistical Society</em>, Series A (159): 289-299.</p>
<p></p>
</div>
<div id="references" class="section level2">
<h2 class="hasAnchor">
<a href="#references" class="anchor"></a>References</h2>
<p>Bates, D., M"{a}chler, M., Bolker, B., &amp; Walker, S. (2015). Fitting linear mixed-effects models using lme4. <em>Journal of Statistical Software</em> 67(1)</p>
<p>Bates, D., M"{a}chler, M., &amp; Bolker, B. (2014). mlmRev: Examples from multilevel modelling software review. <em>R package version, 1</em>.</p>
<p>Agresti, A. (2015). <em>Foundations of linear and generalized linear models</em>. John Wiley &amp; Sons. ISBN: 978-1-118-73003-4</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>
