medical-model-comparison.ipynb

Kernel: R (R-Project)

We consider the eczema medical trial data set again. This time we will compare which of 2 models explain the observed data best.

• Model 1: All studies have the same probability of success.

• Model 2: A hierarchical model where the probability of success in each study is drawn from a beta prior distribution with unknown $\alpha$ and $\beta$ parameters.

Study	Treatment group	Control group
Di Rienzo 2014	20 / 23	9 / 15
Galli 1994	10 / 16	11 / 18
Kaufman 1974	13 / 16	4 / 10
Qin 2014	35 / 45	21 / 39
Sanchez 2012	22 / 31	12 / 29
Silny 2006	7 / 10	0 / 10
Totals	107 / 141	57 / 121

Model 1:

• For each group (treatment and control), all 6 studies have the same fixed, but unknown, probability of success, $\theta_t,\theta_c\in[0,1]$.

• The data follow a binomial distribution in each study, conditioned on the probability of success — $\theta_t$ for treatment or $\theta_c$ for control.

• The priors over $\theta_t$ and $\theta_c$ are uniform.

These assumptions lead to the following model.

• Likelihood: $\prod_{i=1}^6 \text{Binomial}(s_i\,|\,\theta, n_i)$, where $s_i$ is the number of successful recoveries, $f_i$ is the number of failures (did not recover), and $n_i=s_i+f_i$ the number of patients.

• Prior: $\text{Beta}(\theta\,|\,1,1)$ for both $\theta_t$ and $\theta_c$.

• Posterior for treatment group: $\text{Beta}(\theta_t\,|\,108, 35)$.

• Posterior for control group: $\text{Beta}(\theta_c\,|\,58, 65)$.

Since we have closed-form solutions for the posteriors, we can calculate the marginal likelihood by rearranging Bayes' equation: (marginal likelihood) = (likelihood) x (prior) / (posterior).

$$P(\text{data}) = \left[\prod_{i=1}^6 \text{Binomial}(s_i\,|\,\theta, n_i) \right] \text{Beta}(\theta\,|\,\alpha_0,\beta_0) \,/\, \text{Beta}(\theta\,|\,\alpha_1,\beta_1)$$

where $\alpha_0=1$ and $\beta_0=1$ are the parameters of the prior, and $\alpha_1$ and $\beta_1$ are the parameters of the posterior beta distribution.

Since all factors involving $\theta$ cancel out, we are just left with the normalization constants of the likelihood, the prior and the posterior:

$$\begin{align} P(\text{data}) &= \left[ \prod_{i=1}^6 \left(\begin{array}{c}s_i+f_i \\ s_i\end{array}\right) \right] \frac{\text{B}(\alpha_1,\beta_1)}{\text{B}(\alpha_0,\beta_0)} \\ &= \left[\prod_{i=1}^6 \frac{1}{(s_i+f_i+1)\text{B}(s_i+1,f_i+1)}\right] \frac{\text{B}(\alpha_1,\beta_1)}{\text{B}(\alpha_0,\beta_0)} \end{align}$$

We usually compute the log of the marginal likelihood since the results can vary over many orders of magnitude.

Task 1:

• Take the log of the marginal likelihood above.

• Complete the R function below to calculate the log marginal likelihood.

• You can use the built-in function lbeta(a,b) to compute $\log \text{B}(a,b)$.

log_beta_binomial_marginal_likelihood <- function(alpha0, beta0, data) {
    # Compute the log marginal likelihood of a beta-binomial model.
    # alpha0, beta0: prior beta distribution parameters.
    # data: 2 x n matrix, where n is the number of samples from the
    #       binomial distribution. This means the first row contains
    #       the success counts and the second row the failure counts
    #       of the samples.

    alpha1 = alpha0 + sum(data[1,])  # Compute posterior parameters
    beta1 = beta0 + sum(data[2,])
    log_evidence = (
        lbeta(alpha1, beta1) -
        lbeta(alpha0, beta0) -
        sum(lbeta(data[1,] + 1, data[2,] + 1)) -
        sum(log(data[1,] + data[2,] + 1)))
    return(log_evidence)
}

#Vy's
vys_log_beta_binomial_marginal_likelihood <- function(alpha0, beta0, data) {
    alpha1<-0
    beta1<-0
    for (i in 1:length(data)){
        si<-data[i][1]
        fi<-data[i][2]
        sum<-sum+log(1/(si+fi+1))+log(1)-lbeta(si+1,fi+1)
        alpha1<-alpha1+si
    }
    sum<-sum+lbeta(alpha1,beta1)-lbeta(alpha0,beta0) 
    return(sum)
}

log_beta_binomial_marginal_likelihood <- function(alpha0, beta0, data) {
    # Compute the log marginal likelihood of the beta-binomial model for the eczema data set.
    si_treat = c(20,10,13,35,22,7)
    fi_treat = c(23,16,16,45,31,10)
    si_control = c(9,11,4,21,12,0)
    fi_control = c(15,18,10,39,29,10)
    
    pis_treat <- c()
    pis_control <- c()
    
    for s,f in si_treat,fi_treat {
        pi_i = - log(s+f+1) - (lbeta(s+1,f+1)))
        append(pis_treat,pi_i)
    }
    for s,f in si_control,fi_control {
        pi_i = - log(s+f+1) - (lbeta(s+1,f+1)))
        append(pis_control,pi_i)
    }
    
    p_data_treat = sum(pis_treat) + lbeta(108,35) - lbeta(1,1)
    p_data_control = sum(pis_control) + lbeta(58,65) - lbeta(1,1)
    
    # alpha0, beta0: prior beta distribution parameters.
    alpha0 = 1
    beta0 = 1
    # data: 2 x n matrix with number of successes in the first row and number of failures in the second row.
    
}

Error in parse(text = x, srcfile = src): <text>:11:9: unexpected symbol
10:     
11:     for s
            ^
Traceback:

Model 2:

• For each group (intervention and control), each of the 6 studies has a different probability of success.

• Each probability of success is drawn from a beta prior with unknown parameters $\alpha$ and $\beta$.

• Since $\alpha$ and $\beta$ are unknown, we put a broad hyperprior on them — we choose the Gamma(2, 0.5) distribution, which is shown below.

curve(dgamma(x, 2, 0.5), from=0, to=15)

These assumptions lead to the following model:

• Likelihood: $\prod_{i=1}^6 \text{Binomial}(s_i\,|\,\theta_i, n_i)$, where $s_i$ is the number of successful recoveries, $f_i$ is the number of failures (did not recover), and $n_i=s_i+f_i$ the number of patients. Note that each study has its own $\theta_i$, whereas Model 1 had the same $\theta$ for all 6 studies.

• Prior: $\prod_{i=1}^6 \text{Beta}(\theta_i\,|\,\alpha,\beta)$.

• Hyperprior: $P(\alpha,\beta) = \text{Gamma}(\alpha\,|\,2,0.5)\,\text{Gamma}(\beta\,|\,2,0.5)$.

This model has 8 parameters (for each of the treatment and control groups), namely $\theta_1, \ldots, \theta_6$, $\alpha$, and $\beta$.

Since the posterior does not have a closed-form analytical solution, we have to calculate the marginal likelihood by integrating out all of the parameters in the model.

$$P(\text{data}) = \int_0^{\infty} \int_0^{\infty} \int_0^1\cdots\int_0^1 \left[\prod_{i=1}^6 \text{Binomial}(s_i\,|\,\theta_i,n_i)\,\text{Beta}(\theta_i\,|\,\alpha,\beta) \right] P(\alpha,\beta)\ \text{d}\theta_6\cdots\text{d}\theta_1\,\text{d}\beta\,\text{d}\alpha$$

This looks like a crazy 8-dimensional integral, but we can actually integrate out the $\theta_i$ analytically, leaving a 2-dimensional integral over $\alpha$ and $\beta$.

First, note that $P(\alpha,\beta)$ does not contain $\theta_i$, so we can move it outside of the $\theta_i$ integrals.

ParseError: KaTeX parse error: Invalid color: '[' at position 61: …eta) \textcolor[̲rgb]{0.00,0.00,…

Next, since there are no factors containing two different $\theta_i$ variables (go to explanation), we can rearrange the integrals and the products (the blue part) like this:

ParseError: KaTeX parse error: Invalid color: '[' at position 61: …eta) \textcolor[̲rgb]{0.00,0.00,…

Note that we cannot always swap products and integrals.

Since the beta distribution is conjugate to the binomial, the blue integrals above can be evaluated analytically (much like we did for Model 1), to get

$$= \int_0^{\infty}\int_0^{\infty} P(\alpha,\beta) \left[\prod_{i=1}^6 \frac{1}{(s_i+f_i+1)\,\text{B}(s_i+1,f_i+1)}\,\frac{\text{B}(\alpha+s_i, \beta+f_i)}{\text{B}(\alpha,\beta)}\right]\,\text{d}\beta\,\text{d}\alpha$$

Finally, move all the factors that do not depend on $\alpha$ or $\beta$ out of the integrals.

$$= \left[\prod_{i=1}^6 (s_i+f_i+1)\,\text{B}(s_i+1,f_i+1) \right]^{-1} \int_0^{\infty}\int_0^{\infty} P(\alpha,\beta)\, \text{B}(\alpha,\beta)^{-6} \prod_{i=1}^6 \text{B}(\alpha+s_i, \beta+f_i)\ \text{d}\beta\,\text{d}\alpha$$

Unfortunately we cannot evaluate the remaining integrals analytically, so we resort to a numerical calculation.

Task 2:

• Read up about the adaptIntegrate(f, lowerLimit, upperLimit) function in the cubature package — see the function documentation below.

• How would you define a function $f$ so that adaptIntegrate(f, c(0, 0), c(Inf, Inf)) evaluates the 2-dimensional integral over $\alpha$ and $\beta$ above?

library(cubature)
?adaptIntegrate