##############################################
###       BINARY STUDENT REGRESSION        ###
##############################################
### tested on R version 2.5.1 (2007-06-27) ###
##############################################

# Package by M. Plummer, N. Best, K. Cowles and K. Vines: output analysis and diagnostics for MCMC
library(coda)
# Package by J. Albert: functions for learning Bayesian inference
library(LearnBayes)
# note that package 'coda' requires package 'lattice'

######## Data analyzed by Albert (2007, 10.3) ########
dataname <- "Donner survival data"
data(donner)
### maximum likelihood fit for a model of interest ###
mlfit<-glm(survival ~ age + male, data = donner, family = binomial(link = probit), x = TRUE, y = TRUE)
############### prior hyper-parameters ###############
 g <- 45
T0 <- (1/g) * t(mlfit$x) %*% mlfit$x
b0 <- rep(0, ncol(mlfit$x))
nusupp <- c(4, 8, 12, 16, 20)
########### length of Monte Carlo samples ############
mclen <- 10000
######################################################

# Bayesian t (link binary) regression
# MCMC algorithm based on Albert & Chib (1993)
# (log) marginal likelihood computation based on Chib (1995) 
btreg <- function(
                    y, # vector of responses
                    X, # design matrix
                    m, # number of simulations desired
                    b0 = rep(0, ncol(X)), # prior mean vector of regression coefficients
                    T0 = diag(0, ncol(X)), # prior precision matrix of regression coefficients
                    nusupp = c(4, 8, 16, 32) # support values for the degrees of freedom
                 )
{
    # number of regression coefficients
    k <- ncol(X)
    # number of observations
    n <- nrow(X)
    # empty chains
    beta <- matrix(NA, k, m)
    rownames(beta) <- colnames(X)
    z <- matrix(NA, n, m)
    lambda <- matrix(NA, n, m)
    rownames(lambda) <- rownames(X)
    nu <- rep(NA, m)
    # initial values
    fit <- glm(y ~ 0 + X, family = binomial(link = probit))
    beta[,1] <- fit$coef
    z[,1] <- X %*% beta[,1]
    lambda[,1] <- rep(1, n)
    nu[1] <- max(nusupp)
    # Gibbs sampling
    for (i in 2:m)
    {
        # regression coefficients
        B1 <- chol2inv(chol(T0 + t(X) %*% diag(lambda[,i-1]) %*% X))
        b1 <- B1 %*% (T0 %*% b0 + t(X) %*% diag(lambda[,i-1]) %*% z[,i-1])
        beta[,i] <- rmnorm(1, mean = b1, varcov = B1)
        # latent scores
        mu <- X %*% beta[,i]
        aux <- pnorm(-mu, sd = 1 / sqrt(lambda[,i-1]))
        z[,i] <- qnorm((aux * (1 - y) + (1 - aux) * y) * runif(n) + aux * y, sd = 1 / sqrt(lambda[,i-1])) + mu
        # latent precisions
        e <- z[,i] - X %*% beta[,i]
        for(j in 1:n) lambda[j,i] <- rgamma(1, shape = (nu[i-1] + 1) / 2, rate = (nu[i-1] + e[j]^2) / 2)
        # degrees of freedom
        lognuprobs <- n * (0.5 * nusupp * log(0.5 * nusupp) - lgamma(0.5 * nusupp)) +
                      0.5 * nusupp * sum(log(lambda[,i])) - 0.5 * nusupp * sum(lambda[,i])
        lognuprobs <- lognuprobs - mean(lognuprobs)
        nu[i] <- sample(nusupp, 1, prob = exp(lognuprobs))
    }
    # marginal likelihood from the Gibbs output
    betaStar <- apply(beta, 1, mean)
    nuStarPos <- which.max(table(nu))
    psiStar <- pt(X %*% betaStar, nusupp[nuStarPos])
    # full conditional for selected parameter value
    fullCondNuStar <- rep(NA, m)
    for(i in 1:m)
    {
        lognuprobs <- n * (0.5 * nusupp * log(0.5 * nusupp) - lgamma(0.5 * nusupp)) +
                      0.5 * nusupp * sum(log(lambda[,i])) - 0.5 * nusupp * sum(lambda[,i])
        lognuprobs <- lognuprobs - mean(lognuprobs)
        fullCondNuStar[i] <- exp(lognuprobs)[nuStarPos] / sum(exp(lognuprobs))
    }
    fullCondNuMean <- mean(fullCondNuStar)
    # additional sampling
    extrabeta <- matrix(NA, k, m)
    extrabeta[,1] <- beta[,m]
    extralambda <- matrix(NA, n, m)
    extralambda[,1] <- lambda[,m]
    extrazeta <- matrix(NA, n, m)
    extrazeta[,1] <- z[,m]
    for (i in 2:m)
    {
        # regression coefficients
        B1 <- chol2inv(chol(T0 + t(X) %*% diag(extralambda[,i-1]) %*% X))
        b1 <- B1 %*% (T0 %*% b0 + t(X) %*% diag(extralambda[,i-1]) %*% extrazeta[,i-1])
        extrabeta[,i] <- rmnorm(1, mean = b1, varcov = B1)
        # latent scores
        mu <- X %*% extrabeta[,i]
        aux <- pnorm(-mu, sd = 1 / sqrt(extralambda[,i-1]))
        extrazeta[,i] <- qnorm((aux * (1 - y) + (1 - aux) * y) * runif(n) + aux * y,
                               sd = 1 / sqrt(extralambda[,i-1])) + mu
        # latent precisions
        e <- extrazeta[,i] - X %*% extrabeta[,i]
        for(j in 1:n) extralambda[j,i] <- rgamma(1, shape = (nusupp[nuStarPos] + 1) / 2,
                                                    rate = (nusupp[nuStarPos] + e[j]^2) / 2)
    }
    # second full conditional (for selected parameter value)
    fullCondBetaStar <- rep(NA, m)
    for(i in 1:m)
    {
        B1 <- chol2inv(chol(T0 + t(X) %*% diag(extralambda[,i]) %*% X))
        b1 <- B1 %*% (T0 %*% b0 + t(X) %*% diag(extralambda[,i]) %*% extrazeta[,i])
        fullCondBetaStar[i] <-  dmnorm(betaStar, mean = b1, varcov = B1)
    }
    fullCondBetaMean <- mean(fullCondBetaStar)
    # marginal likelihood on the log scale
    logmarlik <- dmnorm(betaStar, mean = b0, varcov = chol2inv(chol(T0)), log = TRUE) - log(length(nusupp))
    logmarlik <- logmarlik - log(fullCondNuMean) - log(fullCondBetaMean)
    for(i in 1:n)
        logmarlik <- logmarlik + dbinom(y[i], size = 1, prob = psiStar[i], log = TRUE)
    names(logmarlik) <- NULL
    # numerical standard error
    logmarSE <- sqrt(var(fullCondNuStar) / (effectiveSize(fullCondNuStar) * fullCondNuMean^2) +
                     var(fullCondBetaStar) / (effectiveSize(fullCondBetaStar) * fullCondBetaMean^2))
    names(logmarSE) <- NULL
    # return the beta and nu chain, the medians of the lambda chains,
    # and the log marginal likelihood (with standard error)
    return(list(beta = t(beta), nu = nu, lambda = apply(lambda, 1, median), logmarlik = logmarlik, logmarSE = logmarSE))
}

# Bayesian fit
cat(system.time(bfit <- btreg(mlfit$y, mlfit$x, mclen, b0 = b0, T0 = T0, nusupp = nusupp)), fill = TRUE)

# parameter summaries
parsum <- matrix(NA, ncol(mlfit$x), 8)
rownames(parsum) <- colnames(mlfit$x)
colnames(parsum) <- c("PriorMean", "PriorSD","PostMean", "PostSD","Post0.025", "Post0.975", "IF(AT)","MLE")
# prior values
if(1/g == 0)
{
    parsum[,2] <- Inf # reference analysis
}else
{
    parsum[,1] <- b0 # prior mean
    parsum[,2] <- sqrt(diag(chol2inv(chol(T0)))) # prior std dev
}
# posterior values
parsum[,3] <- apply(bfit$beta, 2, mean)
parsum[,4] <- apply(bfit$beta, 2, sd)
parsum[,c(5,6)] <- t(apply(bfit$beta, 2, quantile, c(0.025, 0.975)))
# inefficiency factor (ratio of actual to effective sample size)
parsum[,7] <- mclen / apply(bfit$beta, 2, effectiveSize)
# maximum likelihood estimate
parsum[,8] <- coef(mlfit)

# text output to file
sink("binStuReg.txt")
cat("\n")
cat("Bayesian analysis of\n   ", dataname, "\n")
cat("using model\n    ")
show(mlfit$formula)
cat("and Monte Carlo sample size\n   ", mclen, "\n")
cat("Regression coefficients:\n")
show(parsum)
cat("Degrees of freedom:")
show(table(bfit$nu) / sum(table(bfit$nu)))
cat("Log marginal likelihood:", bfit$logmarlik, "with numerical standard error", bfit$logmarSE)
cat("\n")
sink()

# outlier plot to file
pdf("binStuReg.pdf", height = 12, width = 12)
dotchart(bfit$lambda, main = "Outlier plot", xlab = "Posterior median precision",
         bg = "red", lcolor = "black", pch = 22, xlim = c(0, max(bfit$lambda)))
dev.off()