inst/doc/womblR-example.R

In womblR: Spatiotemporal Boundary Detection Model for Areal Unit Data

## ---- echo = FALSE------------------------------------------------------------
###Start with a clean space
# rm(list = ls())

###Take care of some stuff that I don't want the user to see...
# path.package <- "/Users/Sam/Documents/Sam/School/Dissertation/Packages/womblR/"
# suppressMessages(devtools::load_all(path.package)) #loads scripts
# suppressMessages(devtools::document(path.package)) #creates documentation
###Make sure to remove devtools from Suggests line in DESCRIPTION before submission

## -----------------------------------------------------------------------------
library(womblR)

## -----------------------------------------------------------------------------
head(VFSeries)

## ---- fig.align="center", fig.width = 5.5, fig.height = 5.5-------------------
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)",
                 line.col = 1, line.type = 1, line.reg = FALSE)

## -----------------------------------------------------------------------------
blind_spot <- c(26, 35) # define blind spot
VFSeries <- VFSeries[order(VFSeries$Location), ] # sort by location
VFSeries <- VFSeries[order(VFSeries$Visit), ] # sort by visit
VFSeries <- VFSeries[!VFSeries$Location %in% blind_spot, ] # remove blind spot locations
Y <- VFSeries$DLS # define observed outcome data

## -----------------------------------------------------------------------------
Time <- unique(VFSeries$Time) / 365 # years since baseline visit
print(Time)

## -----------------------------------------------------------------------------
W <- HFAII_Queen[-blind_spot, -blind_spot] # visual field adjacency matrix

## -----------------------------------------------------------------------------
DM <- GarwayHeath[-blind_spot] # Garway-Heath angles

## ---- fig.align="center", fig.width = 5.5, fig.height = 5.5-------------------
PlotAdjacency(W = W, DM = DM, zlim = c(0, 180), Visit = NA, 
              main = "Garway-Heath dissimilarity metric\n across the visual field")

## -----------------------------------------------------------------------------
TimeDist <- abs(outer(Time, Time, "-"))
TimeDistVec <- TimeDist[lower.tri(TimeDist)]
minDiff <- min(TimeDistVec)
maxDiff <- max(TimeDistVec)
PhiUpper <- -log(0.01) / minDiff # shortest diff goes down to 1%
PhiLower <- -log(0.95) / maxDiff # longest diff goes up to 95%

## -----------------------------------------------------------------------------
Hypers <- list(Delta = list(MuDelta = c(3, 0, 0), OmegaDelta = diag(c(1000, 1000, 1))),
               T = list(Xi = 4, Psi = diag(3)),
               Phi = list(APhi = PhiLower, BPhi = PhiUpper))

## -----------------------------------------------------------------------------
Starting <- list(Delta = c(3, 0, 0), T = diag(3), Phi = 0.5)

## -----------------------------------------------------------------------------
Nu <- length(Time) # calculate number of visits
Tuning <- list(Theta2 = rep(1, Nu), Theta3 = rep(1, Nu), Phi = 1)

## -----------------------------------------------------------------------------
MCMC <- list(NBurn = 10000, NSims = 10000, NThin = 10, NPilot = 20)

## ---- include = FALSE---------------------------------------------------------
reg.STBDwDM <- STBDwDM(Y = Y, DM = DM, W = W, Time = Time, Starting = Starting, Hypers = Hypers, Tuning = Tuning, MCMC = MCMC)

## ---- eval = FALSE------------------------------------------------------------
#  reg.STBDwDM <- STBDwDM(Y = Y, DM = DM, W = W, Time = Time,
#                         Starting = Starting, Hypers = Hypers, Tuning = Tuning, MCMC = MCMC,
#                         Family = "tobit",
#                         TemporalStructure = "exponential",
#                         Distance = "circumference",
#                         Weights = "continuous",
#                         Rho = 0.99,
#                         ScaleY = 10,
#                         ScaleDM = 100,
#                         Seed = 54)
#  ## Burn-in progress:  |*************************************************|
#  ## Sampler progress:  0%..  10%..  20%..  30%..  40%..  50%..  60%..  70%..  80%..  90%..  100%..

## -----------------------------------------------------------------------------
names(reg.STBDwDM)

## -----------------------------------------------------------------------------
library(coda)

## -----------------------------------------------------------------------------
Mu <- as.mcmc(reg.STBDwDM$mu)
Tau2 <- as.mcmc(reg.STBDwDM$tau2)
Alpha <- as.mcmc(reg.STBDwDM$alpha)
Delta <- as.mcmc(reg.STBDwDM$delta)
T <- as.mcmc(reg.STBDwDM$T)
Phi <- as.mcmc(reg.STBDwDM$phi)

## ---- fig.width = 7.18, fig.height = 5.2, echo = FALSE------------------------
par(mfrow = c(2, 3))
traceplot(Mu[, 1], ylab = expression(mu[1]), main = expression(paste("Posterior of " ~ mu[1])))
traceplot(Tau2[, 1], ylab = expression(tau[1]^2), main = expression(paste("Posterior of " ~ tau[1]^2)))
traceplot(Alpha[, 1], ylab = expression(alpha[1]), main = expression(paste("Posterior of " ~ alpha[1])))
traceplot(Delta[, 1], ylab = expression(delta[1]), main = expression(paste("Posterior of " ~ delta[1])))
traceplot(T[, 1], ylab = expression(paste(T['1,1'])), main = expression(paste("Posterior of " ~ T['1,1'])))
traceplot(Phi, ylab = expression(phi), main = expression(paste("Posterior" ~ phi)))

## ---- echo = FALSE------------------------------------------------------------
c(geweke.diag(Mu)$z[1], geweke.diag(Tau2)$z[1], geweke.diag(Alpha)$z[1],geweke.diag(Delta)$z[1], geweke.diag(T)$z[1],geweke.diag(Phi)$z)

## ---- echo = FALSE, fig.width = 7.18, fig.height = 2.6------------------------
Xlab = "Days from baseline visit"
X = Time * 365 
par(mfcol = c(1, 3))
plot(X, apply(Mu, 2, mean), main = expression(paste("Posterior mean of ",mu," over time")), ylab = expression(mu), xlab = Xlab)
plot(X, apply(Tau2, 2, mean), main = expression(paste("Posterior mean of ",tau^2," over time")), ylab = expression(tau^2), xlab = Xlab)
plot(X, apply(Alpha, 2, mean), main = expression(paste("Posterior mean of ",alpha," over time")), ylab = expression(alpha), xlab = Xlab) 

## -----------------------------------------------------------------------------
CVAlpha <- apply(Alpha, 1, cv <- function(x) sd(x) / mean(x))
STCV <- c(mean(CVAlpha), sd(CVAlpha), quantile(CVAlpha, probs = c(0.025, 0.975)))
names(STCV)[1:2] <- c("Mean", "SD")
print(STCV)

## -----------------------------------------------------------------------------
Wij <- PosteriorAdj(reg.STBDwDM)

## -----------------------------------------------------------------------------
Wij[1:6, 1:7]

## ---- fig.align="center", fig.width = 4.5, fig.height = 4.5-------------------
ColorScheme1 <- c("Black", "#636363", "#bdbdbd", "#f0f0f0", "White")
PlotAdjacency(Wij, Visit = 3, stat = "mean", 
              main = "Posterior mean adjacencies at \n visit 3 across the visual field", 
              color.scheme = ColorScheme1)

## ---- fig.align="center", fig.width = 4.5, fig.height = 4.5-------------------
ColorScheme2 <- rev(ColorScheme1)
zlimSD <- quantile(Wij[,c(5,7,9,11,13,15,17,19,21)], probs = c(0, 1))
PlotAdjacency(Wij, Visit = 4, stat = "sd", 
              main = "Posterior SD of adjacencies at \n visit 4 across the visual field", 
              zlim = zlimSD, color.scheme = ColorScheme2)

## -----------------------------------------------------------------------------
Diags <- diagnostics(reg.STBDwDM, diags = c("dic", "dinf", "waic"), keepDeviance = TRUE)

## ---- fig.align = 'center', fig.width = 4, fig.height = 3.3-------------------
Deviance <- as.mcmc(Diags$deviance)
traceplot(Deviance, ylab = "Deviance", main = "Posterior Deviance")

## ---- eval = FALSE------------------------------------------------------------
#  print(Diags)

## ---- echo = FALSE------------------------------------------------------------
unlist(Diags$dic)
unlist(Diags$dinf)
unlist(Diags$waic)

## -----------------------------------------------------------------------------
NewTimes <- Time[Nu] + c(50, 100) / 365

## -----------------------------------------------------------------------------
Predictions <- predict(reg.STBDwDM, NewTimes)

## -----------------------------------------------------------------------------
names(Predictions)

## -----------------------------------------------------------------------------
names(Predictions$MuTauAlpha)
head(Predictions$MuTauAlpha$alpha)

## -----------------------------------------------------------------------------
names(Predictions$Y)

## ---- fig.align = 'center', fig.width = 5.5, fig.height = 5.5-----------------
PlotSensitivity(Y = apply(Predictions$Y$y10, 2, median),
                main = "Sensitivity estimate (dB) at each \n location on visual field",
                legend.lab = "DLS (dB)", legend.round = 2,
                bins = 250, border = FALSE)