#' linear regression model Baysian residual plot
#'
#' @author Xiaofeng Wang, \email{wangx6@ccf.org}
#' @param inla.obj an object of class "inla".
#' @param covariate the covariate of interest for the residual plot.
#' @param m the size of posterior samples to be generated from the posterior distribution.
#' @param pmedian if pmedian = TRUE, posterior medians of Bayesian residuals will be calculated.
#' @param smooth if smooth = TRUE, a nonparametric smooth curve will be added.
#' @param xlab a title for the x axis.
#' @param cex.lab The magnification to be used for x and y labels relative to the current setting of cex.
#' @param cex.axis The magnification to be used for x and y labels relative to the current setting of cex.
#' @param ylab a title for the y axis.
#' @param main an overall title for the plot.
#' @param ...
#'
#' @details generate a Bayesian residual index plot or a plot of the residual versus the predictor for linear regression using INLA.
#' @return a list consisting of residuals (resid), the covariate used for plot (covariate), the response variable (y)
#' @export
bri.lmresid.plot <- function(inla.obj, covariate = NULL, m = 1000, pmedian = FALSE, smooth = FALSE, xlab = NULL, cex.lab = 1.25, cex.axis = 1.25, ylab = "Bayesian residual", main = "",...){
if(class(inla.obj) != "inla") stop("an 'inla' object is needed!")
if(inla.obj$.args$family != "gaussian") stop("the function only supports Gaussian linear regression model!")
# Get a single random draw of posterior distribution of parameters
p1 <- length(names(inla.obj$marginals.fixed))
beta <- matrix(0, nrow = m, ncol = p1)
for (j in 1:p1){
beta[, j] <- inla.rmarginal(m, marg = inla.obj$marginals.fixed[[j]])
}
yhat <- inla.obj$model.matrix %*% t(beta)
org.formula <- inla.obj$.args$formula
yname <- as.character(org.formula)[2]
y <- inla.obj$.args$data[, yname]
ymat <- matrix(rep(y, times = m), ncol = m)
resid.sample <- ymat - yhat
if (pmedian) {
resid <- apply(resid.sample, 1, median)
} else {
resid <- apply(resid.sample, 1, mean)
}
if (is.null(covariate)) {
covariate <- seq(1: length(resid))
xlab <- "Index"
} else {
if (length(resid) != length(covariate))
stop("the 'covariate' variable does not match with the residual object!")
if (is.null(xlab)) xlab <- "Covariate"
}
if (smooth){
res.dat <- data.frame(covariate = covariate, resid = resid)
res.smooth.inla <- inla(resid ~ -1 + f(covariate, model = 'rw2', constr = FALSE), data = res.dat)
bri.band.plot(res.smooth.inla, name = 'covariate', alpha = 0.05, xlab = xlab, ylab = ylab, main = main, cex.lab = cex.lab, cex.axis = cex.axis, type = 'random', ylim= c(min(resid), max(resid)))
points(res.dat$covariate, res.dat$resid)
} else{
plot(covariate, resid, cex.lab = cex.lab, cex.axis = cex.axis, xlab = xlab, main = main, ylab = ylab, ylim= c(min(resid), max(resid)), ...)
}
return(invisible(list(resid = resid, covariate = covariate, y = y)))
}
#' Baysian Pearson residuals for Poisson regression using INLA
#'
#' @author Xiaofeng Wang, \email{wangx6@ccf.org}
#' @param inla.obj an object of class "inla".
#' @param covariate the covariate of interest for the residual plot.
#' @param m the size of posterior samples to be generated from the posterior distribution.
#' @param pmedian if pmedian = TRUE, posterior medians of Bayesian residuals will be calculated.
#' @param plot if plot = TRUE, a Bayesian residual plot will be generated
#' @param smooth if smooth = TRUE, a nonparametric smooth curve will be added.
#' @param xlab a title for the x axis.
#' @param cex.lab The magnification to be used for x and y labels relative to the current setting of cex.
#' @param cex.axis The magnification to be used for x and y labels relative to the current setting of cex.
#' @param ylab a title for the y axis.
#' @param main an overall title for the plot.
#' @param ...
#' @details compute Bayesian Pearson residuals, and generate a index plot or a plot of the residual versus the predictor for Bayesian Possion regression using INLA.
#' @return a list consisting of residuals (resid), the covariate used for plot (covariate), the response variable (y)
#' @export
bri.Pois.resid <- function(inla.obj, covariate = NULL, m = 1000, pmedian = FALSE, plot = FALSE, smooth = FALSE, xlab = NULL, cex.lab = 1.25, cex.axis = 1.25, ylab = "Bayesian residual", ylim = NULL, main = "",...){
if(class(inla.obj) != "inla") stop("an 'inla' object is needed!")
if(inla.obj$.args$family != "poisson") stop("the function only supports Poisson GLM!")
# Get a single random draw of posterior distribution of parameters
p1 <- length(names(inla.obj$marginals.fixed))
beta <- matrix(0, nrow = m, ncol = p1)
for (j in 1:p1){
beta[, j] <- inla.rmarginal(m, marg = inla.obj$marginals.fixed[[j]])
}
yhat <- exp(inla.obj$model.matrix %*% t(beta))
org.formula <- inla.obj$.args$formula
yname <- as.character(org.formula)[2]
y <- inla.obj$.args$data[, yname]
ymat <- matrix(rep(y, times = m), ncol = m)
if (is.null(inla.obj$.args$E)) {
resid.sample <- (ymat - yhat)/sqrt(yhat)
} else {
Emat <- matrix(rep(inla.obj$.args$E, times = m), ncol = m)
yhat.E <- yhat*Emat
resid.sample <- (ymat - yhat.E)/sqrt(yhat.E)
}
if (pmedian) {
resid <- apply(resid.sample, 1, median)
} else {
resid <- apply(resid.sample, 1, mean)
}
if (plot) {
if (is.null(covariate)) {
covariate <- seq(1: length(resid))
xlab <- "Index"
} else {
if (length(resid) != length(covariate))
stop("the 'covariate' variable does not match with the residual object!")
if (is.null(xlab)) xlab <- "Covariate"
}
if (is.null(ylim)) {ylim <- c(min(resid), max(resid))}
if (smooth){
res.dat <- data.frame(covariate = covariate, resid = resid)
res.smooth.inla <- inla(resid ~ -1 + f(covariate, model = 'rw2', constr = FALSE), data = res.dat)
bri.band.plot(res.smooth.inla, name = 'covariate', alpha = 0.05, xlab = xlab, ylab = ylab, main = main, cex.lab = cex.lab, cex.axis = cex.axis, type = 'random', ylim = ylim)
points(res.dat$covariate, res.dat$resid)
} else{
plot(covariate, resid, cex.lab = cex.lab, cex.axis = cex.axis, xlab = xlab, main = main, ylab = ylab, ylim = ylim, ...)
}
}
return(invisible(list(resid = resid, covariate = covariate, y = y)))
}
#' Baysian residuals for beta regression using INLA
#'
#' @author Xiaofeng Wang, \email{wangx6@ccf.org}
#' @param inla.obj an object of class "inla".
#' @param covariate the covariate of interest for the residual plot.
#' @param m the size of posterior samples to be generated from the posterior distribution.
#' @param pmedian if pmedian = TRUE, posterior medians of Bayesian residuals will be calculated.
#' @param plot if plot = TRUE, a Bayesian residual plot will be generated
#' @param smooth if smooth = TRUE, a nonparametric smooth curve will be added.
#' @param xlab a title for the x axis.
#' @param cex.lab The magnification to be used for x and y labels relative to the current setting of cex.
#' @param cex.axis The magnification to be used for x and y labels relative to the current setting of cex.
#' @param ylab a title for the y axis.
#' @param main an overall title for the plot.
#' @param ...
#' @details compute Bayesian residuals, and generate a index plot or a plot of the residual versus the predictor for Bayesian beta regression using INLA.
#' @return a list consisting of residuals (resid), the covariate used for plot (covariate), the response variable (y)
#' @export
bri.beta.resid <- function(inla.obj, covariate = NULL, m = 1000, pmedian = FALSE, plot = FALSE, smooth = FALSE, xlab = NULL, cex.lab = 1.25, cex.axis = 1.25, ylab = "Bayesian residual", ylim = NULL, main = "",...){
if(class(inla.obj) != "inla") stop("an 'inla' object is needed!")
if(inla.obj$.args$family != "beta") stop("the function only supports beta GLM!")
# Get a single random draw of posterior distribution of parameters
p1 <- length(names(inla.obj$marginals.fixed))
beta <- matrix(0, nrow = m, ncol = p1)
for (j in 1:p1){
beta[, j] <- inla.rmarginal(m, marg = inla.obj$marginals.fixed[[j]])
}
yhat <- exp(inla.obj$model.matrix %*% t(beta))/(1+exp(inla.obj$model.matrix %*% t(beta)))
phihat <- inla.rmarginal(m, marg = inla.obj$marginals.hyperpar[[1]])
phihatmat <- t(matrix(rep(phihat, times = nrow(yhat)), ncol = nrow(yhat)))
varhat <- yhat*(1-yhat)/(1+phihatmat)
org.formula <- inla.obj$.args$formula
yname <- as.character(org.formula)[2]
y <- inla.obj$.args$data[, yname]
ymat <- matrix(rep(y, times = m), ncol = m)
resid.sample <- (ymat - yhat)/sqrt(varhat)
if (pmedian) {
resid <- apply(resid.sample, 1, median)
} else {
resid <- apply(resid.sample, 1, mean)
}
if (plot) {
if (is.null(covariate)) {
covariate <- seq(1: length(resid))
xlab <- "Index"
} else {
if (length(resid) != length(covariate))
stop("the 'covariate' variable does not match with the residual object!")
if (is.null(xlab)) xlab <- "Covariate"
}
if (is.null(ylim)) {ylim <- c(min(resid), max(resid))}
if (smooth){
res.dat <- data.frame(covariate = covariate, resid = resid)
res.smooth.inla <- inla(resid ~ -1 + f(covariate, model = 'rw2', constr = FALSE), data = res.dat)
bri.band.plot(res.smooth.inla, name = 'covariate', alpha = 0.05, xlab = xlab, ylab = ylab, main = main, cex.lab = cex.lab, cex.axis = cex.axis, type = 'random', ylim= ylim)
points(res.dat$covariate, res.dat$resid)
} else{
plot(covariate, resid, cex.lab = cex.lab, cex.axis = cex.axis, xlab = xlab, main = main, ylab = ylab, ylim = ylim, ...)
}
}
return(invisible(list(resid = resid, covariate = covariate, y = y)))
}
#' Plot the baseline hazard functions for survival models using INLA
#' @param inla.obj an object of class "inla".
#' @param plot if plot = TRUE, the plot of the baseline hazard functions will be generated.
#' @param cex.lab The magnification to be used for x and y labels relative to the current setting of cex.
#' @param cex.axis The magnification to be used for x and y labels relative to the current setting of cex.
#' @param ...
#'
#' @details The function supports Weibull model, Expontenial model and Cox proportional hazards model for right censored data.
#' @return A list of evaluated time points (time), and the values of corresponding base line function (basehaz).
#' @export
bri.basehaz.plot <- function(inla.obj, plot = TRUE, cex.lab = 1.25, cex.axis = 1.25, ...){
if(class(inla.obj) != "inla") stop("an 'inla' object is needed!")
family <- inla.obj$.args$family
switch(tolower(family),
poisson = {
if (length(inla.obj$.arg$data$baseline.hazard.values) == nrow(inla.obj$summary.random$baseline.hazard)) {
eval.point <- inla.obj$.arg$data$baseline.hazard.values
basehaz <- inla.obj$summary.random$baseline.hazard$mean + exp(inla.obj$summary.fixed[1,1])
sfun <- stepfun(eval.point[-1], basehaz, f=0)
if (plot) {plot(sfun, xval = eval.point, do.points = F, xaxs = "i", xlim = c(0, max(eval.point)), main="", xlab = "Time", ylab = "Baseline hazard function")}
} else {
eval.point <- inla.obj$.arg$data$baseline.hazard.values
basehaz <- matrix(inla.obj$summary.random$baseline.hazard$mean + exp(inla.obj$summary.fixed[1,1]), nrow = length(eval.point))
maxbh <- max(basehaz)
minbh <- min(basehaz)
if (plot){
if (ncol(basehaz) %in% 1:2) {a <- 1; b <- 2
} else {if (ncol(basehaz) %in% 3:4) {a <- 2; b <- 2}
else {if (ncol(basehaz) %in% 5:6) {a <- 2; b <- 3}
else {if (ncol(basehaz) %in% 7:9) {a <- 3; b <- 3}
else {a <- b <- ceiling(sqrt(ncol(basehaz)))}
}}}
par(mfrow=c(a,b))
for (i in 1:ncol(basehaz)){
sfun <- stepfun(eval.point[-1], basehaz[,i], f=0)
plot(sfun, xval = eval.point, do.points = F, xaxs = "i", xlim = c(0, max(eval.point)), main = paste("Stratified group", i), ylim=c(minbh, maxbh), xlab = "Time", ylab = "Baseline hazard function", cex.lab = cex.lab, cex.axis = cex.axis, ...)
}
}
}
},
weibull = {
org.formula <- inla.obj$.args$formula
org.x <- as.character(org.formula)[2]
splitstr1 <- strsplit(org.x, "(", fixed = TRUE)[[1]]
splitstr2 <- strsplit(splitstr1[2], ",", fixed = TRUE)[[1]]
time <- inla.obj$.args$data[,paste(splitstr2[1])]
alpha <- as.numeric(inla.obj$summary.hyperpar[1])
sigma <- 1/alpha
mu0 <- -1*inla.obj$summary.fixed[1,1]
lambda <- exp(-mu0/sigma)
eval.point <- seq(min(time), max(time), length.out = 101)
h0 <- lambda*alpha*eval.point^(alpha-1)
if (plot){
plot(eval.point, h0, type = "l", xlab = "Time", ylab = "Baseline hazard function", cex.lab = cex.lab, cex.axis = cex.axis, ...)
}
},
exponential = {
org.formula <- inla.obj$.args$formula
org.x <- as.character(org.formula)[2]
splitstr1 <- strsplit(org.x, "(", fixed = TRUE)[[1]]
splitstr2 <- strsplit(splitstr1[2], ",", fixed = TRUE)[[1]]
time <- inla.obj$.args$data[,paste(splitstr2[1])]
mu0 <- -1*inla.obj$summary.fixed[1,1]
lambda <- exp(-mu0)
eval.point <- seq(min(time), max(time), length.out = 101)
h0 <- rep(lambda, length.out = 101)
if (plot){
plot(eval.point, h0, type = "l", xlab = "Time", ylab = "Baseline hazard function", cex.lab = cex.lab, cex.axis = cex.axis, ...)
}
},
stop("The function is only support Weibull, Exponential, and Cox proportional hazards models!!!")
)
return(invisible(list(time = eval.point, basehaz = basehaz)))
}
#' Compute different Bayesian residuals for survival models using INLA
#' @param inla.obj an object of class "inla".
#' @param time the follow up time for right censored data,
#' @param event the status indicator, 1=observed event, 0=right censored event
#'
#' @details The function supports Weibull model, Expontenial model and Cox proportional hazards model for right censored data. Three types of residuals are obtained: Cox-Snell residual; Martingale residuals; Deviance residuals.
#' @return A list of different residuals, and time, event, family information.
#' @export
bri.surv.resid <- function(inla.obj, time, event){
if(class(inla.obj) != "inla") stop("an 'inla' object is needed!")
family <- inla.obj$.args$family
if (!is.numeric(time))
stop("argument 'time' must be numeric")
time <- as.vector(time)
if (!is.numeric(event))
stop("argument 'event' must be numeric")
event <- as.vector(event)
switch(tolower(family),
weibullsurv = {
if (nrow(inla.obj$summary.fitted.values) != length(time))
stop("the 'time' variable does not match with the inla object!")
if (nrow(inla.obj$summary.fitted.values) != length(event))
stop("the 'event' variable does not match with the inla object!")
alpha <- as.numeric(inla.obj$summary.hyperpar[1])
mu0 <- inla.obj$summary.fixed[1,1]
lambda <- exp(mu0*alpha)
theta <- inla.obj$summary.fixed[-1,1]*alpha
# Cox-Snell residual
cs.resid <- as.vector(lambda*exp(inla.obj$model.matrix[,-1] %*% theta)*(time^alpha))
},
exponentialsurv = {
if (nrow(inla.obj$summary.fitted.values) != length(time))
stop("the 'time' variable does not match with the inla object!")
if (nrow(inla.obj$summary.fitted.values) != length(event))
stop("the 'event' variable does not match with the inla object!")
mu0 <- inla.obj$summary.fixed[1,1]
lambda <- exp(mu0)
theta <- inla.obj$summary.fixed[-1,1]
# Cox-Snell residual
cs.resid <- as.vector(lambda*exp(inla.obj$model.matrix[,-1] %*% theta)*time)
},
poisson = {
if (sum(inla.obj$.args$data$baseline.hazard.idx==1) != length(time))
stop("The 'time' variable does not match with the inla object! The function does not support stratified proportional Hazard Model!")
if (sum(inla.obj$.args$data$baseline.hazard.idx==1) != length(event))
stop("The 'event' variable does not match with the inla object! The function does not support stratified proportional Hazard Model!")
h0 <- approxfun(inla.obj$summary.random$baseline.hazard$ID, inla.obj$summary.random$baseline.hazard$mean + exp(inla.obj$summary.fixed[1,1]))
H0.est <- unlist(sapply(time, function(x) integrate(h0, low=0, upper = x))[1,])
idx1 <- which(inla.obj$.arg$data$baseline.hazard.idx == 1)
mm <- inla.obj$model.matrix[idx1,]
beta <- inla.obj$summary.fixed[-1,1]
# Cox-Snell residual
cs.resid <- as.vector(H0.est*exp(mm[,-1] %*% beta))
},
stop("The function is only support Weibull, Exponential, and Cox proportional hazards models!!!")
)
martingale.resid <- event - cs.resid
deviance.resid <- c(sign(martingale.resid)* sqrt(-2* (martingale.resid+ event * log(cs.resid))))
return(list(cs = cs.resid, martingale = martingale.resid, deviance = deviance.resid, time = time, event = event, family = family))
}
#' Cox-Snell residual plot
#'
#' @param resid.obj object from bri.surv.resid function
#' @param lwd the line width, a positive number, defaulting to 2.
#' @param xlab a title for the x axis.
#' @param ylab a title for the y axis.
#' @param main an overall title for the plot.
#' @param cex.lab The magnification to be used for x and y labels relative to the current setting of cex.
#' @param cex.axis The magnification to be used for x and y labels relative to the current setting of cex.
#' @param ...
#'
#' @return
#' @export
bri.csresid.plot <- function(resid.obj, lwd = 2, xlab="Cox-Snell Residual", ylab="Estimated Cumulative Hazard Rates", main = "", cex.lab = 1.25, cex.axis = 1.25, ...){
require(survival)
cs.resid <- resid.obj$cs
event <- resid.obj$event
s.cs.res <- survfit(Surv(cs.resid, event) ~ 1, type="fleming-harrington")
H.est <- cumsum(s.cs.res$n.event/s.cs.res$n.risk)
xy.max <- max(c(s.cs.res$time, H.est))
plot(s.cs.res$time, H.est,type='s',col='black', xlab=xlab, ylab=ylab, lwd=lwd, xlim=c(0,xy.max), ylim=c(0,xy.max), main = main, cex.lab = cex.lab, cex.axis = cex.axis, ...)
abline(0, 1, lty=2)
}
#' Deviance residual plot
#'
#' @param resid.obj object from bri.surv.resid function
#' @param covariate the covariate of interest for the residual plot.
#' @param smooth if smooth = TRUE, a nonparametric smooth curve will be added.
#' @param xlab a title for the x axis.
#' @param ylab a title for the y axis.
#' @param main an overall title for the plot.
#' @param cex.lab The magnification to be used for x and y labels relative to the current setting of cex.
#' @param cex.axis The magnification to be used for x and y labels relative to the current setting of cex.
#' @param ...
#'
#' @return
#' @export
bri.dresid.plot <- function(resid.obj, covariate = NULL, smooth = FALSE, xlab = NULL, ylab = "Deviance residual", main = "", cex.lab = 1.25, cex.axis = 1.25, ...){
d.resid <- resid.obj$deviance
if (is.null(covariate)) {
covariate <- seq(1: length(d.resid))
xlab <- "Index"
} else {
if (length(d.resid) != length(covariate))
stop("the 'covariate' variable does not match with the residual object!")
if (is.null(xlab)) xlab <- "Covariate"
}
if (smooth){
dev.dat <- data.frame(covariate = covariate, deviance = d.resid)
dev.smooth.inla <- inla(deviance ~ -1 + f(covariate, model = 'rw2', constr = FALSE), data = dev.dat)
bri.band.plot(dev.smooth.inla, name = 'covariate', alpha = 0.05, xlab = xlab, ylab = ylab, main = main, type = 'random', ylim= c(min(d.resid), max(d.resid)), cex.lab = cex.lab, cex.axis = cex.axis, ...)
points(dev.dat$covariate, dev.dat$deviance)
} else{
plot(covariate, d.resid, xlab = xlab, main = main, ylab = ylab, ylim= c(min(d.resid), max(d.resid)), cex.lab = cex.lab, cex.axis = cex.axis, ...)
}
}
#' Martingale Residual plot
#'
#' @param resid.obj object from bri.surv.resid function
#' @param covariate the covariate of interest for the residual plot.
#' @param smooth if smooth = TRUE, a nonparametric smooth curve will be added.
#' @param xlab a title for the x axis.
#' @param ylab a title for the y axis.
#' @param main an overall title for the plot.
#' @param cex.lab The magnification to be used for x and y labels relative to the current setting of cex.
#' @param cex.axis The magnification to be used for x and y labels relative to the current setting of cex.
#' @param ...
#'
#' @return
#' @export
bri.mresid.plot <- function(resid.obj, covariate = NULL, smooth = FALSE, xlab = NULL, ylab = "Martingale residual", main = "", cex.lab = 1.25, cex.axis = 1.25, ...){
m.resid <- resid.obj$martingale
if (is.null(covariate)) {
covariate <- seq(1: length(m.resid))
xlab <- "Index"
} else {
if (length(m.resid) != length(covariate))
stop("the 'covariate' variable does not match with the residual object!")
if (is.null(xlab)) xlab <- "Covariate"
}
if (smooth){
mart.dat <- data.frame(covariate = covariate, martingale = m.resid)
mart.smooth.inla <- inla(martingale ~ -1 + f(covariate, model = 'rw2', constr = FALSE), data = mart.dat)
bri.band.plot(mart.smooth.inla, name = 'covariate', alpha = 0.05, xlab = xlab, ylab = ylab, cex.lab = cex.lab, cex.axis = cex.axis, main = main, type = 'random', ylim= c(min(m.resid), max(m.resid)), ...)
points(mart.dat$covariate, mart.dat$martingale)
} else{
plot(covariate, m.resid, xlab = xlab, main = main, ylab = ylab, cex.lab = cex.lab, cex.axis = cex.axis, ylim= c(min(m.resid), max(m.resid)), ...)
}
}
#' Density estimation using INLA
#'
#' @param x the data from which the estimate is to be computed.
#' @param m the number of equally spaced points at which the density is to be estimated.
#' @param from,to the left and right-most points of the grid at which the density is to be estimated.
#' @param cut by default, the values of from and to are cut * diff(range(x)) beyond the extremes of the data.
#' @param diagonal An extra constant added to the diagonal of the precision matrix in INLA.
#' @param constr A boolean variable indicating whater to set a sum to 0 constraint on the term.
#' @param ...
#'
#' @return a list containing the following components: x - the n coordinates of the points where the density is estimated; y - the estimated density values; y.lower - the estimated 2.5 percentile of density; y.lower - the estimated 97.5 percentile of density.
#' @export
bri.density <- function(x, m = 101, from, to, cut = 0.1, diagonal = 1e-03, constr = T, ...){
if (any(is.na(x))) stop("'x' contains missing values!")
if (missing(from)) from <- min(x) - cut * diff(range(x))
if (missing(to)) to <- max(x) + cut * diff(range(x))
Gmatrix <- function(x, sparse = TRUE){
if (any(is.na(x))) stop("'x' contains missing values!")
x <- sort(x)
n <- length(x)
d <- diff(x)
d <- c(Inf, Inf, d, Inf, Inf)
k <- 3:(n + 2)
g <- 2/((d[k - 1]^2) * (d[k - 2] + d[k - 1])) + 2/(d[k - 1]*d[k]) * (1/d[k - 1] + 1/d[k]) + 2/((d[k]^2) * (d[k] + d[k + 1]))
k <- 4:(n + 2)
g1 <- -2/(d[k - 1]^2) * (1/d[k - 2] + 1/d[k])
k <- 5:(n + 2)
g2 <- 2/(d[k - 2] * d[k - 1] * (d[k - 2]+d[k - 1]))
G <- diag(g)
G[row(G) == col(G) + 1] <- g1
G[col(G) == row(G) + 1] <- g1
G[row(G) == col(G) + 2] <- g2
G[col(G) == row(G) + 2] <- g2
if (sparse == TRUE) G <- as(G, "dgTMatrix")
return(G)
}
bins=seq(from, to, length.out = m)
x.bins <- hist(x, breaks = bins, plot=FALSE)
x.bins.root <- sqrt(x.bins$counts+1/4)
idx <- 1:length(x.bins.root)
Q <- Gmatrix(idx)
inla.fit <- inla(x.bins.root ~ f(idx, model = "generic0", Cmatrix = Q,
diagonal = diagonal, constr = constr),
data = as.data.frame(list(x.bins.root = x.bins.root, idx = idx)),
control.predictor = list(compute = T))
inla.est <- inla.fit$summary.linear.predictor[,1]
inla.lower <- inla.fit$summary.linear.predictor[,3]
inla.upper <- inla.fit$summary.linear.predictor[,5]
SimpsonInt <- function (x, f, subdivisions = 256){
ap <- approx(x, f, n = 2 * subdivisions + 1)
integral <- diff(ap$x)[1] * (ap$y[2 * (1:subdivisions) - 1]
+ 4 * ap$y[2 * (1:subdivisions)] + ap$y[2 * (1:subdivisions) + 1])/3
return(sum(integral))
}
normalized <- SimpsonInt(x.bins$mids, inla.est^2)
f <- inla.est^2/normalized
f.lower <- inla.lower^2/normalized
f.upper <- inla.upper^2/normalized
return(structure(list(x = x.bins$mids, y = f, y.lower= f.lower, y.upper=f.upper)))
}
#' Bayesian Regression with Noninformative Priors
#'
#' @param lmfit lm output
#' @param B samples
#'
#' @return posteriors and summary
#' @export
#'
BayesLM.nprior <- function(lmfit, B){
QR <- lmfit$qr
df.res <- lmfit$df.residual
V <- qr.R(QR)
coef <- lmfit$coef
Vb<- chol2inv(V)
s2 <- (t(lmfit$residuals)%*%lmfit$residuals)
s2<- s2[1,1]/df.res
sigma2 <- df.res * s2/rchisq(B,df.res)
coef.post <- data.frame(t(sapply(sigma2,function(x) mvrnorm(1,coef,Vb*x))))
post.dist <- data.frame(coef.post, sigma = sqrt(sigma2))
names(post.dist) <- c(names(lmfit$coef), "sigma")
summary.post <- t(apply(post.dist, 2, function(x)
{
c("mean"=mean(x),
"se"=sd(x),
"0.025quant" = unname(quantile(x,prob=0.025)),
"median"=median(x),
"0.975quant" = unname(quantile(x,prob=0.975)))
}))
return(list(post.dist=post.dist, summary.stat = summary.post))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.