| Title: | Association Among Disease Counts and Socio-Environmental Factors |
|---|---|
| Description: | Estimation of association between disease or death counts (e.g. COVID-19) and socio-environmental risk factors using a zero-inflated Bayesian spatiotemporal model. Non-spatiotemporal models and/or models without zero-inflation are also included for comparison. Functions to produce corresponding maps are also included. See Chakraborty et al. (2022) <doi:10.1007/s13253-022-00487-1> for more details on the method. |
| Authors: | Suman Majumder [cre, aut, cph], Yoon-Bae Jun [aut, cph], Sounak Chakraborty [ctb], Chae-Young Lim [ctb], Tanujit Dey [ctb] |
| Maintainer: | Suman Majumder <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 2.0.0 |
| Built: | 2025-01-30 14:40:22 UTC |
| Source: | https://github.com/sumanm47/bstzinb |
Generate posterior samples for the parameters in a Bayesian Negative Binomial Model
BNB(y, X, A, nchain=3, niter=100, nburn=20, nthin=1)BNB(y, X, A, nchain=3, niter=100, nburn=20, nthin=1)
y |
vector of counts, must be non-negative |
X |
matrix of covariates, numeric |
A |
adjacency matrix, numeric |
nchain |
positive integer, number of MCMC chains to be run |
niter |
positive integer, number of iterations in each chain |
nburn |
non-negative integer, number of iterations to be discarded as burn-in samples |
nthin |
positive integer, thinning interval |
list of posterior samples of the parameters of the model
data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res0 <- BNB(y, X, A, nchain=2, niter=100, nburn=20, nthin=1)data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res0 <- BNB(y, X, A, nchain=2, niter=100, nburn=20, nthin=1)
Generate posterior samples for the parameters in a Bayesian Spatiotemporal Negative Binomial Model
BSTNB(y,X,A, nchain=3,niter=100,nburn=20,nthin=1)BSTNB(y,X,A, nchain=3,niter=100,nburn=20,nthin=1)
y |
vector of counts, must be non-negative |
X |
matrix of covariates, numeric |
A |
adjacency matrix, numeric |
nchain |
positive integer, number of MCMC chains to be run |
niter |
positive integer, number of iterations in each chain |
nburn |
non-negative integer, number of iterations to be discarded as burn-in samples |
nthin |
positive integer, thinning interval |
list of posterior samples of the parameters of the model
data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res2 <- BSTNB(y, X, A, nchain=2, niter=100, nburn=20, nthin=1)data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res2 <- BSTNB(y, X, A, nchain=2, niter=100, nburn=20, nthin=1)
Generate posterior samples for the parameters in a Bayesian Spatiotemporal Zero Inflated Negative Binomial Model
BSTZINB(y,X,A,LinearT = TRUE, nchain=3,niter=100,nburn=20,nthin=1)BSTZINB(y,X,A,LinearT = TRUE, nchain=3,niter=100,nburn=20,nthin=1)
y |
vector of counts, must be non-negative |
X |
matrix of covariates, numeric |
A |
adjacency matrix, numeric |
LinearT |
logical, whether to fit a linear or non-linear temporal trend |
nchain |
positive integer, number of MCMC chains to be run |
niter |
positive integer, number of iterations in each chain |
nburn |
non-negative integer, number of iterations to be discarded as burn-in samples |
nthin |
positive integer, thinning interval |
list of posterior samples of the parameters of the model
data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=2, niter=100, nburn=20, nthin=1)data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=2, niter=100, nburn=20, nthin=1)
Generate posterior samples for the parameters in a Bayesian Zero Inflated Negative Binomial Model
BZINB(y,X,A, nchain=3,niter=100,nburn=20,nthin=1)BZINB(y,X,A, nchain=3,niter=100,nburn=20,nthin=1)
y |
vector of counts, must be non-negative |
X |
matrix of covariates, numeric |
A |
adjacency matrix, numeric |
nchain |
positive integer, number of MCMC chains to be run |
niter |
positive integer, number of iterations in each chain |
nburn |
non-negative integer, number of iterations to be discarded as burn-in samples |
nthin |
positive integer, thinning interval |
list of posterior samples of the parameters of the model
data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res1 <- BSTZINB(y, X, A, nchain=2, niter=100, nburn=20, nthin=1)data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res1 <- BSTZINB(y, X, A, nchain=2, niter=100, nburn=20, nthin=1)
Computes DIC for a BSTNB or BNB fitted object
compute_NB_DIC(y,bstfit,lastit,nchain)compute_NB_DIC(y,bstfit,lastit,nchain)
y |
vector of counts, must be non-negative, the response used for fitting a BSTNB or BSTP model |
bstfit |
BSTNB or BNB fitted object |
lastit |
positive integer, size of the chain used to fit BSTZINB |
nchain |
positive integer, number of chains used to fit BSTZINB |
DIC value
data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res2 <- BSTNB(y, X, A, nchain=3, niter=100, nburn=20, nthin=1) compute_NB_DIC(y,res2,lastit=(100-20)/1,nchain=3)data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res2 <- BSTNB(y, X, A, nchain=3, niter=100, nburn=20, nthin=1) compute_NB_DIC(y,res2,lastit=(100-20)/1,nchain=3)
Computes DIC for a BSTZINB fitted object
compute_ZINB_DIC(y,bstfit,lastit,nchain)compute_ZINB_DIC(y,bstfit,lastit,nchain)
y |
vector of counts, must be non-negative, the response used for fitting a BSTZINB model |
bstfit |
BSTZINB fitted object |
lastit |
positive integer, size of the chain used to fit BSTZINB |
nchain |
positive integer, number of chains used to fit BSTZINB |
DIC value
data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1) compute_ZINB_DIC(y,res3,lastit=(100-20)/1,nchain=3)data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1) compute_ZINB_DIC(y,res3,lastit=(100-20)/1,nchain=3)
Conducts a test of convergence for a given parameter in the fitted objects using the posterior samples for the said parameter
conv.test(params,nchain=3,thshold=1.96)conv.test(params,nchain=3,thshold=1.96)
params |
numeric matrix of dimension 2 (iterations x number of parameters, single chain) or 3 (iterations x number of parameters x chain, multiple chains) of posterior samples |
nchain |
positive integer, number of chains used to fit BSTZINB, BSTNB or BSTP |
thshold |
positive scalar, the threshold for testing the convergence. Defaults to 1.96 |
logical vector indicating whether convergence was achieved or not
data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1) conv.test(res3$Alpha,nchain=3)data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1) conv.test(res3$Alpha,nchain=3)
Data set containing neighborhood information for counties in the US, to be used to create adjacency matrices
county.adjacencycounty.adjacency
county.adjacencyA dataframe with 22200 rows and 4 columns
Creates the adjacency matrix for the supplied counties within the United States using the available neighborhood information
get_adj_mat(county.adjacency,Countyvec,Statevec)get_adj_mat(county.adjacency,Countyvec,Statevec)
county.adjacency |
data frame containing the neighborhood information for the counties of the entire US |
Countyvec |
character vector containing the names of the counties for which the adjacency matrix is to be computed |
Statevec |
character vector containing the names of the states the supplied counties belong to |
the corresponding adjacency matrix
data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA"))data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA"))
Produce a descending order of bar plot for time-averaged log-q estimates over quantile-representative counties
qRankPar(state.set,cname,bstfit,vn=12, cex.title=18, cex.lab=18, cex.legend=18)qRankPar(state.set,cname,bstfit,vn=12, cex.title=18, cex.lab=18, cex.legend=18)
state.set |
character vector of set of states on which the the graphics is to be made |
cname |
character vector of the names of the counties |
bstfit |
the fitted data for BSTP, BSTNB or BSTZINB |
vn |
positive integer, number of sample counties to display |
cex.title |
Positive number to control the size of the text of the main title. Defaults to 18. |
cex.lab |
Positive number to control the size of the text in the axes labels. Defaults to 18. |
cex.legend |
Positive number to control the size of the text in the legend. Defaults to 18. |
bar graph
data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1) qRankPar(state.set=c("IA"),cname=countyname,bstfit=res3,vn=12, cex.title=18, cex.lab=12, cex.legend=12)data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1) qRankPar(state.set=c("IA"),cname=countyname,bstfit=res3,vn=12, cex.title=18, cex.lab=12, cex.legend=12)
Produce a descending order of bar plot for time-averaged log-q estimates over top ranking counties
qRankParTop(state.set,cname,bstfit,vn=12, cex.title=18, cex.lab=18, cex.legend=18)qRankParTop(state.set,cname,bstfit,vn=12, cex.title=18, cex.lab=18, cex.legend=18)
state.set |
character vector of set of states on which the the graphics is to be made |
cname |
character vector of the names of the counties |
bstfit |
the fitted data for BSTP, BSTNB or BSTZINB |
vn |
positive integer, number of sample counties to display |
cex.title |
Positive number to control the size of the text of the main title. Defaults to 18. |
cex.lab |
Positive number to control the size of the text in the axes labels. Defaults to 18. |
cex.legend |
Positive number to control the size of the text in the legend. Defaults to 18. |
bar graph
data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1) qRankParTop(state.set=c("IA"),cname=countyname,bstfit=res3,vn=12, cex.title=18, cex.lab=12, cex.legend=12)data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1) qRankParTop(state.set=c("IA"),cname=countyname,bstfit=res3,vn=12, cex.title=18, cex.lab=12, cex.legend=12)
Generates a short summary table for a fitted object using BSTP, BSTNB or BSTZINB function
ResultTableSummary(bstfit)ResultTableSummary(bstfit)
bstfit |
fitted object using the function BSTP, BSTNB or BSTZINB |
summary table
data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1) ResultTableSummary(res3)data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1) ResultTableSummary(res3)
Fits BSTP, BSTNB and BSTZINB (with linear or non-linear temporal trend) to a given data and summarizes the results in a table
ResultTableSummary2(y,X,A,LinearT=FALSE, nchain=3,niter=100,nburn=20,nthin=1)ResultTableSummary2(y,X,A,LinearT=FALSE, nchain=3,niter=100,nburn=20,nthin=1)
y |
vector of counts, must be non-negative |
X |
matrix of covariates, numeric |
A |
adjacency matrix, numeric |
LinearT |
logical, whether to fit a linear or non-linear temporal trend |
nchain |
positive integer, number of MCMC chains to be run |
niter |
positive integer, number of iterations in each chain |
nburn |
non-negative integer, number of iterations to be discarded as burn-in samples |
nthin |
positive integer, thinning interval |
summary tables for the different methods
data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) ResultTableSummary2(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1)data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) ResultTableSummary2(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1)
Synthetic dataframe to be used for examples and trial runs
simdatsimdat
simdatA dataframe with 2376 rows and 5 columns: sid (region ID), tid (timepoint), y (count response), V1 (intercept), and x (covariate).
Produce a time-trend curve over the study time domain for counties in the US
TimetrendCurve(bstfit,cname,vn=5,smooth.mode=TRUE, cex.title=18, cex.lab=18, cex.legend=18)TimetrendCurve(bstfit,cname,vn=5,smooth.mode=TRUE, cex.title=18, cex.lab=18, cex.legend=18)
bstfit |
fitted object from BSTP, BSTNB or BSTZINB |
cname |
character vector of county names to use |
vn |
positive integer, number of sample counties to use |
smooth.mode |
logical, should splines be fitted to make it smooth |
cex.title |
Positive number to control the size of the text of the main title. Defaults to 18. |
cex.lab |
Positive number to control the size of the text in the axes labels. Defaults to 18. |
cex.legend |
Positive number to control the size of the text in the legend. Defaults to 18. |
time-trend curves
data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1) TimetrendCurve(res3,cname=countyname,vn=5,smooth.mode=TRUE,cex.title=18, cex.lab=12, cex.legend=12)data(simdat) y <- simdat$y X <- cbind(simdat$V1,simdat$x) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) A <- get_adj_mat(county.adjacency,countyname,c("IA")) res3 <- BSTZINB(y, X, A, LinearT=TRUE, nchain=3, niter=100, nburn=20, nthin=1) TimetrendCurve(res3,cname=countyname,vn=5,smooth.mode=TRUE,cex.title=18, cex.lab=12, cex.legend=12)
Dataframe to be used internally to make maps and get county information
USAcitiesUSAcities
USAcitiesA dataframe with 3232 rows and 4 columns: state_id (State abbreviation), county_name (County name), county_fips (FIPS codes for the counties) and population (County population).
Creates a map of any given quantity (at a selected time or averaged over time) for regions in the US specified by state and county
USDmapCount(state.sel,dat,scol,tcol=NULL,tsel=NULL,cname,uplim=NULL)USDmapCount(state.sel,dat,scol,tcol=NULL,tsel=NULL,cname,uplim=NULL)
state.sel |
character vector giving the selected states |
dat |
data frame having named components: y - the necessary quantity (numeric), sid - the region indices, tid - the time indices |
scol |
column index of the spatial regions |
tcol |
(optional) column index of the time points |
tsel |
(optional) selected time point |
cname |
character vector of county names, must match those in USAcities |
uplim |
(optional) numeric, upper limit for the given quantity |
spatial map of the required quantity over the specified region
data(simdat) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) USDmapCount(state.sel="IA",dat=simdat,scol=1,tcol=2,tsel=150,cname=countyname)data(simdat) data(county.adjacency) data(USAcities) IAcities <- subset(USAcities,state_id=="IA") countyname <- unique(IAcities$county_name) USDmapCount(state.sel="IA",dat=simdat,scol=1,tcol=2,tsel=150,cname=countyname)