R/RcppExports.R

In bayesefa: Bayesian Exploratory Factor

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' @title Sample Uniquenesses.
#' @description Sample the unique factor variances from the posterior distribution.
#' @param Y A N by J matrix of responses or augmented data.
#' @param F A N by M matrix of factor scores.
#' @param Lambda A J by M matrix of loadings.
#' @param psis_inv A J vector of inverse uniquenessess.
#' @author Albert X Man
#' 
#' @noRd
update_uniquenesses <- function(Y, F, Lambda, psis_inv) {
    invisible(.Call(`_bayesefa_update_uniquenesses`, Y, F, Lambda, psis_inv))
}

#' @title Exploratory Factor Analysis of Continuous Response Data
#' @description Implement the Man and Culpepper (2020) mode-jumping algorithm to factor analyze continuous response data. 
#' @param Y A N by J matrix of mean-centered, continuous variables. 
#' @param M An integer specifying the number of factors. 
#' @param gamma The value of the mode-jumping tuning parameter. Man and Culpepper (2020) used gamma = 0.5. 
#' @param burnin Number of burn-in iterations to discard.
#' @param chain_length The total number of iterations (burn-in + post-burn-in).
#' 
#' @return A list that contains nsamples = chain_length - burnin array draws from the posterior distribution:
#' 
#' - `LAMBDA`: A J by M by nsamples array of sampled loading matrices.
#' - `PSIs`: A J by nsamples matrix of vector of variable uniquenesses.
#' - `ROW_OUT`: A matrix of sampled row indices of founding variables for mode-jumping algorithm.
#' - `F_OUT`: An array of sampled factor scores. 
#' - `ACCEPTED`: Acceptance rates for mode-jumping Metropolis-Hastings (MH) steps.
#' 
#' @author Steven Andrew Culpepper, Albert Man
#' 
#' @references
#' 
#' Man, A. & Culpepper, S. A. (2020). A mode-jumping algorithm for Bayesian factor analysis. Journal of the American Statistical Association, doi:10.1080/01621459.2020.1773833.
#' 
#' @export
#' 
#' @examples
#' 
#' data(exchangerate)
#' 
#' #Retain complete cases and drop Month_Year column
#' X<-exchangerate[complete.cases(exchangerate),-1]
#' X<-apply(X,2, diff)
#' X<-as.matrix(scale(X))
#' 
#' #Specify the number of factors
#' my_M<-2
#' 
#' #Run the mode-jumping EFA algorithm
#' burn_in<-150
#' chain_length<-300
#' out <- EFA_Mode_Jumper(X,my_M,gamma=0.5,burnin=burn_in,chain_length)
#'   
#' #Rotate all samples to permutation positive lower triangular (PPLT) 
#' #structure with USD and FRF as factor founding variables
#'   my_lambda_sample = out$LAMBDA
#'     for (j in 1:dim(my_lambda_sample)[3]) {
#'       my_rotate = proposal2(c(1,4),my_lambda_sample[,,j],out$F_OUT[,,j])
#'       my_lambda_sample[,,j] = my_rotate$lambda
#'     }
#'     
#' #compute posterior mean of PPLT loadings
#' mLambda<-apply(my_lambda_sample,c(1,2),mean)
#'       
EFA_Mode_Jumper <- function(Y, M, gamma, burnin, chain_length = 10000L) {
    .Call(`_bayesefa_EFA_Mode_Jumper`, Y, M, gamma, burnin, chain_length)
}

#' @title Sample Bounded Normal Random Variables.
#' @description Sample truncated normal random variables.
#' @param mean location parameter.
#' @param sd scale parameter.
#' @param upper upper bound.
#' @param bound lower bound.
#' @author Steven A Culpepper
#' 
#' @noRd
rTruncNorm_bounds <- function(mean, sd, upper, bound) {
    .Call(`_bayesefa_rTruncNorm_bounds`, mean, sd, upper, bound)
}

#' @title Sample Unique Factor Variances for Mixed-Type Variables.
#' @description Sample unique factor variables for more general ordinal and bounded model.
#' @param Y A N by J matrix of responses or augmented data.
#' @param Ms model indicator where 0 = "bounded", 1 = "continuous", 2 = "binary", >2 = "ordinal".
#' @param F A N by M matrix of factor scores.
#' @param Lambda A J by M matrix of loadings.
#' @param psis_inv A J vector of inverse uniquenessess.
#' @param continuous_indicator A J vector indicated whether a variable is continuous.
#' @author Albert X Man
#' 
#' @noRd
update_uniquenesses_mixed <- function(Y, Ms, F, Lambda, psis_inv, continuous_indicator) {
    invisible(.Call(`_bayesefa_update_uniquenesses_mixed`, Y, Ms, F, Lambda, psis_inv, continuous_indicator))
}

#' @title Sample augmented data for mixed-type model
#' @description Sample augmented data for binary, ordinal, and bounded response data 
#' @param Y A N by J matrix of ordinal responses.
#' @param MISS A N by J matrix of missing data indicators.
#' @param Z A N by J matrix of augmented data.
#' @param as A J by M matrix of loadings.
#' @param bs A J vector of intercepts.
#' @param theta A N by M matrix of factor scores.
#' @param Ms A vector of the number of score categories. Note 2 = two parameter normal ogive, >2 ordinal normal ogive.
#' @param Kaps A matrix of category thresholds.
#' @param sdMH A vector of tuning parameters of length J for implementing the Cowles (1996) Metropolis-Hastings threshold sampler.
#' @param psis_inv A J vector of inverse uniquenessess.
#' @param bounds A J by 2 matrix denoting the min and max variable values. Note that bounds are only used for variable j if element j of Ms is zero.
#' @author Steven A Culpepper
#' @return A list containing:
#' 
#' - `Kaps_new`: An updated matrix of category thresholds.
#' - `MHaccept`: A binary vector indicating whether an MH proposal for thresholds was accepted. 
#'  
#' @noRd
update_WKappaZ_NA_mixed <- function(Y, MISS, Z, as, bs, theta, Ms, Kaps, sdMH, psis_inv, bounds) {
    .Call(`_bayesefa_update_WKappaZ_NA_mixed`, Y, MISS, Z, as, bs, theta, Ms, Kaps, sdMH, psis_inv, bounds)
}

#' @title Sample Intercept Parameters for Mixed-Type Model
#' @description Sample intercept parameters from the posterior distribution.
#' @param intercept A J vector of variable intercepts.
#' @param Z A N by J matrix of augmented data.
#' @param F A N by M matrix of factor scores.
#' @param Lambda A J by M matrix of loadings.
#' @param psis_inv A J vector of inverse-uniquenesses.
#' @param intercept_var A hyperparameter for the scale of item intercepts.
#' @author Albert X Man
#' 
#' @noRd
update_intercept_mixed <- function(intercept, Z, F, Lambda, psis_inv, intercept_var) {
    invisible(.Call(`_bayesefa_update_intercept_mixed`, intercept, Z, F, Lambda, psis_inv, intercept_var))
}

#' @title Exploratory Factor Analysis of Mixed-Type Responses
#' @description Implement the Man and Culpepper (2020) mode-jumping algorithm to factor analyze mixed-type response data. Missing values should be specified as a non-numeric value such as NA.
#' @param Y A N by J matrix of mixed-type item responses. 
#' @param M An interger specifying the number of factors. 
#' @param gamma The value of the mode-jumping tuning parameter. Man and Culpepper (2020) used gamma = 0.5. 
#' @param Ms model indicator where 0 = "bounded", 1 = "continuous", 2 = "binary", >2 = "ordinal".
#' @param sdMH A J vector of tuning parameters for the Cowles (1996) Metropolis-Hastings sampler for ordinal data latent thresholds.
#' @param bounds A J by 2 matrix denoting the min and max variable values. Note that bounds are only used for variable j if element j of Ms is zero.
#' @param burnin Number of burn-in iterations to discard.
#' @param chain_length The total number of iterations (burn-in + post-burn-in).
#' @return A list that contains nsamples = chain_length - burnin array draws from the posterior distribution:
#' 
#' - `LAMBDA`: A J by M by nsamples array of sampled loading matrices on the standardized metric.
#' - `PSIs`: A J by nsamples matrix of vector of variable uniquenesses on the standardized metric.
#' - `ROW_OUT`: A matrix of sampled row indices of founding variables for mode-jumping algorithm.
#' - `THRESHOLDS`: An array of sampled thresholds. 
#' - `INTERCEPTS`: Sampled variable thresholds on the standardized metric.
#' - `ACCEPTED`: Acceptance rates for mode-jumping Metropolis-Hastings (MH) steps.
#' - `MHACCEPT`: A J vector of acceptance rates for item threshold parameters. Note that binary items have an acceptance rate of zero, because MH steps are never performed.
#' - `LAMBDA_unst`: An array of unstandardized loadings.
#' - `PSIs_inv_unst`: A matrix of unstandardized uniquenesses.
#' - `THRESHOLDS_unst`: Unstandardized thresholds.
#' - `INTERCEPTS_unst`: Unstandardized intercepts.
#' 
#' @author Albert X Man, Steven Andrew Culpepper
#' @references
#' Cowles, M. K. (1996), Accelerating Monte Carlo Markov chain convergence for cumulative link generalized linear models," Statistics and Computing, 6, 101-111.
#' 
#' Man, A. & Culpepper, S. A. (2020). A mode-jumping algorithm for Bayesian factor analysis. Journal of the American Statistical Association, doi:10.1080/01621459.2020.1773833.
#' @export
IFA_Mode_Jumper_MixedResponses <- function(Y, M, gamma, Ms, sdMH, bounds, burnin, chain_length = 10000L) {
    .Call(`_bayesefa_IFA_Mode_Jumper_MixedResponses`, Y, M, gamma, Ms, sdMH, bounds, burnin, chain_length)
}

#' @title Sample Observed Data Given Augmented Data
#' @description Draw observed data given augmented values for the general mixed-type response data for binary, ordinal, and bounded variables.
#' @param threshold A vector of category thresholds.
#' @param Msj The model indicator variable.
#' @param MISS A N by J matrix of missing data indicators.
#' @param Z A scalar augmented value.
#' @param bounds The lower and upper bounds (only used if Msj = 0).
#' @author Steven Andrew Culpepper
#' 
#' @noRd
sampleY_given_Z <- function(threshold, Msj, Z, bounds) {
    .Call(`_bayesefa_sampleY_given_Z`, threshold, Msj, Z, bounds)
}

#' @title Posterior Predictions of Item Factor Analysis with Mixed Response Types
#' @description Generate posterior predictions for new variables using posterior samples.
#' @param OUTPUT A list of output from IFA_Mode_Jumper_MixedResponses.
#' @param Y A N by J matrix of item responses for predictions. Variables to predict are indicated in Y by NAs.
#' @param Ms model indicator where 0 = "bounded", 1 = "continuous", 2 = "binary", >2 = "ordinal".
#' @param variable_predict_flag A J vector. 0 = do not predict the variable; 1 = predict the variable.
#' @param bounds A J by 2 matrix denoting the min and max variable values. Note that bounds are only used for variable j if element j of Ms is zero.
#' @param n_mcmc_iterations The number of Gibbs iterations for sampling posterior predictions for factor scores and missing data. The default is 10 iterations.
#' @return array of predictions for all posterior samples provided in OUTPUT.
#' @author Steven Andrew Culpepper
#' @export
mixedresponse_posterior_prediction <- function(OUTPUT, Y, Ms, variable_predict_flag, bounds, n_mcmc_iterations = 10L) {
    .Call(`_bayesefa_mixedresponse_posterior_prediction`, OUTPUT, Y, Ms, variable_predict_flag, bounds, n_mcmc_iterations)
}

#' @title Sample truncated normal/ 
#' @description Sample a truncated normal random variable that is bounded below.
#' @param mean Location parameter for truncated normal.
#' @param sd Square root of truncated normal scale parameter.
#' @param b_lb The lower bound.
#' @author Albert X Man, Steven Andrew Culpepper
#' @return A value from the truncated normal distribution.
#' @noRd
rTruncNorm_lb <- function(mean, sd, b_lb) {
    .Call(`_bayesefa_rTruncNorm_lb`, mean, sd, b_lb)
}

#' @title Sample augmented data
#' @description Sample augmented data for binary and ordinal response data 
#' @param Y A N by J matrix of ordinal responses.
#' @param MISS A N by J matrix of missing data indicators.
#' @param Z A N by J matrix of augmented data.
#' @param as A J by M matrix of loadings.
#' @param bs A J vector of intercepts.
#' @param theta A N by M matrix of factor scores.
#' @param Ms A vector of the number of score categories. Note 2 = two parameter normal ogive, >2 ordinal normal ogive.
#' @param Kaps A matrix of category thresholds.
#' @param sdMH A vector of tuning parameters of length J for implementing the Cowles (1996) Metropolis-Hastings threshold sampler.
#' @author Steven Andrew Culpepper
#' @return A list containing:
#' 
#' - `Kaps_new`: An updated matrix of category thresholds.
#' - `MHaccept`: A binary vector indicating whether an MH proposal for thresholds was accepted. 
#'  
#' @noRd
update_WKappaZ_NA <- function(Y, MISS, Z, as, bs, theta, Ms, Kaps, sdMH) {
    .Call(`_bayesefa_update_WKappaZ_NA`, Y, MISS, Z, as, bs, theta, Ms, Kaps, sdMH)
}

#' @title Sample Intercept Parameters
#' @description Sample intercept parameters from the posterior distribution.
#' @param intercept A J vector of variable intercepts.
#' @param Z A N by J matrix of augmented data.
#' @param F A N by M matrix of factor scores.
#' @param psis_inv A J vector of inverse-uniquenesses.
#' @param intercept_var A hyperparameter for the scale of item intercepts.
#' @author Albert X Man
#' 
#' @noRd
update_intercept <- function(intercept, Z, F, Lambda, psis_inv, intercept_var) {
    invisible(.Call(`_bayesefa_update_intercept`, intercept, Z, F, Lambda, psis_inv, intercept_var))
}

#' @title Compute Deviance for Ordinal Probit Factor Model
#' @description Internal function to compute -2LL using unstandardized parameters.
#' @param N The number of observations.  (> 0)
#' @param J The number of items.  (> 0)
#' @param Y A N by J matrix of item responses.
#' @param MISS A N by J matrix of missing data indicators.
#' @param as A matrix of item loadings.
#' @param bs A vector of threshold parameters.
#' @param theta A matrix of factor scores.
#' @param Ms A vector of the number of score categories. Note 2 = two parameter normal ogive, >2 ordinal normal ogive.
#' @param Kaps A matrix of category thresholds.
#' @return -2LL.
#' @author Steven Andrew Culpepper
#' @export
min2LL_ono <- function(N, J, Y, MISS, as, bs, theta, Ms, Kaps) {
    .Call(`_bayesefa_min2LL_ono`, N, J, Y, MISS, as, bs, theta, Ms, Kaps)
}

#' @title Exploratory Item Factor Analysis for Ordinal Response Data
#' @description Implement the Man and Culpepper (2020) mode-jumping algorithm to factor analyze ordinal response data. Missing values should be specified as a non-numeric value such as NA.
#' @param Y A N by J matrix of item responses.
#' @param M The number of factors.
#' @param gamma Mode-jumping tuning parameter. Man & Culpepper used 0.5.
#' @param Ms A vector of the number of score categories. Note 2 = two parameter normal ogive, >2 ordinal normal ogive.
#' @param sdMH A vector of tuning parameters of length J for implementing the Cowles (1996) Metropolis-Hastings threshold sampler.
#' @param burnin Number of burn-in iterations to discard.
#' @param chain_length The total number of iterations (burn-in + post-burn-in).
#' @return A list that contains nsamples = chain_length - burnin array draws from the posterior distribution:
#' 
#' - `LAMBDA`: A J by M by nsamples array of sampled loading matrices on the standardized metric.
#' - `PSIs`: A J by nsamples matrix of vector of variable uniquenesses on the standardized metric.
#' - `ROW_OUT`: A matrix of sampled row indices of founding variables for mode-jumping algorithm.
#' - `THRESHOLDS`: An array of sampled thresholds. 
#' - `INTERCEPTS`: Sampled variable thresholds on the standardized metric.
#' - `ACCEPTED`: Acceptance rates for mode-jumping Metropolis-Hastings (MH) steps.
#' - `MHACCEPT`: A J vector of acceptance rates for item threshold parameters. Note that binary items have an acceptance rate of zero, because MH steps are never performed.
#' 
#' @author Albert X Man, Steven Andrew Culpepper
#' @references
#' Cowles, M. K. (1996), Accelerating Monte Carlo Markov chain convergence for cumulative link generalized linear models," Statistics and Computing, 6, 101-111.
#' 
#' Man, A. & Culpepper, S. A. (2020). A mode-jumping algorithm for Bayesian factor analysis. Journal of the American Statistical Association, doi:10.1080/01621459.2020.1773833.
#' 
#' @examples
#' 
#' #Load the psych package to apply the algorithm with the bfi dataset
#' library(psych)
#' data(bfi)
#' 
#' #subtract 1 from all items so scores start at 0
#' ydat<-as.matrix(bfi[,1:25])-1
#' J<-ncol(ydat)
#' N<-nrow(ydat)
#' 
#' #compute the number of item categories
#' Ms<-apply(ydat,2,max)+1
#' 
#' #specify burn-in and chain length
#' #in full application set these to 10000 and 20000
#' burnin = 10 
#' chain_length=20
#' 
#' out5<-IFA_Mode_Jumper(ydat,M=5,gamma=.5,Ms,sdMH=rep(.025,J),burnin,chain_length)
#' 
#' #check the acceptance rates for the threshold tuning parameter fall between .2 and .5
#' mh_acceptance_rates<-t(out5$MHACCEPT)
#' 
#' #compute mean thresholds
#' mthresholds<-apply(out5$THRESHOLDS,c(1,2),mean)
#' 
#' #compute mean intercepts
#' mintercepts<-apply(out5$INTERCEPTS,1,mean)
#' 
#' #rotate loadings to PLT
#' my_lambda_sample = out5$LAMBDA
#' my_M<-5
#' for (j in 1:dim(my_lambda_sample)[3]) {
#'   my_rotate = proposal2(1:my_M,my_lambda_sample[,,j],matrix(0,nrow=N,ncol = my_M))
#'   my_lambda_sample[,,j] = my_rotate$lambda
#' }
#' 
#' #compute mean of PLT loadings
#' mLambda<-apply(my_lambda_sample,c(1,2),mean)
#' 
#' #find promax rotation of posterior mean
#' rotatedLambda<-promax(mLambda)$loadings[1:J,1:my_M]
#' 
#' #save promax rotation matrix
#' promaxrotation<-promax(mLambda)$rotmat
#' 
#' #compute the factor correlation matrix
#' phi<-solve(t(promaxrotation)%*%promaxrotation)
#' 
#' @export
IFA_Mode_Jumper <- function(Y, M, gamma, Ms, sdMH, burnin, chain_length = 10000L) {
    .Call(`_bayesefa_IFA_Mode_Jumper`, Y, M, gamma, Ms, sdMH, burnin, chain_length)
}

#' @title Rotate loadings and factor scores to a permuted positive lower triangle
#' @description Rotates loading matrix according to specified founding variable row indices. 
#' @param new_r_idx vector of row indices of length equal to the number of factors.
#' @param lambda Loading matrix.
#' @param factors A n x m matrix of factor scores that is rotated. 
#' @author Albert X Man
#' @return A list of rotated loadings and factors.
#'
#' @export
proposal2 <- function(new_r_idx, lambda, factors) {
    .Call(`_bayesefa_proposal2`, new_r_idx, lambda, factors)
}

#' @title Propose Mode-Jumping Rotation
#' @description Function to propose a new rotation 
#' @param X A N by J matrix of responses.
#' @param lambda_mean A J by M matrix of loadings.
#' @param f_mean A N by M matrix of factor scores.
#' @param invClam The loading precision hyper-parameter.
#' @param sigma This does not appear to be used.
#' @param r_idx A M vector of permuted positive lower triangular (PPLT) row indices.
#' @param my_gamma The mode-jumping tuning parameter.
#' @author Albert X Man
#' @return A list containing:
#' 
#' - `lambda`: Rotated loading matrix. 
#' - `r_idx`: PPLT row indices.
#' - `f_mat`: Rotated factor scores.
#' - `accepted`: An indicator denoting whether a rotated candidate was accepted.
#' 
#' @noRd
mode_jump <- function(X, lambda_mean, f_mean, invClam, sigma, r_idx, my_gamma) {
    .Call(`_bayesefa_mode_jump`, X, lambda_mean, f_mean, invClam, sigma, r_idx, my_gamma)
}

#' @title Generate Random Multivariate Normal Distribution
#' @description Creates a random Multivariate Normal when given number of obs, mean, and sigma. 
#' @param n An \code{int}, which gives the number of observations.  (> 0)
#' @param mu A \code{vector} length m that represents the means of the normals.
#' @param S A \code{matrix} with dimensions m x m that provides Sigma, the covariance matrix. 
#' @return A \code{matrix} that is a Multivariate Normal distribution
#' @author James J Balamuta
#' @noRd
#' 
#' @examples 
#' #Call with the following data:
#' rmvnorm(2, c(0,0), diag(2))
#' 
rmvnorm <- function(n, mu, S) {
    .Call(`_bayesefa_rmvnorm`, n, mu, S)
}

#' @title Simulate a Gamma-Type Random Variable
#' @description Simulates a gamma-type random variable with the slice sampler. 
#' @param x The value of the given parameter.
#' @param alph The gamma-type shape parameter.
#' @param g The value of the g parameter.
#' @param d The value of the d parameter
#' @author Albert X Man
#' @return A scalar that is a gamma-type distribution.
#' @noRd
#' 
sim_gamma_type <- function(x, alph, g, d) {
    .Call(`_bayesefa_sim_gamma_type`, x, alph, g, d)
}

#' @title Sample Factor Scores
#' @description Sample factor scores from the posterior distribution. 
#' @param Y A N by J matrix of responses or augmented data.
#' @param I_K A M by M identity matrix.
#' @param F A N by M matrix of factor scores.
#' @param Lambda A J by M matrix of loadings.
#' @param psis_inv A J vector of inverse uniquenessess.
#' @author Albert X Man, Steven A Culpepper
#' 
#' @noRd
update_F_matrix <- function(Y, I_K, F, Lambda, psis_inv) {
    invisible(.Call(`_bayesefa_update_F_matrix`, Y, I_K, F, Lambda, psis_inv))
}

#' @title Sample the Loading Precision Hyperparameter.
#' @description Sample the loading precision parameter from the posterior distribution.
#' @param Lambda A J by M matrix of loadings.
#' @param invClam The loading precision hyper-parameter.
#' @author Albert X Max
#' 
#' @noRd
update_invClam <- function(Lambda, invClam) {
    invisible(.Call(`_bayesefa_update_invClam`, Lambda, invClam))
}

#' @title Set difference function
#' @description Find indices in x that are not included in y. 
#' @param x The first set.
#' @param y The second set
#' @return A vec of elements that are in x and not y.
#' @author Albert X Man
#' @noRd
#' 
my_setdiff <- function(x, y) {
    .Call(`_bayesefa_my_setdiff`, x, y)
}

#' @title Sample the Loadings.
#' @description Sample the loadings from the posterior distribution.
#' @param Y A N by J matrix of responses or augmented data.
#' @param r_idx A M vector of permuted positive lower triangular (PPLT) row indices.
#' @param F A N by M matrix of factor scores.
#' @param Lambda A J by M matrix of loadings.
#' @param psis_inv A J vector of inverse uniquenessess.
#' @param invClam The loading precision hyper-parameter.
#' @param my_gamma The mode-jumping tuning parameter.
#' @author Albert X Man
#' 
#' @noRd
update_Lambda_loadings_hard_zero <- function(Y, r_idx, F, Lambda, psis_inv, invClam, my_gamma) {
    invisible(.Call(`_bayesefa_update_Lambda_loadings_hard_zero`, Y, r_idx, F, Lambda, psis_inv, invClam, my_gamma))
}

#' @title Initialize Thresholds.
#' @description Initialize category threshold parameters.
#' @param Ms A J vector indicating the number of categories.
#' @author Steven A Culpepper
#' @noRd
kappa_initialize <- function(Ms) {
    .Call(`_bayesefa_kappa_initialize`, Ms)
}