This repository includes code for BTST (Bayesian Tensor Skew-t) model with tensor spike-and-slab lasso/tensor normal prior described in the paper "A New Class of Tensor Skew Distributions", by Lee et al (2023+).
The BTST package can be installed with the devtools package:
library(devtools)
devtools::install_github(repo='InkooLee/BTST')
Vignette is available at /InkooLee/BTST/blob/master/vignettes/BTST_vignette.html
The BTST package provides posterior sampling algorithm for tensor skew-t model.
## Load library and required package
library(BTST)
library(abind)
## load GAAD data
data <- read.csv("PDdata.csv")
## load required functions
source("./functions_tensor.R")
source("./BTST_TSSL.R")
source("./BTST_TN.R")
#'@ Inputs
#' Y \in R^{t x s x b x n} three way tensor of responses with n subjects
#' X \in R^{p} vector of subject-level covariates (with intercept, (p+1) x 1)
t = 28 # maximum number of teeth per one subject
s = 6 # maximum number of surfaces per each tooth
b = 2 # two bio-markers (i.e., PPD and CAL)
n = 290 # number of subjects
# Three way tensor response
YP <- array(data$PPD, dim = c(s,t,n)) # PPD
YC <- array(data$CAL, dim = c(s,t,n)) # CAL
YP <- aperm(YP, c(2,1,3)) # t x s matrices for PPD
YC <- aperm(YC, c(2,1,3)) # t x s matrices for CAL
Y <- abind(YP, YC) # t x s x nl
Y <- array(Y, dim = c(t,s,b,n)) # t x s x b x n
#'@PPD/CAL
#'Y[,,,1] # n subjects for PPD
#'Y[,,,2] # n subjects for CAL
## covariates
X <- data[,7:11]
X <- X[seq(1, nrow(X), 168), ]
rownames(X) <- NULL
X <- t(as.matrix(X)) # matrix of covariates (Age,Gender,BMI,Smoker,HbA1c) x number of subjects
X <- rbind(t(matrix(1,n)),X) # including intercept term
p <- dim(X)[1]
# Missing Responses
vecy <- t(mat(Y,4))
vecym <- matrix(NA, t*s*b, n)
vecym[which(!is.na(vecy))] <- 0; vecym[which(is.na(vecy))] <- 1; delta_p <- vecym; # delta_p is t*s*b by n matrix with missing = 1, observed = 0
We can fit the model as
## Fit the model
result_TSSL <- BTST_TSSL(Y,X,vecy, n.burn = 100, n.save = 1000, thin = 5)
result_TN <- BTST_TN(Y,X,vecy, n.burn = 100, n.save = 1000, thin = 5)