---
title: "`MCMCpack` Software Tutorial"
author: "Huiling Liu"
output:
  ioslides_presentation:
    smaller: yes
---

## MCMCpack Introduction
- An R package that contains functions to perform Bayesian inference using posterior simulation.
- Aimed primarily at social scientists.
- The implementation of MCMC algorithms are model-specific.



## Functions in MCMCpack

Statistical Models:

- Linear regression models (linear regression with Gaussian errors, a singular value decomposition regression, and regression for a censored dependent variable)
- Discrete choice models (logistic/multinomial logistic/ordinal probit/probit regression)
- Measurement models (one-dimensional IRT,k-dimensional IRT, k-dimensional ordinal factor, k-dimensional linear factor, k-dimensional mixed factor, and k-dimensional robust IRT)
- A model for count data (Poisson regression model)
- Models for ecological inference (hierarchical ecological inference model and dynamic ecological inference model)
- Time-series models for change-point problems (binary/ probit/ ordinal probit/ Poisson change-point model).

Also contains some useful utility functions, including some additional density functions and pseudo-random number generators for statistical distributions, a general purpose Metropolis sampling algorithm, and tools for visualization.


## Installation of MCMCpack
```{r,eval=FALSE}
install.packages(MCMCpack)
```

```{r,message=FALSE}
library(MCMCpack)
```

## Ex1: Linear Regression Model
The model takes the following form:

$$y_i = x_i^{'}\beta + \varepsilon_{i}, \quad  \varepsilon_{i} \sim \mathcal{N}(0, 1/\tau)$$

We assume standard, semi-conjugate priors:

$\beta \sim \mathcal{N}(b_0,B_0^{-1})$ and $\tau \sim \mathcal{G}amma(c_0/2, d_0/2)$

where $\beta$ and $\tau$ are assumed a priori independent. 

## 
```{r,eval=FALSE}
MCMCregress(formula, data, burnin = 1000, mcmc = 10000, b0 = 0, B0 = 0,
            
            marginal.likelihood = c("none", "Laplace", "Chib95"), ...)
```

burnin:	The number of burn-in iterations for the sampler.

mcmc:	The number of MCMC iterations after burnin.

b0:	The prior mean of $\beta$. 

B0:	The prior precision of $\beta$.Default value of 0 is equivalent to an improper uniform prior for beta.

marginal.likelihood:  "none" in which case the marginal likelihood will not be calculated, "Laplace" in which case the Laplace approximation is used, and "Chib95" in which case the method of Chib (1995) is used.


## Birthwt example
```{r}
data(birthwt)
post.lm1 <- MCMCregress(bwt~age+lwt+as.factor(race) + smoke + ht,
                       data=birthwt,mcmc=50000,
                       b0=c(2700, 0, 0, -500, -500, -500, -500),
                       B0=c(1e-6, .01, .01, 1.6e-5, 1.6e-5, 1.6e-5,
                            1.6e-5),  marginal.likelihood="Chib95")
```

##
```{r}
summary(post.lm1)
```

## Model Comparison
Suppose that the observed data y could have been generated under one of two models $M_1$ and $M_2$.

A natural thing to ask from the Bayesian perspective is, âwhat is the posterior probability that $M_1$ is true (assuming either $M_1$ or $M_2$ is true)?â

Using Bayes theorem we can write:

$$Pr(M_k|y)=\frac{p(y|M_k)Pr(M_k)}{p(y|M_1)Pr(M_1)+p(y|M_1)Pr(M_1)}, \quad k=1,2$$

$p(y|M_1)$ and $p(y|M_2)$ here are marginal liklihoods.

## Model Comparison
It is instructive to look at the posterior odds in favor of one model (say $M_1$):
        $$\frac{Pr(M_1|y)}{Pr(M_2|y)}=\frac{p(y|M_1)}{p(y|M_2)} \times \frac{Pr(M_1)}{Pr(M_2)}$$                                                                        
Define Bayes factor for $M_1$ relative to $M_2$ as

$$B_{12}=\frac{p(y|M_1)}{p(y|M_2)} $$
                                                                    
Since the posterior odds equal the Bayes factor when the models are equally likely a priori, the Bayes factor is a measure of how much support is available in the data for one model relative to another.


## Model Comparison Example
PostProbMod: calculate the posterior probability that each model under study is correct given that one of the models under study is correct.
```{r}
post.lm1 <- MCMCregress(bwt~age+lwt+as.factor(race) + smoke + ht,
                       data=birthwt,mcmc=50000,
                       b0=c(2700, 0, 0, -500, -500, -500, -500),
                       B0=c(1e-6, .01, .01, 1.6e-5, 1.6e-5, 1.6e-5,
                            1.6e-5),  marginal.likelihood="Chib95")
post.lm2 <- MCMCregress(bwt~age+lwt+as.factor(race) + smoke,
                       data=birthwt,mcmc=50000,
                       b0=c(2700, 0, 0, -500, -500, -500),
                       B0=c(1e-6, .01, .01, 1.6e-5, 1.6e-5, 1.6e-5),                                               marginal.likelihood="Chib95")
BF <- BayesFactor(post.lm1,post.lm2)
mod.probs <- PostProbMod(BF)
```

## Model Comparison Example
```{r}
print(mod.probs)
```

## Ex2: Item Response Theory (IRT) Model

- $y_{ij}$ represents the observed choice by subject $i$ on item $j$.

- each subject has an ability (ideal point) dentoed by $\theta_{i(K \times 1)}$.

- each item has a difficulty parameter $\alpha_j$ and discrimination parameter $\beta_{j(K \times 1)}$. 

Consider 

$$y_{ij} \sim Bernoulli(\pi_{ij}),\quad  i=1,...,I \quad  j=1,...,J$$
$$\pi_{ij}=\Phi (-\alpha_{j}+\beta_j^T\theta_i)$$

where $\Phi$ is the the CDF of the standard normal distribution.

Prior distribution for the model parameters are:
$$(\alpha_j,\beta_j)^T \sim N(a_0,A_0^{-1}) \quad j=1,...,J$$
and $$\theta_i \sim N(0,1) \quad i=1,...,I$$

## 
```{r,eval=FALSE}
MCMCirt1d(datamatrix, theta.constraints = list(), burnin = 1000, mcmc = 20000)
```
- datamatrix:	The matrix of data. Must be 0, 1, or missing values. The rows of datamatrix correspond to subjects and the columns correspond to items.

- theta.constraints: A list specifying possible simple equality or inequality constraints on the ability parameters. A typical entry in the list has one of three forms: 

(1). varname=c which will constrain the ability parameter for the subject named varname to be equal to c, 

(2). varname="+" which will constrain the ability parameter for the subject named varname to be positive, 

(3). varname="-" which will constrain the ability parameter for the subject named varname to be negative.


## SupremeCourt Example

- Data for 9 justices Rehnquist, Stevens, O'Connor, Scalia, Kennedy, Souter, Thomas, Ginsburg, and Breyer for the U.S. Supreme Court from 43 non-unanimous cases.
- The votes are coded liberal (1), conservative (0) or missing values.
```{r,message=FALSE}
data(SupremeCourt)
post.irt <- MCMCirt1d(t(SupremeCourt),
  theta.constraints=list(Stevens="-",
  Scalia="+"), burnin=5000, mcmc=100000)
```


##
```{r}
summary(post.irt)
```

## 
```{r,echo=0}
library(grDevices)
theta.post <- as.matrix(post.irt)
set.seed(1)
x <- rnorm(20000,0,1)
plot(density(x),ylim=c(0,2),xlim=c(-4,4),col=gray(6/8),lwd=2,xlab="Ideal Point",ylab="f",main="")
cl <- colors()
col <- c()
r <- c(2,9,8,6,3,5,1,7,4)
for(i in 1:9){
  lines(density(theta.post[,r[i]]),col=rgb(red=31*(i-1),blue=248-31*(i-1),green=0,maxColorValue = 248),lwd=2)
  col[i] <- rgb(red=31*(i-1),blue=248-31*(i-1),green=0,maxColorValue = 248)
}
legend("topleft",col=c(col,gray(6/8)),c("stevens","breyer","ginsburg","souter","o'connor","kennedy","rehnquist","thomas","scalia","prior"),lty=c(1,1,1,1,1,1,1,1,1,1),lwd=2,cex=0.6)
```

## EX3:Metropolis Sampling for a User-Defined Model
MCMCpack also has some facilities for fitting user-specified models.

The MCMCmetrop1R() function uses a random walk Metropolis algorithm to
sample from a user-defined log-posterior density.

## 
```{r,eval=FALSE}
MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000,
  tune = 1, logfun = TRUE, force.samp = FALSE,
  V = NULL, optim.method = "BFGS", ...)
```

fun: The unnormalized (log)density of the distribution from which to take a sample.

theta.init: Starting values for the sampling. Must be of the appropriate dimension.
It must also be the case that fun(theta.init, ...) is greater than -Inf if fun() is a logdensity or greater than 0 if fun() is a density.

tune: The tuning parameter for the Metropolis sampling. Can be either a positive scalar or a k-vector, where k is the length of Î¸.

## Example for a User-Defined Model
Suppose one is interested in fitting a logistic regression with an improper uniform prior. 

One could do the following:
```{r}
 logitfun <- function(beta, y, X){
      eta <- X %*% beta
      p <- 1.0/(1.0+exp(-eta))
      sum( y * log(p) + (1-y)*log(1-p) )
    }

x1 <- rnorm(1000)
x2 <- rnorm(1000)
Xdata <- cbind(1,x1,x2)
p <- exp(.5 - x1 + x2)/(1+exp(.5 - x1 + x2))
yvector <- rbinom(1000, 1, p)
```

```{r}
post.samp <- MCMCmetrop1R(logitfun, theta.init=c(0,0,0),
                              X=Xdata, y=yvector,
                              mcmc=40000, burnin=500,
                              tune=c(1.5, 1.5, 1.5),
                              logfun=TRUE)
```

## 
```{r}
plot(post.samp)
```

## 
```{r}
summary(post.samp)
```
