Pajor vs Glmer Simulation Results.Rmd

---
title: "Pajor vs Glmer Simulation Results"
author: "Hillary Heiling"
date: "October 11, 2019"
output: html_document
---

## Description of Pajor Method

Paper: 

Pajor, Anna. Estimating the Marginal Likelihood Using the Arithmetic Mean Identity. Bayesian Anal. 12 (2017), no. 1, 261--287. doi:10.1214/16-BA1001. https://projecteuclid.org/euclid.ba/1459772735

### Notation

$\theta$ = random effects

$f(y|\beta, \gamma, X, \theta)$ = likelihood of y given fixed and random effect parameters, the covariates, and the draws of the random effects

$p(\theta)$ = prior, standard MVN

$s(\theta)$ = importance function

$A$ = subset of the parameter space of $\theta$

$P(A|y)$ = posterior probability of being in subset A

Suggestion by author: Let $ P(A|y) = 1 $ such that we restrict ourselves to the range (min, max) of the posterior samples for each element of the posterior.

### Regular Corrected Arithmetic Mean (CAM) Formula

$$ f(y) = \frac{1}{P(A|y)} \int_\Theta f(y|\beta,\gamma, X, \theta) I_A(\theta) p(\theta) d\theta $$
Draws of $\theta$ from prior $(\theta)$. The $I_A(\theta)$ function is the indicator function of the draws of $\theta$ falling into the region of A

### Importance Sample version of CAM (IS CAM) Formula 

$$ f(y|\beta, \gamma, X) = \frac{1}{P(A|y)} \int_\Theta f(y|\beta, \gamma, X, \theta) I_A(\theta) \frac{p(\theta)}{s(\theta)} s(\theta) d\theta $$
Draws of $\theta$ from importance function $s(\theta)$

### IS CAME Estimate

$$ f(y | \beta, \gamma, X) \approx \prod_{k=1}^d \left ( \frac{1}{M} \sum_{m=1}^M \left ( \prod_{i \in k} f(y|\beta, \gamma, X, \theta_{km}) \right ) I_A(\theta_{km}) \frac{p(\theta_{km})}{s_k(\theta_{km})} \right ) $$

$$ log(f(y | \beta, \gamma, X)) \approx \sum_{k=1}^d log \left ( \frac{1}{M} \sum_{m=1}^M \left ( \prod_{i \in k} f(y|\beta, \gamma, X, \theta_{km}) \right ) I_A(\theta_{km}) \frac{p(\theta_{km})}{s_k(\theta_{km})} \right ) $$

Draws of $\theta_{km}$ from the importance sampling function: group-specific, with mean = posterior mean for group k, and covariance = posterior covariance for group k.

When posterior from gibbs sampling, thinned samples for calculation of posterior covariance.

## Pajor vs Glmer Simulation Results

Simulated Data:

* N = 500
* K = 5
* Random effect sd = {0.5, 1.0, 2.0}
* Random intercept + 2 random slopes
* Each sd situation repeated 100 times

Number of draws from the posterior: 10^4

Number of draws from the importance function for each group: M = {10^4, 10^5, 3*10^5}

Value calculated: (ll_Pajor - ll_glmer) / ll_glmer

```{r, echo=FALSE}
load("sim_table2_Pajor_output.RData")

out_1 = output
s_l = 100
```

Boxplots of results:

```{r, echo=FALSE}
par(mfrow = c(1,2))
boxplot(out_1[1:s_l,1], out_1[1:s_l,2], out_1[1:s_l,3], names = c("M 10^4", "M 10^5", "M 3x10^5"),
        main = "Results for SD = 0.5",
        ylab = "(ll_Pajor - ll_glmer)/ll_glmer", ylim = c(0, 1.5))
boxplot(out_1[(s_l+1):200,1], out_1[(s_l+1):200,2], out_1[(s_l+1):200,3], 
        names = c("M 10^4", "M 10^5", "M 3x10^5"),
        main = "Results for SD = 1.0",
        ylab = "(ll_Pajor - ll_glmer)/ll_glmer", ylim = c(0, 1.5))



```

Further examination of SD = 0.5 situation:

```{r, echo=FALSE}
boxplot(out_1[1:s_l,1], out_1[1:s_l,2], out_1[1:s_l,3], names = c("M 10^4", "M 10^5", "M 3x10^5"),
        main = "Results for SD = 0.5",
        ylab = "(ll_Pajor - ll_glmer)/ll_glmer", ylim = c(0, 0.25))
```


## Diagnostic Output

```{r, echo=FALSE}
load("glmer_Pajor_diagnostics.RData")

diagnostic = output

# d_out = diagnotic output
# sd_ranef = 0.5
d_out1 = diagnostic[1:100]

# sd_ranef = 1.0
d_out2 = diagnostic[101:200]

# Parameters
d = 5
q = 3
```

### Acceptance Rates - All

Acceptance rates for sd_ranef = 0.5

```{r, echo=FALSE}
library(stringr)
# Acceptance rates for sd_ranef = 0.5
acc_rate1 = matrix(0, nrow = 100, ncol = d*q)

vars_q = rep(c("(Intercept)","X1","X2"), each = d)
grp = rep(1:5, times = q)
colnames(acc_rate1) = str_c(vars_q, ":", grp)

for(i in 1:100){
  gibbs_acc_rate = d_out1[[i]]$gibbs_accept_rate
  rate_vec = c(gibbs_acc_rate[,1], gibbs_acc_rate[,2], gibbs_acc_rate[,3])
  acc_rate1[i,] = matrix(rate_vec, nrow = 1)
}

```

```{r, echo=FALSE}
par(mfrow = c(1,3))
boxplot(acc_rate1[,1], acc_rate1[,2], acc_rate1[,3], acc_rate1[,4], acc_rate1[,5], 
        names = colnames(acc_rate1)[1:5],
        main = "Intercept: Accept Rates for SD = 0.5",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
boxplot(acc_rate1[,6], acc_rate1[,7], acc_rate1[,8], acc_rate1[,9], acc_rate1[,10], 
        names = colnames(acc_rate1)[6:10],
        main = "X1: Accept Rates for SD = 0.5",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
boxplot(acc_rate1[,11], acc_rate1[,12], acc_rate1[,13], acc_rate1[,14], acc_rate1[,15], 
        names = colnames(acc_rate1)[11:15],
        main = "X2: Accept Rates for SD = 0.5",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))

```

Acceptance rates for sd_ranef = 1.0

```{r, echo=FALSE}

# Acceptance rates for sd_ranef = 1.0
acc_rate2 = matrix(0, nrow = 100, ncol = d*q)

vars_q = rep(c("(Intercept)","X1","X2"), each = d)
grp = rep(1:5, times = q)
colnames(acc_rate2) = str_c(vars_q, ":", grp)

for(i in 1:100){
  gibbs_acc_rate = d_out2[[i]]$gibbs_accept_rate
  rate_vec = c(gibbs_acc_rate[,1], gibbs_acc_rate[,2], gibbs_acc_rate[,3])
  acc_rate2[i,] = matrix(rate_vec, nrow = 1)
}

full_accept_rates_sd1 = acc_rate2

```

```{r, echo=FALSE}
par(mfrow = c(1,3))
boxplot(acc_rate2[,1], acc_rate2[,2], acc_rate2[,3], acc_rate2[,4], acc_rate2[,5], 
        names = colnames(acc_rate2)[1:5],
        main = "Intercept: Accept Rates for SD = 1.0",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
boxplot(acc_rate2[,6], acc_rate2[,7], acc_rate2[,8], acc_rate2[,9], acc_rate2[,10], 
        names = colnames(acc_rate2)[6:10],
        main = "X1: Accept Rates for SD = 1.0",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
boxplot(acc_rate2[,11], acc_rate2[,12], acc_rate2[,13], acc_rate2[,14], acc_rate2[,15], 
        names = colnames(acc_rate2)[11:15],
        main = "X2: Accept Rates for SD = 1.0",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))

```

### Acceptance Rates - Percent Difference for ll_Pajor vs ll_glmer > 10%

Acceptance rates for sd_ranef = 0.5

```{r, echo=FALSE}

# Acceptance rates for sd_ranef = 0.5
acc_rate1 = matrix(0, nrow = 100, ncol = d*q)

vars_q = rep(c("(Intercept)","X1","X2"), each = d)
grp = rep(1:5, times = q)
colnames(acc_rate1) = str_c(vars_q, ":", grp)

for(i in 1:100){
  gibbs_acc_rate = d_out1[[i]]$gibbs_accept_rate
  ll_pdif = d_out1[[i]]$ll_pdif
  
  rate_vec = c(gibbs_acc_rate[,1], gibbs_acc_rate[,2], gibbs_acc_rate[,3])
  if(any(ll_pdif > 0.1)){
    rate_vec2 = rate_vec
  }else{
    rate_vec2 = rep(NA, times = d*q)
  }
  acc_rate1[i,] = matrix(rate_vec2, nrow = 1)
}

```

```{r, echo=FALSE}
par(mfrow = c(1,3))
boxplot(acc_rate1[,1], acc_rate1[,2], acc_rate1[,3], acc_rate1[,4], acc_rate1[,5], 
        names = colnames(acc_rate1)[1:5],
        main = "Intercept: Accept Rates for SD = 0.5",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
boxplot(acc_rate1[,6], acc_rate1[,7], acc_rate1[,8], acc_rate1[,9], acc_rate1[,10], 
        names = colnames(acc_rate1)[6:10],
        main = "X1: Accept Rates for SD = 0.5",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
boxplot(acc_rate1[,11], acc_rate1[,12], acc_rate1[,13], acc_rate1[,14], acc_rate1[,15], 
        names = colnames(acc_rate1)[11:15],
        main = "X2: Accept Rates for SD = 0.5",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))

```

Acceptance rates for sd_ranef = 1.0

```{r, echo=FALSE}

# Acceptance rates for sd_ranef = 1.0
acc_rate2 = matrix(0, nrow = 100, ncol = d*q)

vars_q = rep(c("(Intercept)","X1","X2"), each = d)
grp = rep(1:5, times = q)
colnames(acc_rate2) = str_c(vars_q, ":", grp)

for(i in 1:100){
  gibbs_acc_rate = d_out2[[i]]$gibbs_accept_rate
  ll_pdif = d_out2[[i]]$ll_pdif
  
  rate_vec = c(gibbs_acc_rate[,1], gibbs_acc_rate[,2], gibbs_acc_rate[,3])
  if(any(ll_pdif > 0.1)){
    rate_vec2 = rate_vec
  }else{
    rate_vec2 = rep(NA, times = d*q)
  }
  
  acc_rate2[i,] = matrix(rate_vec2, nrow = 1)
}

```

```{r, echo=FALSE}
par(mfrow = c(1,3))
boxplot(acc_rate2[,1], acc_rate2[,2], acc_rate2[,3], acc_rate2[,4], acc_rate2[,5], 
        names = colnames(acc_rate2)[1:5],
        main = "Intercept: Accept Rates for SD = 1.0",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
boxplot(acc_rate2[,6], acc_rate2[,7], acc_rate2[,8], acc_rate2[,9], acc_rate2[,10], 
        names = colnames(acc_rate2)[6:10],
        main = "X1: Accept Rates for SD = 1.0",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
boxplot(acc_rate2[,11], acc_rate2[,12], acc_rate2[,13], acc_rate2[,14], acc_rate2[,15], 
        names = colnames(acc_rate2)[11:15],
        main = "X2: Accept Rates for SD = 1.0",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
```

### Acceptance Rates - Percent Difference for ll_Pajor vs ll_glmer > 30%

Acceptance rates for sd_ranef = 1.0

```{r, echo=FALSE}

# Acceptance rates for sd_ranef = 1.0
acc_rate2 = matrix(0, nrow = 100, ncol = d*q)

vars_q = rep(c("(Intercept)","X1","X2"), each = d)
grp = rep(1:5, times = q)
colnames(acc_rate2) = str_c(vars_q, ":", grp)

for(i in 1:100){
  gibbs_acc_rate = d_out2[[i]]$gibbs_accept_rate
  ll_pdif = d_out2[[i]]$ll_pdif
  
  rate_vec = c(gibbs_acc_rate[,1], gibbs_acc_rate[,2], gibbs_acc_rate[,3])
  if(any(ll_pdif > 0.3)){
    rate_vec2 = rate_vec
  }else{
    rate_vec2 = rep(NA, times = d*q)
  }
  
  acc_rate2[i,] = matrix(rate_vec2, nrow = 1)
}

```

```{r, echo=FALSE}
par(mfrow = c(1,3))
boxplot(acc_rate2[,1], acc_rate2[,2], acc_rate2[,3], acc_rate2[,4], acc_rate2[,5], 
        names = colnames(acc_rate2)[1:5],
        main = "Intercept: Accept Rates for SD = 1.0",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
boxplot(acc_rate2[,6], acc_rate2[,7], acc_rate2[,8], acc_rate2[,9], acc_rate2[,10], 
        names = colnames(acc_rate2)[6:10],
        main = "X1: Accept Rates for SD = 1.0",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
boxplot(acc_rate2[,11], acc_rate2[,12], acc_rate2[,13], acc_rate2[,14], acc_rate2[,15], 
        names = colnames(acc_rate2)[11:15],
        main = "X2: Accept Rates for SD = 1.0",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
```


### Acceptance Rates - Percent Difference for ll_Pajor vs ll_glmer < 5%

Acceptance rates for sd_ranef = 1.0

```{r, echo=FALSE}

# Acceptance rates for sd_ranef = 1.0
acc_rate2 = matrix(0, nrow = 100, ncol = d*q)

vars_q = rep(c("(Intercept)","X1","X2"), each = d)
grp = rep(1:5, times = q)
colnames(acc_rate2) = str_c(vars_q, ":", grp)

for(i in 1:100){
  gibbs_acc_rate = d_out2[[i]]$gibbs_accept_rate
  ll_pdif = d_out2[[i]]$ll_pdif
  
  rate_vec = c(gibbs_acc_rate[,1], gibbs_acc_rate[,2], gibbs_acc_rate[,3])
  if(any(ll_pdif < 0.05)){
    rate_vec2 = rate_vec
  }else{
    rate_vec2 = rep(NA, times = d*q)
  }
  
  acc_rate2[i,] = matrix(rate_vec2, nrow = 1)
}

```

```{r, echo=FALSE}
par(mfrow = c(1,3))
boxplot(acc_rate2[,1], acc_rate2[,2], acc_rate2[,3], acc_rate2[,4], acc_rate2[,5], 
        names = colnames(acc_rate2)[1:5],
        main = "Intercept: Accept Rates for SD = 1.0",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
boxplot(acc_rate2[,6], acc_rate2[,7], acc_rate2[,8], acc_rate2[,9], acc_rate2[,10], 
        names = colnames(acc_rate2)[6:10],
        main = "X1: Accept Rates for SD = 1.0",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
boxplot(acc_rate2[,11], acc_rate2[,12], acc_rate2[,13], acc_rate2[,14], acc_rate2[,15], 
        names = colnames(acc_rate2)[11:15],
        main = "X2: Accept Rates for SD = 1.0",
        ylab = "gibbs accept_rate", ylim = c(0, 1.0))
```

The End
The End