<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/estimating-treatment-effects-and-iccs-from-glmms-on-the-observed-scale-using-bayes-part-1-lognormal-models/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://rpsychologist.com/GLMM-part1-lognormal"> R Psychologist - R</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>

<p>When a <strong>multilevel</strong> model includes either a non-linear transformation (such as the log-transformation) of the response variable, or of the expectations via a GLM link-function, then the interpretation of the results will be different compared to a standard Gaussian multilevel model; specifically, the estimates will be on a transformed scale and not in the original units, and the effects will no longer refer to the average effect in the population, instead they are conditional/cluster-specific. In this post, I will deal with linear mixed-effects models (LMM) that use a log-transformed outcome variable. This will be the first part of a three-part tutorial on some of the finer details of (G)LMMs, and how Bayes can make your (frequentist) life easier. In each post I will focus on:</p>
<ul>
<li><strong>log-transformation</strong> of the dependent variable in a <strong>multilevel</strong> model</li>
<li>Calculating both <strong>multiplicative effects</strong> (% change) and differences on the <strong>untransformed scale</strong>.</li>
<li><strong>Subject-specific</strong>/cluster-specific versus <strong>population-average effects</strong> (conditional versus marginal effects).</li>
<li><strong>Intraclass correlations</strong> on the transformed/link/latent scale or on the response/data/original scale.</li>
</ul>
<p>In part 2 I will cover a GLMM with a binary outcome, and part 3 will focus on semicontinuous (hurdle/two-part) models when the outcome is a skewed continuous variable that include zeros.</p>
<p>The aim of this post is also to show:</p>
<ul>
<li>How simulation-based approaches like MCMC make it much easier to make inferences about transformed parameters.</li>
<li>How the estimates from a multilevel model can be transformed to answer the same questions as population-average models or fixed effects models. </li>
</ul>
<h2>Data-generating model: a lognormal LMM</h2>
<p>As an example, we will use a simple hierarchical design with clusters nested within either a treatment or a control condition. This is a classic 2-level multilevel model.</p>
<div>$$ \text{log}(y_{ij} \mid u_i) \sim \mathcal{N}(\beta_0 + \beta_1 TX_{i} + u_i, \sigma_e^2) $$</div>
<div>$$ u_i \sim \mathcal{N}(0, \sigma^2_u) $$</div>
<p>
In this model, the natural logarithm of <span>\(y\)</span> is normally distributed conditional on the cluster-specific effect <span>\(u_i\)</span>. Let’s assume <span>\(y\)</span> is the average money (in USD) lost gambling per week (since I do research on gambling disorder). Gambling expenditure is sometimes quite well described by a lognormal distribution, or a (generalized) Gamma distribution. Importantly, the methods in this post can quite easily be tweaked to work with other response distributions. Each treatment arm has <span>\(n_2\)</span> clusters, and each cluster has <span>\(n_1\)</span> observations. In this example clusters refer to e.g. schools, centers, groups, and the first level are subjects. However, it would be easy to use the contents of this post for a longitudinal model where clusters refer to subjects, and the subjects have repeated observations.</p>
<h2>Simulate data</h2>
<p>Let’s create a simple function that generates data from this model</p>
<div>
<pre><span></span><span>#' Generate lognormal data with a random intercept</span>
<span>#'</span>
<span>#' @param n1 patients per cluster</span>
<span>#' @param n2 clusters per treatment</span>
<span>#' @param B0 log intercept</span>
<span>#' @param B1 log treatment effect</span>
<span>#' @param sd_log log sd</span>
<span>#' @param u0 SD of log intercepts (random intercept)</span>
<span>#'</span>
<span>#' @return a data.frame</span>
gen_data <span>&lt;-</span> <span>function</span><span>(</span>n1<span>,</span> n2<span>,</span> B0<span>,</span> B1<span>,</span> sd_log<span>,</span> u0<span>)</span> <span>{</span>

  cluster <span>&lt;-</span> <span>rep</span><span>(</span><span>1</span><span>:</span><span>(</span><span>2</span> <span>*</span> n2<span>),</span> each <span>=</span> n1<span>)</span>
  TX <span>&lt;-</span> <span>rep</span><span>(</span><span>c</span><span>(</span><span>0</span><span>,</span> <span>1</span><span>),</span> each <span>=</span> n1 <span>*</span> n2<span>)</span>
  u0 <span>&lt;-</span> rnorm<span>(</span><span>2</span> <span>*</span> n2<span>,</span> sd <span>=</span> u0<span>)[</span>cluster<span>]</span>

  mulog <span>&lt;-</span> <span>(</span>B0 <span>+</span> B1 <span>*</span> TX <span>+</span> u0<span>)</span>
  y <span>&lt;-</span> rlnorm<span>(</span><span>2</span> <span>*</span> n1 <span>*</span> n2<span>,</span> meanlog <span>=</span> mulog<span>,</span> sdlog <span>=</span> sd_log<span>)</span>

  d <span>&lt;-</span> <span>data.frame</span><span>(</span>cluster<span>,</span>
                  TX<span>,</span>
                  y<span>)</span>
  d
<span>}</span>
</pre>
</div>
<p>The model parameters are</p>
<div>
<pre><span></span><span>set.seed</span><span>(</span><span>4445</span><span>)</span>
pars <span>&lt;-</span> <span>list</span><span>(</span><span>&quot;n1&quot;</span> <span>=</span> <span>30</span><span>,</span> <span># observations per cluster</span>
             <span>&quot;n2&quot;</span> <span>=</span> <span>10</span><span>,</span> <span># clusters per treatment</span>
             <span>&quot;B0&quot;</span> <span>=</span> <span>log</span><span>(</span><span>500</span><span>),</span>
             <span>&quot;B1&quot;</span> <span>=</span> <span>log</span><span>(</span><span>0.5</span><span>),</span>
             <span>&quot;sd_log&quot;</span> <span>=</span> <span>0.5</span><span>,</span>
             <span>&quot;u0&quot;</span> <span>=</span> <span>0.5</span><span>)</span>

d <span>&lt;-</span> <span>do.call</span><span>(</span>gen_data<span>,</span>
             pars<span>)</span>
</pre>
</div>
<h2>Fit the model</h2>
<p>We can fit this model using <code>brms</code> (or rstanarm, lme4, etc).</p>
<div>
<pre><span></span>fit <span>&lt;-</span> brm<span>(</span>y <span>~</span> <span>1</span> <span>+</span> TX <span>+</span> <span>(</span><span>1</span> <span>|</span> cluster<span>),</span> 
           family <span>=</span> lognormal<span>(),</span> 
           data <span>=</span> d<span>,</span>
           cores <span>=</span> <span>4</span><span>,</span> 
           seed <span>=</span> <span>4445</span><span>)</span>
fit
</pre>
</div>
<div>
<pre><span></span><span>##</span>  <span>Family</span><span>:</span> <span>lognormal</span> 
<span>##</span>   <span>Links</span><span>:</span> <span>mu</span> <span>=</span> <span>identity</span><span>;</span> <span>sigma</span> <span>=</span> <span>identity</span> 
<span>##</span> <span>Formula</span><span>:</span> <span>y</span> <span>~</span> <span>1</span> <span>+</span> <span>TX</span> <span>+</span> <span>(</span><span>1</span> <span>|</span> <span>cluster</span><span>)</span> 
<span>##</span>    <span>Data</span><span>:</span> <span>d</span> <span>(</span><span>Number</span> <span>of</span> <span>observations</span><span>:</span> <span>600</span><span>)</span> 
<span>##</span> <span>Samples</span><span>:</span> <span>4</span> <span>chains</span><span>,</span> <span>each</span> <span>with</span> <span>iter</span> <span>=</span> <span>2000</span><span>;</span> <span>warmup</span> <span>=</span> <span>1000</span><span>;</span> <span>thin</span> <span>=</span> <span>1</span><span>;</span>
<span>##</span>          <span>total</span> <span>post-warmup</span> <span>samples</span> <span>=</span> <span>4000</span>
<span>##</span> 
<span>##</span> <span>Group-Level</span> <span>Effects</span><span>:</span> 
<span>##</span> <span>~</span><span>cluster</span> <span>(</span><span>Number</span> <span>of</span> <span>levels</span><span>:</span> <span>20</span><span>)</span> 
<span>##</span>               <span>Estimate</span> <span>Est</span><span>.</span><span>Error</span> <span>l-95</span><span>%</span> <span>CI</span> <span>u-95</span><span>%</span> <span>CI</span> <span>Eff</span><span>.</span><span>Sample</span> <span>Rhat</span>
<span>##</span> <span>sd</span><span>(</span><span>Intercept</span><span>)</span>     <span>0</span><span>.</span><span>64</span>      <span>0</span><span>.</span><span>12</span>     <span>0</span><span>.</span><span>45</span>     <span>0</span><span>.</span><span>91</span>        <span>640</span> <span>1</span><span>.</span><span>01</span>
<span>##</span> 
<span>##</span> <span>Population-Level</span> <span>Effects</span><span>:</span> 
<span>##</span>           <span>Estimate</span> <span>Est</span><span>.</span><span>Error</span> <span>l-95</span><span>%</span> <span>CI</span> <span>u-95</span><span>%</span> <span>CI</span> <span>Eff</span><span>.</span><span>Sample</span> <span>Rhat</span>
<span>##</span> <span>Intercept</span>     <span>6</span><span>.</span><span>21</span>      <span>0</span><span>.</span><span>21</span>     <span>5</span><span>.</span><span>80</span>     <span>6</span><span>.</span><span>61</span>        <span>583</span> <span>1</span><span>.</span><span>00</span>
<span>##</span> <span>TX</span>           <span>-0</span><span>.</span><span>70</span>      <span>0</span><span>.</span><span>30</span>    <span>-1</span><span>.</span><span>25</span>    <span>-0</span><span>.</span><span>09</span>        <span>548</span> <span>1</span><span>.</span><span>00</span>
<span>##</span> 
<span>##</span> <span>Family</span> <span>Specific</span> <span>Parameters</span><span>:</span> 
<span>##</span>       <span>Estimate</span> <span>Est</span><span>.</span><span>Error</span> <span>l-95</span><span>%</span> <span>CI</span> <span>u-95</span><span>%</span> <span>CI</span> <span>Eff</span><span>.</span><span>Sample</span> <span>Rhat</span>
<span>##</span> <span>sigma</span>     <span>0</span><span>.</span><span>51</span>      <span>0</span><span>.</span><span>01</span>     <span>0</span><span>.</span><span>48</span>     <span>0</span><span>.</span><span>54</span>       <span>1834</span> <span>1</span><span>.</span><span>00</span>
<span>##</span> 
<span>##</span> <span>Samples</span> <span>were</span> <span>drawn</span> <span>using</span> <span>sampling</span><span>(</span><span>NUTS</span><span>).</span> <span>For</span> <span>each</span> <span>parameter</span><span>,</span> <span>Eff</span><span>.</span><span>Sample</span> 
<span>##</span> <span>is</span> <span>a</span> <span>crude</span> <span>measure</span> <span>of</span> <span>effective</span> <span>sample</span> <span>size</span><span>,</span> <span>and</span> <span>Rhat</span> <span>is</span> <span>the</span> <span>potential</span> 
<span>##</span> <span>scale</span> <span>reduction</span> <span>factor</span> <span>on</span> <span>split</span> <span>chains</span> <span>(</span><span>at</span> <span>convergence</span><span>,</span> <span>Rhat</span> <span>=</span> <span>1</span><span>).</span>
</pre>
</div>
<p>Using the notation of our 2-level model, we get these results</p>
<table>
<thead>
<tr>
<th>Parameter</th>
<th>Estimate</th>
</tr>
</thead>
<tbody>
<tr>
<td><span>\(\hat\beta_0\)</span></td>
<td>6.21</td>
</tr>
<tr>
<td><span>\(\hat\beta_1\)</span></td>
<td>-0.70</td>
</tr>
<tr>
<td><span>\(\hat\sigma_u\)</span></td>
<td>0.64</td>
</tr>
<tr>
<td><span>\(\hat\sigma_e\)</span></td>
<td>0.51</td>
</tr>
</tbody>
</table>
<h2>Interpreting the parameters</h2>
<p>All of the estimates are on the log scale, and they do not have an immediate interpretation on the natural scale. To further complicate things we should also remind ourselves that the effects from (G)LMM have cluster-specific interpretations (“subject-specific” in a longitudinal model when clusters == subjects). Researchers who work mostly with LMMs sometimes tend to forget this, since for LMMs the cluster-specific effects are equal to the population-average (marginal) effects. The non-linear transformation of either the data (like the log-transformation) or via the link function in (G)LMMs both means that back-transformation of the model’s estimates from e.g. the log scale does not give the expected values on the natural scale (cf. Jensen’s inequality). </p>
<p>In the figure below I try to illustrate both the cluster-specific nature of the model, and the transformation between the log scale and the natural data scale, without being too technical. </p>
<p><img alt="lognormal multilevel LMM visualizing cluster-specific and population-average values" src="https://i0.wp.com/rpsychologist.com/img/GLMM-part1-lognormal/interpret-lognormal-LMM.png?w=456" data-recalc-dims="1"></p>
<h3>How to calculate marginal effects on the data scale</h3>
<p>So how do we derive the numbers presented in the figure above? Well, for the lognormal distribution with the log mean <span>\(\mu\)</span> and the log SD <span>\(\sigma_e\)</span> we have:</p>
<ul>
<li>mean: <span>\(\text{exp}(\mu + \frac{\sigma_e^2}{2})\)</span></li>
<li>median: <span>\(\text{exp}(\mu)\)</span></li>
</ul>
<p>If we want to get the marginal mean we need to average this back-transformation over the random effects distribution. Since, as noted earlier, simply back-transforming the log coefficient (the expected mean on the log scale) will not give the expected mean on the original scale (money lost). </p>
<div>$$ \text{exp}(\beta_0 + \frac{\sigma^2_e}{2}) \ne \text{E} (Y_{ij} \mid TX = 0)$$</div>
<p>The correct expectation is</p>
<div>$$ \text{E}(y_{ij} \mid TX = 0) = \int \text{exp}(x + \sigma^2_e/2) f_{\mathcal{N}}(x, \beta_{0}, \sigma_{u}^2)dx$$</div>
<p>
This integral is easily solved using numerical integration. </p>
<div>
<pre><span></span>mu <span>&lt;-</span> integrate<span>(</span>
  f <span>=</span> <span>function</span><span>(</span>x<span>)</span> <span>{</span>
    <span>exp</span><span>(</span>x <span>+</span> sd_log <span>^</span> <span>2</span> <span>/</span> <span>2</span><span>)</span> <span>*</span> dnorm<span>(</span>x<span>,</span> B0<span>,</span> sd <span>=</span> u0<span>)</span>
  <span>},</span>
  lower <span>=</span> B0 <span>-</span> <span>10</span> <span>*</span> u0<span>,</span>
  upper <span>=</span> B0 <span>+</span> <span>10</span> <span>*</span> u0
<span>)</span><span>$</span>value
mu
</pre>
</div>
<div>
<pre><span></span>## [1] 642.0127
</pre>
</div>
<p>In this case, we can also derive an exact and simple solution using the moment-generating function of the lognormal distribution (see Proof 1: left as an exercise to the reader <img alt="🤓" src="https://s.w.org/images/core/emoji/2.4/72x72/1f913.png" class="wp-smiley">)</p>
<div>
<pre><span></span><span>exp</span><span>(</span>B0 <span>+</span> <span>(</span>u0<span>^</span><span>2</span> <span>+</span> sd_log<span>^</span><span>2</span><span>)</span><span>/</span><span>2</span><span>)</span>
</pre>
</div>
<div>
<pre><span></span>## [1] 642.0127
</pre>
</div>
<p>A third option is Monte Carlo integration. Were we simply simulate a large number of cluster outcomes on the log scale, back-transform to the natural scale, and then get the average. For multivariate integrals (e.g. models with correlated random effects) this is often the simplest solution that also works well in practice, although the package <code>cubature</code> offers a fast solution using adaptive numerical integration. </p>
<div>
<pre><span></span>u_i <span>&lt;-</span> rnorm<span>(</span><span>1e5</span><span>,</span> <span>0</span><span>,</span> u0<span>)</span>
y <span>&lt;-</span> <span>exp</span><span>(</span>B0 <span>+</span> u_i <span>+</span> sd_log<span>^</span><span>2</span><span>/</span><span>2</span><span>)</span>
<span>mean</span><span>(</span>y<span>)</span>
</pre>
</div>
<div>
<pre><span></span>## [1] 643.239
</pre>
</div>
<p>We see that all three methods give the same result (within Monte Carlo error).</p>
<h3>Can we interpret the cluster-specific values?</h3>
<p>So, if the expected value for the outcome group is 642 USD, and <span>\(\text{exp}(\beta_0 + \frac{\sigma^2_e}{2})\)</span> = 567 USD is there a meaningful interpretation of 567? As shown in the figure above, this is the expected value for the cluster at the center of the distribution of cluster means on the <em>log scale</em>. Since percentiles are invariant under a monotonic transformation (and the mean is equal to the median for the normal distribution), this cluster (with <span>\(u_i = 0\)</span>) is also at the median on the natural scale. </p>
<div>
<pre><span></span><span>exp</span><span>(</span>B0 <span>+</span> sd_log<span>^</span><span>2</span><span>/</span><span>2</span><span>)</span>
</pre>
</div>
<div>
<pre><span></span>## [1] 566.5742
</pre>
</div>
<p>So to summarize, the cluster-specific effects are,</p>
<div>$$ \text{E}(y_{ij} \mid u_i, TX_i) = \text{exp}(\beta_0 + \beta_1TX_i + u_i + \frac{\sigma^2_e}{2}),$$</div>
<p>and the expected marginal mean in each treatment group is,</p>
<div>$$ \text{E}(y_{ij} | TX_i) = \text{exp}(\beta_0 + \beta_1TX_i + \frac{\sigma^2_u}{2} + \frac{\sigma^2_e}{2})$$</div>
<h3>Average treatment effects on the data scale</h3>
<p>Using the results above, the cluster-specific average treatment effect on the original scale is,</p>
<div>$$ \text{E}(Y_{ij} \mid u_i, \text{TX} = 1) – \text{E}(Y_{ij} \mid u_i, \text{TX} = 0) = \\\<br>
\text{exp}(\beta_0 + \beta_1 + u_i + \frac{\sigma^2_e}{2}) – \text{exp}(\beta_0 + u_i + \frac{\sigma^2_e}{2})$$</div>
<p>which is (for <span>\(u_i = 0\)</span>),</p>
<div>
<pre><span></span>## [1] -283.2871
</pre>
</div>
<p>And the marginal treatment effect on the original scale is,</p>
<div>$$ \text{E}(Y_{ij} \mid \text{TX} = 1) – \text{E}(Y_{ij} \mid \text{TX} = 0) =<br>
\text{exp}(\beta_0 +\beta_1 + \frac{\sigma^2_e + \sigma^2_u}{2}) – \text{exp}(\beta_0 + \frac{\sigma^2_e + \sigma^2_u}{2})$$</div>
<div>
<pre><span></span>## [1] -321.0064
</pre>
</div>
<h3>Multiplicate effects (% change)</h3>
<p>Most know that if we take the exponent of the log differences we get the % change. Generally, these multiplicative effects are also cluster-specific. However, for our simple model <span>\(\text{exp}(\beta_1)\)</span> gives both the cluster-specific and population-average effect. We can see this if we look at the ratio of the expected conditional outcome in the treatment and control group,</p>
<div>$$<br>
\begin{align*}<br>
\frac{\text{E}(Y_{ij} \mid u_i, \text{TX} = 1)}{\text{E}(Y_{ij} \mid u_i, \text{TX} = 0)} &amp;= \frac{\text{exp}(\beta_0 + \beta_1 +u_i + \frac{\sigma^2_e}{2})}<br>
{\text{exp}(\beta_0 + u_i + \frac{\sigma^2_e}{2})} \\\ &amp;= \text{exp}(\beta_1)<br>
\end{align*}<br>
$$</div>
<p>and similarly for the marginal outcomes
</p>
<div>$$<br>
\begin{align*}<br>
\frac{\text{E}(Y_{ij} \mid \text{TX} = 1)}{\text{E}(Y_{ij} \mid \text{TX} = 0)} &amp;= \frac{\text{exp}(\beta_0 + \beta_1 + \frac{\sigma^2_u + \sigma^2_e}{2})}<br>
{\text{exp}(\beta_0 + \frac{\sigma^2_u + \sigma^2_e}{2})} \\\ &amp;= \text{exp}(\beta_1).<br>
\end{align*}<br>
$$</div>
<p>The cluster-specific and population-average effect is the same. However, this is not true in general for nonlinear transformation or GLMM link functions (which we shall see in part 2 on binary outcomes).</p>
<h3>95 % CIs for the backtransformed marginal effects</h3>
<p>Getting back to our fitted model, let’s see how we can transform the estimates into the posterior distribution of e.g. the marginal mean outcome in the control group. We do this by applying the back-transformation and marginalization using all the posterior samples of <span>\(\hat\beta_0\)</span>, <span>\(\hat\sigma_u\)</span>, and <span>\(\hat\sigma_e\)</span>. </p>
<div>
<pre><span></span>marg_f <span>&lt;-</span> <span>function</span><span>(</span>B0<span>,</span> u0<span>,</span> sd_log<span>)</span> <span>{</span>
  <span>data.frame</span><span>(</span><span>&quot;marginal&quot;</span> <span>=</span> <span>exp</span><span>(</span>B0 <span>+</span> <span>(</span>u0<span>^</span><span>2</span> <span>+</span> sd_log<span>^</span><span>2</span><span>)</span><span>/</span><span>2</span><span>),</span>
             <span>&quot;conditional&quot;</span> <span>=</span> <span>exp</span><span>(</span>B0 <span>+</span> <span>(</span>sd_log<span>^</span><span>2</span><span>)</span><span>/</span><span>2</span><span>)</span>
             <span>)</span>
<span>}</span>
ss <span>&lt;-</span> posterior_samples<span>(</span>fit<span>)</span>

<span># mapply would also work</span>
mu <span>&lt;-</span> pmap_df<span>(</span>
    <span>list</span><span>(</span>
      B0 <span>=</span> ss<span>$</span>b_Intercept<span>,</span>
      u0 <span>=</span> ss<span>$</span>sd_cluster__Intercept<span>,</span>
      sd_log <span>=</span> ss<span>$</span>sigma
    <span>),</span>
    <span>.</span>f <span>=</span> marg_f
  <span>)</span>

map<span>(</span>mu<span>,</span> quantile<span>,</span> probs <span>=</span> <span>c</span><span>(</span><span>0.025</span><span>,</span> <span>0.5</span><span>,</span> <span>0.975</span><span>))</span> <span>%>%</span> 
    <span>do.call</span><span>(</span><span>rbind</span><span>,</span> <span>.</span><span>)</span> 
</pre>
</div>
<div>
<pre><span></span>##                 2.5%      50%     97.5%
## marginal    454.4591 697.9116 1123.3794
## conditional 374.4601 570.5069  847.8781
</pre>
</div>
<h3>Summarizing the model</h3>
<p>Let’s summarize all the results from the model, both the conditional and marginal, absolute and multiplicative treatment effects. </p>
<p>

</p>
<div>
<div>
<pre><span></span>ss <span>&lt;-</span> posterior_samples<span>(</span>fit<span>,</span> pars <span>=</span> <span>c</span><span>(</span><span>&quot;b_&quot;</span><span>,</span> 
                                      <span>&quot;sd&quot;</span><span>,</span> 
                                      <span>&quot;sigma&quot;</span><span>))</span>


res <span>&lt;-</span> ss <span>%>%</span> 
  select<span>(</span>B0 <span>=</span> b_Intercept<span>,</span> B1 <span>=</span> b_TX<span>,</span> u0 <span>=</span> sd_cluster__Intercept<span>,</span> sigma<span>)</span> <span>%>%</span> 
  mutate<span>(</span>y_cond_cc <span>=</span> <span>exp</span><span>(</span>B0 <span>+</span> sigma<span>^</span><span>2</span><span>/</span><span>2</span><span>),</span>
         y_cond_tx <span>=</span> <span>exp</span><span>(</span>B0 <span>+</span> B1 <span>+</span> sigma<span>^</span><span>2</span><span>/</span><span>2</span><span>),</span>
         y_marg_cc <span>=</span> <span>exp</span><span>(</span>B0 <span>+</span> <span>(</span>u0<span>^</span><span>2</span> <span>+</span> sigma<span>^</span><span>2</span><span>)</span><span>/</span><span>2</span><span>),</span>
         y_marg_tx <span>=</span> <span>exp</span><span>(</span>B0 <span>+</span> B1 <span>+</span> <span>(</span>u0<span>^</span><span>2</span> <span>+</span> sigma<span>^</span><span>2</span><span>)</span><span>/</span><span>2</span><span>),</span>
         TE_cond_ATE <span>=</span> y_cond_tx <span>-</span> y_cond_cc<span>,</span>
         TE_marg_ATE <span>=</span> y_marg_tx <span>-</span> y_marg_cc<span>,</span>
         TE_cond_multi <span>=</span> y_cond_tx<span>/</span>y_cond_cc<span>,</span>
         TE_marg_multi <span>=</span> y_marg_tx<span>/</span>y_marg_cc<span>,</span>
         TE_multi <span>=</span> <span>exp</span><span>(</span>B1<span>)</span>
         <span>)</span> <span>%>%</span> 
  select<span>(</span><span>-</span>B0<span>,</span> <span>-</span>B1<span>,</span> <span>-</span>u0<span>,</span> <span>-</span>sigma<span>)</span>

tab <span>&lt;-</span> res <span>%>%</span> 
  gather<span>(</span>key<span>,</span> val<span>)</span> <span>%>%</span> 
  mutate<span>(</span>key <span>=</span> <span>factor</span><span>(</span>key<span>,</span> levels <span>=</span> <span>unique</span><span>(</span>key<span>)))</span> <span>%>%</span> 
  group_by<span>(</span>key<span>)</span> <span>%>%</span> 
  do<span>(</span>  <span>data.frame</span><span>(</span>posterior_median <span>=</span> median<span>(</span><span>.</span><span>$</span>val<span>),</span> 
                  l_95_CI <span>=</span> quantile<span>(</span><span>.</span><span>$</span>val<span>,</span> <span>0.025</span><span>),</span> 
                  u_95_CI <span>=</span> quantile<span>(</span><span>.</span><span>$</span>val<span>,</span> <span>0.975</span><span>),</span>
                  stringsAsFactors <span>=</span> <span>FALSE</span><span>))</span> <span>%>%</span> 
  ungroup<span>()</span> <span>%>%</span> 
  mutate<span>(</span>key <span>=</span> <span>c</span><span>(</span><span>&quot;$\\text{E}(Y_{ij} \\mid u_i = 0, \\text{TX} = 0)$&quot;</span><span>,</span>
                 <span>&quot;$\\text{E}(Y_{ij} \\mid u_i = 0, \\text{TX} = 1)$&quot;</span><span>,</span>
                 <span>&quot;$\\text{E}(Y_{ij} \\mid \\text{TX} = 0)$&quot;</span><span>,</span>
                 <span>&quot;$\\text{E}(Y_{ij} \\mid \\text{TX} = 1)$&quot;</span><span>,</span>
                 <span>&quot;$\\text{E}(Y_{ij} \\mid u_i = 0, \\text{TX} = 1) - \\\\ \\text{E}(Y_{ij} \\mid u_i = 0, \\text{TX} = 0)$&quot;</span><span>,</span>
                 <span>&quot;$\\text{E}(Y_{ij} \\mid \\text{TX} = 1) - \\\\ \\text{E}(Y_{ij} \\mid \\text{TX} = 0)$&quot;</span><span>,</span>
                 <span>&quot;$\\frac{\\text{E}(Y_{ij} \\mid u_i = 0, \\text{TX} = 1)}{\\text{E}(Y_{ij} \\mid u_i = 0, \\text{TX} = 0)}$&quot;</span><span>,</span>
                 <span>&quot;$\\frac{\\text{E}(Y_{ij} \\mid \\text{TX} = 1)}{\\text{E}(Y_{ij} \\mid \\text{TX} = 0)}$&quot;</span><span>,</span>
                 <span>&quot;$\\text{exp}(\\beta_1)$&quot;</span><span>),</span>
         Effect <span>=</span> <span>c</span><span>(</span><span>&quot;Average outcome for a median cluster in the control group&quot;</span><span>,</span>
                     <span>&quot;Average outcome for a median cluster in the treatment group&quot;</span><span>,</span>
                     <span>&quot;Average outcome in the control group&quot;</span><span>,</span>
                     <span>&quot;Average outcome in the treatment group&quot;</span><span>,</span>
                     <span>&quot;Difference between the average outcome in a median cluster in the treatment versus control group&quot;</span><span>,</span>
                     <span>&quot;Difference between the average outcome in the treatment versus control group&quot;</span><span>,</span>
                     <span>&quot;Ratio of the average outcome in a median cluster in the treatment and the control group&quot;</span><span>,</span>
                     <span>&quot;Ratio of the average outcome in the treatment and the control group&quot;</span><span>,</span>
                     <span>&quot;Ratio of the average outcome in the treatment and the control group&quot;</span><span>))</span> <span>%>%</span> 
  select<span>(</span>Effect<span>,</span> Formula <span>=</span> key<span>,</span> posterior_median<span>,</span> l_95_CI<span>,</span> u_95_CI<span>)</span>
</pre>
</div>
</div>
<div>
<table>
<thead>
<tr>
<th align="left">Effect</th>
<th align="left">Formula</th>
<th align="right">posterior_median</th>
<th align="right">l_95_CI</th>
<th align="right">u_95_CI</th>
</tr>
</thead>
<tbody>
<tr>
<td align="left">Average outcome for a median cluster in the control group</td>
<td align="left"><span>\(\text{E}(Y_{ij} \mid u_i = 0, \text{TX} = 0)\)</span></td>
<td align="right">570.5</td>
<td align="right">374.5</td>
<td align="right">847.9</td>
</tr>
<tr>
<td align="left">Average outcome for a median cluster in the treatment group</td>
<td align="left"><span>\(\text{E}(Y_{ij} \mid u_i = 0, \text{TX} = 1)\)</span></td>
<td align="right">281.8</td>
<td align="right">186.7</td>
<td align="right">435.2</td>
</tr>
<tr>
<td align="left">Average outcome in the control group</td>
<td align="left"><span>\(\text{E}(Y_{ij} \mid \text{TX} = 0)\)</span></td>
<td align="right">697.9</td>
<td align="right">454.5</td>
<td align="right">1123.4</td>
</tr>
<tr>
<td align="left">Average outcome in the treatment group</td>
<td align="left"><span>\(\text{E}(Y_{ij} \mid \text{TX} = 1)\)</span></td>
<td align="right">342.7</td>
<td align="right">231.3</td>
<td align="right">575.5</td>
</tr>
<tr>
<td align="left">Difference between the average outcome in a median cluster in the treatment versus control group</td>
<td align="left"><span>\(\text{E}(Y_{ij} \mid u_i = 0, \text{TX} = 1) – \\ \text{E}(Y_{ij} \mid u_i = 0, \text{TX} = 0)\)</span></td>
<td align="right">-286.1</td>
<td align="right">-580.2</td>
<td align="right">-37.9</td>
</tr>
<tr>
<td align="left">Difference between the average outcome in the treatment versus control group</td>
<td align="left"><span>\(\text{E}(Y_{ij} \mid \text{TX} = 1) – \\ \text{E}(Y_{ij} \mid \text{TX} = 0)\)</span></td>
<td align="right">-349.0</td>
<td align="right">-747.5</td>
<td align="right">-47.4</td>
</tr>
<tr>
<td align="left">Ratio of the average outcome in a median cluster in the treatment and the control group</td>
<td align="left"><span>\(\frac{\text{E}(Y_{ij} \mid u_i = 0, \text{TX} = 1)}{\text{E}(Y_{ij} \mid u_i = 0, \text{TX} = 0)}\)</span></td>
<td align="right">0.5</td>
<td align="right">0.3</td>
<td align="right">0.9</td>
</tr>
<tr>
<td align="left">Ratio of the average outcome in the treatment and the control group</td>
<td align="left"><span>\(\frac{\text{E}(Y_{ij} \mid \text{TX} = 1)}{\text{E}(Y_{ij} \mid \text{TX} = 0)}\)</span></td>
<td align="right">0.5</td>
<td align="right">0.3</td>
<td align="right">0.9</td>
</tr>
<tr>
<td align="left">Ratio of the average outcome in the treatment and the control group</td>
<td align="left"><span>\(\text{exp}(\beta_1)\)</span></td>
<td align="right">0.5</td>
<td align="right">0.3</td>
<td align="right">0.9</td>
</tr>
</tbody>
</table>
</div>
<h3>Using brms::fitted</h3>
<p>Let’s compare our results to predicted means we get with <code>brms::fitted</code> and <code>brms::marginal_effects</code></p>
<div>
<pre><span></span>##      Estimate Est.Error     Q2.5    Q97.5
## [1,] 570.5069  114.8815 374.4601 847.8781
</pre>
</div>
<div>
<pre><span></span>##   TX        y cluster cond__ estimate__     se__  lower__  upper__
## 1  0 481.4254      NA      1   570.5069 114.8815 374.4601 847.8781
</pre>
</div>
<p>We see that this is exactly the same as the median outcome (cluster-specific outcome) in the control group. I’ve seen some people confused about this since the function is called <code>marginal_effects</code>. </p>
<h3>Limiting inference to only the included clusters</h3>
<p>We can also use <code>fitted</code> to average over the included clusters, this is very similar to fitting the clusters as “fixed effects”. So let’s compare averaging over the included clusters in the multilevel model to fitting a fixed effects (no pooling) model. </p>
<div>
<pre><span></span>fit_fixed <span>&lt;-</span> brm<span>(</span>y <span>~</span> <span>0</span> <span>+</span> <span>factor</span><span>(</span>cluster<span>),</span> 
                    family <span>=</span> lognormal<span>(),</span> 
                    data <span>=</span> d<span>,</span>
                    cores <span>=</span> <span>4</span><span>,</span> 
                    iter <span>=</span> <span>2000</span><span>)</span>
</pre>
</div>
<p>Let’s look at the predicted mean for the control group using the fixed effects model.</p>
<div>
<pre><span></span>cluster_fixed <span>&lt;-</span> fitted<span>(</span>fit_fixed<span>,</span> 
                         newdata <span>=</span> <span>expand.grid</span><span>(</span>TX <span>=</span> <span>0</span><span>,</span> 
                                               cluster <span>=</span> <span>1</span><span>:</span><span>10</span><span>),</span> 
                         re_formula <span>=</span> <span>NULL</span><span>,</span>
                         summary <span>=</span> <span>FALSE</span><span>)</span>
cluster_fixed <span>&lt;-</span> <span>apply</span><span>(</span>cluster_fixed<span>,</span> <span>1</span><span>,</span> <span>mean</span><span>)</span>
quantile<span>(</span>cluster_fixed<span>,</span> <span>c</span><span>(</span><span>0.025</span><span>,</span> <span>0.5</span><span>,</span> <span>0.975</span><span>))</span>
</pre>
</div>
<div>
<pre><span></span>##     2.5%      50%    97.5% 
## 600.4927 638.8134 682.9103
</pre>
</div>
<p>Which is really similar to the control groups sample mean which is 636. Now let’s compare the posterior distribution from the fixed effects model, to using our multilevel model to average over only the included clusters (and not the estimated distribution of clusters).</p>
<p>

</p>
<div>
<div>
<pre><span></span>get_fixed_cluster_marginal <span>&lt;-</span> <span>function</span><span>(</span>fit<span>,</span> newdat<span>,</span> label<span>)</span> <span>{</span>
  cluster_fixed_ <span>&lt;-</span> fitted<span>(</span>fit<span>,</span> 
                             newdata <span>=</span> newdat<span>,</span> 
                             re_formula <span>=</span> <span>NULL</span><span>,</span>
                             summary <span>=</span> <span>FALSE</span><span>)</span>
  cluster_fixed <span>&lt;-</span> <span>apply</span><span>(</span>cluster_fixed_<span>,</span> <span>1</span><span>,</span> <span>mean</span><span>)</span>

  <span># posterior samples</span>
  <span>list</span><span>(</span>ss_fixed <span>=</span> cluster_fixed_<span>,</span> 
       ss_marg <span>=</span> cluster_fixed<span>,</span>
       label <span>=</span> label<span>)</span>
<span>}</span>

<span># Get prediction for each treatment group</span>
fixed_cluster_C_RE <span>&lt;-</span> get_fixed_cluster_marginal<span>(</span>fit <span>=</span> fit<span>,</span> 
                                              newdat <span>=</span> <span>expand.grid</span><span>(</span>TX <span>=</span> <span>0</span><span>,</span> 
                                                                   cluster <span>=</span> <span>1</span><span>:</span><span>10</span><span>),</span>
                                              label <span>=</span> <span>&quot;C_RE&quot;</span><span>)</span> 
fixed_cluster_TX_RE <span>&lt;-</span> get_fixed_cluster_marginal<span>(</span>fit <span>=</span> fit<span>,</span> 
                                              newdat <span>=</span> <span>expand.grid</span><span>(</span>TX <span>=</span> <span>1</span><span>,</span> 
                                                                   cluster <span>=</span> <span>11</span><span>:</span><span>20</span><span>),</span>
                                              label <span>=</span> <span>&quot;TX_RE&quot;</span><span>)</span> 



fixed_cluster_C_FE <span>&lt;-</span> get_fixed_cluster_marginal<span>(</span>fit <span>=</span> fit_fixed<span>,</span> 
                                                 newdat <span>=</span> <span>expand.grid</span><span>(</span>TX <span>=</span> <span>0</span><span>,</span> 
                                                                      cluster <span>=</span> <span>1</span><span>:</span><span>10</span><span>),</span>
                                                 label <span>=</span> <span>&quot;C_FE&quot;</span><span>)</span> 
fixed_cluster_TX_FE <span>&lt;-</span> get_fixed_cluster_marginal<span>(</span>fit <span>=</span> fit_fixed<span>,</span> 
                                                 newdat <span>=</span> <span>expand.grid</span><span>(</span>TX <span>=</span> <span>1</span><span>,</span> 
                                                                      cluster <span>=</span> <span>11</span><span>:</span><span>20</span><span>),</span>
                                                 label <span>=</span> <span>&quot;TX_FE&quot;</span><span>)</span> 

<span># Summarise posterior median and 95 % CI</span>
sums <span>&lt;-</span> <span>data.frame</span><span>(</span><span>&quot;model&quot;</span> <span>=</span> <span>rep</span><span>(</span><span>c</span><span>(</span><span>&quot;RE&quot;</span><span>,</span> 
                                   <span>&quot;FE&quot;</span><span>),</span> <span>3</span><span>),</span>
                   <span>&quot;variable&quot;</span> <span>=</span> <span>rep</span><span>(</span><span>c</span><span>(</span><span>&quot;Control&quot;</span><span>,</span> 
                                      <span>&quot;Treatment&quot;</span><span>,</span> 
                                      <span>&quot;Difference&quot;</span><span>),</span> each <span>=</span> <span>2</span><span>))</span>  
qs <span>&lt;-</span> <span>rbind</span><span>(</span>quantile<span>(</span>fixed_cluster_C_RE<span>$</span>ss_marg<span>,</span> <span>c</span><span>(</span><span>0.025</span><span>,</span> <span>0.5</span><span>,</span> <span>0.975</span><span>)),</span>
            quantile<span>(</span>fixed_cluster_C_FE<span>$</span>ss_marg<span>,</span> <span>c</span><span>(</span><span>0.025</span><span>,</span> <span>0.5</span><span>,</span> <span>0.975</span><span>)),</span>
            quantile<span>(</span>fixed_cluster_TX_RE<span>$</span>ss_marg<span>,</span> <span>c</span><span>(</span><span>0.025</span><span>,</span> <span>0.5</span><span>,</span> <span>0.975</span><span>)),</span>
            quantile<span>(</span>fixed_cluster_TX_FE<span>$</span>ss_marg<span>,</span> <span>c</span><span>(</span><span>0.025</span><span>,</span> <span>0.5</span><span>,</span> <span>0.975</span><span>)),</span>
            quantile<span>(</span>fixed_cluster_TX_RE<span>$</span>ss_marg <span>-</span> fixed_cluster_C_RE<span>$</span>ss_marg<span>,</span> <span>c</span><span>(</span><span>0.025</span><span>,</span> <span>0.5</span><span>,</span> <span>0.975</span><span>)),</span>
            quantile<span>(</span>fixed_cluster_TX_FE<span>$</span>ss_marg <span>-</span> fixed_cluster_C_FE<span>$</span>ss_marg<span>,</span> <span>c</span><span>(</span><span>0.025</span><span>,</span> <span>0.5</span><span>,</span> <span>0.975</span><span>)))</span>
<span>colnames</span><span>(</span>qs<span>)</span> <span>&lt;-</span> <span>c</span><span>(</span><span>&quot;l_95_CI&quot;</span><span>,</span> <span>&quot;median&quot;</span><span>,</span> <span>&quot;u_95_CI&quot;</span><span>)</span>

sums <span>&lt;-</span> <span>cbind</span><span>(</span>sums<span>,</span> qs<span>)</span>


ps <span>&lt;-</span> <span>rbind</span><span>(</span><span>data.frame</span><span>(</span><span>&quot;model&quot;</span> <span>=</span> <span>&quot;RE&quot;</span><span>,</span> variable <span>=</span> <span>&quot;Control&quot;</span><span>,</span> samples <span>=</span> fixed_cluster_C_RE<span>$</span>ss_marg<span>),</span>
            <span>data.frame</span><span>(</span><span>&quot;model&quot;</span> <span>=</span> <span>&quot;RE&quot;</span><span>,</span> variable <span>=</span> <span>&quot;Treatment&quot;</span><span>,</span> samples <span>=</span> fixed_cluster_TX_RE<span>$</span>ss_marg<span>),</span>
            <span>data.frame</span><span>(</span><span>&quot;model&quot;</span> <span>=</span> <span>&quot;RE&quot;</span><span>,</span> variable <span>=</span> <span>&quot;Difference&quot;</span><span>,</span> samples <span>=</span> fixed_cluster_TX_RE<span>$</span>ss_marg <span>-</span> fixed_cluster_C_RE<span>$</span>ss_marg<span>),</span>
            <span>data.frame</span><span>(</span><span>&quot;model&quot;</span> <span>=</span> <span>&quot;FE&quot;</span><span>,</span> variable <span>=</span> <span>&quot;Control&quot;</span><span>,</span> samples <span>=</span> fixed_cluster_C_FE<span>$</span>ss_marg<span>),</span>
            <span>data.frame</span><span>(</span><span>&quot;model&quot;</span> <span>=</span> <span>&quot;FE&quot;</span><span>,</span> variable <span>=</span> <span>&quot;Treatment&quot;</span><span>,</span> samples <span>=</span> fixed_cluster_TX_FE<span>$</span>ss_marg<span>),</span>
            <span>data.frame</span><span>(</span><span>&quot;model&quot;</span> <span>=</span> <span>&quot;FE&quot;</span><span>,</span> variable <span>=</span> <span>&quot;Difference&quot;</span><span>,</span> samples <span>=</span> fixed_cluster_TX_FE<span>$</span>ss_marg <span>-</span> fixed_cluster_C_FE<span>$</span>ss_marg<span>))</span>

<span># plot</span>
p <span>&lt;-</span> ggplot<span>(</span>ps<span>,</span> aes<span>(</span>x <span>=</span> samples<span>,</span> 
               color <span>=</span> model<span>,</span> 
               fill <span>=</span> model<span>))</span> <span>+</span> 
  geom_density<span>(</span>alpha <span>=</span> <span>0.25</span><span>,</span> aes<span>(</span>y <span>=</span> <span>..</span>scaled..<span>))</span> <span>+</span> 
  geom_segment<span>(</span>data <span>=</span> sums<span>,</span>
               aes<span>(</span>x <span>=</span> l_95_CI<span>,</span>
                   xend <span>=</span> u_95_CI<span>,</span>
                   y <span>=</span> <span>0</span><span>,</span> 
                   yend <span>=</span> <span>0</span><span>),</span> 
               position <span>=</span> position_nudge<span>(</span>y <span>=</span> <span>c</span><span>(</span><span>-0.025</span><span>,</span> <span>-0.05</span><span>)))</span> <span>+</span>
  geom_point<span>(</span>data <span>=</span> sums<span>,</span>
             aes<span>(</span>x <span>=</span> median<span>,</span>
                 y <span>=</span> <span>0</span><span>),</span> 
             position <span>=</span> position_nudge<span>(</span>y <span>=</span> <span>c</span><span>(</span><span>-0.025</span><span>,</span> <span>-0.05</span><span>)))</span> <span>+</span>
  facet_wrap<span>(</span><span>~</span>variable<span>,</span> scales <span>=</span> <span>&quot;free_x&quot;</span><span>)</span> <span>+</span>
  labs<span>(</span>x <span>=</span> <span>&quot;Posterior samples&quot;</span><span>,</span>
       y <span>=</span> <span>&quot;Cluster mean ($)&quot;</span><span>,</span>
       title <span>=</span> <span>&quot;Treating clusters as \&quot;fixed\&quot;&quot;</span><span>,</span>
       subtitle <span>=</span> <span>&quot;Comparing no pooling (FE) vs partial pooling (RE)\nwhen both models limit inferences to the included clusters only&quot;</span>
       <span>)</span> <span>+</span>
  theme_minimal<span>()</span> <span>+</span>
  theme<span>(</span>legend.position <span>=</span> <span>&quot;bottom&quot;</span><span>)</span>
</pre>
</div>
</div>
<p><img alt="Treating clusters as 'fixed' in a 'fixed effects' model versus a multilevel model" title="Treating clusters as 'fixed' in a 'fixed effects' model versus a multilevel model" src="https://i1.wp.com/rpsychologist.com/img/GLMM-part1-log-transform/fixed-effects-cluster-multilevel-1.png?w=456" data-recalc-dims="1"></p>
<p>We see that the results are highly similar. Although, the models behind the results do differ. If we look at the estimates for the specific clusters we can see the—in this model very weak—effect of partial pooling in the multilevel model.</p>
<p><img alt="No pooling (fixed effects) versus partial pooling (multilevel) for the cluster means" title="No pooling (fixed effects) versus partial pooling (multilevel) for the cluster means" src="https://i1.wp.com/rpsychologist.com/img/GLMM-part1-log-transform/partial-vs-no-pooling-1.png?w=456" data-recalc-dims="1"></p>
<h2>Intraclass correlations (ICC)</h2>
<p>Since the variance components are on the log scale ICCs calculated using them will also be on the log scale. This is not always what we want, for instance, if we are calculating agreement, reliability, or heritability, it can make a lot of sense to calculate ICCs on the data scale. If I am calculating ICCs for family members’ reports on gambling losses, then I am interested in agreement on the data scale (USD) and not on the log scale.</p>
<p>ICCs on the data scale can be calculated using the same techniques we used for the marginalization of the outcomes. You can read more about this in the context of GLMMs in Villemereuil et al. (2016; and the related R package <code>QCglmm</code>). Here is an example, first we calculate the marginal mean</p>
<div>$$ \text{E}(y_{i} \mid TX_i) = \bar\mu_i = \int \text{exp}(x + \sigma_e^2/2) f_{\mathcal{N}}(x, \mu_{i}, \sigma_{u_0}^2)dx$$</div>
<p>
then we use the marginal mean to calculate the expected variance between clusters on the data scale,
</p>
<div>$$ \sigma_B^2 = \int (e^x – \bar\mu)^2 f_{\mathcal{N}}(x, \mu_{i}, \sigma_{u_0^2})dx$$</div>
<p>and lastly we calculate the expected variance within clusters on the data scale,</p>
<div>$$ \sigma_W^2 = \int \left(( \text{exp}(\sigma_e^2) – 1 ) \text{exp}(2x + \sigma_e^2) \right)f_{\mathcal{N}}(x, \mu_i, \sigma_{u_0^2})dx.$$</div>
<p>Where <span>\((\text{exp}(\sigma_e^2) – 1 ) \text{exp}(2x + \sigma_e^2)\)</span> is the variance of the lognormal distribution evaluated at <span>\(x\)</span>. </p>
<p>ICC on response scale in then
</p>
<div>$$ \text{Cor}(y_{ij}, y_{ij’}) = \frac{\sigma_B^2}{\sigma_B^2 + \sigma_W^2} $$</div>
<p>
Whereas, the log scale ICC is,
</p>
<div>$$ \text{Cor}\left(\text{log}(y_{ij}), \text{log}(y_{ij’})\right) = \frac{\sigma_u^2}{\sigma_u^2 + \sigma_e^2} $$</div>
<p>
Since this is also a one-dimensional integral we can use <code>integrate</code>,</p>
<p>

</p>

<div>
<div>
<pre><span></span>get_ICC_ln <span>&lt;-</span> <span>function</span><span>(</span>B0<span>,</span> u0<span>,</span> sd_log<span>,</span> <span>...</span><span>)</span> <span>{</span>

  mu <span>&lt;-</span>  integrate<span>(</span>
    f <span>=</span> <span>function</span><span>(</span>x<span>)</span> <span>{</span>
      <span>exp</span><span>(</span>x <span>+</span> sd_log<span>^</span><span>2</span> <span>/</span> <span>2</span><span>)</span> <span>*</span> dnorm<span>(</span>x<span>,</span> B0<span>,</span> sd <span>=</span> u0<span>)</span>
    <span>},</span>
    lower <span>=</span> B0 <span>-</span> <span>10</span> <span>*</span> u0<span>,</span>
    upper <span>=</span> B0 <span>+</span> <span>10</span> <span>*</span> u0
  <span>)</span><span>$</span>value

  sigma2_bw <span>&lt;-</span>  integrate<span>(</span>
    f <span>=</span> <span>function</span><span>(</span>x<span>)</span> <span>{</span>
      <span>(</span><span>exp</span><span>(</span>x <span>+</span> sd_log<span>^</span><span>2</span> <span>/</span> <span>2</span><span>)</span> <span>-</span> mu<span>)</span><span>^</span><span>2</span>  <span>*</span> dnorm<span>(</span>x<span>,</span> B0<span>,</span> sd <span>=</span> u0<span>)</span>
    <span>},</span>
    lower <span>=</span> B0 <span>-</span> <span>10</span> <span>*</span> u0<span>,</span>
    upper <span>=</span> B0 <span>+</span> <span>10</span> <span>*</span> u0
  <span>)</span><span>$</span>value

  sigma2_wi <span>&lt;-</span>  integrate<span>(</span>
    f <span>=</span> <span>function</span><span>(</span>x<span>)</span> <span>{</span>
      <span>(</span><span>exp</span><span>(</span>sd_log<span>^</span><span>2</span><span>)</span> <span>-</span> <span>1</span><span>)</span> <span>*</span> <span>exp</span><span>(</span><span>2</span> <span>*</span> x <span>+</span> sd_log<span>^</span><span>2</span><span>)</span>  <span>*</span> dnorm<span>(</span>x<span>,</span> B0<span>,</span> sd <span>=</span> u0<span>)</span>
    <span>},</span>
    lower <span>=</span> B0 <span>-</span> <span>10</span> <span>*</span> u0<span>,</span>
    upper <span>=</span> B0 <span>+</span> <span>10</span> <span>*</span> u0
  <span>)</span><span>$</span>value

  ICC <span>&lt;-</span>  sigma2_bw <span>/</span> <span>(</span>sigma2_bw <span>+</span> sigma2_wi<span>)</span>
  <span>list</span><span>(</span>
    <span>&quot;ICC&quot;</span> <span>=</span> ICC<span>,</span>
    <span>&quot;WI_sd&quot;</span> <span>=</span> <span>sqrt</span><span>(</span>sigma2_wi<span>),</span>
    <span>&quot;BW_sd&quot;</span> <span>=</span> <span>sqrt</span><span>(</span>sigma2_bw<span>),</span>
    <span>&quot;mean&quot;</span> <span>=</span> mu
  <span>)</span>
<span>}</span>
</pre>
</div>
</div>
<p>If you find code easier to read, the Monte Carlo version would be</p>
<p>

</p>

<div>
<div>
<pre><span></span>get_ICC_ln_sim <span>&lt;-</span> <span>function</span><span>(</span>B0<span>,</span> u0<span>,</span> sd_log<span>,</span> R <span>=</span> <span>1e4</span><span>)</span> <span>{</span>
  RE <span>&lt;-</span> B0 <span>+</span> rnorm<span>(</span>R<span>,</span> <span>0</span><span>,</span> u0<span>)</span>
  y <span>&lt;-</span> <span>exp</span><span>(</span>RE <span>+</span> sd_log<span>^</span><span>2</span><span>/</span><span>2</span><span>)</span>

  sigma2_bw <span>&lt;-</span> var<span>(</span>y<span>)</span> 
  sigma2_wi <span>&lt;-</span> <span>mean</span><span>(</span> <span>(</span><span>exp</span><span>(</span>sd_log<span>^</span><span>2</span><span>)</span> <span>-</span> <span>1</span><span>)</span> <span>*</span> <span>exp</span><span>(</span><span>2</span> <span>*</span> RE <span>+</span> sd_log<span>^</span><span>2</span><span>))</span>

  ICC <span>&lt;-</span> sigma2_bw <span>/</span> <span>(</span>sigma2_bw <span>+</span> sigma2_wi<span>)</span>

  <span>list</span><span>(</span><span>&quot;ICC&quot;</span> <span>=</span> ICC<span>,</span>
       <span>&quot;WI_sd&quot;</span> <span>=</span> <span>sqrt</span><span>(</span>sigma2_wi<span>),</span>
       <span>&quot;BW_sd&quot;</span> <span>=</span> <span>sqrt</span><span>(</span>sigma2_bw<span>),</span>
       <span>&quot;mean&quot;</span> <span>=</span> <span>mean</span><span>(</span>y<span>))</span>
<span>}</span>
</pre>
</div>
</div>
<p>Let’s see if the two methods agree.</p>
<p>

</p>

<div>
<div>
<pre><span></span>ICCS <span>&lt;-</span> <span>list</span><span>(</span>get_ICC_ln<span>(</span>B0<span>,</span> 
                        u0<span>,</span> 
                        sd_log<span>),</span> 
             get_ICC_ln_sim<span>(</span>B0<span>,</span>
                            u0<span>,</span> 
                            sd_log<span>,</span>
                            R <span>=</span> <span>1e6</span><span>),</span>
             <span>data.frame</span><span>(</span>ICC <span>=</span> u0<span>^</span><span>2</span><span>/</span><span>(</span>u0<span>^</span><span>2</span> <span>+</span> sd_log<span>^</span><span>2</span><span>),</span>
                        WI_sd <span>=</span> <span>NA</span><span>,</span> 
                        BW_sd <span>=</span> <span>NA</span><span>,</span> 
                        mean <span>=</span> <span>NA</span><span>)</span>
             <span>)</span>

res <span>&lt;-</span> <span>lapply</span><span>(</span>ICCS<span>,</span> <span>data.frame</span><span>)</span> <span>%>%</span> 
  <span>do.call</span><span>(</span><span>rbind</span><span>,</span> <span>.</span><span>)</span>
</pre>
</div>
</div>
<table>
<thead>
<tr>
<th align="left"></th>
<th align="right">Data scale (Integrate)</th>
<th align="right">Data scale (Simulate)</th>
<th align="right">Log scale</th>
</tr>
</thead>
<tbody>
<tr>
<td align="left">ICC</td>
<td align="right">0.44</td>
<td align="right">0.44</td>
<td align="right">0.5</td>
</tr>
<tr>
<td align="left">WI_sd</td>
<td align="right">387.71</td>
<td align="right">387.69</td>
<td align="right">NA</td>
</tr>
<tr>
<td align="left">BW_sd</td>
<td align="right">342.15</td>
<td align="right">342.22</td>
<td align="right">NA</td>
</tr>
<tr>
<td align="left">mean</td>
<td align="right">642.01</td>
<td align="right">641.93</td>
<td align="right">NA</td>
</tr>
</tbody>
</table>
<p>Great! We also see that by calculating variances on the data scale the ICCs are lower, in this case, the difference is small, but if the log variance is large this difference can be <em>very</em> large. So when you see ICCs from lognormal models or GLMMs you should ask if they are on the transformed (latent) scale or on the observed (data) scale, and if the scale is approriate for the research question.</p>
<h2>Simulation study</h2>
<p>I began by promising that the methods in this post could make your <em>frequentist</em> life easier, and then continued to only fit Bayesian models. So let’s do a small simulation study to evaluate the frequentist properties of this model and all the transformations and marginalizations. I will try two different models, shown in the table below, each model will be evaluated using 5000 replications, and I will fit both the correct model using a lognormal response distribution and a standard LMM using the Gaussian response distribution.</p>
<table>
<thead>
<tr>
<th>Parameter</th>
<th>Sim1: less skew</th>
<th>Sim2: more skew</th>
</tr>
</thead>
<tbody>
<tr>
<td><span>\(\beta_0\)</span></td>
<td>log(500)</td>
<td>log(500)</td>
</tr>
<tr>
<td><span>\(\beta_1\)</span></td>
<td>log(0.5)</td>
<td>log(0.5)</td>
</tr>
<tr>
<td><span>\(\sigma_u\)</span></td>
<td>0.5</td>
<td>1</td>
</tr>
<tr>
<td><span>\(\sigma_e\)</span></td>
<td>0.5</td>
<td>1</td>
</tr>
<tr>
<td><span>\(n_1\)</span></td>
<td>30</td>
<td>5</td>
</tr>
<tr>
<td><span>\(n_2\)</span></td>
<td>10</td>
<td>75</td>
</tr>
</tbody>
</table>
<h3>Simulation results</h3>
<p>In the figures below we see that the posterior median show mostly negligible bias for both the original parameters and the transformed marginal effects. The Gaussian model has the most trouble with the variance components, which show substantial bias.<br>
<img alt="Simulation: Relative bias for models with skewed data fit using either a normal or lognormal multilevel model" title="Simulation: Relative bias for models with skewed data fit using either a normal or lognormal multilevel model" src="https://i1.wp.com/rpsychologist.com/img/GLMM-part1-log-transform/sim-skewed-data-normal-vs-lognormal-multilevel-rel-bias-1.png?w=456" data-recalc-dims="1"></p>
<p>If we look at the coverage probabilities we see that the lognormal model’s posterior intervals have nominal coverage rates for all parameters. The Gaussian model has slightly worse coverage for the treatment effects, but not that bad. However, the CI’s from the Gaussian model does a poor job of capturing the true ICC and variance components. </p>
<p><img alt="Simulation: 95 % CI coverage for models with skewed data fit using either a normal or lognormal multilevel model" title="Simulation: 95 % CI coverage for models with skewed data fit using either a normal or lognormal multilevel model" src="https://i1.wp.com/rpsychologist.com/img/GLMM-part1-log-transform/sim-skewed-data-normal-vs-lognormal-multilevel-CI-coverage-1.png?w=456" data-recalc-dims="1"></p>
<p>

</p>
<div>
Normal model
<div>
<table>
<thead>
<tr>
<th> var </th>
<th> model </th>
<th> est_mean </th>
<th> est_median </th>
<th> theta </th>
<th> est_rel_bias </th>
<th> est_median_rel_bias </th>
<th> se </th>
<th> coverage </th>
<th> power </th>
</tr>
</thead>
<tbody>
<tr>
<td> ICC </td>
<td> gaussian </td>
<td> 0.4 </td>
<td> 0.4 </td>
<td> 0.4 </td>
<td> 0.0 </td>
<td> 0.0 </td>
<td> 0.1 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> marg_TE </td>
<td> gaussian </td>
<td> -323.5 </td>
<td> -323.4 </td>
<td> NA </td>
<td> NA </td>
<td> NA </td>
<td> 125.0 </td>
<td> NA </td>
<td> 0.8 </td>
</tr>
<tr>
<td> marg_mean_C </td>
<td> gaussian </td>
<td> 635.6 </td>
<td> 635.9 </td>
<td> 642.0 </td>
<td> 0.0 </td>
<td> 0.0 </td>
<td> 105.2 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> marg_mean_TX </td>
<td> gaussian </td>
<td> 312.1 </td>
<td> 312.4 </td>
<td> 321.0 </td>
<td> 0.0 </td>
<td> 0.0 </td>
<td> 53.6 </td>
<td> 1.0 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> multiplicative_TE_ratio_marg </td>
<td> gaussian </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.0 </td>
<td> 0.0 </td>
<td> 0.1 </td>
<td> 1.0 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> sd_BW_C </td>
<td> gaussian </td>
<td> 267.3 </td>
<td> 261.3 </td>
<td> 342.2 </td>
<td> -0.2 </td>
<td> -0.2 </td>
<td> 84.8 </td>
<td> 0.5 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> sd_BW_TX </td>
<td> gaussian </td>
<td> 267.3 </td>
<td> 261.3 </td>
<td> 171.1 </td>
<td> 0.6 </td>
<td> 0.5 </td>
<td> 84.8 </td>
<td> 0.4 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> sd_WI_C </td>
<td> gaussian </td>
<td> 302.7 </td>
<td> 302.4 </td>
<td> 387.7 </td>
<td> -0.2 </td>
<td> -0.2 </td>
<td> 57.0 </td>
<td> 0.1 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> sd_WI_TX </td>
<td> gaussian </td>
<td> 302.7 </td>
<td> 302.4 </td>
<td> 193.9 </td>
<td> 0.6 </td>
<td> 0.6 </td>
<td> 57.0 </td>
<td> 0.0 </td>
<td> 1.0 </td>
</tr>
</tbody>
</table>
</div>
<p>Lognormal model</p>
<div>
<table>
<thead>
<tr>
<th> var </th>
<th> model </th>
<th> est_mean </th>
<th> est_median </th>
<th> theta </th>
<th> est_rel_bias </th>
<th> est_median_rel_bias </th>
<th> se </th>
<th> coverage </th>
<th> power </th>
</tr>
</thead>
<tbody>
<tr>
<td> b0 </td>
<td> LN </td>
<td> 6.2 </td>
<td> 6.2 </td>
<td> 6.2 </td>
<td> 0.0 </td>
<td> 0.0 </td>
<td> 0.2 </td>
<td> 1.0 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> b1 </td>
<td> LN </td>
<td> -0.7 </td>
<td> -0.7 </td>
<td> -0.7 </td>
<td> 0.0 </td>
<td> 0.0 </td>
<td> 0.2 </td>
<td> 0.9 </td>
<td> 0.8 </td>
</tr>
<tr>
<td> ICC </td>
<td> LN </td>
<td> 0.4 </td>
<td> 0.4 </td>
<td> 0.4 </td>
<td> 0.0 </td>
<td> 0.0 </td>
<td> 0.1 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> ICC_log </td>
<td> LN </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.0 </td>
<td> 0.0 </td>
<td> 0.1 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> marg_TE </td>
<td> LN </td>
<td> -342.1 </td>
<td> -330.4 </td>
<td> NA </td>
<td> NA </td>
<td> NA </td>
<td> 126.7 </td>
<td> NA </td>
<td> 0.8 </td>
</tr>
<tr>
<td> marg_mean_C </td>
<td> LN </td>
<td> 683.1 </td>
<td> 664.0 </td>
<td> 642.0 </td>
<td> 0.1 </td>
<td> 0.0 </td>
<td> 118.8 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> marg_mean_TX </td>
<td> LN </td>
<td> 341.0 </td>
<td> 331.5 </td>
<td> 321.0 </td>
<td> 0.1 </td>
<td> 0.0 </td>
<td> 58.0 </td>
<td> 1.0 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> median_clust_TE </td>
<td> LN </td>
<td> -258.5 </td>
<td> -253.1 </td>
<td> NA </td>
<td> NA </td>
<td> NA </td>
<td> 94.2 </td>
<td> NA </td>
<td> 0.8 </td>
</tr>
<tr>
<td> median_clust_C </td>
<td> LN </td>
<td> 516.2 </td>
<td> 508.3 </td>
<td> 500.0 </td>
<td> 0.0 </td>
<td> 0.0 </td>
<td> 84.4 </td>
<td> 1.0 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> median_clust_TX </td>
<td> LN </td>
<td> 257.7 </td>
<td> 253.7 </td>
<td> 250.0 </td>
<td> 0.0 </td>
<td> 0.0 </td>
<td> 41.3 </td>
<td> 1.0 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> multiplicative_TE_exp_beta1 </td>
<td> LN </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.1 </td>
<td> 0.0 </td>
<td> 0.1 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> multiplicative_TE_ratio_marg </td>
<td> LN </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.1 </td>
<td> 0.0 </td>
<td> 0.1 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> multiplicative_TE_ratio_medians </td>
<td> LN </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.1 </td>
<td> 0.0 </td>
<td> 0.1 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> sd_BW_C </td>
<td> LN </td>
<td> 409.0 </td>
<td> 370.9 </td>
<td> 342.2 </td>
<td> 0.2 </td>
<td> 0.1 </td>
<td> 131.4 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> sd_BW_TX </td>
<td> LN </td>
<td> 204.1 </td>
<td> 185.1 </td>
<td> 171.1 </td>
<td> 0.2 </td>
<td> 0.1 </td>
<td> 65.0 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> sd_log </td>
<td> LN </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.0 </td>
<td> 0.0 </td>
<td> 0.0 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> sd_u0_log </td>
<td> LN </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.5 </td>
<td> 0.1 </td>
<td> 0.0 </td>
<td> 0.1 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> sd_WI_C </td>
<td> LN </td>
<td> 429.6 </td>
<td> 409.3 </td>
<td> 387.7 </td>
<td> 0.1 </td>
<td> 0.1 </td>
<td> 88.8 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
<tr>
<td> sd_WI_TX </td>
<td> LN </td>
<td> 214.4 </td>
<td> 204.3 </td>
<td> 193.9 </td>
<td> 0.1 </td>
<td> 0.1 </td>
<td> 43.5 </td>
<td> 0.9 </td>
<td> 1.0 </td>
</tr>
</tbody>
</table>
</div>
</div>
<h3>Power</h3>
<p>Lastly, let’s take a closer look at the marginalized treatment effect on the data scale. The figure below shows a sample of the confidence intervals for the difference in the average money lost gambling from the two conditions, from Model 2 with more skewed data. We can see that the, although, the normal model is quite robust and the estimated effect is unbiased, it is less precise which leads to reduced power and CIs with worse coverage compared to the true lognormal model.</p>
<p><img alt="lognormal multilevel LMM compared to a normal LMM with skewed data, average treatment effect shown, the effect of skewness on precision and CI coverage" src="https://i0.wp.com/rpsychologist.com/img/GLMM-part1-lognormal/illustrate_CI_ATE.png?w=456" data-recalc-dims="1"></p>
<h1>Summary</h1>
<p>There is a lot of information in this post, so let’s summarize the key take-home messages: </p>
<ul>
<li>Both cluster-specific (conditional) and population-average (marginal) treatment effects are useful—they just answer different questions and we should be careful to not mixed them up.</li>
<li>The multilevel model can be used to give both cluster-specific (conditional) and population-average effects on both the original and transformed scale while retaining all of the advantages of likelihood-based/Bayesian methods.</li>
<li>There is no free lunch. We saw that the calculations of the marginal effects on the data scale use estimates of both the lognormal variance, as well as the random intercept Thus, model misspecification will have a large impact on the appropriateness of the back-transformation. In that sense, directly modeling the population-average effects using e.g. GEE and applying nonparametric retransformations (e.g. smearing) is more robust in some aspects. However, once we add complexities like MNAR missing data in combination with weak adherence (both by patients and therapists…) to complex interventions, then it can be argued that a carefully specified generative model that combines both statistical and clinical expertise will be more likely to give useful answers regarding the treatment effect. In my experience, such models will often have multilevel and cluster-specific components. </li>
<li>For outcomes that are highly meaningful on the original scale (such as expenditure), absolute effects on the data scale can be very meaningful. </li>
<li>Using the default priors in brms resulted in frequentist CI’s with nominal coverage probabilities.</li>
<li>It is useful to understand how the conditional and marginal effects relate to each other, to avoid misinterpretations of the default cluster-specific effects (very common in the clinical studies I read).</li>
</ul>
<h1>Further reading</h1>
<p>The issues covered in this post is dealt with in many articles, here is a selection of relevant articles:</p>
<ul>
<li>Gromping, U. (1996). A Note on Fitting a Marginal Model to Mixed Effects Log-Linear Regression Data via GEE. Biometrics, 52(1), 280–285. https://doi.org/10.2307/2533162</li>
<li>Heagerty, P. J., &amp; Zeger, S. L. (2000). Marginalized multilevel models and likelihood inference (with comments and a rejoinder by the authors). Statistical Science, 15(1), 1–26. https://doi.org/10.1214/ss/1009212671</li>
<li>Hedeker, D., Toit, S. H. C. du, Demirtas, H., &amp; Gibbons, R. D. (2018). A note on marginalization of regression parameters from mixed models of binary outcomes. Biometrics, 74(1), 354–361. https://doi.org/10.1111/biom.12707</li>
<li>Villemereuil, P. de, Schielzeth, H., Nakagawa, S., &amp; Morrissey, M. (2016). General Methods for Evolutionary Quantitative Genetic Inference from Generalized Mixed Models. Genetics, 204(3), 1281–1294. https://doi.org/10.1534/genetics.115.186536</li>
</ul>
<h1>Simulation code</h1>
<p><code>run_sim.R</code>: perform the simulation</p>
<p>

</p>

<div>
<div>
<pre><span></span><span>library</span><span>(</span>tidyr<span>)</span>
<span>library</span><span>(</span>dplyr<span>)</span>
<span>library</span><span>(</span>brms<span>)</span>
<span>library</span><span>(</span>future.apply<span>)</span>
<span>library</span><span>(</span>purrr<span>)</span>

<span>source</span><span>(</span><span>&quot;sim/functions.R&quot;</span><span>)</span>

<span># Parallelization setup</span>
cores <span>&lt;-</span> <span>16</span>
nsim <span>&lt;-</span> <span>5000</span>
plan<span>(</span>tweak<span>(</span>multiprocess<span>,</span> 
           workers <span>=</span> cores<span>))</span> 

<span># ---- model-setup ----</span>
<span>set.seed</span><span>(</span><span>4445</span><span>)</span>
pars <span>&lt;-</span> <span>list</span><span>(</span><span>&quot;n1&quot;</span> <span>=</span> <span>30</span><span>,</span> <span># observations per cluster</span>
             <span>&quot;n2&quot;</span> <span>=</span> <span>10</span><span>,</span> <span># clusters per treatment</span>
             <span>&quot;B0&quot;</span> <span>=</span> <span>log</span><span>(</span><span>500</span><span>),</span>
             <span>&quot;B1&quot;</span> <span>=</span> <span>log</span><span>(</span><span>0.5</span><span>),</span>
             <span>&quot;sd_log&quot;</span> <span>=</span> <span>0.5</span><span>,</span>
             <span>&quot;u0&quot;</span> <span>=</span> <span>0.5</span><span>)</span>

d <span>&lt;-</span> <span>do.call</span><span>(</span>gen_data<span>,</span>
             pars<span>)</span>

<span># ---- compile models ----</span>
fit_LN <span>&lt;-</span> brm<span>(</span>y <span>~</span> <span>1</span> <span>+</span> TX <span>+</span> <span>(</span><span>1</span> <span>|</span> cluster<span>),</span> 
              family <span>=</span> lognormal<span>(),</span> 
              data <span>=</span> d<span>,</span>
              cores <span>=</span> <span>1</span><span>,</span> 
              iter <span>=</span> <span>1</span><span>)</span>

fit_gauss <span>&lt;-</span> brm<span>(</span>y <span>~</span> <span>1</span> <span>+</span> TX <span>+</span> <span>(</span><span>1</span> <span>|</span> cluster<span>),</span> 
                 family <span>=</span> gaussian<span>(),</span> 
                 data <span>=</span> d<span>,</span>
                 cores <span>=</span> <span>1</span><span>,</span> 
                 iter <span>=</span> <span>1</span><span>)</span>

<span># Run ---------------------------------------------------------------------</span>
res1 <span>&lt;-</span> run_sim<span>(</span>nsim <span>=</span> nsim<span>,</span> 
               pars <span>=</span> pars<span>,</span> 
               fit_LN <span>=</span> fit_LN<span>,</span>
               fit_gauss <span>=</span> fit_gauss<span>)</span>

pars <span>&lt;-</span> <span>list</span><span>(</span><span>&quot;n1&quot;</span> <span>=</span> <span>5</span><span>,</span> 
             <span>&quot;n2&quot;</span> <span>=</span> <span>75</span><span>,</span>
             <span>&quot;B0&quot;</span> <span>=</span> <span>log</span><span>(</span><span>500</span><span>),</span>
             <span>&quot;B1&quot;</span> <span>=</span> <span>log</span><span>(</span><span>0.5</span><span>),</span>
             <span>&quot;sd_log&quot;</span> <span>=</span> <span>1</span><span>,</span>
             <span>&quot;u0&quot;</span> <span>=</span> <span>1</span><span>)</span>

res2 <span>&lt;-</span> run_sim<span>(</span>nsim <span>=</span> nsim<span>,</span> 
               pars <span>=</span> pars<span>,</span> 
               fit_LN <span>=</span> fit_LN<span>,</span>
               fit_gauss <span>=</span> fit_gauss<span>)</span>

<span>saveRDS</span><span>(</span>res1<span>,</span> file <span>=</span> <span>&quot;sim_res_GLMM_LN1.rds&quot;</span><span>)</span>
<span>saveRDS</span><span>(</span>res2<span>,</span> file <span>=</span> <span>&quot;sim_res_GLMM_LN2.rds&quot;</span><span>)</span>
</pre>
</div>
</div>
<p><code>summarise_sim.R</code>: Summarize results</p>
<p>

</p>

<div>
<div>
<pre><span></span><span>library</span><span>(</span>tidyr<span>)</span>
<span>library</span><span>(</span>dplyr<span>)</span>
<span>library</span><span>(</span>ggplot2<span>)</span>
<span>source</span><span>(</span><span>&quot;sim/functions.R&quot;</span><span>)</span>

res1 <span>&lt;-</span> <span>readRDS</span><span>(</span><span>&quot;sim_res_GLMM_LN1.rds&quot;</span><span>)</span>
res2 <span>&lt;-</span> <span>readRDS</span><span>(</span><span>&quot;sim_res_GLMM_LN2.rds&quot;</span><span>)</span>

<span># summary</span>
res_sum1 <span>&lt;-</span> summarise_sim<span>(</span>res1<span>)</span>
res_sum2 <span>&lt;-</span> summarise_sim<span>(</span>res2<span>)</span>
</pre>
</div>
</div>
<p><code>functions.R</code>: Functions used</p>
<p>

</p>

<div>
<div>
<pre><span></span><span># ---- gen-data ----</span>

<span>#' Generate lognormal data with a random intercept</span>
<span>#'</span>
<span>#' @param n1 patients per cluster</span>
<span>#' @param n2 clusters per treatment</span>
<span>#' @param B0 log intercept</span>
<span>#' @param B1 log treatment effect</span>
<span>#' @param sd_log log sd</span>
<span>#' @param u0 SD of log intercepts (random intercept)</span>
<span>#'</span>
<span>#' @return a data.frame</span>
gen_data <span>&lt;-</span> <span>function</span><span>(</span>n1<span>,</span> n2<span>,</span> B0<span>,</span> B1<span>,</span> sd_log<span>,</span> u0<span>)</span> <span>{</span>

  cluster <span>&lt;-</span> <span>rep</span><span>(</span><span>1</span><span>:</span><span>(</span><span>2</span> <span>*</span> n2<span>),</span> each <span>=</span> n1<span>)</span>
  TX <span>&lt;-</span> <span>rep</span><span>(</span><span>c</span><span>(</span><span>0</span><span>,</span> <span>1</span><span>),</span> each <span>=</span> n1 <span>*</span> n2<span>)</span>
  u0 <span>&lt;-</span> rnorm<span>(</span><span>2</span> <span>*</span> n2<span>,</span> sd <span>=</span> u0<span>)[</span>cluster<span>]</span>

  mulog <span>&lt;-</span> <span>(</span>B0 <span>+</span> B1 <span>*</span> TX <span>+</span> u0<span>)</span>
  y <span>&lt;-</span> rlnorm<span>(</span><span>2</span> <span>*</span> n1 <span>*</span> n2<span>,</span> meanlog <span>=</span> mulog<span>,</span> sdlog <span>=</span> sd_log<span>)</span>

  d <span>&lt;-</span> <span>data.frame</span><span>(</span>cluster<span>,</span>
                  TX<span>,</span>
                  y<span>)</span>
  d
<span>}</span>


<span># ---- icc-integrate ----</span>
get_ICC_ln <span>&lt;-</span> <span>function</span><span>(</span>B0<span>,</span> u0<span>,</span> sd_log<span>,</span> <span>...</span><span>)</span> <span>{</span>

  mu <span>&lt;-</span>  integrate<span>(</span>
    f <span>=</span> <span>function</span><span>(</span>x<span>)</span> <span>{</span>
      <span>exp</span><span>(</span>x <span>+</span> sd_log<span>^</span><span>2</span> <span>/</span> <span>2</span><span>)</span> <span>*</span> dnorm<span>(</span>x<span>,</span> B0<span>,</span> sd <span>=</span> u0<span>)</span>
    <span>},</span>
    lower <span>=</span> B0 <span>-</span> <span>10</span> <span>*</span> u0<span>,</span>
    upper <span>=</span> B0 <span>+</span> <span>10</span> <span>*</span> u0
  <span>)</span><span>$</span>value

  sigma2_bw <span>&lt;-</span>  integrate<span>(</span>
    f <span>=</span> <span>function</span><span>(</span>x<span>)</span> <span>{</span>
      <span>(</span><span>exp</span><span>(</span>x <span>+</span> sd_log<span>^</span><span>2</span> <span>/</span> <span>2</span><span>)</span> <span>-</span> mu<span>)</span><span>^</span><span>2</span>  <span>*</span> dnorm<span>(</span>x<span>,</span> B0<span>,</span> sd <span>=</span> u0<span>)</span>
    <span>},</span>
    lower <span>=</span> B0 <span>-</span> <span>10</span> <span>*</span> u0<span>,</span>
    upper <span>=</span> B0 <span>+</span> <span>10</span> <span>*</span> u0
  <span>)</span><span>$</span>value

  sigma2_wi <span>&lt;-</span>  integrate<span>(</span>
    f <span>=</span> <span>function</span><span>(</span>x<span>)</span> <span>{</span>
      <span>(</span><span>exp</span><span>(</span>sd_log<span>^</span><span>2</span><span>)</span> <span>-</span> <span>1</span><span>)</span> <span>*</span> <span>exp</span><span>(</span><span>2</span> <span>*</span> x <span>+</span> sd_log<span>^</span><span>2</span><span>)</span>  <span>*</span> dnorm<span>(</span>x<span>,</span> B0<span>,</span> sd <span>=</span> u0<span>)</span>
    <span>},</span>
    lower <span>=</span> B0 <span>-</span> <span>10</span> <span>*</span> u0<span>,</span>
    upper <span>=</span> B0 <span>+</span> <span>10</span> <span>*</span> u0
  <span>)</span><span>$</span>value

  ICC <span>&lt;-</span>  sigma2_bw <span>/</span> <span>(</span>sigma2_bw <span>+</span> sigma2_wi<span>)</span>
  <span>list</span><span>(</span>
    <span>&quot;ICC&quot;</span> <span>=</span> ICC<span>,</span>
    <span>&quot;WI_sd&quot;</span> <span>=</span> <span>sqrt</span><span>(</span>sigma2_wi<span>),</span>
    <span>&quot;BW_sd&quot;</span> <span>=</span> <span>sqrt</span><span>(</span>sigma2_bw<span>),</span>
    <span>&quot;mean&quot;</span> <span>=</span> mu
  <span>)</span>
<span>}</span>

<span># ---- icc-sim ----</span>
get_ICC_ln_sim <span>&lt;-</span> <span>function</span><span>(</span>B0<span>,</span> u0<span>,</span> sd_log<span>,</span> R <span>=</span> <span>1e4</span><span>)</span> <span>{</span>
  RE <span>&lt;-</span> B0 <span>+</span> rnorm<span>(</span>R<span>,</span> <span>0</span><span>,</span> u0<span>)</span>
  y <span>&lt;-</span> <span>exp</span><span>(</span>RE <span>+</span> sd_log<span>^</span><span>2</span><span>/</span><span>2</span><span>)</span>

  sigma2_bw <span>&lt;-</span> var<span>(</span>y<span>)</span> 
  sigma2_wi <span>&lt;-</span> <span>mean</span><span>(</span> <span>(</span><span>exp</span><span>(</span>sd_log<span>^</span><span>2</span><span>)</span> <span>-</span> <span>1</span><span>)</span> <span>*</span> <span>exp</span><span>(</span><span>2</span> <span>*</span> RE <span>+</span> sd_log<span>^</span><span>2</span><span>))</span>

  ICC <span>&lt;-</span> sigma2_bw <span>/</span> <span>(</span>sigma2_bw <span>+</span> sigma2_wi<span>)</span>

  <span>list</span><span>(</span><span>&quot;ICC&quot;</span> <span>=</span> ICC<span>,</span>
       <span>&quot;WI_sd&quot;</span> <span>=</span> <span>sqrt</span><span>(</span>sigma2_wi<span>),</span>
       <span>&quot;BW_sd&quot;</span> <span>=</span> <span>sqrt</span><span>(</span>sigma2_bw<span>),</span>
       <span>&quot;mean&quot;</span> <span>=</span> <span>mean</span><span>(</span>y<span>))</span>
<span>}</span>


<span># ---- sim-function ----</span>
summarise_posterior <span>&lt;-</span> <span>function</span><span>(</span>samples<span>,</span> name<span>)</span> <span>{</span>
  <span>data.frame</span><span>(</span>
    <span>&quot;var&quot;</span> <span>=</span> name<span>,</span>
    <span>&quot;est&quot;</span> <span>=</span> <span>mean</span><span>(</span>samples<span>),</span>
    <span>&quot;est_median&quot;</span> <span>=</span> median<span>(</span>samples<span>),</span>
    <span>&quot;est_se&quot;</span> <span>=</span> sd<span>(</span>samples<span>),</span>
    <span>&quot;lwr&quot;</span> <span>=</span> quantile<span>(</span>samples<span>,</span> <span>0.025</span><span>),</span>
    <span>&quot;upr&quot;</span> <span>=</span> quantile<span>(</span>samples<span>,</span> <span>0.975</span><span>),</span>
    row.names <span>=</span> <span>NULL</span><span>,</span>
    stringsAsFactors <span>=</span> <span>FALSE</span>
  <span>)</span>
<span>}</span>
add_thetas <span>&lt;-</span> <span>function</span><span>(</span>x<span>,</span> pars<span>)</span> <span>{</span>
  theta_C <span>&lt;-</span> <span>do.call</span><span>(</span>get_ICC_ln<span>,</span> pars<span>)</span>
  theta_TX <span>&lt;-</span> get_ICC_ln<span>(</span>B0 <span>=</span> pars<span>$</span>B0 <span>+</span> pars<span>$</span>B1<span>,</span>
                             u0 <span>=</span> pars<span>$</span>u0<span>,</span> 
                             sd_log <span>=</span> pars<span>$</span>sd_log<span>)</span>

    mutate<span>(</span>x<span>,</span>
      theta <span>=</span> case_when<span>(</span>
        var <span>==</span> <span>&quot;b0&quot;</span> <span>~</span> pars<span>$</span>B0<span>,</span>
        var <span>==</span> <span>&quot;b1&quot;</span> <span>~</span> pars<span>$</span>B1<span>,</span>
        var <span>==</span> <span>&quot;sd_u0_log&quot;</span> <span>~</span> pars<span>$</span>u0<span>,</span>
        var <span>==</span> <span>&quot;sd_BW_C&quot;</span> <span>~</span> theta_C<span>$</span>BW_sd<span>,</span>
        var <span>==</span> <span>&quot;sd_BW_TX&quot;</span> <span>~</span> theta_TX<span>$</span>BW_sd<span>,</span>
        var <span>==</span> <span>&quot;sd_WI_C&quot;</span> <span>~</span> theta_C<span>$</span>WI_sd<span>,</span>
        var <span>==</span> <span>&quot;sd_WI_TX&quot;</span> <span>~</span> theta_TX<span>$</span>WI_sd<span>,</span>
        var <span>==</span> <span>&quot;marg_mean_C&quot;</span> <span>~</span> theta_C<span>$</span><span>mean</span><span>,</span>
        var <span>==</span> <span>&quot;marg_mean_TX&quot;</span> <span>~</span> theta_TX<span>$</span><span>mean</span><span>,</span>
        var <span>==</span> <span>&quot;marg_TE&quot;</span> <span>~</span> theta_TX<span>$</span>mean <span>-</span> theta_C<span>$</span><span>mean</span><span>,</span>
        var <span>==</span> <span>&quot;multiplicative_TE_ratio_marg&quot;</span> <span>~</span> theta_TX<span>$</span><span>mean</span><span>/</span>theta_C<span>$</span><span>mean</span><span>,</span>
        var <span>==</span> <span>&quot;multiplicative_TE_ratio_medians&quot;</span> <span>~</span> <span>with</span><span>(</span>pars<span>,</span> <span>exp</span><span>(</span>B0 <span>+</span> B1<span>)</span><span>/</span><span>exp</span><span>(</span>B0<span>)),</span>
        var <span>==</span> <span>&quot;multiplicative_TE_exp_beta1&quot;</span> <span>~</span> <span>exp</span><span>(</span>pars<span>$</span>B1<span>),</span>
        var <span>==</span> <span>&quot;median_clust_C&quot;</span> <span>~</span> <span>with</span><span>(</span>pars<span>,</span> <span>exp</span><span>(</span>B0<span>)),</span>
        var <span>==</span> <span>&quot;median_clust_TX&quot;</span> <span>~</span> <span>with</span><span>(</span>pars<span>,</span> <span>exp</span><span>(</span>B0 <span>+</span> B1<span>)),</span>
        var <span>==</span> <span>&quot;median_clust_TE&quot;</span> <span>~</span> <span>with</span><span>(</span>pars<span>,</span> <span>exp</span><span>(</span>B0 <span>+</span> B1<span>)</span> <span>-</span> <span>exp</span><span>(</span>B0<span>)),</span>
        var <span>==</span> <span>&quot;sd_log&quot;</span> <span>~</span> pars<span>$</span>sd_log<span>,</span>
        var <span>==</span> <span>&quot;ICC&quot;</span> <span>~</span> theta_C<span>$</span>ICC<span>,</span>
        var <span>==</span> <span>&quot;ICC_log&quot;</span> <span>~</span> <span>with</span><span>(</span>pars<span>,</span> u0<span>^</span><span>2</span><span>/</span><span>(</span>u0<span>^</span><span>2</span> <span>+</span> sd_log<span>^</span><span>2</span><span>)),</span>
        <span>TRUE</span> <span>~</span> <span>NA_real_</span>
      <span>)</span>
    <span>)</span>
<span>}</span>


<span># extract results</span>

extract_results_LN <span>&lt;-</span> <span>function</span><span>(</span>fit_LN<span>)</span> <span>{</span>
  ss <span>&lt;-</span> posterior_samples<span>(</span>fit_LN<span>)</span>

  icc_C <span>&lt;-</span> pmap_df<span>(</span>
    <span>list</span><span>(</span>
      B0 <span>=</span> ss<span>$</span>b_Intercept<span>,</span>
      u0 <span>=</span> ss<span>$</span>sd_cluster__Intercept<span>,</span>
      sd_log <span>=</span> ss<span>$</span>sigma
    <span>),</span>
    <span>.</span>f <span>=</span> get_ICC_ln
  <span>)</span>

  icc_TX <span>&lt;-</span> pmap_df<span>(</span>
    <span>list</span><span>(</span>
      B0 <span>=</span> ss<span>$</span>b_Intercept <span>+</span> ss<span>$</span>b_TX<span>,</span>
      u0 <span>=</span> ss<span>$</span>sd_cluster__Intercept<span>,</span>
      sd_log <span>=</span> ss<span>$</span>sigma
    <span>),</span>
    <span>.</span>f <span>=</span> get_ICC_ln
  <span>)</span>

 <span># posterior quantities</span>
  median_clust_C <span>&lt;-</span> <span>exp</span><span>(</span>ss<span>$</span>b_Intercept<span>)</span>
  median_clust_TX <span>&lt;-</span> <span>exp</span><span>(</span>ss<span>$</span>b_Intercept <span>+</span> ss<span>$</span>b_TX<span>)</span>

  vars <span>&lt;-</span> <span>list</span><span>(</span>
    <span>&quot;b0&quot;</span> <span>=</span> ss<span>$</span>b_Intercept<span>,</span>
    <span>&quot;b1&quot;</span> <span>=</span> ss<span>$</span>b_TX<span>,</span>
    <span>&quot;sd_u0_log&quot;</span> <span>=</span> ss<span>$</span>sd_cluster__Intercept<span>,</span>
    <span>&quot;sd_BW_C&quot;</span> <span>=</span> icc_C<span>$</span>BW_sd<span>,</span>
    <span>&quot;sd_BW_TX&quot;</span> <span>=</span> icc_TX<span>$</span>BW_sd<span>,</span>
    <span>&quot;sd_WI_C&quot;</span> <span>=</span> icc_C<span>$</span>WI_sd<span>,</span>
    <span>&quot;sd_WI_TX&quot;</span> <span>=</span> icc_TX<span>$</span>WI_sd<span>,</span>
    <span>&quot;marg_mean_C&quot;</span> <span>=</span> icc_C<span>$</span><span>mean</span><span>,</span>
    <span>&quot;marg_mean_TX&quot;</span> <span>=</span> icc_TX<span>$</span><span>mean</span><span>,</span>
    <span>&quot;marg_TE&quot;</span> <span>=</span> icc_TX<span>$</span>mean <span>-</span> icc_C<span>$</span><span>mean</span><span>,</span>
    <span>&quot;median_clust_C&quot;</span> <span>=</span> median_clust_C<span>,</span>
    <span>&quot;median_clust_TX&quot;</span> <span>=</span> median_clust_TX<span>,</span>
    <span>&quot;median_clust_TE&quot;</span> <span>=</span> median_clust_TX <span>-</span> median_clust_C<span>,</span>
    <span>&quot;multiplicative_TE_ratio_marg&quot;</span> <span>=</span> icc_TX<span>$</span><span>mean</span><span>/</span>icc_C<span>$</span><span>mean</span><span>,</span>
    <span>&quot;multiplicative_TE_ratio_medians&quot;</span> <span>=</span> median_clust_TX<span>/</span>median_clust_C<span>,</span>
    <span>&quot;multiplicative_TE_exp_beta1&quot;</span> <span>=</span> <span>exp</span><span>(</span>ss<span>$</span>b_TX<span>),</span>
    <span>&quot;sd_log&quot;</span> <span>=</span> ss<span>$</span>sigma<span>,</span>
    <span>&quot;ICC&quot;</span> <span>=</span> icc_C<span>$</span>ICC<span>,</span>
    <span>&quot;ICC_log&quot;</span> <span>=</span> ss<span>$</span>sd_cluster__Intercept<span>^</span><span>2</span><span>/</span><span>(</span>ss<span>$</span>sd_cluster__Intercept<span>^</span><span>2</span> <span>+</span> ss<span>$</span>sigma<span>^</span><span>2</span><span>)</span>
  <span>)</span>
  res_LN <span>&lt;-</span> map2_df<span>(</span>vars<span>,</span> 
                    <span>names</span><span>(</span>vars<span>),</span> 
                    summarise_posterior<span>)</span>
  res_LN<span>$</span>model <span>&lt;-</span> <span>&quot;LN&quot;</span>

  res_LN
<span>}</span>

<span># gaussian</span>
extract_results_gauss <span>&lt;-</span> <span>function</span><span>(</span>fit_gauss<span>)</span> <span>{</span>
  ss <span>&lt;-</span> posterior_samples<span>(</span>fit_gauss<span>)</span>

  <span># posterior quantities</span>
  vars <span>&lt;-</span> <span>list</span><span>(</span>
    <span>&quot;sd_BW_C&quot;</span> <span>=</span> ss<span>$</span>sd_cluster__Intercept<span>,</span>
    <span>&quot;sd_BW_TX&quot;</span> <span>=</span> ss<span>$</span>sd_cluster__Intercept<span>,</span>
    <span>&quot;sd_WI_C&quot;</span> <span>=</span> ss<span>$</span>sigma<span>,</span>
    <span>&quot;sd_WI_TX&quot;</span> <span>=</span> ss<span>$</span>sigma<span>,</span>
    <span>&quot;marg_mean_C&quot;</span> <span>=</span> ss<span>$</span>b_Intercept<span>,</span>
    <span>&quot;marg_mean_TX&quot;</span> <span>=</span> ss<span>$</span>b_Intercept <span>+</span> ss<span>$</span>b_TX<span>,</span>
    <span>&quot;marg_TE&quot;</span> <span>=</span> ss<span>$</span>b_TX<span>,</span>
    <span>&quot;multiplicative_TE_ratio_marg&quot;</span> <span>=</span> <span>(</span>ss<span>$</span>b_Intercept <span>+</span> ss<span>$</span>b_TX<span>)</span><span>/</span>ss<span>$</span>b_Intercept<span>,</span>
    <span>&quot;ICC&quot;</span> <span>=</span> ss<span>$</span>sd_cluster__Intercept<span>^</span><span>2</span><span>/</span><span>(</span>ss<span>$</span>sd_cluster__Intercept<span>^</span><span>2</span> <span>+</span> ss<span>$</span>sigma<span>^</span><span>2</span><span>)</span>
  <span>)</span>
  res_gauss <span>&lt;-</span> map2_df<span>(</span>vars<span>,</span> 
                       <span>names</span><span>(</span>vars<span>),</span> 
                       summarise_posterior<span>)</span>
  res_gauss<span>$</span>model <span>&lt;-</span> <span>&quot;gaussian&quot;</span>

  res_gauss
<span>}</span>

<span># performs 1 sim replications</span>
refit_mod <span>&lt;-</span> <span>function</span><span>(</span>i<span>,</span> fit_LN<span>,</span> fit_gauss<span>,</span> pars<span>)</span> <span>{</span>
  d <span>&lt;-</span> <span>do.call</span><span>(</span>gen_data<span>,</span> pars<span>)</span>

  <span># fit models</span>
  fit_LN <span>&lt;-</span> update<span>(</span>
    fit_LN<span>,</span>
    newdata <span>=</span> d<span>,</span>
    recompile <span>=</span> <span>FALSE</span><span>,</span>
    chains <span>=</span> <span>1</span><span>,</span>
    iter <span>=</span> <span>2000</span>
  <span>)</span>
  fit_gauss <span>&lt;-</span> update<span>(</span>
    fit_gauss<span>,</span>
    newdata <span>=</span> d<span>,</span>
    recompile <span>=</span> <span>FALSE</span><span>,</span>
    chains <span>=</span> <span>1</span><span>,</span>
    iter <span>=</span> <span>2000</span>
  <span>)</span>

  <span># save results</span>
  res_LN <span>&lt;-</span> extract_results_LN<span>(</span>fit_LN<span>)</span>
  res_gauss <span>&lt;-</span> extract_results_gauss<span>(</span>fit_gauss<span>)</span>

  <span># Combine</span>
  res <span>&lt;-</span> <span>rbind</span><span>(</span>res_LN<span>,</span> 
               res_gauss<span>)</span>
  res<span>$</span>sim <span>&lt;-</span> i

  res
<span>}</span>



<span># Run simulation ----------------------------------------------------------</span>
run_sim <span>&lt;-</span> <span>function</span><span>(</span>nsim<span>,</span> pars<span>,</span> fit_LN<span>,</span> fit_gauss<span>)</span> <span>{</span>
  res <span>&lt;-</span> future.apply<span>::</span>future_lapply<span>(</span><span>1</span><span>:</span>nsim<span>,</span> 
                       refit_mod<span>,</span>
                       fit <span>=</span> fit_LN<span>,</span> 
                       fit_gauss <span>=</span> fit_gauss<span>,</span> 
                       pars <span>=</span> pars<span>)</span>

  res <span>&lt;-</span> <span>do.call</span><span>(</span><span>rbind</span><span>,</span> res<span>)</span>
  <span>list</span><span>(</span>res <span>=</span> res<span>,</span> 
       pars <span>=</span> pars<span>)</span>
<span>}</span> 

<span># Summarise simulation ----------------------------------------------------</span>
summarise_sim <span>&lt;-</span> <span>function</span><span>(</span>x<span>)</span> <span>{</span>
  x<span>$</span>res <span>%>%</span>
    add_thetas<span>(</span>x<span>$</span>pars<span>)</span> <span>%>%</span>
    group_by<span>(</span>var<span>,</span> model<span>)</span> <span>%>%</span>
    summarise<span>(</span>
      est_mean <span>=</span> <span>mean</span><span>(</span>est<span>),</span>
      est_median <span>=</span> <span>mean</span><span>(</span>est_median<span>),</span>
      theta <span>=</span> <span>unique</span><span>(</span>theta<span>),</span>
      est_rel_bias <span>=</span> <span>(</span>est_mean <span>-</span> theta<span>)</span> <span>/</span> theta<span>,</span>
      est_median_rel_bias <span>=</span> <span>(</span>est_median <span>-</span> theta<span>)</span> <span>/</span> theta<span>,</span>
      se <span>=</span> sd<span>(</span>est<span>),</span>
      coverage <span>=</span> <span>mean</span><span>(</span>lwr <span></span> theta <span>&amp;</span> theta <span></span> upr<span>),</span>
      power <span>=</span> <span>mean</span><span>(</span><span>sign</span><span>(</span>lwr<span>)</span> <span>==</span> <span>sign</span><span>(</span>upr<span>))</span>
    <span>)</span>
<span>}</span>
</pre>
</div>
</div>
<p></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="http://rpsychologist.com/GLMM-part1-lognormal"> R Psychologist - R</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/PqZxUjNQ1Pg" height="1"></a></div>