Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Bayesian normal meta-analysis with descriptives
---
title: "Normal meta-analysis"
author: "Richard D. Morey"
date: "04/07/2017"
output: html_document
---
```{r}
# data here
mu.true = 100
tau.true = 1/15^2
N = rpois(30, 20) + 100
M = length(N)
y = rnorm(M, mu.true, 1 / sqrt(tau.true * N))
s2 = rgamma(M, (N - 1)/2, (N - 1)*tau.true/2)
```
## Fixed effects
Let $y_i$, $s^2$, and $N_i$ be the observed mean, variance, and sample size for the $i$th study ($i = 1,\ldots,M$). Assume that the population is normal with mean $\mu$ and precision $\tau$, and hence,
\[
y_i \stackrel{iid}{\sim} \mbox{Normal}(\mu, N_i \tau)
\]
and
\[
s^2_i \stackrel{iid}{\sim} \Gamma( (N_i - 1)/2, (N_i - 1)\tau/2 )
\]
(note the parametrization of the normal distribution in terms of precision, and the parametrization of the gamma distribution in terms of rate). $y_i$ and $s^2_i$ are independent thanks to the assumption of normality.
For convenience let
\[
\begin{eqnarray*}
N_t &=& \sum N_i\\
\bar{y} &=& \frac{\sum_{i=1}^M N_iy_i}{N_t}\\
S_y &=& \sum_{i=1}^MN_i(y_i - \bar{y})^2\\
S_w &=& \sum_{i=1}^M (N_i-1)s^2_i
\end{eqnarray*}
\]
The joint likelihood is thus proportional to
\[
p(\boldsymbol y,\boldsymbol {s^2}\mid \mu, \tau) \propto\tau^{M/2}\exp\left\{-\frac{\tau}{2}\sum_{i=1}^MN_i\left(y_i - \mu\right)^2\right\}\exp\left\{-\frac{\tau}{2}\sum_{i=1}^M(N_i - 1)q_i\right\}\prod_{i=1}^M\tau^{(N_i - 1)/2}
\]
which can be simplified to
\[
p(\boldsymbol y,\boldsymbol {s^2}\mid \mu, \tau) \propto\tau^{N_t/2}\exp\left\{-\frac{\tau}{2}\left(S_y + S_w + N_t(\mu - \bar{y})^2\right)\right\}
\]
### Prior 1: "Non-informative"
We might choose the prior
\[
p(\mu, \tau) \propto \tau^{-1}.
\]
In this case, the joint posterior is proportional to
\[
p(\mu, \tau\mid \boldsymbol y,\boldsymbol {s^2}) \propto\tau^{N_t/2 - 1}\exp\left\{-\frac{\tau}{2}\left(S_y + S_w + N_t(\mu - \bar{y})^2\right)\right\}
\]
It is easy to show that the marginal posterior for $\tau$ is
\[
\tau\mid\boldsymbol y,\boldsymbol {s^2}\sim \Gamma(a, b)
\]
where
\[
\begin{eqnarray}
a &=& (N_t - 1)/2\\
b &=& \frac{S_y + S_w}{2}
\end{eqnarray}
\]
and the marginal posterior for $\mu$ is
\[
\mu\mid\boldsymbol y,\boldsymbol {s^2}\sim \mbox{Shifted, scaled }t_{N_t-1}(\bar{y},\sigma_\mu)
\]
where
\[
\begin{eqnarray*}
\sigma_\mu&=& \sqrt{\frac{S_y + S_w}{N_t(N_t - 1)}}
\end{eqnarray*}
\]
The quantity $\sigma_\mu$ can be interpreted similarly to the way people often interpret standard errors.
```{r}
Nt = sum(N)
ybar = sum(y * N)/Nt
Sy = sum( N * (y - ybar)^2 )
Sw = sum((N-1)*s2)
nu_mu = Nt - 1
scale_mu = sqrt( (Sy + Sw) / (Nt * nu_mu) )
location_mu = ybar
shape_tau = (Nt - 1)/2
rate_tau = (Sy + Sw)/2
noninformed = c(
scale_mu = scale_mu,
location_mu = location_mu,
shape_tau = shape_tau,
rate_tau = rate_tau,
mean_tau = shape_tau / rate_tau,
sd_tau = sqrt(shape_tau) / rate_tau
)
```
### Prior 2: Conjugate
\[
\begin{eqnarray*}
\mu \mid \tau \sim \mbox{Normal}(\mu_0, \tau N_0)\\
\tau \sim \Gamma\left(\nu_0/2, \nu_0/(2\tau_0)\right)
\end{eqnarray*}
\]
The prior parameter $N_0$ can be interpreted on the sample size scale; it is the weight given to the prior mean $\mu_0$, in equivalent sample size units.
The prior parameter $\tau_0$ is the prior mean of the parameter $\tau$; the prior parameter $\nu_0$ can be interpreted as the weight given to the prior mean $\tau_0$, in equivalent sample size units.
The joint prior is
\[
p(\mu, \tau) \propto \tau^{1/2}\exp\left\{-\frac{\tau N_0}{2}(\mu - \mu_0)^2\right\} \tau^{\nu_0/2 - 1} \exp\left\{-\frac{\tau\nu_0}{2\tau_0}\right\}
\]
The joint posterior is proportional to
\[
p(\mu, \tau\mid \boldsymbol y,\boldsymbol {s^2}) \propto\tau^{\frac{N_t + \nu_0 - 1}{2} - 1} \exp\left\{-\frac{\tau}{2}\left(\nu_0/\tau_0 + S_y + S_w + N_t(\mu - \bar{y})^2 + N_0(\mu - \mu_0)^2\right)\right\}
\]
which simplifies to
\[
p(\mu, \tau\mid \boldsymbol y,\boldsymbol {s^2}) \propto\tau^{\frac{N_t + \nu_0 - 1}{2} - 1} \exp\left\{-\frac{\tau}{2}\left(\nu_0/\tau_0 + S_0 + S_y + S_w + N_1(\mu - \bar{y}_1)^2\right)\right\}
\]
where
\[
\begin{eqnarray*}
N_1 &=& N_t + N_0\\
\bar{y}_1 &=&\frac{N_t\bar{y} + N_0\mu_0}{N_1}\\
S_0 &=& N_t\bar{y}^2 + N_0\mu_0^2 - N_1\bar{y}_1^2
\end{eqnarray*}
\]
This is a similar form as before. The marginal posterior for $\tau$ is
\[
\tau\mid\boldsymbol y,\boldsymbol {s^2}\sim \Gamma(a, b)
\]
where
\[
\begin{eqnarray}
a &=& (N_t + \nu_0 - 1)/2\\
b &=& \frac{\nu_0/\tau_0 + S_0 + S_y + S_w}{2}
\end{eqnarray}
\]
and the marginal posterior for $\mu$ is
\[
\mu\mid\boldsymbol y,\boldsymbol {s^2}\sim \mbox{Shifted, scaled }t_{\nu_1}(\bar{y}_1,\sigma^{'}_\mu)
\]
where
\[
\begin{eqnarray*}
\nu_1 &=& N_t + \nu_0 - 2\\
\sigma^{'}_\mu&=& \sqrt{\frac{\nu_0/\tau_0 + S_0 + S_y + S_w}{N_1\nu_1}}
\end{eqnarray*}
\]
```{r}
# Prior parameters
N0 = 1
mu0 = 90
nu0 = 1
tau0 = 10
### end prior parameters
N1 = Nt + N0
ybar1 = (Nt * ybar + N0 * mu0) / N1
S0 = Nt*ybar^2 + N0*mu0^2 - N1*ybar1^2
nu1 = Nt + nu0 - 2
scale_mu = sqrt( (Sy + Sw + S0 + nu0/tau0) / (N1 * nu1) )
location_mu = ybar1
shape_tau = (nu1 + 1)/2
rate_tau = (Sy + Sw + S0 + nu0/tau0)/2
conjugate = c(
scale_mu = scale_mu,
location_mu = location_mu,
shape_tau = shape_tau,
rate_tau = rate_tau,
mean_tau = shape_tau / rate_tau,
sd_tau = sqrt(shape_tau) / rate_tau
)
```
```{r}
knitr::kable(cbind(noninformed = noninformed, conjugate = conjugate))
```
### Check with JAGS
```{r}
mod = "
model{
for(i in 1:length(y)){
y[i] ~ dnorm(mu, tau * N[i])
s2[i] ~ dgamma( (N[i] - 1)/2, tau * (N[i] - 1) / 2 )
}
mu ~ dnorm(mu0, N0 * tau)
tau ~ dgamma( nu0 / 2, nu0 / (tau0 * 2) )
}
"
forJAGS = list(y = y,
s2 = s2,
N = N,
N0 = N0,
mu0 = mu0,
nu0 = nu0,
tau0 = tau0)
library(rjags)
comp.mod = jags.model(file = textConnection(mod),
data = forJAGS)
samps = coda.samples(model = comp.mod, variable.names = c("mu", "tau"), n.iter = 100000)[[1]]
plot(density(samps[,"mu"]), main = "mu (mean)")
curve( dt( (x - location_mu)/scale_mu, nu1 )/scale_mu,
from = par()$usr[1],
to = par()$usr[2],
col = "red", add = TRUE)
plot(density(samps[,"tau"]), main = "tau (precision)")
curve( dgamma( x, shape_tau, rate = rate_tau),
from = par()$usr[1],
to = par()$usr[2],
col = "red", add = TRUE)
plot(density(1/sqrt(samps[,"tau"])), main = "sigma (standard deviation")
curve( 2*dgamma( 1/x^2, shape_tau, rate = rate_tau)/x^3,
from = par()$usr[1],
to = par()$usr[2],
col = "red", add = TRUE)
```
## Random effects
Random effects are easier to do with JAGS.
```{r}
mod = "
model{
for(i in 1:length(y)){
y[i] ~ dnorm(mu[i], N[i] / sig2[i])
s2[i] ~ dgamma( (N[i] - 1)/2, (N[i] - 1) / (2 * sig2[i]) )
mu[i] ~ dnorm(mu_mu, prec_mu)
# model on log(sigma)
sig2[i] <- exp(2*logsig[i])
logsig[i] ~ dnorm(mu_logsig, prec_logsig)
}
# These priors are only meant for demonstration
mu_mu ~ dnorm(0, 1/10000)
mu_logsig ~ dnorm(0, 1/10000)
prec_mu <- pow(sigma_mu, -2)
prec_logsig <- pow(sigma_logsig, -2)
# Half Cauchy (to allow heavy pooling)
sigma_mu ~ dt(0, 1, 1) I(0,)
sigma_logsig ~ dt(0, 1, 1) I(0,)
}
"
forJAGS = list(y = y,
s2 = s2,
N = N)
library(rjags)
comp.mod = jags.model(file = textConnection(mod),
data = forJAGS)
samps = coda.samples(model = comp.mod, variable.names = c("mu", "sig2", "mu_mu", "mu_logsig", "sigma_mu", "sigma_logsig"), n.iter = 100000)[[1]]
mu.cols = grepl(pattern = "mu[", colnames(samps), fixed = TRUE)
mu.est = colMeans(samps[,mu.cols])
# Show pooling
plot(y, mu.est, asp = TRUE, pch = 19)
abline(0,1, lty=2)
sig2.cols = grepl(pattern = "sig2[", colnames(samps), fixed = TRUE)
sig2.est = colMeans(samps[,sig2.cols])
# Show pooling
plot(s2, sig2.est, asp = TRUE, pch = 19)
abline(0,1, lty=2)
```
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.