spBFA (version 1.5.0)

bfa_sp: Spatial factor analysis using a Bayesian hierarchical model.

Description

bfa_sp is a Markov chain Monte Carlo (MCMC) sampler for a Bayesian spatial factor analysis model. The spatial component is introduced using a Probit stick-breaking process prior on the factor loadings. The model is implemented using a Bayesian hierarchical framework.

Usage

bfa_sp(
  formula,
  data,
  dist,
  time,
  K,
  L = Inf,
  trials = NULL,
  family = "normal",
  temporal.structure = "exponential",
  spatial.structure = "discrete",
  starting = NULL,
  hypers = NULL,
  tuning = NULL,
  mcmc = NULL,
  seed = 54,
  gamma.shrinkage = TRUE,
  include.space = TRUE,
  clustering = TRUE
)

Value

bfa_sp returns a list containing the following objects

lambda

NKeep x (M x O x K) matrix of posterior samples for factor loadings matrix lambda. The labels for each column are Lambda_O_M_K.

eta

NKeep x (Nu x K) matrix of posterior samples for the latent factors eta. The labels for each column are Eta_Nu_K.

beta

NKeep x P matrix of posterior samples for beta.

sigma2

NKeep x (M * (O - C)) matrix of posterior samples for the variances sigma2. The labels for each column are Sigma2_O_M.

kappa

NKeep x ((O * (O + 1)) / 2) matrix of posterior samples for kappa. The columns have names that describe the samples within them. The row is listed first, e.g., Kappa3_2 refers to the entry in row 3, column 2.

delta

NKeep x K matrix of posterior samples for delta.

tau

NKeep x K matrix of posterior samples for tau.

upsilon

NKeep x ((K * (K + 1)) / 2) matrix of posterior samples for Upsilon. The columns have names that describe the samples within them. The row is listed first, e.g., Upsilon3_2 refers to the entry in row 3, column 2.

psi

NKeep x 1 matrix of posterior samples for psi.

xi

NKeep x (M x O x K) matrix of posterior samples for factor loadings cluster labels xi. The labels for each column are Xi_O_M_K.

rho

NKeep x 1 matrix of posterior samples for rho.

metropolis

2 (or 1) x 3 matrix of metropolis acceptance rates, updated tuners, and original tuners that result from the pilot adaptation.

runtime

A character string giving the runtime of the MCMC sampler.

datobj

A list of data objects that are used in future bfa_sp functions and should be ignored by the user.

dataug

A list of data augmentation objects that are used in future bfa_sp functions and should be ignored by the user.

Arguments

formula

A formula object, corresponding to the spatial factor analysis model. The response must be on the left of a ~ operator, and the terms on the right must indicate the covariates to be included in the fixed effects. If no covariates are desired a zero should be used, ~ 0.

data

A required data.frame containing the variables in the model. The data frame must contain M x O x Nu rows. Here, M represents the number of spatial locations, O the number of different observation types and Nu the number of temporal visits. The observations must be first be ordered spatially, second by observation type and then temporally. This means that the first M x O observations come from the first time point and the first M observations come the first spatial observation type.

dist

A M x M dimensional distance matrix. For a discrete spatial process the matrix contains binary adjacencies that dictate the spatial neighborhood structure and for continuous spatial processes the matrix should be a continuous distance matrix (e.g., Euclidean).

time

A Nu dimensional vector containing the observed time points in increasing order.

K

A scalar that indicates the dimension (i.e., quantity) of latent factors.

L

The number of latent clusters. If finite, a scalar indicating the number of clusters for each column of the factor loadings matrix. By default L is set at Inf so that the Probit stick-breaking process becomes an infinite mixture model.

trials

A variable in data that contains the number of trials for each of the binomial observations. If there is no count data, trials should be left missing.

family

Character string indicating the distribution of the observed data. Options include: "normal", "probit", "tobit", and "binomial". family must have either O or 1 dimension(s) (the one populates the rest). Any combination of likelihoods can be used.

temporal.structure

Character string indicating the temporal kernel. Options include: "exponential" and "ar1".

spatial.structure

Character string indicating the type of spatial process. Options include: "continuous" (i.e., Gaussian process with exponential kernel) and "discrete" (i.e., proper CAR).

starting

Either NULL or a list containing starting values to be specified for the MCMC sampler. If NULL is not chosen then none, some or all of the starting values may be specified.

When NULL is chosen then default starting values are automatically generated. Otherwise a list must be provided with names Beta, Delta, Sigma2, Kappa, Rho, Upsilon or Psi containing appropriate objects. Beta (or Delta) must either be a P (or K) dimensional vector or a scalar (the scalar populates the entire vector). Sigma2 must be either a M x (O - C) matrix or a scalar. Kappa must be a O x O dimensional matrix, Rho a scalar, Upsilon a K x K matrix, and Psi a scalar.

hypers

Either NULL or a list containing hyperparameter values to be specified for the MCMC sampler. If NULL is not chosen then none, some or all of the hyperparameter values may be specified.

When NULL is chosen then default hyperparameter values are automatically generated. These default hyperparameters are described in detail in (Berchuck et al.). Otherwise a list must be provided with names Beta, Delta, Sigma2, Kappa, Rho, Upsilon or Psi containing further hyperparameter information. These objects are themselves lists and may be constructed as follows.

Beta is a list with two objects, MuBeta and SigmaBeta. These values represent the prior mean and variance parameters for the multivariate normal prior.

Delta is a list with two objects, A1 and A2. These values represent the prior shape parameters for the multiplicative Gamma shrinkage prior.

Sigma2 is a list with two objects, A and B. These values represent the shape and scale for the variance parameters.

Kappa is a list with two objects, SmallUpsilon and BigTheta. SmallUpsilon represents the degrees of freedom parameter for the inverse-Wishart hyperprior and must be a real number scalar, while BigTheta represents the scale matrix and must be a O x O dimensional positive definite matrix.

Rho is a list with two objects, ARho and BRho. ARho represents the lower bound for the uniform hyperprior, while BRho represents the upper bound. The bounds must be specified carefully. This is only specified for continuous spatial processes.

Upsilon is a list with two objects, Zeta and Omega. Zeta represents the degrees of freedom parameter for the inverse-Wishart hyperprior and must be a real number scalar, while Omega represents the scale matrix and must be a K x K dimensional positive definite matrix.

Psi is a list with two objects, dependent on if the temporal kernel is exponential or ar1. For exponential, the two objects are APsi and BPsi. APsi represents the lower bound for the uniform hyperprior, while BPsi represents the upper bound. The bounds must be specified carefully. For ar1, the two objects are Beta and Gamma, which are the two shape parameters of a Beta distribution shifted to have domain in (-1, 1).

tuning

Either NULL or a list containing tuning values to be specified for the MCMC Metropolis steps. If NULL is not chosen then all of the tuning values must be specified.

When NULL is chosen then default tuning values are automatically generated to 1. Otherwise a list must be provided with names Psi, or Rho. Each of these entries must be scalars containing tuning variances for their corresponding Metropolis updates. Rho is only specified for continuous spatial processes.

mcmc

Either NULL or a list containing input values to be used for implementing the MCMC sampler. If NULL is not chosen then all of the MCMC input values must be specified.

NBurn: The number of sampler scans included in the burn-in phase. (default = 10,000)

NSims: The number of post-burn-in scans for which to perform the sampler. (default = 10,000)

NThin: Value such that during the post-burn-in phase, only every NThin-th scan is recorded for use in posterior inference (For return values we define, NKeep = NSims / NThin (default = 1).

NPilot: The number of times during the burn-in phase that pilot adaptation is performed (default = 20)

seed

An integer value used to set the seed for the random number generator (default = 54).

gamma.shrinkage

A logical indicating whether a gamma shrinkage process prior is used for the variances of the factor loadings columns. If FALSE, the hyperparameters (A1 and A2) indicate the shape and rate for a gamma prior on the precisions. Default is TRUE.

include.space

A logical indicating whether a spatial process should be included. Default is TRUE, however if FALSE the spatial correlation matrix is fixed as an identity matrix. This specification overrides the spatial.structure input.

clustering

A logical indicating whether the Bayesian non-parametric process should be used, default is TRUE. If FALSE is specified each column is instead modeled with an independent spatial process.

Details

Details of the underlying statistical model proposed by Berchuck et al. 2019. are forthcoming.

References

Reference for Berchuck et al. 2019 is forthcoming.

Examples

# \donttest{
###Load womblR for example visual field data
library(womblR)

###Format data for MCMC sampler
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
dat <- data.frame(Y = VFSeries$DLS / 10) # create data frame with scaled data
Time <- unique(VFSeries$Time) / 365 # years since baseline visit
W <- HFAII_Queen[-blind_spot, -blind_spot] # visual field adjacency matrix (data object from womblR)
M <- dim(W)[1] # number of locations

###Prior bounds for psi
TimeDist <- as.matrix(dist(Time))
BPsi <- log(0.025) / -min(TimeDist[TimeDist > 0])
APsi <- log(0.975) / -max(TimeDist)

###MCMC options
K <- 10 # number of latent factors
O <- 1 # number of spatial observation types
Hypers <- list(Sigma2 = list(A = 0.001, B = 0.001),
               Kappa = list(SmallUpsilon = O + 1, BigTheta = diag(O)),
               Delta = list(A1 = 1, A2 = 20),
               Psi = list(APsi = APsi, BPsi = BPsi),
               Upsilon = list(Zeta = K + 1, Omega = diag(K)))
Starting <- list(Sigma2 = 1,
                 Kappa = diag(O),
                 Delta = 2 * (1:K),
                 Psi = (APsi + BPsi) / 2,
                 Upsilon = diag(K))
Tuning <- list(Psi = 1)
MCMC <- list(NBurn = 1000, NSims = 1000, NThin = 2, NPilot = 5)

###Fit MCMC Sampler
reg.bfa_sp <- bfa_sp(Y ~ 0, data = dat, dist = W, time = Time,  K = 10, 
                     starting = Starting, hypers = Hypers, tuning = Tuning, mcmc = MCMC,
                     L = Inf,
                     family = "tobit",
                     trials = NULL,
                     temporal.structure = "exponential",
                     spatial.structure = "discrete",
                     seed = 54, 
                     gamma.shrinkage = TRUE,
                     include.space = TRUE,
                     clustering = TRUE)

###Note that this code produces the pre-computed data object "reg.bfa_sp"
###More details can be found in the vignette.

# }