##############################################
###       STUDENT LINEAR 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, 9.2) ######
dataname <- "Bird measurements from British islands"
data(birdextinct)
obsnames <- birdextinct$species
### least square fit for a model of interest ###
birdextinct$logtime <- log(birdextinct$time)
lsfit<-lm(logtime ~ nesting + size + status, data = birdextinct, x = TRUE, y = TRUE)
############ prior hyper-parameters ############
 g <- 62
T0 <- (1/g) * t(lsfit$x) %*% lsfit$x
b0 <- rep(0, ncol(lsfit$x))
n0 <- 3
s0 <- 1
nu <- c(40, 4)
######### length of Monte Carlo samples ########
mclen <- 5000
################################################

# Bayesian linear model with Student errors
# (log) marginal likelihood computation based on Chib (1995) 
bstulm <- 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
                    n0 = 0, # prior sample size
                    s0 = 1, # prior guess of error standard deviation
                    nu = 4 # Student degrees of freedom
                  )
{
    # number of regression coefficients
    p <- ncol(X)
    # number of observations
    n <- nrow(X)
    # empty chains
    beta <- matrix(NA, p, m)
    rownames(beta) <- colnames(X)
    sigma <- rep(NA, m)
    lambda <- matrix(NA, n, m)
    # initial parameter values
    fit <- lm(y ~ 0 + X)
    beta[,1] <- fit$coef
    sigma[1] <- sqrt(sum(fit$residuals^2)/fit$df.residual)
    lambda[,1] <- rep(1, n)
    # posterior sample size
    n1 <- n0 + n
    # Gibbs sampling
    for (i in 2:m)
    {
        # regression coefficients
        tau <- (1 / sigma[i-1]^2)
        B1 <- chol2inv(chol(T0 + tau * t(X) %*% diag(lambda[,i-1]) %*% X))
        b1 <- B1 %*% (T0 %*% b0 + tau * t(X) %*% diag(lambda[,i-1]) %*% y)
        beta[,i] <- rmnorm(1, mean = b1, varcov = B1)
        # latent precisions
        e <- y - X %*% beta[,i]
        for(j in 1:n) lambda[j,i] <- rgamma(1, shape = (nu + 1) / 2, rate = (nu + tau * e[j]^2) / 2)
        # error standard deviation
        shape <- n1 / 2
        rate <- (n0 * s0^2 + t(e) %*% diag(lambda[,i]) %*% e) / 2
        sigma[i] <- sqrt(1 / rgamma(1, shape = shape, rate = rate))
    }
    # marginal likelihood from the Gibbs output
    betaStar <- apply(beta, 1, mean)
    tauStar <- mean(1 / sigma^2)
    muStar <- X %*% betaStar
    # full conditional for selected parameter value
    fullCondTauStar <- rep(NA, m)
    for(i in 1:m)
    {
        e <- y - X %*% beta[,i]
        fullCondTauStar[i] <-  dgamma(tauStar, shape = n1 / 2, rate = (n0 * s0^2 + t(e) %*% diag(lambda[,i]) %*% e) / 2)
    }
    fullCondTauMean <- mean(fullCondTauStar)
    # additional sampling
    extrabeta <- matrix(NA, p, m)
    extrabeta[,1] <- beta[,m]
    extralambda <- matrix(NA, n, m)
    extralambda[,1] <- lambda[,m]
    for (i in 2:m)
    {
        # regression coefficients
        B1 <- chol2inv(chol(T0 + tauStar * t(X) %*% diag(extralambda[,i-1]) %*% X))
        b1 <- B1 %*% (T0 %*% b0 + tauStar * t(X) %*% diag(extralambda[,i-1]) %*% y)
        extrabeta[,i] <- rmnorm(1, mean = b1, varcov = B1)
        # latent precisions
        e <- y - X %*% extrabeta[,i]
        for(j in 1:n) extralambda[j,i] <- rgamma(1, shape = (nu + 1) / 2, rate = (nu + tauStar * e[j]^2) / 2)
    }
    # second full conditional (for selected parameter value)
    fullCondBetaStar <- rep(NA, m)
    for(i in 1:m)
    {
        B1 <- chol2inv(chol(T0 + tauStar * t(X) %*% diag(extralambda[,i]) %*% X))
        b1 <- B1 %*% (T0 %*% b0 + tauStar * t(X) %*% diag(extralambda[,i]) %*% y)
        fullCondBetaStar[i] <-  dmnorm(betaStar, mean = b1, varcov = B1)
    }
    fullCondBetaMean <- mean(fullCondBetaStar)
    # marginal likelihood on the log scale
    logmarlik <- dgamma(tauStar, shape = n0 / 2, rate = n0 * s0^2 / 2, log = TRUE) +
                 dmnorm(betaStar, mean = b0, varcov = chol2inv(chol(T0)), log = TRUE)
    logmarlik <- logmarlik - log(fullCondTauMean) - log(fullCondBetaMean)
    for(i in 1:length(y))
        logmarlik <- logmarlik + dmt(y[i], mean = muStar[i], S = 1 / tauStar, df = nu, log = TRUE)
    # numerical standard error
        logmarSE <- sqrt(var(fullCondTauStar) / (effectiveSize(fullCondTauStar) * fullCondTauMean^2) +
                         var(fullCondBetaStar) / (effectiveSize(fullCondBetaStar) * fullCondBetaMean^2))
    # return the beta and sigma chains, the medians of the lambda chains,
    # and the log marginal likelihood (with standard error)
    return(list(beta = t(beta), sigma = sigma, lambda = apply(lambda, 1, median),
           logmarlik = logmarlik, logmarSE = logmarSE))
}

# Bayesian fit and parameter summaries: prior summaries, posterior summaries,
# latent precision (posterior) medians and log marginal likelihood
priorSum <- matrix(NA, ncol(lsfit$x)+1, 3)
rownames(priorSum) <- c(colnames(lsfit$x),"Sigma")
colnames(priorSum) <- c("PriorMean", "PriorSD", "LSE")
#
postSum <- array(NA, c(ncol(lsfit$x)+1, 5, length(nu)))
dimnames(postSum) <- list(c(colnames(lsfit$x),"Sigma"),
                         c("PostMean", "PostSD", "Post0.025", "Post0.975", "IF(AT)"),
                         paste("df =", nu))
#
latpremed <- matrix(NA, nrow(lsfit$x), length(nu))
rownames(latpremed) <- obsnames
colnames(latpremed) <- paste("df =", nu)
#
logmarlik <- matrix(NA, 2, length(nu))
rownames(logmarlik) <- c("Estimate", "NumStdErr")
colnames(logmarlik) <- paste("df =", nu)
# prior values for regression coefficient
if(1/g == 0)
{
    priorSum[1:ncol(lsfit$x),2] <- Inf # reference analysis
}else
{
    priorSum[1:ncol(lsfit$x),1] <- b0 # prior mean
    priorSum[1:ncol(lsfit$x),2] <- sqrt(diag(chol2inv(chol(T0)))) # prior std dev
}
# prior values for error standard deviation
if(n0 <= 1)
{
    priorSum[nrow(priorSum),2] <- Inf # prior mean undefined and infinite prior variance (reference analysis)
}else
{
    auxsgm<-sqrt(1 / rgamma(mclen, shape = n0/2, rate = n0 * s0^2 / 2))
    priorSum[nrow(priorSum),1] <- mean(auxsgm) # prior mean
    if(n0 <= 2)
    {
        priorSum[nrow(priorSum),2] <- Inf # prior mean defined, but infinite prior variance
    }else
    {
        priorSum[nrow(priorSum),2] <- sd(auxsgm) # prior std dev
    }
    rm(auxsgm)
}
# least square estimates
priorSum[1:ncol(lsfit$x),3] <- coef(lsfit)
priorSum[nrow(priorSum),3] <- sqrt(sum(lsfit$residuals^2) / lsfit$df.residual)
# values from MCMC output
for(i in 1:length(nu))
{
    cat(system.time(bfit <- bstulm(lsfit$y, lsfit$x, mclen, b0 = b0, T0 = T0, n0 = n0, s0 = s0, nu = nu[i])),
        fill = TRUE)
    # posterior means and std devs
    postSum[1:ncol(lsfit$x),1,i] <- apply(bfit$beta, 2, mean)
    postSum[nrow(postSum),1,i] <- mean(bfit$sigma)
    postSum[1:ncol(lsfit$x),2,i] <- apply(bfit$beta, 2, sd)
    postSum[nrow(postSum),2,i] <- sd(bfit$sigma)
    # posterior quantiles
    postSum[1:ncol(lsfit$x),c(3,4),i] <- t(apply(bfit$beta, 2, quantile, c(0.025, 0.975)))
    postSum[nrow(postSum),c(3,4),i] <- quantile(bfit$sigma, c(0.025, 0.975))
    # inefficiency factor (ratio of actual to effective sample size)
    postSum[1:ncol(lsfit$x),5,i] <- mclen / apply(bfit$beta, 2, effectiveSize)
    postSum[nrow(postSum),5,i] <- mclen / effectiveSize(bfit$sigma)
    # log marginal likelihood
    logmarlik[1,i] <- bfit$logmarlik
    logmarlik[2,i] <- bfit$logmarSE
    # latent precision medians
    latpremed[,i] <- bfit$lambda
}
rm(i)

# text output to file
sink("stuLinReg.txt")
cat("\n")
cat("Bayesian analysis of\n   ", dataname, "\n")
cat("using model\n    ")
show(formula(terms(lsfit)))
cat("and Monte Carlo sample size\n   ", mclen, "\n\n")
cat("Prior summaries:\n")
show(priorSum)
cat("\n")
cat("Posterior summaries:\n")
show(postSum)
cat("Log marginal likelihood:\n")
show(logmarlik)
cat("\n")
cat("Posterior medians of latent precisions:\n")
show(round(latpremed, 2))
sink()

# outlier plot to file
pdf("stuLinReg.pdf", height = 12, width = 12)
mycolors <- heat.colors(length(nu))
dotchart(latpremed[,1], main = "Outlier plot", xlab = "Posterior median precision",
         bg = mycolors[1], lcolor = "black", pch = 22, xlim = c(0, max(latpremed)))
if(length(nu) > 1) for (i in 2:length(nu)) points(latpremed[,i], 1:nrow(latpremed), pch = 22, bg = mycolors[i])
legend("topleft", legend = nu, fill = mycolors, title = "df")
rm(i, mycolors)
dev.off()