# Brms meta-analysis - how to include standard errors of zero?

**URL:** <https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477>\
**Category:** Modeling\
**Tags:** meta-analysis, brms\
**Created:** [September 17, 2021, 7:34pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477 "2021-09-17T19:34:44Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![xiazhu111](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/xiazhu111/32/12772_2.png) [@xiazhu111](https://discourse.mc-stan.org/u/xiazhu111)\
**Post date:** [September 17, 2021, 7:34pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/1 "2021-09-17T19:34:44Z")

</div>

Hi everyone,

I am building a model an additive model that incorporates standard errors (SE) and regression weight, which are the sample sizes. The data I am using come from the literature, so the approach is meta-analytic. Some authors published their raw data, in order words, individual observations with sample size = 1 and standard error = 0. Some authors reported average +/- SE with sample size \> 1.

The problem is that the family of distributions I am using is gaussian with a log link, and an error arises when I input SEs of 0. Here is the structure of my model and the output when I try to run it:

`fit1 <- brm(response|weights(Weight) + se(SE) ~ X1 + X2 + X3..., family = gaussian(link="log"), data = Data)`

Output:

```
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: normal_lpdf: Scale parameter is 0, but must be > 0! (in 'model125839f2ad0_c2a2bf402f16cab2ec43c0c861aee5d3' at line 40)

```

I am wondering if there is some way to combine averages and individual samples in brms, or if anyone has suggestions on how to move forward?

Thank you for your time!

---

<div class="post-metadata">

**Author:** ![xiazhu111](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/xiazhu111/32/12772_2.png) [@xiazhu111](https://discourse.mc-stan.org/u/xiazhu111)\
**Post date:** [September 17, 2021, 7:42pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/2 "2021-09-17T19:42:17Z")

</div>

The data looks something like this, hypothetically:

```
Response SE Weight
10 0 1
6 0 1
15 7 5

```

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [September 24, 2021, 4:55pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/3 "2021-09-24T16:55:27Z")

</div>

Unfortunately, I don’t think combining individual and study-level data can be achieved with `brms` right now - at least not easily. I think you might be able to use the `subset` addition term to provide a different formula (without `se`) for the individual observations. However that would IMHO make the model coefficients differ between the two formulas. You might be able to go around this by rewriting your model via the non-linear syntax ([Estimating Non-Linear Models with brms • brms](https://paul-buerkner.github.io/brms/articles/brms_nonlinear.html)), which let’s you have the same parameter in multiple formulas, but makes the whole thing much more complex. At this point, you are somewhat fighting against `brms` and I am not sure if implementing the model directly in Stan wouldn’t be easier.

I’ll tag @paul.buerkner for second opinion (as he’s also AFAIK working on metanalyses in `brms`).

Best of luck with your project!

---

<div class="post-metadata">

**Author:** ![xiazhu111](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/xiazhu111/32/12772_2.png) [@xiazhu111](https://discourse.mc-stan.org/u/xiazhu111)\
**Post date:** [September 24, 2021, 9:28pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/4 "2021-09-24T21:28:21Z")

</div>

Thank you for your input, Mr. Modrak! Unless Mr. Burkner has other suggestions, we might end up doing it in a roundabout way I believe similar to what you proposed, or averaging the individual observations while making sure sample sizes are as small as can be.

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [September 25, 2021, 9:09am UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/5 "2021-09-25T09:09:18Z")

</div>

```
fit1 <- brm(response|weights(Weight) + se(SE, sigma = TRUE) ~ X1 + X2 + X3..., family = gaussian(link="log"), data = Data)

```

could work. that way the scale will not be zero as you still have a residual standard deviation parameter (which would probably want anyway). Note that, In classical meta analysis with one data row (effect size) per study, `sigma = TRUE` would replace the observation-level random effect as both measure the residual SD and would not be jointly identified.

---

<div class="post-metadata">

**Author:** ![xiazhu111](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/xiazhu111/32/12772_2.png) [@xiazhu111](https://discourse.mc-stan.org/u/xiazhu111)\
**Post date:** [September 25, 2021, 7:34pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/6 "2021-09-25T19:34:46Z")

</div>

Thank you for your response, Dr. Buerkner. I tried the code above and obtained the following output:  
“The model does not contain posterior samples.”  
I am wondering how to resolve this issue?

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [September 26, 2021, 9:18pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/7 "2021-09-26T21:18:00Z")

</div>

This does not really tell much. Can you try with one chain and show all the messages that (r)stan prints out?

---

<div class="post-metadata">

**Author:** ![xiazhu111](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/xiazhu111/32/12772_2.png) [@xiazhu111](https://discourse.mc-stan.org/u/xiazhu111)\
**Post date:** [September 27, 2021, 5:03pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/8 "2021-09-27T17:03:37Z")

</div>

Hi Dr. Burkner, here is the format of the code I ran:

```
fit1 <- brm(response|weights(Weight) + se(SE, sigma=TRUE) ~ 
X1 + X2 + X3 + ..., family = gaussian(link="log"),chains=1,iter=4000,warmup=500,data = Data, control = list(adapt_delta = 0.99, max_treedepth = 17))

```

And the output I got was:

```
Warning messages:
1: There were 500 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 17. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded 
2: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low 
3: Examine the pairs() plot to diagnose sampling problems
4: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess 
5: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess

```

Seems like it should work if I just increase max\_treedepth? This output is different from last time - previously, when I ran with `chains = 4`, I didn’t get any warning messages at all but also no posterior samples were drawn. I think the model might be okay now!

---

<div class="post-metadata">

**Author:** ![xiazhu111](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/xiazhu111/32/12772_2.png) [@xiazhu111](https://discourse.mc-stan.org/u/xiazhu111)\
**Post date:** [September 27, 2021, 5:18pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/9 "2021-09-27T17:18:48Z")

</div>

Good news, Dr. Burkner: I re-ran with `max_treedepth = 18` and obtained no warning messages. Also the model sampled successfully and all Rhats were 1.00. Then my question is why does the model not work for `chains = 4`?

I tried running using chains = 3 and got this warning message:

```
Warning messages:
1: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low 
2: Examine the pairs() plot to diagnose sampling problems

```

But besides this, it seems the model worked. The next step is to use this code on my real data, and I hope it goes well. For now, I think this issue is resolved! Thank you again for your help.

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [September 28, 2021, 6:37am UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/10 "2021-09-28T06:37:26Z")

</div>

You should still look at Rhat and ESS etc. and perhaps trace plots (via plot()) or similar but it seems the model at least runs now :-)

---

<div class="post-metadata">

**Author:** ![xiazhu111](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/xiazhu111/32/12772_2.png) [@xiazhu111](https://discourse.mc-stan.org/u/xiazhu111)\
**Post date:** [September 28, 2021, 2:57pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/11 "2021-09-28T14:57:13Z")

</div>

Hi Dr. Burkner, I just want to confirm one more thing: setting `sigma = TRUE` does not replace the SEs we give the model? In other words, it is there to help the model run, but it does not replace the information we give the model? @paul.buerkner

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [September 28, 2021, 3:34pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/12 "2021-09-28T15:34:07Z")

</div>

it just estimates a residual standard deviation parameter in addition to the known SEs. something you would want in your model anyway.

---

<div class="post-metadata">

**Author:** ![xiazhu111](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/xiazhu111/32/12772_2.png) [@xiazhu111](https://discourse.mc-stan.org/u/xiazhu111)\
**Post date:** [September 28, 2021, 3:48pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/13 "2021-09-28T15:48:53Z")

</div>

Great, thank you!

---

<div class="post-metadata">

**Author:** ![drjhw](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/drjhw/32/9524_2.png) [@drjhw](https://discourse.mc-stan.org/u/drjhw)\
**Post date:** [May 23, 2023, 2:43pm UTC](https://discourse.mc-stan.org/t/brms-meta-analysis-how-to-include-standard-errors-of-zero/24477/14 "2023-05-23T14:43:25Z")

</div>

I have what I think is a parallel problem, that I think is very prevalent in biomedical research.

I want to analyse a multi-level data-set that includes both nested and crossed random effects. The data are from an experiment to measure the intensity of a marker in individual living cells. We observe cells occurs in “subplots” (microscope fields, where each field may show several cells), which are samples of plots (culture dishes where the cells grow, that usually contain hundreds of cells). We measure the intensity of the marker (in individual cells, nested in subplots, within plots) at different times. The cells are from different people (blocks) and we sometimes measure people’s cells at different times - but there is partial overlap between people and times. The cells that we measure for each dish and person are not the same from time to time.

It’s a complicated and unbalanced design - but it reflects the realities of when cells become available and how many of them grow at different times in different culture dishes and different microscope fields (which we do not control). As I say, complex data of this kind are very common in biomedical research. We ignore differences between microscope fields (subplots) for now.

There are very many cells in the total study (~200k) - and some individual microscope fields have no cells, so that estimation of the full model for all the data, using the zero-inflated negative binomial family is practically impossible (it takes several hours to become reasonably stable using variational Bayes, but even if I use the estimates from the vB output as starting values for estimation using NUTS, the warm-up phase ran for 72 hours without reporting any progress at all).

Since I cannot analyse the raw data, I wondered about using a meta-analysis approach, in line with the method suggested previously in this post. Specifically, I obtained the mean and standard deviation of the intensity in each culture dish at each time, along with the number of cells where we measured the intensity. So, after taking the mean intensity per cell and its standard deviation for each dish, the data look something like this:-

Person culture\_dish time Total\_N\_of\_cells mean\_intensity\_per\_cell sd\_of\_intensity\_per\_cell  
1 1 1 123 3.5 1.1  
1 2 2 322 4.2 1.5  
2 1 2 126 2.1 0.6  
2 2 1 139 1.4 0.4  
etc. etc. (there are roughly 200 cells per person - so the total data-set has about 1000 observations).

I want to estimate the mean\_intensity\_per\_cell for each person, taking account - _as far as possible_ - the variation between culture dishes and times. I take it that the problem shows formal similarity to Xiazhu111’s question. That is, it seems reasonable (to me) to conceptualise the design as a meta-analysis, where each person constitutes an ‘experiment’. The main difference from XiaZhu111’s original post is that there are no culture dishes with no cells, so every culture dish has a real/positive standard deviation for the mean\_intensity\_per\_cell. Given that there are no standard deviations of zero, following Paul Buerker’s post (above), I wrote the following brmsformula to separate the effects of individual people from the other design factors:-

mean\_intensity\_per\_cell|se(sd\_of\_intensity\_per\_cell/sqrt(Total\_N\_of\_cells+1))+weights(Total\_N\_of\_cells)~1+(1|Person)+(1|culture\_dish)+(1|time), family=gaussian(link=‘log’)

(I use family=gaussian(link=‘log’) for 3 reasons: (1) se() is not available with family=lognormal; (2) the distribution of mean\_intensity\_per\_cell is actually right-skewed; (3) Paul Buerkner used it in his previous post).

So, now I would like to ask:-

1. Does this formula seem appropriate to take into account both the variation in the mean\_intensity\_per\_cell and the Total\_N\_of\_cells between culture dishes?
2. is it reasonable to use family=gaussian(link=‘log’) to account for the skewness of the data (mean\_intensity\_per\_cell), or would it be better to take log(mean\_intensity\_per\_cell) and use family=gaussian, or family=student?
3. If I use log(mean\_intensity\_per\_cell), then should I also log-transform its standard error and the Total\_N\_of\_cells?

I think it likely that there are many biomedical researchers whose data are similar to those that I describe here. Therefore, figuring out the best way to analyse these complex, unbalanced and numerous data (initial N~200k, or more) could be very helpful to many other researchers.

With many thanks, in anticipation of your help