# Spatial factor analysis using a Bayesian hierarchical model. `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. ```r 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 `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. ## Returns `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. ## Details Details of the underlying statistical model proposed by Berchuck et al. 2019. are forthcoming. ## Examples ```r ###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. ``` ## References Reference for Berchuck et al. 2019 is forthcoming.

bfa_sp function