In InkooLee/BSTN: Bayesian Regression Analayis of Skewed Tensor Responses
knitr::opts_chunk$set(echo = TRUE)
library(kableExtra)
options(knitr.table.format = 'markdown')
Inkoo Lee
Last Updated - August 10, 2022
Paper
- Lee, I., Sinha, D., Mai, Q., Zhang, X. & Bandyopadhyay, D. (2022+). Bayesian Regression Analysis of Skewed Tensor Responses. (Accepted), Biometrics
Skewed tensor response model
For an order-$K$ tensor response $\mathcal{Y}i \in \mathbb{R}^{d_1 \times \cdots \times d_K}$ with a vector of covariates $\textbf{x}_i \in \mathbb{R}^p$, we consider a tensor response regression model
\begin{equation}\label{eq1}
\mathcal{Y}{i} = \mathcal{B}\bar{\times}{(K+1)} \textbf{x}{i} + \mathcal{E}{i},~~for~~ i=1,\ldots,n,~
\end{equation}
where $\mathcal{B} \in \mathbb{R}^{d_1 \times \cdots \times d_K \times p}$ is an order-$(K+1)$ unknown tensor of regression coefficients, $\bar{\times}{(K+1)}$ is the $(K+1)$-mode vector product, and the error $\mathcal{E}{i} \in \mathbb{R}^{d_1 \times \cdots \times d_K}$ is an order-$K$ tensor. We model the skewness in
the distribution of $\mathcal{Y}$ via
\begin{equation}\label{eq2}
\mathcal{E}{i} = |\mathcal{Z}{2i}|\times{K} \pmb{\Lambda} + \mathcal{Z}{1i},
\end{equation}
where $\pmb{\Lambda} = diag(\lambda_1,\ldots,\lambda{d_K}) \in \mathbb{R}^{d_K \times d_K}$ with skewness parameters $\pmb{\lambda}=(\lambda_1,\ldots,\lambda_{d_K})$ and $|\mathbf{M}|$ denotes the
matrix whose each element of is the absolute value of the corresponding element of matrix $\mathbf{M}$.
Tensor Spike-and-Slab prior
\begin{equation}\label{eq:IdentifyingCoef}
\begin{aligned}
&\mathcal{B}{i_1 i_2 i_3 j} = \eta{i_3j} + \gamma_{i_1i_2i_3 j},~~ \eta_{i_3j} \sim N(0, \tau^2)~~ for~~ \tau > 0
\end{aligned}
\end{equation}
and introduce sparsity via spike-and-slab priors \citep{ishwaran2005spike} as
\begin{equation}\label{eq:SparsePrior}
\gamma_{i_1i_2i_3j} \overset{iid}{\sim} (1 - \omega_{i_1
i_2 i_3 j})\delta_0(\gamma_{i_1i_2i_3 j}) + \omega_{i_1i_2i_3 j}N(0, \nu^2)~~ for~~ \nu > 0,
\end{equation}
where, $\delta_0(\cdot)$ is the indicator function at 0 and $\omega_{i_1i_2i_3 j} \overset{iid}{\sim} Bernoulli(\psi)$ with hyperprior $\psi \sim Beta(a,b)$. For the GAAD study, we choose $a = b = 0.1$, representing a vague hyperprior for $\psi$.
Since sparsity inducing prior on $\gamma_{i_1i_2i_3j}$ tends to concentrate around zero, $\mathcal{B}$ in \eqref{eq:IdentifyingCoef} is mostly determined by $\eta_{i_3j}$ allowing parsimonious interpretation of covariate effects, while estimates of $\mathcal{B}$ using Tensor Normal prior (with large-scale covariances) are much more different from each other.
GAAD data analysis using BSTN package
## Load package
library(BSTN)
We fit the BSTN (Bayesian skewed tensor normal) model with tensor SAS (spike-and-slab) prior.
# n.burn = burn-in period, n.save = number of posterior samples, thin = thinning size
BSTN_SAS <- BSTN_SAS(Y,X,vecy, n.burn = 100, n.save = 1000, thin = 5)
D-Statistic
$D = d_1/|d_2|$ for
each regression effects $\mathcal{B}{i_1i_2i_3j}$, where $d_1$ is distance of the Bayesian estimate (posterior median) from null value (0 here) and $d_2$ is the distance of the Bayes estimate from the nearest endpoint of the 95\% CI of the parameter.
Here, $D > 1.5$ $(D < -1.5)$ implies strong posterior evidence of positive (negative) effect of covariate $x_j$ on the biomarker $i_3$ at site $i_2$ of tooth $i_1$.
Similarly, the range $-0.8 < D < 0.8$ suggests moderate evidence
for null hypothesis, and $0.8 \leq D \leq 1.5$ represents moderate positive effect, and $-1.5 \leq D \leq -0.8$ denotes moderate negative effect of $x_j$ on $\mathcal{Y}{i_1i_2i_3}$.
The following plot (Figure 1 in the main text) shows D-statistic map of tensor regression coefficients with sparsity inducing prior.
load("C:/Users/inkku/OneDrive/Desktop/BSTN/GAAD_result/PPD.rda")
load("C:/Users/inkku/OneDrive/Desktop/BSTN/GAAD_result/CAL.rda")
mat<-function(A,j)
{
Aj<-t(apply(A,j,"c"))
if(nrow(Aj)!=dim(A)[j]) { Aj<-t(Aj) }
Aj
}
PPD_B.est <- mat(PPD3$B.est,1)[2:6,] # exclude intercept
CAL_B.est <- mat(CAL$B.est,1)[2:6,] # exclude intercept
PPD2init <- as.matrix(PPD_B.est[1,])
PPD3init <- as.matrix(PPD_B.est[2,])
PPD4init <- as.matrix(PPD_B.est[3,])
PPD5init <- as.matrix(PPD_B.est[4,])
PPD6init <- as.matrix(PPD_B.est[5,])
CAL2init <- as.matrix(CAL_B.est[1,])
CAL3init <- as.matrix(CAL_B.est[2,])
CAL4init <- as.matrix(CAL_B.est[3,])
CAL5init <- as.matrix(CAL_B.est[4,])
CAL6init <- as.matrix(CAL_B.est[5,])
library(latticeExtra)
library(RColorBrewer)
library(lattice)
library(rasterVis)
library(circlize)
t=28;s=6;
#Acol.regions <- c("#FF80FFFF", "#FFFFFF", "#1E90FF")
Acol.regions <- c("#FF80FFFF","#FFB6C1", "#FFFFFF", "#BFFFFFFF", "#1E90FF")
n.iter = 1000
## PPD AGE ##
PPD2 <- matrix(NA,168,n.iter)
for (i in 1:n.iter){
PPD2[,i] <- PPD2init[(168*(i-1) + 1):(168*i),]
}
CIPPD2 <- matrix(NA,168,3)
for (i in 1:168){
CIPPD2[i,] <- quantile(PPD2[i,], c(0.025,0.5,0.975))
}
PPD2_d1 <- as.matrix(CIPPD2[,2]); PPD2_d2 <- matrix(NA, 168, 1)
for (i in 1:168){
if (PPD2_d1[i,] > 0){
PPD2_d2[i,] <- PPD2_d1[i,] - CIPPD2[i,1]
}
else {PPD2_d2[i,] <- CIPPD2[i,3] - PPD2_d1[i,]}
}
PPD2map <- PPD2_d1/abs(PPD2_d2)
PPD2map <- matrix(PPD2map, t,s)
PPD2map <- levelplot(t(PPD2map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4),
colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)),
margin=F, scales = list(draw = FALSE),main = list(label = "PPD_Age", cex = 1))
print(PPD2map, split = c(1, 1, 5, 2), more = TRUE)
## PPD Gender ##
PPD3 <- matrix(NA,168,n.iter)
for (i in 1:n.iter){
PPD3[,i] <- PPD3init[(168*(i-1) + 1):(168*i),]
}
CIPPD3 <- matrix(NA,168,3)
for (i in 1:168){
CIPPD3[i,] <- quantile(PPD3[i,], c(0.025,0.5,0.975))
}
PPD3_d1 <- as.matrix(CIPPD3[,2]); PPD3_d2 <- matrix(NA, 168, 1)
for (i in 1:168){
if (PPD3_d1[i,] > 0){
PPD3_d2[i,] <- PPD3_d1[i,] - CIPPD3[i,1]
}
else {PPD3_d2[i,] <- CIPPD3[i,3] - PPD3_d1[i,]}
}
PPD3map <- PPD3_d1/abs(PPD3_d2)
PPD3map <- matrix(PPD3map, t,s)
PPD3map <- levelplot(t(PPD3map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4),
colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)),
margin=F, scales = list(draw = FALSE),main = list(label = "PPD_Gender", cex = 1))
print(PPD3map, split = c(2, 1, 5, 2), more = TRUE)
## PPD BMI ##
PPD4 <- matrix(NA,168,n.iter)
for (i in 1:n.iter){
PPD4[,i] <- PPD4init[(168*(i-1) + 1):(168*i),]
}
CIPPD4 <- matrix(NA,168,3)
for (i in 1:168){
CIPPD4[i,] <- quantile(PPD4[i,], c(0.025,0.5,0.975))
}
PPD4_d1 <- as.matrix(CIPPD4[,2]); PPD4_d2 <- matrix(NA, 168, 1)
for (i in 1:168){
if (PPD4_d1[i,] > 0){
PPD4_d2[i,] <- PPD4_d1[i,] - CIPPD4[i,1]
}
else {PPD4_d2[i,] <- CIPPD4[i,3] - PPD4_d1[i,]}
}
PPD4map <- PPD4_d1/abs(PPD4_d2)
PPD4map <- matrix(PPD4map, t,s)
PPD4map <- levelplot(t(PPD4map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4),
colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)),
margin=F, scales = list(draw = FALSE),main = list(label = "PPD_BMI", cex = 1))
print(PPD4map, split = c(3, 1, 5, 2), more = TRUE)
## PPD Smoker ##
PPD5 <- matrix(NA,168,n.iter)
for (i in 1:n.iter){
PPD5[,i] <- PPD5init[(168*(i-1) + 1):(168*i),]
}
CIPPD5 <- matrix(NA,168,3)
for (i in 1:168){
CIPPD5[i,] <- quantile(PPD5[i,], c(0.025,0.5,0.975))
}
PPD5_d1 <- as.matrix(CIPPD5[,2]); PPD5_d2 <- matrix(NA, 168, 1)
for (i in 1:168){
if (PPD5_d1[i,] > 0){
PPD5_d2[i,] <- PPD5_d1[i,] - CIPPD5[i,1]
}
else {PPD5_d2[i,] <- CIPPD5[i,3] - PPD5_d1[i,]}
}
PPD5map <- PPD5_d1/abs(PPD5_d2)
PPD5map <- matrix(PPD5map, t,s)
PPD5map <- levelplot(t(PPD5map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4),
colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)),
margin=F, scales = list(draw = FALSE),main = list(label = "PPD_Smoker", cex = 1))
print(PPD5map, split = c(4, 1, 5, 2), more = TRUE)
## PPD HbA1c ##
PPD6 <- matrix(NA,168,n.iter)
for (i in 1:n.iter){
PPD6[,i] <- PPD6init[(168*(i-1) + 1):(168*i),]
}
CIPPD6 <- matrix(NA,168,3)
for (i in 1:168){
CIPPD6[i,] <- quantile(PPD6[i,], c(0.025,0.5,0.975))
}
PPD6_d1 <- as.matrix(CIPPD6[,2]); PPD6_d2 <- matrix(NA, 168, 1)
for (i in 1:168){
if (PPD6_d1[i,] > 0){
PPD6_d2[i,] <- PPD6_d1[i,] - CIPPD6[i,1]
}
else {PPD6_d2[i,] <- CIPPD6[i,3] - PPD6_d1[i,]}
}
PPD6map <- PPD6_d1/abs(PPD6_d2)
PPD6map <- matrix(PPD6map, t,s)
PPD6map <- levelplot(t(PPD6map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4),
colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)),
margin=F, scales = list(draw = FALSE),main = list(label = "PPD_HbA1c", cex = 1))
print(PPD6map, split = c(5, 1, 5, 2), more = TRUE)
## CAL AGE ##
CAL2 <- matrix(NA,168,n.iter)
for (i in 1:n.iter){
CAL2[,i] <- CAL2init[(168*(i-1) + 1):(168*i),]
}
CICAL2 <- matrix(NA,168,3)
for (i in 1:168){
CICAL2[i,] <- quantile(CAL2[i,], c(0.025,0.5,0.975))
}
CAL2_d1 <- as.matrix(CICAL2[,2]); CAL2_d2 <- matrix(NA, 168, 1)
for (i in 1:168){
if (CAL2_d1[i,] > 0){
CAL2_d2[i,] <- CAL2_d1[i,] - CICAL2[i,1]
}
else {CAL2_d2[i,] <- CICAL2[i,3] - CAL2_d1[i,]}
}
CAL2map <- CAL2_d1/abs(CAL2_d2)
CAL2map <- matrix(CAL2map, t,s)
CAL2map <- levelplot(t(CAL2map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4),
colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)),
margin=F, scales = list(draw = FALSE),main = list(label = "CAL_Age", cex = 1))
print(CAL2map, split = c(1, 2, 5, 2), more = TRUE)
## CAL Gender ##
CAL3 <- matrix(NA,168,n.iter)
for (i in 1:n.iter){
CAL3[,i] <- CAL3init[(168*(i-1) + 1):(168*i),]
}
CICAL3 <- matrix(NA,168,3)
for (i in 1:168){
CICAL3[i,] <- quantile(CAL3[i,], c(0.025,0.5,0.975))
}
CAL3_d1 <- as.matrix(CICAL3[,2]); CAL3_d2 <- matrix(NA, 168, 1)
for (i in 1:168){
if (CAL3_d1[i,] > 0){
CAL3_d2[i,] <- CAL3_d1[i,] - CICAL3[i,1]
}
else {CAL3_d2[i,] <- CICAL3[i,3] - CAL3_d1[i,]}
}
CAL3map <- CAL3_d1/abs(CAL3_d2)
CAL3map <- matrix(CAL3map, t,s)
CAL3map <- levelplot(t(CAL3map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4),
colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)),
margin=F, scales = list(draw = FALSE),main = list(label = "CAL_Gender", cex = 1))
print(CAL3map, split = c(2, 2, 5, 2), more = TRUE)
## CAL BMI ##
CAL4 <- matrix(NA,168,n.iter)
for (i in 1:n.iter){
CAL4[,i] <- CAL4init[(168*(i-1) + 1):(168*i),]
}
CICAL4 <- matrix(NA,168,3)
for (i in 1:168){
CICAL4[i,] <- quantile(CAL4[i,], c(0.025,0.5,0.975))
}
CAL4_d1 <- as.matrix(CICAL4[,2]); CAL4_d2 <- matrix(NA, 168, 1)
for (i in 1:168){
if (CAL4_d1[i,] > 0){
CAL4_d2[i,] <- CAL4_d1[i,] - CICAL4[i,1]
}
else {CAL4_d2[i,] <- CICAL4[i,3] - CAL4_d1[i,]}
}
CAL4map <- CAL4_d1/abs(CAL4_d2)
CAL4map <- matrix(CAL4map, t,s)
CAL4map <- levelplot(t(CAL4map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4),
colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)),
margin=F, scales = list(draw = FALSE),main = list(label = "CAL_BMI", cex = 1))
print(CAL4map, split = c(3, 2, 5, 2), more = TRUE)
## CAL Smoker ##
CAL5 <- matrix(NA,168,n.iter)
for (i in 1:n.iter){
CAL5[,i] <- CAL5init[(168*(i-1) + 1):(168*i),]
}
CICAL5 <- matrix(NA,168,3)
for (i in 1:168){
CICAL5[i,] <- quantile(CAL5[i,], c(0.025,0.5,0.975))
}
CAL5_d1 <- as.matrix(CICAL5[,2]); CAL5_d2 <- matrix(NA, 168, 1)
for (i in 1:168){
if (CAL5_d1[i,] > 0){
CAL5_d2[i,] <- CAL5_d1[i,] - CICAL5[i,1]
}
else {CAL5_d2[i,] <- CICAL5[i,3] - CAL5_d1[i,]}
}
CAL5map <- CAL5_d1/abs(CAL5_d2)
CAL5map <- matrix(CAL5map, t,s)
CAL5map <- levelplot(t(CAL5map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4),
colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)),
margin=F, scales = list(draw = FALSE),main = list(label = "CAL_Smoker", cex = 1))
print(CAL5map, split = c(4, 2, 5, 2), more = TRUE)
## CAL HbA1c ##
CAL6 <- matrix(NA,168,n.iter)
for (i in 1:n.iter){
CAL6[,i] <- CAL6init[(168*(i-1) + 1):(168*i),]
}
CICAL6 <- matrix(NA,168,3)
for (i in 1:168){
CICAL6[i,] <- quantile(CAL6[i,], c(0.025,0.5,0.975))
}
CAL6_d1 <- as.matrix(CICAL5[,2]); CAL6_d2 <- matrix(NA, 168, 1)
for (i in 1:168){
if (CAL6_d1[i,] > 0){
CAL6_d2[i,] <- CAL6_d1[i,] - CICAL6[i,1]
}
else {CAL6_d2[i,] <- CICAL6[i,3] - CAL6_d1[i,]}
}
CAL6map <- CAL6_d1/abs(CAL6_d2)
CAL6map <- matrix(CAL6map, t,s)
CAL6map <- levelplot(t(CAL6map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4),
colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)),
margin=F, scales = list(draw = FALSE),main = list(label = "CAL_HbA1c", cex = 1))
print(CAL6map, split = c(5, 2, 5, 2), more = TRUE)
Now, the following plots (Figure A.3 in the Supplementary materials) represent residuals (real values - fitted values) for PPD & CAL of randomly picked patient.
setwd('C:/Users/inkku/Downloads')
load("resd_PPD.rda"); load("resd_CAL.rda")
library(lattice);library(viridisLite)
t = 28; s = 6; b = 2; n = 290;
PPD_resd <- array(resd_PPD[,38], dim = c(t,s))
CAL_resd <- array(resd_CAL[,38], dim = c(t,s))
col <- viridis(100)
PPD <- array(t(PPD_resd), dim = c(6,28,1))
dimnames(PPD) <- list(c(1:6),c(1:28),c("PPD residual"))
levelplot(PPD,par.strip.text=list(cex=0.55), xlab="Sites",ylab="Teeth", scales = list(draw = FALSE), col.regions = col)
CAL <- array(t(CAL_resd), dim = c(6,28,1))
dimnames(CAL) <- list(c(1:6),c(1:28),c("CAL residual"))
levelplot(CAL,par.strip.text=list(cex=0.55), xlab="Sites",ylab="Teeth", scales = list(draw = FALSE), col.regions = col)
InkooLee/BSTN documentation built on Aug. 13, 2022, 5:42 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE) library(kableExtra) options(knitr.table.format = 'markdown')
Inkoo Lee
Last Updated - August 10, 2022
Paper
- Lee, I., Sinha, D., Mai, Q., Zhang, X. & Bandyopadhyay, D. (2022+). Bayesian Regression Analysis of Skewed Tensor Responses. (Accepted), Biometrics
Skewed tensor response model
For an order-$K$ tensor response $\mathcal{Y}i \in \mathbb{R}^{d_1 \times \cdots \times d_K}$ with a vector of covariates $\textbf{x}_i \in \mathbb{R}^p$, we consider a tensor response regression model \begin{equation}\label{eq1} \mathcal{Y}{i} = \mathcal{B}\bar{\times}{(K+1)} \textbf{x}{i} + \mathcal{E}{i},~~for~~ i=1,\ldots,n,~ \end{equation} where $\mathcal{B} \in \mathbb{R}^{d_1 \times \cdots \times d_K \times p}$ is an order-$(K+1)$ unknown tensor of regression coefficients, $\bar{\times}{(K+1)}$ is the $(K+1)$-mode vector product, and the error $\mathcal{E}{i} \in \mathbb{R}^{d_1 \times \cdots \times d_K}$ is an order-$K$ tensor. We model the skewness in the distribution of $\mathcal{Y}$ via \begin{equation}\label{eq2} \mathcal{E}{i} = |\mathcal{Z}{2i}|\times{K} \pmb{\Lambda} + \mathcal{Z}{1i}, \end{equation} where $\pmb{\Lambda} = diag(\lambda_1,\ldots,\lambda{d_K}) \in \mathbb{R}^{d_K \times d_K}$ with skewness parameters $\pmb{\lambda}=(\lambda_1,\ldots,\lambda_{d_K})$ and $|\mathbf{M}|$ denotes the matrix whose each element of is the absolute value of the corresponding element of matrix $\mathbf{M}$.
Tensor Spike-and-Slab prior
\begin{equation}\label{eq:IdentifyingCoef} \begin{aligned} &\mathcal{B}{i_1 i_2 i_3 j} = \eta{i_3j} + \gamma_{i_1i_2i_3 j},~~ \eta_{i_3j} \sim N(0, \tau^2)~~ for~~ \tau > 0 \end{aligned} \end{equation} and introduce sparsity via spike-and-slab priors \citep{ishwaran2005spike} as \begin{equation}\label{eq:SparsePrior} \gamma_{i_1i_2i_3j} \overset{iid}{\sim} (1 - \omega_{i_1 i_2 i_3 j})\delta_0(\gamma_{i_1i_2i_3 j}) + \omega_{i_1i_2i_3 j}N(0, \nu^2)~~ for~~ \nu > 0, \end{equation} where, $\delta_0(\cdot)$ is the indicator function at 0 and $\omega_{i_1i_2i_3 j} \overset{iid}{\sim} Bernoulli(\psi)$ with hyperprior $\psi \sim Beta(a,b)$. For the GAAD study, we choose $a = b = 0.1$, representing a vague hyperprior for $\psi$.
Since sparsity inducing prior on $\gamma_{i_1i_2i_3j}$ tends to concentrate around zero, $\mathcal{B}$ in \eqref{eq:IdentifyingCoef} is mostly determined by $\eta_{i_3j}$ allowing parsimonious interpretation of covariate effects, while estimates of $\mathcal{B}$ using Tensor Normal prior (with large-scale covariances) are much more different from each other.
GAAD data analysis using BSTN package
## Load package library(BSTN)
We fit the BSTN (Bayesian skewed tensor normal) model with tensor SAS (spike-and-slab) prior.
# n.burn = burn-in period, n.save = number of posterior samples, thin = thinning size BSTN_SAS <- BSTN_SAS(Y,X,vecy, n.burn = 100, n.save = 1000, thin = 5)
D-Statistic
$D = d_1/|d_2|$ for each regression effects $\mathcal{B}{i_1i_2i_3j}$, where $d_1$ is distance of the Bayesian estimate (posterior median) from null value (0 here) and $d_2$ is the distance of the Bayes estimate from the nearest endpoint of the 95\% CI of the parameter. Here, $D > 1.5$ $(D < -1.5)$ implies strong posterior evidence of positive (negative) effect of covariate $x_j$ on the biomarker $i_3$ at site $i_2$ of tooth $i_1$. Similarly, the range $-0.8 < D < 0.8$ suggests moderate evidence for null hypothesis, and $0.8 \leq D \leq 1.5$ represents moderate positive effect, and $-1.5 \leq D \leq -0.8$ denotes moderate negative effect of $x_j$ on $\mathcal{Y}{i_1i_2i_3}$.
The following plot (Figure 1 in the main text) shows D-statistic map of tensor regression coefficients with sparsity inducing prior.
load("C:/Users/inkku/OneDrive/Desktop/BSTN/GAAD_result/PPD.rda") load("C:/Users/inkku/OneDrive/Desktop/BSTN/GAAD_result/CAL.rda") mat<-function(A,j) { Aj<-t(apply(A,j,"c")) if(nrow(Aj)!=dim(A)[j]) { Aj<-t(Aj) } Aj } PPD_B.est <- mat(PPD3$B.est,1)[2:6,] # exclude intercept CAL_B.est <- mat(CAL$B.est,1)[2:6,] # exclude intercept PPD2init <- as.matrix(PPD_B.est[1,]) PPD3init <- as.matrix(PPD_B.est[2,]) PPD4init <- as.matrix(PPD_B.est[3,]) PPD5init <- as.matrix(PPD_B.est[4,]) PPD6init <- as.matrix(PPD_B.est[5,]) CAL2init <- as.matrix(CAL_B.est[1,]) CAL3init <- as.matrix(CAL_B.est[2,]) CAL4init <- as.matrix(CAL_B.est[3,]) CAL5init <- as.matrix(CAL_B.est[4,]) CAL6init <- as.matrix(CAL_B.est[5,]) library(latticeExtra) library(RColorBrewer) library(lattice) library(rasterVis) library(circlize) t=28;s=6; #Acol.regions <- c("#FF80FFFF", "#FFFFFF", "#1E90FF") Acol.regions <- c("#FF80FFFF","#FFB6C1", "#FFFFFF", "#BFFFFFFF", "#1E90FF") n.iter = 1000 ## PPD AGE ## PPD2 <- matrix(NA,168,n.iter) for (i in 1:n.iter){ PPD2[,i] <- PPD2init[(168*(i-1) + 1):(168*i),] } CIPPD2 <- matrix(NA,168,3) for (i in 1:168){ CIPPD2[i,] <- quantile(PPD2[i,], c(0.025,0.5,0.975)) } PPD2_d1 <- as.matrix(CIPPD2[,2]); PPD2_d2 <- matrix(NA, 168, 1) for (i in 1:168){ if (PPD2_d1[i,] > 0){ PPD2_d2[i,] <- PPD2_d1[i,] - CIPPD2[i,1] } else {PPD2_d2[i,] <- CIPPD2[i,3] - PPD2_d1[i,]} } PPD2map <- PPD2_d1/abs(PPD2_d2) PPD2map <- matrix(PPD2map, t,s) PPD2map <- levelplot(t(PPD2map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4), colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)), margin=F, scales = list(draw = FALSE),main = list(label = "PPD_Age", cex = 1)) print(PPD2map, split = c(1, 1, 5, 2), more = TRUE) ## PPD Gender ## PPD3 <- matrix(NA,168,n.iter) for (i in 1:n.iter){ PPD3[,i] <- PPD3init[(168*(i-1) + 1):(168*i),] } CIPPD3 <- matrix(NA,168,3) for (i in 1:168){ CIPPD3[i,] <- quantile(PPD3[i,], c(0.025,0.5,0.975)) } PPD3_d1 <- as.matrix(CIPPD3[,2]); PPD3_d2 <- matrix(NA, 168, 1) for (i in 1:168){ if (PPD3_d1[i,] > 0){ PPD3_d2[i,] <- PPD3_d1[i,] - CIPPD3[i,1] } else {PPD3_d2[i,] <- CIPPD3[i,3] - PPD3_d1[i,]} } PPD3map <- PPD3_d1/abs(PPD3_d2) PPD3map <- matrix(PPD3map, t,s) PPD3map <- levelplot(t(PPD3map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4), colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)), margin=F, scales = list(draw = FALSE),main = list(label = "PPD_Gender", cex = 1)) print(PPD3map, split = c(2, 1, 5, 2), more = TRUE) ## PPD BMI ## PPD4 <- matrix(NA,168,n.iter) for (i in 1:n.iter){ PPD4[,i] <- PPD4init[(168*(i-1) + 1):(168*i),] } CIPPD4 <- matrix(NA,168,3) for (i in 1:168){ CIPPD4[i,] <- quantile(PPD4[i,], c(0.025,0.5,0.975)) } PPD4_d1 <- as.matrix(CIPPD4[,2]); PPD4_d2 <- matrix(NA, 168, 1) for (i in 1:168){ if (PPD4_d1[i,] > 0){ PPD4_d2[i,] <- PPD4_d1[i,] - CIPPD4[i,1] } else {PPD4_d2[i,] <- CIPPD4[i,3] - PPD4_d1[i,]} } PPD4map <- PPD4_d1/abs(PPD4_d2) PPD4map <- matrix(PPD4map, t,s) PPD4map <- levelplot(t(PPD4map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4), colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)), margin=F, scales = list(draw = FALSE),main = list(label = "PPD_BMI", cex = 1)) print(PPD4map, split = c(3, 1, 5, 2), more = TRUE) ## PPD Smoker ## PPD5 <- matrix(NA,168,n.iter) for (i in 1:n.iter){ PPD5[,i] <- PPD5init[(168*(i-1) + 1):(168*i),] } CIPPD5 <- matrix(NA,168,3) for (i in 1:168){ CIPPD5[i,] <- quantile(PPD5[i,], c(0.025,0.5,0.975)) } PPD5_d1 <- as.matrix(CIPPD5[,2]); PPD5_d2 <- matrix(NA, 168, 1) for (i in 1:168){ if (PPD5_d1[i,] > 0){ PPD5_d2[i,] <- PPD5_d1[i,] - CIPPD5[i,1] } else {PPD5_d2[i,] <- CIPPD5[i,3] - PPD5_d1[i,]} } PPD5map <- PPD5_d1/abs(PPD5_d2) PPD5map <- matrix(PPD5map, t,s) PPD5map <- levelplot(t(PPD5map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4), colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)), margin=F, scales = list(draw = FALSE),main = list(label = "PPD_Smoker", cex = 1)) print(PPD5map, split = c(4, 1, 5, 2), more = TRUE) ## PPD HbA1c ## PPD6 <- matrix(NA,168,n.iter) for (i in 1:n.iter){ PPD6[,i] <- PPD6init[(168*(i-1) + 1):(168*i),] } CIPPD6 <- matrix(NA,168,3) for (i in 1:168){ CIPPD6[i,] <- quantile(PPD6[i,], c(0.025,0.5,0.975)) } PPD6_d1 <- as.matrix(CIPPD6[,2]); PPD6_d2 <- matrix(NA, 168, 1) for (i in 1:168){ if (PPD6_d1[i,] > 0){ PPD6_d2[i,] <- PPD6_d1[i,] - CIPPD6[i,1] } else {PPD6_d2[i,] <- CIPPD6[i,3] - PPD6_d1[i,]} } PPD6map <- PPD6_d1/abs(PPD6_d2) PPD6map <- matrix(PPD6map, t,s) PPD6map <- levelplot(t(PPD6map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4), colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)), margin=F, scales = list(draw = FALSE),main = list(label = "PPD_HbA1c", cex = 1)) print(PPD6map, split = c(5, 1, 5, 2), more = TRUE) ## CAL AGE ## CAL2 <- matrix(NA,168,n.iter) for (i in 1:n.iter){ CAL2[,i] <- CAL2init[(168*(i-1) + 1):(168*i),] } CICAL2 <- matrix(NA,168,3) for (i in 1:168){ CICAL2[i,] <- quantile(CAL2[i,], c(0.025,0.5,0.975)) } CAL2_d1 <- as.matrix(CICAL2[,2]); CAL2_d2 <- matrix(NA, 168, 1) for (i in 1:168){ if (CAL2_d1[i,] > 0){ CAL2_d2[i,] <- CAL2_d1[i,] - CICAL2[i,1] } else {CAL2_d2[i,] <- CICAL2[i,3] - CAL2_d1[i,]} } CAL2map <- CAL2_d1/abs(CAL2_d2) CAL2map <- matrix(CAL2map, t,s) CAL2map <- levelplot(t(CAL2map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4), colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)), margin=F, scales = list(draw = FALSE),main = list(label = "CAL_Age", cex = 1)) print(CAL2map, split = c(1, 2, 5, 2), more = TRUE) ## CAL Gender ## CAL3 <- matrix(NA,168,n.iter) for (i in 1:n.iter){ CAL3[,i] <- CAL3init[(168*(i-1) + 1):(168*i),] } CICAL3 <- matrix(NA,168,3) for (i in 1:168){ CICAL3[i,] <- quantile(CAL3[i,], c(0.025,0.5,0.975)) } CAL3_d1 <- as.matrix(CICAL3[,2]); CAL3_d2 <- matrix(NA, 168, 1) for (i in 1:168){ if (CAL3_d1[i,] > 0){ CAL3_d2[i,] <- CAL3_d1[i,] - CICAL3[i,1] } else {CAL3_d2[i,] <- CICAL3[i,3] - CAL3_d1[i,]} } CAL3map <- CAL3_d1/abs(CAL3_d2) CAL3map <- matrix(CAL3map, t,s) CAL3map <- levelplot(t(CAL3map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4), colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)), margin=F, scales = list(draw = FALSE),main = list(label = "CAL_Gender", cex = 1)) print(CAL3map, split = c(2, 2, 5, 2), more = TRUE) ## CAL BMI ## CAL4 <- matrix(NA,168,n.iter) for (i in 1:n.iter){ CAL4[,i] <- CAL4init[(168*(i-1) + 1):(168*i),] } CICAL4 <- matrix(NA,168,3) for (i in 1:168){ CICAL4[i,] <- quantile(CAL4[i,], c(0.025,0.5,0.975)) } CAL4_d1 <- as.matrix(CICAL4[,2]); CAL4_d2 <- matrix(NA, 168, 1) for (i in 1:168){ if (CAL4_d1[i,] > 0){ CAL4_d2[i,] <- CAL4_d1[i,] - CICAL4[i,1] } else {CAL4_d2[i,] <- CICAL4[i,3] - CAL4_d1[i,]} } CAL4map <- CAL4_d1/abs(CAL4_d2) CAL4map <- matrix(CAL4map, t,s) CAL4map <- levelplot(t(CAL4map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4), colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)), margin=F, scales = list(draw = FALSE),main = list(label = "CAL_BMI", cex = 1)) print(CAL4map, split = c(3, 2, 5, 2), more = TRUE) ## CAL Smoker ## CAL5 <- matrix(NA,168,n.iter) for (i in 1:n.iter){ CAL5[,i] <- CAL5init[(168*(i-1) + 1):(168*i),] } CICAL5 <- matrix(NA,168,3) for (i in 1:168){ CICAL5[i,] <- quantile(CAL5[i,], c(0.025,0.5,0.975)) } CAL5_d1 <- as.matrix(CICAL5[,2]); CAL5_d2 <- matrix(NA, 168, 1) for (i in 1:168){ if (CAL5_d1[i,] > 0){ CAL5_d2[i,] <- CAL5_d1[i,] - CICAL5[i,1] } else {CAL5_d2[i,] <- CICAL5[i,3] - CAL5_d1[i,]} } CAL5map <- CAL5_d1/abs(CAL5_d2) CAL5map <- matrix(CAL5map, t,s) CAL5map <- levelplot(t(CAL5map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4), colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)), margin=F, scales = list(draw = FALSE),main = list(label = "CAL_Smoker", cex = 1)) print(CAL5map, split = c(4, 2, 5, 2), more = TRUE) ## CAL HbA1c ## CAL6 <- matrix(NA,168,n.iter) for (i in 1:n.iter){ CAL6[,i] <- CAL6init[(168*(i-1) + 1):(168*i),] } CICAL6 <- matrix(NA,168,3) for (i in 1:168){ CICAL6[i,] <- quantile(CAL6[i,], c(0.025,0.5,0.975)) } CAL6_d1 <- as.matrix(CICAL5[,2]); CAL6_d2 <- matrix(NA, 168, 1) for (i in 1:168){ if (CAL6_d1[i,] > 0){ CAL6_d2[i,] <- CAL6_d1[i,] - CICAL6[i,1] } else {CAL6_d2[i,] <- CICAL6[i,3] - CAL6_d1[i,]} } CAL6map <- CAL6_d1/abs(CAL6_d2) CAL6map <- matrix(CAL6map, t,s) CAL6map <- levelplot(t(CAL6map),col.regions = Acol.regions,xlab="sites",ylab="teeth", at = c(-4,-1.5,-0.8,0.8,1.5,4), colorkey = FALSE, labels=as.character(c(-4,-1.5,-0.8,0.8,1.5,4)), margin=F, scales = list(draw = FALSE),main = list(label = "CAL_HbA1c", cex = 1)) print(CAL6map, split = c(5, 2, 5, 2), more = TRUE)
Now, the following plots (Figure A.3 in the Supplementary materials) represent residuals (real values - fitted values) for PPD & CAL of randomly picked patient.
setwd('C:/Users/inkku/Downloads') load("resd_PPD.rda"); load("resd_CAL.rda")
library(lattice);library(viridisLite) t = 28; s = 6; b = 2; n = 290; PPD_resd <- array(resd_PPD[,38], dim = c(t,s)) CAL_resd <- array(resd_CAL[,38], dim = c(t,s)) col <- viridis(100) PPD <- array(t(PPD_resd), dim = c(6,28,1)) dimnames(PPD) <- list(c(1:6),c(1:28),c("PPD residual")) levelplot(PPD,par.strip.text=list(cex=0.55), xlab="Sites",ylab="Teeth", scales = list(draw = FALSE), col.regions = col) CAL <- array(t(CAL_resd), dim = c(6,28,1)) dimnames(CAL) <- list(c(1:6),c(1:28),c("CAL residual")) levelplot(CAL,par.strip.text=list(cex=0.55), xlab="Sites",ylab="Teeth", scales = list(draw = FALSE), col.regions = col)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.