bfa_sp {spBFA} | R Documentation |
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
)
Arguments
formula |
A |
data |
A required |
dist |
A |
time |
A |
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 |
trials |
A variable in |
family |
Character string indicating the distribution of the observed data. Options
include: |
temporal.structure |
Character string indicating the temporal kernel. Options include:
|
spatial.structure |
Character string indicating the type of spatial process. Options include:
|
starting |
Either When |
hypers |
Either When
|
tuning |
Either When |
mcmc |
Either
|
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 |
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.
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 matrixlambda
. The labels for each column are Lambda_O_M_K.eta
NKeep x (Nu x K)
matrix
of posterior samples for the latent factorseta
. The labels for each column are Eta_Nu_K.beta
NKeep x P
matrix
of posterior samples forbeta
.sigma2
NKeep x (M * (O - C))
matrix
of posterior samples for the variancessigma2
. The labels for each column are Sigma2_O_M.kappa
NKeep x ((O * (O + 1)) / 2)
matrix
of posterior samples forkappa
. The columns have names that describe the samples within them. The row is listed first, e.g.,Kappa3_2
refers to the entry in row3
, column2
.delta
NKeep x K
matrix
of posterior samples fordelta
.tau
NKeep x K
matrix
of posterior samples fortau
.upsilon
NKeep x ((K * (K + 1)) / 2)
matrix
of posterior samples forUpsilon
. The columns have names that describe the samples within them. The row is listed first, e.g.,Upsilon3_2
refers to the entry in row3
, column2
.psi
NKeep x 1
matrix
of posterior samples forpsi
.xi
NKeep x (M x O x K)
matrix
of posterior samples for factor loadings cluster labelsxi
. The labels for each column are Xi_O_M_K.rho
NKeep x 1
matrix
of posterior samples forrho
.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 futurebfa_sp
functions and should be ignored by the user.dataug
A
list
of data augmentation objects that are used in futurebfa_sp
functions and should be ignored by the user.
References
Reference for Berchuck et al. 2019 is forthcoming.
Examples
###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.