Title: | Spatiotemporal Boundary Detection Model for Areal Unit Data |
---|---|
Description: | Implements a spatiotemporal boundary detection model with a dissimilarity metric for areal data with inference in a Bayesian setting using Markov chain Monte Carlo (MCMC). The response variable can be modeled as Gaussian (no nugget), probit or Tobit link and spatial correlation is introduced at each time point through a conditional autoregressive (CAR) prior. Temporal correlation is introduced through a hierarchical structure and can be specified as exponential or first-order autoregressive. Full details of the package can be found in the accompanying vignette. Furthermore, the details of the package can be found in "Diagnosing Glaucoma Progression with Visual Field Data Using a Spatiotemporal Boundary Detection Method", by Berchuck et al (2018), <arXiv:1805.11636>. The paper is in press at the Journal of the American Statistical Association. |
Authors: | Samuel I. Berchuck [aut, cre] |
Maintainer: | Samuel I. Berchuck <[email protected]> |
License: | GPL (>= 2) |
Version: | 1.0.5 |
Built: | 2024-11-19 06:35:39 UTC |
Source: | CRAN |
Calculates diagnostic metrics using output from the STBDwDM
model.
diagnostics( obj, diags = c("dic", "dinf", "waic"), keepDeviance = FALSE, keepPPD = FALSE )
diagnostics( obj, diags = c("dic", "dinf", "waic"), keepDeviance = FALSE, keepPPD = FALSE )
obj |
A |
diags |
A vector of character strings indicating the diagnostics to compute. Options include: Deviance Information Criterion ("dic"), d-infinity ("dinf") and Watanabe-Akaike information criterion ("waic"). At least one option must be included. Note: The probit model cannot compute the DIC or WAIC diagnostics due to computational issues with computing the multivariate normal CDF. |
keepDeviance |
A logical indicating whether the posterior deviance distribution is returned (default = FALSE). |
keepPPD |
A logical indicating whether the posterior predictive distribution at each observed location is returned (default = FALSE). |
To assess model fit, DIC, d-infinity and WAIC are used. DIC is based on the deviance statistic and penalizes for the complexity of a model with an effective number of parameters estimate pD (Spiegelhalter et al 2002). The d-infinity posterior predictive measure is an alternative diagnostic tool to DIC, where d-infinity=P+G. The G term decreases as goodness of fit increases, and P, the penalty term, inflates as the model becomes over-fit, so small values of both of these terms and, thus, small values of d-infinity are desirable (Gelfand and Ghosh 1998). WAIC is invariant to parametrization and is asymptotically equal to Bayesian cross-validation (Watanabe 2010). WAIC = -2 * (lppd - p_waic_2). Where lppd is the log pointwise predictive density and p_waic_2 is the estimated effective number of parameters based on the variance estimator from Vehtari et al. 2016. (p_waic_1 is the mean estimator).
diagnostics
returns a list containing the diagnostics requested and
possibly the deviance and/or posterior predictive distribution objects.
Samuel I. Berchuck
Gelfand, A. E., & Ghosh, S. K. (1998). Model choice: a minimum posterior predictive loss approach. Biometrika, 1-11.
Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & Van Der Linde, A. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64(4), 583-639.
Vehtari, A., Gelman, A., & Gabry, J. (2016). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 1-20.
Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11(Dec), 3571-3594.
These Garway-Heath angles are used as the dissimilarity metric when implementing the boundary detection model for a longitudinal series of visual fields.
data(GarwayHeath)
data(GarwayHeath)
A vector with length 54, where each entry represents the angle (in degrees) that the underlying retinal nerve fiber enters the optic nerve head. The measure ranges from 0-360, where 0 is designated at the 9-o’clock position (right eye) and angles are counted counter clockwise. These angles are estimates for the Humphrey Field Analyzer-II (Carl Zeiss Meditec Inc., Dublin, CA). The 26th and 35th entries are missing as they correspond to a natural blind spot.
Garway-Heath, et al. (2000). Ophthalmology 107:10:1809–1815. (PubMed)
Binary adjacency matrix for the Humphrey Field Analyzer-II (Carl Zeiss Meditec Inc., Dublin, CA)
data(HFAII_Queen)
data(HFAII_Queen)
This adjacency matrix is formated using queen neighbor criteria, meaning two locations on the visual field are only considered neighbors if they share an edge or corner. The adjacency matrix is a 54 x 54 dimensional binary object with zeros on the diagonal and the column and row sums are equal to the number of neighbors.
Binary adjacency matrix for the Humphrey Field Analyzer-II (Carl Zeiss Meditec Inc., Dublin, CA)
data(HFAII_QueenHF)
data(HFAII_QueenHF)
This adjacency matrix is formated using queen neighbor criteria, meaning two locations on the visual field are only considered neighbors if they share an edge or corner. An additional criterion is included that locations are not considered neighbors if they fall within different hemispheres on the visual field. The adjacency matrix is a 54 x 54 dimensional binary object with zeros on the diagonal and the column and row sums are equal to the number of neighbors.
Binary adjacency matrix for the Humphrey Field Analyzer-II (Carl Zeiss Meditec Inc., Dublin, CA)
data(HFAII_Rook)
data(HFAII_Rook)
This adjacency matrix is formated using rook neighbor criteria, meaning two locations on the visual field are only considered neighbors if they share an edge. The adjacency matrix is a 54 x 54 dimensional binary object with zeros on the diagonal and the column and row sums are equal to the number of neighbors.
is.PosteriorAdj
is a general test of an object being interpretable as a
PosteriorAdj
object.
is.PosteriorAdj(x)
is.PosteriorAdj(x)
x |
object to be tested. |
The PosteriorAdj
class is defined as the posterior adjacency
object that results from the PosteriorAdj
function.
is.STBDwDM
is a general test of an object being interpretable as a
STBDwDM
object.
is.STBDwDM(x)
is.STBDwDM(x)
x |
object to be tested. |
The STBDwDM
class is defined as the regression object that
results from the STBDwDM
regression function.
Plots a heat map of the differential light sensitivity on the Humphrey Field Analyzer-II visual field.
PlotAdjacency( Wij, Visit = 1, stat = "mean", main = "Estimated Adjacencies", color.scheme = c("Black", "White"), edgewidth = 2, cornerwidth = 1/4, lwd.border = 3, color.bs = "gray", zlim = c(0, 1), legend = TRUE, DM = NULL, W = NULL )
PlotAdjacency( Wij, Visit = 1, stat = "mean", main = "Estimated Adjacencies", color.scheme = c("Black", "White"), edgewidth = 2, cornerwidth = 1/4, lwd.border = 3, color.bs = "gray", zlim = c(0, 1), legend = TRUE, DM = NULL, W = NULL )
Wij |
a |
Visit |
either an integer |
stat |
either "mean" or "sd" (only used for Visit != NA). |
main |
an overall title for the plot. |
color.scheme |
a vector of colors to be used to show the adjacencies changing. |
edgewidth |
a scalar indicating the width of the edges. |
cornerwidth |
a scalar indicating the width of the corners. |
lwd.border |
a scalar indicating width of the visual field border. |
color.bs |
one color specifying the blind spot. |
zlim |
the limits used for the legend (default are c(0,1)). |
legend |
logical, indicating whether the legend should be present (default = TRUE). |
DM |
a dissimilarity metric to be plotted at each location on the visual field (default = NULL). |
W |
an adjacency matrix that specifies the visual field, required if Wij is not provided (default = NULL). |
PlotAdjacency
is used in the application of glaucoma progression to
plot the posterior mean and standard deviation neighborhood adjacencies across the
visual field.
Samuel I. Berchuck
###Define blind spot locations on the HFA-II blind_spot <- c(26, 35) ###Load visual field adjacency matrix W <- HFAII_Queen[ -blind_spot, -blind_spot] ###Load Garway-Heath angles for dissimiliarity metric DM <- GarwayHeath[-blind_spot] #Uses Garway-Heath angles object "GarwayHeath" ###Adjacency plots PlotAdjacency(W = W, DM = DM, zlim = c(0, 180), Visit = NA, main = "Garway-Heath dissimilarity metric\n across the visual field")
###Define blind spot locations on the HFA-II blind_spot <- c(26, 35) ###Load visual field adjacency matrix W <- HFAII_Queen[ -blind_spot, -blind_spot] ###Load Garway-Heath angles for dissimiliarity metric DM <- GarwayHeath[-blind_spot] #Uses Garway-Heath angles object "GarwayHeath" ###Adjacency plots PlotAdjacency(W = W, DM = DM, zlim = c(0, 180), Visit = NA, main = "Garway-Heath dissimilarity metric\n across the visual field")
Plots a heat map of the differential light sensitivity on the Humphrey Field Analyzer-II visual field.
PlotSensitivity( Y = Y, main = "Sensitivity Estimate (dB) at each \nlocation on visual field", legend.lab = "DLS (dB)", zlim = c(10, 35), bins = 200, border = TRUE, legend = TRUE, color = c("yellow", "orange", "red"), col.bs = "grey", legend.round = 0, legend.vals = 5 )
PlotSensitivity( Y = Y, main = "Sensitivity Estimate (dB) at each \nlocation on visual field", legend.lab = "DLS (dB)", zlim = c(10, 35), bins = 200, border = TRUE, legend = TRUE, color = c("yellow", "orange", "red"), col.bs = "grey", legend.round = 0, legend.vals = 5 )
Y |
variable to be plotted on the visual field (e.g. differential light sensitivity). |
main |
an overall title for the plot. |
legend.lab |
a label for the legend (default = "DLS (dB)"). |
zlim |
the limits used for the legend (default are the minimum and maximum of Y). |
bins |
the number of bins used to refine the color palette for the figure and legend. |
border |
logical, indicating whether there should be a border around the visual field (default = TRUE). |
legend |
logical, indicating whether the legend should be present (default = TRUE). |
color |
a vector of character strings representing the color palette. |
col.bs |
color of the blind spot locations (default = "grey"). |
legend.round |
integer, indicating the digits that the legend labels are rounded to (default = 0). |
legend.vals |
integer, indicating the number of labels values to be included on the legend (default = 5). |
PlotSensitivity
is used in the application of glaucoma progression to
plot a variable across the visual field in the form of a heat map.
Samuel I. Berchuck
data(VFSeries) PlotSensitivity(Y = VFSeries$DLS[VFSeries$Visit == 1], main = "Sensitivity estimate (dB) at each \n location on visual field", legend.lab = "DLS (dB)", zlim = c(10, 35), bins = 250)
data(VFSeries) PlotSensitivity(Y = VFSeries$DLS[VFSeries$Visit == 1], main = "Sensitivity estimate (dB) at each \n location on visual field", legend.lab = "DLS (dB)", zlim = c(10, 35), bins = 250)
Plots a time series at each location of the Humphrey Field Analyzer-II visual field .
PlotVfTimeSeries( Y, Location, Time, main = "Visual field sensitivity time series \n at each location", xlab = "Time from first visit (days)", ylab = "Sensitivity (dB)", line.col = "red", line.reg = TRUE, line.type = 1 )
PlotVfTimeSeries( Y, Location, Time, main = "Visual field sensitivity time series \n at each location", xlab = "Time from first visit (days)", ylab = "Sensitivity (dB)", line.col = "red", line.reg = TRUE, line.type = 1 )
Y |
a time series variable to be plotted. |
Location |
a variable corresponding to the location on the visual field that the time series variable was observed. |
Time |
a variable corresponding to the time that the time series variable was observed. |
main |
an overall title for the plot. |
xlab |
a title for the x axis. |
ylab |
a title for the y axis. |
line.col |
color for the regression line, either character string corresponding to a color or a integer (default = "red"). |
line.reg |
logical, determines if there are regression lines printed (default = TRUE). |
line.type |
integer, specifies the type of regression line printed (default = 1). |
PlotVfTimeSeries
is used in the application of glaucoma progression.
In each cell is the observed DLS at each location over visits, with the red line
representing a linear regression trend.
Samuel I. Berchuck
data(VFSeries) PlotVfTimeSeries(Y = VFSeries$DLS, Location = VFSeries$Location, Time = VFSeries$Time, main = "Visual field sensitivity time series \n at each location", xlab = "Days from baseline visit", ylab = "Differential light sensitivity (dB)")
data(VFSeries) PlotVfTimeSeries(Y = VFSeries$DLS, Location = VFSeries$Location, Time = VFSeries$Time, main = "Visual field sensitivity time series \n at each location", xlab = "Days from baseline visit", ylab = "Differential light sensitivity (dB)")
Calculates the posterior mean and standard deviation for the neighborhood adjacencies
from the STBDwDM
model.
PosteriorAdj(object)
PosteriorAdj(object)
object |
a |
The function PosteriorAdj
calculates the posterior mean and standard
deviation of the neighborhood adjacencies for each pairwise location. The neighborhood
structure used to do this comes from Berchuck et al. 2017.
PosteriorAdj
returns a matrix containing the following columns.
i
Location i
(i.e. which row/column on the adjacency matrix W).
j
Location j
(i.e. which row/column on the adjacency matrix W).
DM
The dissimilarity metric between locations i
and j
.
meant
The posterior mean of the neighborhood adjacency between location
i
and j
at time t, t = 1, ... , Nu
.
sdt
The posterior mean of the neighborhood adjacency between location
i
and j
at time t, t = 1, ... , Nu
.
Samuel I. Berchuck
Predicts future observations from the STBDwDM
model.
## S3 method for class 'STBDwDM' predict(object, NewTimes, ...)
## S3 method for class 'STBDwDM' predict(object, NewTimes, ...)
object |
a |
NewTimes |
a numeric vector including desired time(s) points for prediction. |
... |
other arguments. |
predict.STBDwDM
uses Bayesian krigging to predict vectors at future
time points. The function returns the krigged observed outcomes along with the
observational level parameters (mu
, tau
, and alpha
).
predict.STBDwDM
returns a list containing the following objects.
MuTauAlpha
A list
containing three matrices, mu
,
tau
and alpha
. Each matrix is dimension NKeep x s
, where
s
is the number of new time points. Each matrix contains posterior
samples obtained by Bayesian krigging.
Y
A list
containing s
posterior predictive distribution
matrices. Each matrix is dimension NKeep x s
, where s
is the number of new time points. Each matrix is obtained through Bayesian krigging.
Samuel I. Berchuck
STBDwDM
is a Markov chain Monte Carlo (MCMC) sampler for a spatiotemporal
boundary detection model using the Bayesian hierarchical framework.
STBDwDM( Y, DM, W, Time, Starting = NULL, Hypers = NULL, Tuning = NULL, MCMC = NULL, Family = "tobit", TemporalStructure = "exponential", Distance = "circumference", Weights = "continuous", Rho = 0.99, ScaleY = 10, ScaleDM = 100, Seed = 54 )
STBDwDM( Y, DM, W, Time, Starting = NULL, Hypers = NULL, Tuning = NULL, MCMC = NULL, Family = "tobit", TemporalStructure = "exponential", Distance = "circumference", Weights = "continuous", Rho = 0.99, ScaleY = 10, ScaleDM = 100, Seed = 54 )
Y |
An |
DM |
An |
W |
An |
Time |
A |
Starting |
Either When |
Hypers |
Either When
|
Tuning |
Either When |
MCMC |
Either
|
Family |
Character string indicating the distribution of the observed data. Options
include: |
TemporalStructure |
Character string indicating the temporal structure of the
time observations. Options include: |
Distance |
Character string indicating the distance metric for computing the
dissimilarity metric. Options include: |
Weights |
Character string indicating the type of weight used. Options include:
|
Rho |
A scalar in |
ScaleY |
A positive scalar used for scaling the observed data, |
ScaleDM |
A positive scalar used for scaling the dissimilarity metric distances,
|
Seed |
An integer value used to set the seed for the random number generator (default = 54). |
Details of the underlying statistical model can be found in the article by Berchuck et al. (2018), "Diagnosing Glaucoma Progression with Visual Field Data Using a Spatiotemporal Boundary Detection Method", <arXiv:1805.11636>.
STBDwDM
returns a list containing the following objects
mu
NKeep x Nu
matrix
of posterior samples for mu
. The
t-th column contains posterior samples from the the t-th time point.
tau2
NKeep x Nu
matrix
of posterior samples for tau2
.
The t-th column contains posterior samples from the the t-th time point.
alpha
NKeep x Nu
matrix
of posterior samples for alpha
.
The t-th column contains posterior samples from the the t-th time point.
delta
NKeep x 3
matrix
of posterior samples for delta
.
The columns have names that describe the samples within them.
T
NKeep x 6
matrix
of posterior samples for T
. The
columns have names that describe the samples within them. The row is listed first, e.g.,
t32
refers to the entry in row 3
, column 2
.
phi
NKeep x 1
matrix
of posterior samples for phi
.
metropolis
(2 * Nu + 1) x 2
matrix
of metropolis
acceptance rates and tuners that result from the pilot adaptation. The first Nu
correspond to the Theta2
(i.e. tau2
) parameters, the next Nu
correspond to
the Theta3
(i.e. alpha
) parameters and the last row give the phi
values.
runtime
A character
string giving the runtime of the MCMC sampler.
datobj
A list
of data objects that are used in future STBDwDM
functions
and should be ignored by the user.
dataug
A list
of data augmentation objects that are used in future
STBDwDM
functions and should be ignored by the user.
Samuel I. Berchuck
Berchuck et al. (2018), "Diagnosing Glaucoma Progression with Visual Field Data Using a Spatiotemporal Boundary Detection Method", <arXiv:1805.11636>.
A dataset containing 9 visual field series from a patient of the Vein Pulsation Study Trial in Glaucoma and the Lions Eye Institute trial registry, Perth, Western Australia.
data(VFSeries)
data(VFSeries)
A data frame with 486 rows and 4 variables:
The visual field test visit number, (1, 2, ... , 9).
The observed outcome variable, differential light sensitivity (DLS).
The time of the visual field test (in days from baseline).
The location on the visual field of a Humphrey Field Analyzer-II (Carl Zeiss Meditec Inc., Dublin, CA) (1, 2, ... , 54).
https://anzctr.org.au/Trial/Registration/TrialReview.aspx?ACTRN=12608000274370
This package implements a spatiotemporal boundary detection with a dissimilarity metric for areal data with inference in a Bayesian setting using Markov chain Monte Carlo (MCMC). The response variable can be modeled as Gaussian (no nugget), probit or Tobit link and spatial correlation is introduced at each time point through a conditional autoregressive (CAR) prior. Temporal correlation is introduced through a hierarchical structure and can be specified as exponential or first-order autoregressive. Full details of the the package can be found in the accompanying vignette. Furthermore, the details of the package can be found in "Diagnosing Glaucoma Progression with Visual Field Data Using a Spatiotemporal Boundary Detection Method", by Berchuck et al (2018), <arXiv:1805.11636>. The paper is in press at the Journal of the American Statistical Association.
Samuel I. Berchuck [email protected]