<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/bayesian-nonparametric-models-in-nimble-part-2-nonparametric-random-effects/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://r-nimble.org/bayesian-nonparametric-models-in-nimble-part-2-nonparametric-random-effects"> R – NIMBLE</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>

<p></p>
<h1>Bayesian nonparametrics in NIMBLE: Nonparametric random effects</h1>
<h2>Overview</h2>
<p>NIMBLE is a hierarchical modeling package that uses nearly the same language for model specification as the popular MCMC packages WinBUGS, OpenBUGS and JAGS, while making the modeling language extensible — you can add distributions and functions — and also allowing customization of the algorithms used to estimate the parameters of the model.</p>
<p>Recently, we added support for Markov chain Monte Carlo (MCMC) inference for Bayesian nonparametric (BNP) mixture models to NIMBLE. In particular, starting with version 0.6-11, NIMBLE provides functionality for fitting models involving Dirichlet process priors using either the Chinese Restaurant Process (CRP) or a truncated stick-breaking (SB) representation of the Dirichlet process prior.</p>
<p>We will illustrate NIMBLE’s BNP capabilities using two examples. In a previous post, we showed how to use nonparametric mixture models with different kernels for density estimation. In this post, we will take a parametric generalized linear mixed model and show how to switch to a nonparametric representation of the random effects that avoids the assumption of normally-distributed random effects.</p>
<p>For more detailed information on NIMBLE and Bayesian nonparametrics in NIMBLE, see the <a rel="nofollow" href="https://r-nimble.org/documentation" target="_blank">NIMBLE User Manual</a>.</p>
<h2>Parametric meta analysis of Avandia myocardial infarctions (MIs)</h2>
<p>We will illustrate the use of nonparametric mixture models for modeling random effects distributions in the context of a meta-analysis of the side effects of a formerly very popular drug for diabetes called Avandia. The data we analyze played a role in raising serious questions about the safety of this drug. The question is whether Avandia use increases the risk of myocardial infarction (heart attack). There are 48 studies (the 49th study in the data file is different in some ways and excluded here), each with treatment and control arms.</p>
<div>
<div>
<div>
<pre><span>dat</span> <span>&lt;-</span> <span>read.csv</span><span>(</span><span>'https://rawgit.com/nimble-dev/nimble-demos/master/intro_bnp/avandia.csv'</span><span>)</span>
<span>head</span><span>(dat)</span>
</pre>
</div>
<div>
<pre>##   trial nAvandia avandiaMI nControl controlMI
## 1     1      357         2      176         0
## 2     2      391         2      207         1
## 3     3      774         1      185         1
## 4     4      213         0      109         1
## 5     5      232         1      116         0
## 6     6       43         0       47         1
</pre>
</div>
<div>
<pre><span>dat</span> <span>&lt;-</span> <span>dat[</span><span>-</span><span>49</span><span>, ]</span>
</pre>
</div>
</div>
</div>
<h3>Model formulation</h3>
<p>We begin with a standard generalized linear mixed model (GLMM)-based meta analysis. The vectors <img alt="n" title="n" src="http://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> and <img alt="x" title="x" src="http://s0.wp.com/latex.php?latex=x&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> contain the total number of patients in the control and the number of patients suffering from myocardial infarctions in the control group of each study, respectively. Similarly, the vectors <img alt="m" title="m" src="http://s0.wp.com/latex.php?latex=m&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> and <img alt="y" title="y" src="http://s0.wp.com/latex.php?latex=y&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> contain similar information for patients receiving the drug Avandia. The model takes the form</p>
<p><img alt="  x_{i} mid theta, gamma_i sim mbox{Bin} left(n_i, frac{expleft{ gamma_i right}}{1 + expleft{ gamma_i right}} right) , quadquad y_{i} mid theta, gamma_i sim mbox{Bin} left(m_i, frac{expleft{ theta + gamma_i right}}{1 + expleft{ theta + gamma_i right}} right)  " title="  x_{i} mid theta, gamma_i sim mbox{Bin} left(n_i, frac{expleft{ gamma_i right}}{1 + expleft{ gamma_i right}} right) , quadquad y_{i} mid theta, gamma_i sim mbox{Bin} left(m_i, frac{expleft{ theta + gamma_i right}}{1 + expleft{ theta + gamma_i right}} right)  " src="http://s0.wp.com/latex.php?latex=++x_%7Bi%7D+%5Cmid+%5Ctheta%2C+%5Cgamma_i+%5Csim+%5Cmbox%7BBin%7D+%5Cleft%28n_i%2C+%5Cfrac%7B%5Cexp%5Cleft%5C%7B+%5Cgamma_i+%5Cright%5C%7D%7D%7B1+%2B+%5Cexp%5Cleft%5C%7B+%5Cgamma_i+%5Cright%5C%7D%7D+%5Cright%29+%2C+%5Cquad%5Cquad+y_%7Bi%7D+%5Cmid+%5Ctheta%2C+%5Cgamma_i+%5Csim+%5Cmbox%7BBin%7D+%5Cleft%28m_i%2C+%5Cfrac%7B%5Cexp%5Cleft%5C%7B+%5Ctheta+%2B+%5Cgamma_i+%5Cright%5C%7D%7D%7B1+%2B+%5Cexp%5Cleft%5C%7B+%5Ctheta+%2B+%5Cgamma_i+%5Cright%5C%7D%7D+%5Cright%29++&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"></p>
<p>where the random effects, <img alt="gamma_i" title="gamma_i" src="http://s0.wp.com/latex.php?latex=%5Cgamma_i&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, follow a common normal distribution, <img alt="gamma_i sim mbox{N}(0, tau^2)" title="gamma_i sim mbox{N}(0, tau^2)" src="http://s0.wp.com/latex.php?latex=%5Cgamma_i+%5Csim+%5Cmbox%7BN%7D%280%2C+%5Ctau%5E2%29&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, and the <img alt="theta" title="theta" src="http://s0.wp.com/latex.php?latex=%5Ctheta&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> and <img alt="tau^2" title="tau^2" src="http://s0.wp.com/latex.php?latex=%5Ctau%5E2&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> are given reasonably non-informative priors. The parameter <img alt="theta" title="theta" src="http://s0.wp.com/latex.php?latex=%5Ctheta&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> quantifies the difference in risk between the control and treatment arms, while the <img alt="gamma_i" title="gamma_i" src="http://s0.wp.com/latex.php?latex=%5Cgamma_i&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> quantify study-specific variation.</p>
<p>This model can be specified in NIMBLE using the following code:</p>
<div>
<div>
<div>
<pre><span>x</span> <span>&lt;-</span> <span>dat</span><span>$</span><span>controlMI</span>
<span>n</span> <span>&lt;-</span> <span>dat</span><span>$</span><span>nControl</span>
<span>y</span> <span>&lt;-</span> <span>dat</span><span>$</span><span>avandiaMI</span>
<span>m</span> <span>&lt;-</span> <span>dat</span><span>$</span><span>nAvandia</span>

<span>nStudies</span> <span>&lt;-</span> <span>nrow</span><span>(dat)</span>
<span>data</span> <span>&lt;-</span> <span>list</span><span>(</span><span>x</span> <span>= x,</span> <span>y</span> <span>= y)</span>
<span>constants</span> <span>=</span> <span>list</span><span>(</span><span>n</span> <span>= n,</span> <span>m</span> <span>= m,</span> <span>nStudies</span> <span>= nStudies)</span>

<span>codeParam</span> <span>&lt;-</span> <span>nimbleCode</span><span>({</span>
    <span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>nStudies) {</span>
        <span>y[i]</span> <span>~</span> <span>dbin</span><span>(</span><span>size</span> <span>= m[i],</span> <span>prob</span> <span>= q[i])</span> <span># avandia MIs</span>
        <span>x[i]</span> <span>~</span> <span>dbin</span><span>(</span><span>size</span> <span>= n[i],</span> <span>prob</span> <span>= p[i])</span> <span># control MIs</span>
        <span>q[i]</span> <span>&lt;-</span> <span>expit</span><span>(theta</span> <span>+</span> <span>gamma[i])</span>       <span># Avandia log-odds</span>
        <span>p[i]</span> <span>&lt;-</span> <span>expit</span><span>(gamma[i])</span>               <span># control log-odds</span>
        <span>gamma[i]</span> <span>~</span> <span>dnorm</span><span>(mu,</span> <span>var</span> <span>= tau2)</span>      <span># study effects</span>
    <span>}</span>
    <span>theta</span> <span>~</span> <span>dflat</span><span>()</span>        <span># effect of Avandia</span>
    <span># random effects hyperparameters</span>
    <span>mu</span> <span>~</span> <span>dnorm</span><span>(</span><span>0</span><span>,</span> <span>10</span><span>)</span>
    <span>tau2</span> <span>~</span> <span>dinvgamma</span><span>(</span><span>2</span><span>,</span> <span>1</span><span>)</span>
<span>})</span>
</pre>
</div>
</div>
</div>
<h3>Running the MCMC</h3>
<p>Let’s run a basic MCMC.</p>
<div>
<div>
<div>
<pre><span>set.seed</span><span>(</span><span>9</span><span>)</span>
<span>inits</span> <span>=</span> <span>list</span><span>(</span><span>theta</span> <span>=</span> <span>0</span><span>,</span> <span>mu</span> <span>=</span> <span>0</span><span>,</span> <span>tau2</span> <span>=</span> <span>1</span><span>,</span> <span>gamma</span> <span>=</span> <span>rnorm</span><span>(nStudies))</span>

<span>samples</span> <span>&lt;-</span> <span>nimbleMCMC</span><span>(</span><span>code</span> <span>= codeParam,</span> <span>data</span> <span>= data,</span> <span>inits</span> <span>= inits,</span>
                      <span>constants</span> <span>= constants,</span> <span>monitors</span> <span>=</span> <span>c</span><span>(</span><span>&quot;mu&quot;</span><span>,</span> <span>&quot;tau2&quot;</span><span>,</span> <span>&quot;theta&quot;</span><span>,</span> <span>&quot;gamma&quot;</span><span>),</span>
                      <span>thin</span> <span>=</span> <span>10</span><span>,</span> <span>niter</span> <span>=</span> <span>22000</span><span>,</span> <span>nburnin</span> <span>=</span> <span>2000</span><span>,</span> <span>nchains</span> <span>=</span> <span>1</span><span>,</span> <span>setSeed</span> <span>=</span> <span>TRUE</span><span>)</span>
</pre>
</div>
<div>
<pre>## defining model...
</pre>
</div>
<div>
<pre>## building model...
</pre>
</div>
<div>
<pre>## setting data and initial values...
</pre>
</div>
<div>
<pre>## running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
## checking model sizes and dimensions...
## checking model calculations...
## model building finished.
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
## compilation finished.
## runMCMC's handling of nburnin changed in nimble version 0.6-11. Previously, nburnin samples were discarded *post-thinning*.  Now nburnin samples are discarded *pre-thinning*.  The number of samples returned will be floor((niter-nburnin)/thin).
## running chain 1...
</pre>
</div>
<div>
<pre>## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
</pre>
</div>
<div>
<pre><span>par</span><span>(</span><span>mfrow</span> <span>=</span> <span>c</span><span>(</span><span>1</span><span>,</span> <span>4</span><span>),</span> <span>cex</span> <span>=</span> <span>1.1</span><span>,</span> <span>mgp</span> <span>=</span> <span>c</span><span>(</span><span>1.8</span><span>,</span><span>.7</span><span>,</span><span>0</span><span>))</span>
<span>ts.plot</span><span>(samples[ ,</span> <span>'theta'</span><span>],</span> <span>xlab</span> <span>=</span> <span>'iteration'</span><span>,</span> <span>ylab</span> <span>=</span> <span>expression</span><span>(theta))</span>
<span>hist</span><span>(samples[ ,</span> <span>'theta'</span><span>],</span> <span>xlab</span> <span>=</span> <span>expression</span><span>(theta),</span> <span>main</span> <span>=</span> <span>'effect of Avandia'</span><span>)</span>

<span>gammaCols</span> <span>&lt;-</span> <span>grep</span><span>(</span><span>'gamma'</span><span>,</span> <span>colnames</span><span>(samples))</span>
<span>gammaMn</span> <span>&lt;-</span> <span>colMeans</span><span>(samples[ , gammaCols])</span>
<span>hist</span><span>(gammaMn,</span> <span>xlab</span> <span>=</span> <span>'posterior means'</span><span>,</span> <span>main</span> <span>=</span> <span>'random effects distribution'</span><span>)</span>
<span>hist</span><span>(samples[</span><span>1000</span><span>, gammaCols],</span> <span>xlab</span> <span>=</span> <span>'single draw'</span><span>,</span>
                   <span>main</span> <span>=</span> <span>'random effects distribution'</span><span>)</span>
</pre>
</div>
</div>
<div><img src="https://i2.wp.com/r-nimble.org/wp-content/uploads/2018/12/mcmc-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
</div>
<p>The results suggests there is an overall difference in risk between the control and treatment arms. But what about the normality assumption? Are our conclusions robust to that assumption? Perhaps the random effects distribution are skewed. (And recall that the estimates above of the random effects are generated under the normality assumption, which pushes the estimated effects to look more normal…)</p>
<h2>DP-based random effects modeling for meta analysis</h2>
<h3>Model formulation</h3>
<p>Now, we use a nonparametric distribution for the <img alt="gamma_i" title="gamma_i" src="http://s0.wp.com/latex.php?latex=%5Cgamma_i&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">s. More specifically, we assume that each <img alt="gamma_i" title="gamma_i" src="http://s0.wp.com/latex.php?latex=%5Cgamma_i&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> is generated from a location-scale mixture of normal distributions:</p>
<p><img alt="  gamma_i mid mu_i, tau_i^2 sim mbox{N}(mu_i, tau_i^2), quadquad (mu_i, tau_i^2) mid G sim G, quadquad G sim mbox{DP}(alpha, H),  " title="  gamma_i mid mu_i, tau_i^2 sim mbox{N}(mu_i, tau_i^2), quadquad (mu_i, tau_i^2) mid G sim G, quadquad G sim mbox{DP}(alpha, H),  " src="http://s0.wp.com/latex.php?latex=++%5Cgamma_i+%5Cmid+%5Cmu_i%2C+%5Ctau_i%5E2+%5Csim+%5Cmbox%7BN%7D%28%5Cmu_i%2C+%5Ctau_i%5E2%29%2C+%5Cquad%5Cquad+%28%5Cmu_i%2C+%5Ctau_i%5E2%29+%5Cmid+G+%5Csim+G%2C+%5Cquad%5Cquad+G+%5Csim+%5Cmbox%7BDP%7D%28%5Calpha%2C+H%29%2C++&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"></p>
<p>where <img alt="H" title="H" src="http://s0.wp.com/latex.php?latex=H&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> is a normal-inverse-gamma distribution.</p>
<p>This specification induces clustering among the random effects. As in the case of density estimation problems, the DP prior allows the data to determine the number of components, from as few as one component (i.e., simplifying to the parametric model), to as many as <img alt="n" title="n" src="http://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> components, i.e., one component for each observation. This allows the distribution of the random effects to be multimodal if the data supports such behavior, greatly increasing its flexibility. This model can be specified in NIMBLE using the following code:</p>
<div>
<div>
<div>
<pre><span>codeBNP</span> <span>&lt;-</span> <span>nimbleCode</span><span>({</span>
    <span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>nStudies) {</span>
        <span>y[i]</span> <span>~</span> <span>dbin</span><span>(</span><span>size</span> <span>= m[i],</span> <span>prob</span> <span>= q[i])</span>   <span># avandia MIs</span>
        <span>x[i]</span> <span>~</span> <span>dbin</span><span>(</span><span>size</span> <span>= n[i],</span> <span>prob</span> <span>= p[i])</span>   <span># control MIs</span>
        <span>q[i]</span> <span>&lt;-</span> <span>expit</span><span>(theta</span> <span>+</span> <span>gamma[i])</span>         <span># Avandia log-odds</span>
        <span>p[i]</span> <span>&lt;-</span> <span>expit</span><span>(gamma[i])</span>                 <span># control log-odds</span>
        <span>gamma[i]</span> <span>~</span> <span>dnorm</span><span>(mu[i],</span> <span>var</span> <span>= tau2[i])</span>  <span># random effects from mixture dist.</span>
        <span>mu[i]</span> <span>&lt;-</span> <span>muTilde[xi[i]]</span>                 <span># mean for random effect from cluster xi[i]</span>
        <span>tau2[i]</span> <span>&lt;-</span> <span>tau2Tilde[xi[i]]</span>             <span># var for random effect from cluster xi[i]</span>
    <span>}</span>
    <span># mixture component parameters drawn from base measures</span>
    <span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>nStudies) {</span>
        <span>muTilde[i]</span> <span>~</span> <span>dnorm</span><span>(mu0,</span> <span>var</span> <span>= var0)</span>
        <span>tau2Tilde[i]</span> <span>~</span> <span>dinvgamma</span><span>(a0, b0)</span>
    <span>}</span>
    <span># CRP for clustering studies to mixture components</span>
    <span>xi[</span><span>1</span><span>:</span><span>nStudies]</span> <span>~</span> <span>dCRP</span><span>(alpha,</span> <span>size</span> <span>= nStudies)</span>
    <span># hyperparameters</span>
    <span>alpha</span> <span>~</span> <span>dgamma</span><span>(</span><span>1</span><span>,</span> <span>1</span><span>)</span>
    <span>mu0</span> <span>~</span> <span>dnorm</span><span>(</span><span>0</span><span>,</span> <span>10</span><span>)</span>
    <span>var0</span> <span>~</span> <span>dinvgamma</span><span>(</span><span>2</span><span>,</span> <span>1</span><span>)</span>
    <span>a0</span> <span>~</span> <span>dinvgamma</span><span>(</span><span>2</span><span>,</span> <span>1</span><span>)</span>
    <span>b0</span> <span>~</span> <span>dinvgamma</span><span>(</span><span>2</span><span>,</span> <span>1</span><span>)</span>
    <span>theta</span> <span>~</span> <span>dflat</span><span>()</span>          <span># effect of Avandia</span>
<span>})</span>
</pre>
</div>
</div>
</div>
<h3>Running the MCMC</h3>
<p>The following code compiles the model and runs a collapsed Gibbs sampler for the model</p>
<div>
<div>
<div>
<pre><span>inits</span> <span>&lt;-</span> <span>list</span><span>(</span><span>gamma</span> <span>=</span> <span>rnorm</span><span>(nStudies),</span> <span>xi</span> <span>=</span> <span>sample</span><span>(</span><span>1</span><span>:</span><span>2</span><span>, nStudies,</span> <span>replace</span> <span>=</span> <span>TRUE</span><span>),</span>
              <span>alpha</span> <span>=</span> <span>1</span><span>,</span> <span>mu0</span> <span>=</span> <span>0</span><span>,</span> <span>var0</span> <span>=</span> <span>1</span><span>,</span> <span>a0</span> <span>=</span> <span>1</span><span>,</span> <span>b0</span> <span>=</span> <span>1</span><span>,</span> <span>theta</span> <span>=</span> <span>0</span><span>,</span>
              <span>muTilde</span> <span>=</span> <span>rnorm</span><span>(nStudies),</span> <span>tau2Tilde</span> <span>=</span> <span>rep</span><span>(</span><span>1</span><span>, nStudies))</span>

<span>samplesBNP</span> <span>&lt;-</span> <span>nimbleMCMC</span><span>(</span><span>code</span> <span>= codeBNP,</span> <span>data</span> <span>= data,</span> <span>inits</span> <span>= inits,</span>
               <span>constants</span> <span>= constants,</span>
               <span>monitors</span> <span>=</span> <span>c</span><span>(</span><span>&quot;theta&quot;</span><span>,</span> <span>&quot;gamma&quot;</span><span>,</span> <span>&quot;alpha&quot;</span><span>,</span> <span>&quot;xi&quot;</span><span>,</span> <span>&quot;mu0&quot;</span><span>,</span> <span>&quot;var0&quot;</span><span>,</span> <span>&quot;a0&quot;</span><span>,</span> <span>&quot;b0&quot;</span><span>),</span>
               <span>thin</span> <span>=</span> <span>10</span><span>,</span> <span>niter</span> <span>=</span> <span>22000</span><span>,</span> <span>nburnin</span> <span>=</span> <span>2000</span><span>,</span> <span>nchains</span> <span>=</span> <span>1</span><span>,</span> <span>setSeed</span> <span>=</span> <span>TRUE</span><span>)</span>
</pre>
</div>
<div>
<pre>## defining model...
</pre>
</div>
<div>
<pre>## building model...
</pre>
</div>
<div>
<pre>## setting data and initial values...
</pre>
</div>
<div>
<pre>## running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
## checking model sizes and dimensions...
## checking model calculations...
## model building finished.
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
## compilation finished.
## runMCMC's handling of nburnin changed in nimble version 0.6-11. Previously, nburnin samples were discarded *post-thinning*.  Now nburnin samples are discarded *pre-thinning*.  The number of samples returned will be floor((niter-nburnin)/thin).
## running chain 1...
</pre>
</div>
<div>
<pre>## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
</pre>
</div>
<div>
<pre><span>gammaCols</span> <span>&lt;-</span> <span>grep</span><span>(</span><span>'gamma'</span><span>,</span> <span>colnames</span><span>(samplesBNP))</span>
<span>gammaMn</span> <span>&lt;-</span> <span>colMeans</span><span>(samplesBNP[ , gammaCols])</span>
<span>xiCols</span> <span>&lt;-</span> <span>grep</span><span>(</span><span>'xi'</span><span>,</span> <span>colnames</span><span>(samplesBNP))</span>

<span>par</span><span>(</span><span>mfrow</span> <span>=</span> <span>c</span><span>(</span><span>1</span><span>,</span><span>5</span><span>),</span> <span>cex</span> <span>=</span> <span>1.1</span><span>,</span> <span>mgp</span> <span>=</span> <span>c</span><span>(</span><span>1.8</span><span>,</span><span>.7</span><span>,</span><span>0</span><span>))</span>
<span>ts.plot</span><span>(samplesBNP[ ,</span> <span>'theta'</span><span>],</span> <span>xlab</span> <span>=</span> <span>'iteration'</span><span>,</span> <span>ylab</span> <span>=</span> <span>expression</span><span>(theta),</span>
   <span>main</span> <span>=</span> <span>expression</span><span>(</span><span>paste</span><span>(</span><span>'traceplot for '</span><span>, theta)))</span>
<span>hist</span><span>(samplesBNP[ ,</span> <span>'theta'</span><span>],</span> <span>xlab</span> <span>=</span> <span>expression</span><span>(theta),</span> <span>main</span> <span>=</span> <span>'effect of Avandia'</span><span>)</span>
<span>hist</span><span>(gammaMn,</span> <span>xlab</span> <span>=</span> <span>'posterior means'</span><span>,</span>
              <span>main</span> <span>=</span> <span>&quot;random effects distrib'n&quot;</span><span>)</span>
<span>hist</span><span>(samplesBNP[</span><span>1000</span><span>, gammaCols],</span> <span>xlab</span> <span>=</span> <span>'single draw'</span><span>,</span>
                   <span>main</span> <span>=</span> <span>&quot;random effects distrib'n&quot;</span><span>)</span>

<span># How many mixture components are inferred?</span>
<span>xiRes</span> <span>&lt;-</span> <span>samplesBNP[ , xiCols]</span>
<span>nGrps</span> <span>&lt;-</span> <span>apply</span><span>(xiRes,</span> <span>1</span><span>,</span> <span>function</span><span>(</span><span>x</span><span>)</span> <span>length</span><span>(</span><span>unique</span><span>(x)))</span>
<span>ts.plot</span><span>(nGrps,</span> <span>xlab</span> <span>=</span> <span>'iteration'</span><span>,</span> <span>ylab</span> <span>=</span> <span>'number of components'</span><span>,</span>
   <span>main</span> <span>=</span> <span>'number of components'</span><span>)</span>
</pre>
</div>
</div>
<div><img src="https://i0.wp.com/r-nimble.org/wp-content/uploads/2018/12/DP-MCMC-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
</div>
<p>The primary inference seems robust to the original parametric assumption. This is probably driven by the fact that there is not much evidence of lack of normality in the random effects distribution (as evidenced by the fact that the posterior distribution of the number of mixture components places a large amount of probability on exactly one component).</p>
<h2>More information and future development</h2>
<p>Please see our <a rel="nofollow" href="https://r-nimble.org/documentation" target="_blank">User Manual</a> for more details.</p>
<p>We’re in the midst of improvements to the existing BNP functionality as well as adding additional Bayesian nonparametric models, such as hierarchical Dirichlet processes and Pitman-Yor processes, so please add yourself to our <a rel="nofollow" href="https://groups.google.com/forum/#!forum/nimble-announce" target="_blank">announcement</a> or <a rel="nofollow" href="https://groups.google.com/forum/#!forum/nimble-users" target="_blank">user support/discussion</a> Google groups.</p>
<p> <img width="1" src="https://r-nimble.org/?feed-stats-post-id=902" height="1"></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://r-nimble.org/bayesian-nonparametric-models-in-nimble-part-2-nonparametric-random-effects"> R – NIMBLE</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/sfxYVL8qUZA" height="1"></a></div>