brinla: Bayesian Regression with INLA

library(brinla)

# Splines as a Linear Mixed Model

library(INLA)
data(age.income, package = 'SemiPar')
str(age.income)
library(brinla)
X <- spline.mixed(age.income$age, degree = 2, type = 'TPB')$X
Z <- spline.mixed(age.income$age, degree = 2, type = 'TPB')$Z
age.income$ID <- 1:nrow(age.income)
formula <- log.income ~ -1 + X + f(ID, model = 'z', Z = Z, hyper = list(prec = list(param = c(1e-6, 1e-6))))
result.tpb <- inla(formula, data = age.income, control.predictor = list(compute = TRUE))
result.tpb$summary.fitted.values
bri.band.plot(result.tpb, x = age.income$age, alpha = 0.05, type = 'fitted', xlab = 'Age (years)', ylab = 'Log(annual income)', ylim = range(age.income$log.income))
points(age.income$age, age.income$log.income, cex = 0.5)
result.tpb$summary.random$ID[-(1:nrow(age.income)),]
X <- spline.mixed(age.income$age, type = 'OSS')$X
Z <- spline.mixed(age.income$age, type = 'OSS')$Z
age.income$ID <- 1:nrow(age.income)
formula <- log.income ~ -1 + X + f(ID, model = 'z', Z = Z, hyper = list(prec = list(param = c(1e-6, 1e-6))))
result.oss <- inla(formula, data = age.income, control.predictor = list(compute = TRUE))
B.tmp <- spline.mixed(age.income$age, type = 'OSS')$B
n.row <- nrow(B.tmp)
n.col <- ncol(B.tmp)
attributes(B.tmp) <- NULL
Bmat <- as(matrix(B.tmp, n.row, n.col), 'sparseMatrix')
Q <- as(spline.mixed(age.income$age, type = 'OSS')$Q, 'sparseMatrix')
formula <- y ~ -1 + f(x, model = 'generic0', Cmatrix = Q, hyper = list(prec = list(param = c(1e-6, 1e-6))))
data.inla <- list(y = age.income$log.income, x = 1:ncol(Bmat))
result.oss2 <- inla(formula, data = data.inla, control.predictor = list(A = Bmat, compute = TRUE))
result.oss2$summary.fitted.values[1:n.row, ]
library(brinla)
bri.band.plot(result.oss2, ind = 1:n.row, x = age.income$age, alpha = 0.05, type = 'fitted', xlab = 'Age (years)', ylab = 'Log(annual income)', ylim = range(age.income$log.income))
points(age.income$age, age.income$log.income, cex = 0.5)

# Analysis of Variance for Functional Data

data(DTI, package = 'refund')
str(DTI)
DTI.sub <- DTI[DTI$Nscans==4 & DTI$case == 1, ]
dat.v1 <- DTI.sub[DTI.sub$visit==1,]$rcst
plot(dat.v1[1,], type = "n", ylim = c(0.1, 1))
for(i in 1:dim(dat.v1)[1]) lines(dat.v1[i,])
y <- as.vector(t(DTI.sub$rcst))
ns <- dim(DTI.sub$rcst)[2]	  ## number of locations
ng <- length(unique(DTI.sub$visit))  ## number of visits (groups)
n <- length(unique(DTI.sub$ID))  ## number of subjects

D1 <- Matrix(rep(1,ng*n),ng*n,1)
A.mu <- kronecker(D1, Diagonal(n=ns, x=1))

D1 <- Diagonal(n = ns, x = 1)
D2 <- Diagonal(n = (ng-1), x = 1)
D3 <- Matrix(rep(0, ng-1), 1, ng-1)
D4 <- kronecker(rBind(D2, D3), D1)
A.a <- kronecker(Matrix(rep(1, n), n, 1), D4)

A <- cBind(A.mu, A.a)
mu <- 1:ns
alpha <- rep(1:ns, ng-1)
alpha.rep <- rep(1:(ng-1), each = ns)
mu2 <- c(mu, rep(NA, length(alpha)))
alpha2 <- c(rep(NA, length(mu)), alpha)
alpha2.rep <- c(rep(NA, length(mu)), alpha.rep)
data.inla <- list(y=y, mu=mu2, alpha=alpha2, alpha.rep=alpha2.rep)
formula <- y ~ -1 + f(mu, model = 'rw2', constr = FALSE, scale.model = T) + f(alpha, model = 'rw2', constr = FALSE, scale.model = TRUE, replicate = alpha.rep)
A1.lc <- kronecker(Matrix(rep(1,ng-1),ng-1,1), Diagonal(n=ns, x=1))
A2.lc <- Diagonal(n = (ng - 1)*ns, x = 1)
lc <- inla.make.lincombs(mu = A1.lc, alpha = A2.lc)
result <- inla(formula, data = data.inla, family = 'beta', control.predictor = list(A = A, compute = TRUE), control.compute = list(config = TRUE), lincomb = lc)
bri.band.ggplot(result, ind = 1:ns, type = 'lincomb') #1st visit
bri.band.ggplot(result, ind = 1:ns + ns, type = 'lincomb') #2nd visit
bri.band.ggplot(result, ind = 1:ns + 2*ns, type='lincomb') #3rd visit
bri.band.ggplot(result, name = 'mu', type = 'random') #4th visit
bri.band.ggplot(result, name='alpha', ind=1:ns, type='random')
bri.band.ggplot(result, name='alpha', ind=1:ns+ns, type='random')
bri.band.ggplot(result, name='alpha', ind=1:ns+2*ns, type='random')
res.exc1 <- excursions.brinla(result, name = 'alpha', ind = 1:ns, u = 0, type = '!=', alpha = 0.05, method = 'NIQC')
res.exc2 <- excursions.brinla(result, name = 'alpha', ind = 1:ns + ns, u = 0, type = '!=', alpha = 0.05, method = 'NIQC')
res.exc3 <- excursions.brinla(result, name = 'alpha', ind = 1:ns + 2*ns, u = 0, type = '!=', alpha = 0.05, method = 'NIQC')
bri.excursions.ggplot(res.exc1)
bri.excursions.ggplot(res.exc2)
bri.excursions.ggplot(res.exc3)

# Extreme Values

data(calder, package="brinla")
plot(Flow ~ WaterYear, calder)
library(INLA)
calder$year <- calder$WaterYear-1973
imod <- inla(Flow ~ 1+year, data=calder, family="gev",scale=0.1)
imod$summary.fixed
inla.pmarginal(0, imod$marginals.fixed$year)
library(brinla)
plot(bri.hyper.sd(imod$marginals.hyperpar$`precision for GEV observations`), type="l", xlab="SD", ylab="density")
plot(imod$marginals.hyperpar$`tail parameter for GEV observations`,type="l", xlim=c(-0.2,0.5), xlab="xi", ylab="density")
pgev <- function(y,xi,tau,eta,sigma=1){
  exp(-(1+xi*sqrt(tau*sigma)*(y-eta))^(-1/xi))
}
yr <- 1997-1973
maxflow <- 173.17
eta <- sum(c(1,yr)*imod$summary.fixed$mean)
tau <- imod$summary.hyperpar$mean[1]
xi <- imod$summary.hyperpar$mean[2]
sigma <- 0.1
pless <- pgev(maxflow, xi, tau, eta,sigma)
1-pless
1/(1-pless)
nsamp <- 999
imod <- inla(Flow ~ 1+year, data=calder, family="gev",scale=0.1, control.compute = list(config=TRUE))
postsamp <- inla.posterior.sample(nsamp, imod)
pps <- t(sapply(postsamp, function(x) c(x$hyperpar, x$latent[42:43])))
colnames(pps) <- c("precision","shape","beta0","beta1")
plot(shape*0.01 ~ precision, pps, ylab="shape")
imod <- inla(Flow ~ 1+year, data=calder, family="gev",scale=0.1, control.compute = list(config=TRUE),control.inla=list(int.strategy='grid', dz=0.2))
postsamp <- inla.posterior.sample(nsamp, imod)
pps <- t(sapply(postsamp, function(x) c(x$hyperpar, x$latent[42:43])))
colnames(pps) <- c("precision","shape","beta0","beta1")
plot(shape*0.01 ~ precision, pps, ylab="shape")
sigma <- 0.1
maxflow <- 173.17
retp = numeric(nsamp)
for(i in 1:nsamp){
  eta <- sum(c(1,yr)*pps[i,3:4])
  tau <- pps[i,1]
  xi <- 0.01*pps[i,2]
  pless <-  pgev(maxflow, xi, tau, eta,sigma)
  retp[i] <- 1/(1-pless)
}
quantile(retp, c(0.025, 0.5, 0.975))

# Density Estimation using INLA

library(brinla)
## Example 1: Density estimation Case 1
set.seed(123)
n <- 500
x <- rnorm(n)
x.den1 <- bri.density(x, cut = 0.3)
x.den2 <- density(x, bw = "SJ")

curve(dnorm(x, mean = 0, sd = 1), from = -4, to = 4, lwd = 3, lty = 3, col = "red", xlab = "x", ylab = "f(x)", cex.lab = 1.5, cex.axis = 1.5, ylim=c(0, 0.45))
lines(x.den1, col = "black", lty = 1, lwd = 3, ylim = c(0, 0.25))
lines(x.den1$x, x.den1$y.upper, lty = 2, lwd = 3)
lines(x.den1$x, x.den1$y.lower, lty = 2, lwd = 3)
lines(x.den2, col = "blue", lty = 4, lwd = 3)

## Example 2: Density estimation Case 2
set.seed(123)
n <- 1000
x <- c(rnorm(n/2, mean = -1.5, sd = 1), rnorm(n/2, mean = 2.5, sd = 0.75))
x.den1 <- bri.density(x, cut = 0.3)
x.den2 <- density(x, bw = "SJ")

curve(dnorm(x, mean = -1.5, sd = 1)/2 + dnorm(x, mean = 2.5, sd = 0.75)/2, from = -6, to = 6, lwd = 3, lty = 3, col = "red", xlab = "x", ylab = "f(x)", ylim=c(0, 0.3), cex.lab = 1.5, cex.axis = 1.5)
lines(x.den1, col = "black", lty = 1, lwd = 3)
lines(x.den1$x, x.den1$y.upper, lty = 2, lwd = 3)
lines(x.den1$x, x.den1$y.lower, lty = 2, lwd = 3)
lines(x.den2, col = "blue", lty = 4, lwd = 3)


sessionInfo()