From 17435c1481a4b630e3e6d8298ac92da57ddabf68 Mon Sep 17 00:00:00 2001 From: James Balamuta Date: Thu, 11 Feb 2021 16:58:23 -0600 Subject: [PATCH] Check-in version of code uploaded to CRAN --- DESCRIPTION | 24 + NAMESPACE | 10 + R/RcppExports.R | 476 +++++++++++ R/bayesefa-internal.R | 16 + R/exchangerate-data.R | 31 + data/exchangerate.rda | Bin 0 -> 5027 bytes inst/include/commonEFAfunctions.h | 28 + man/EFA_Mode_Jumper.Rd | 67 ++ man/IFA_Mode_Jumper.Rd | 97 +++ man/IFA_Mode_Jumper_MixedResponses.Rd | 61 ++ man/bayesefa-package.Rd | 22 + man/exchangerate.Rd | 41 + man/min2LL_ono.Rd | 36 + man/mixedresponse_posterior_prediction.Rd | 37 + man/proposal2.Rd | 24 + src/EFA_mode_jumper4.cpp | 204 +++++ ...A_mode_jumper_mixedresponses_v4bounded.cpp | 756 ++++++++++++++++++ src/IFA_mode_jumper_v5.cpp | 615 ++++++++++++++ src/Makevars | 14 + src/Makevars.win | 14 + src/RcppExports.cpp | 385 +++++++++ src/commonEFAfunctions.cpp | 361 +++++++++ 22 files changed, 3319 insertions(+) create mode 100644 DESCRIPTION create mode 100644 NAMESPACE create mode 100644 R/RcppExports.R create mode 100644 R/bayesefa-internal.R create mode 100644 R/exchangerate-data.R create mode 100644 data/exchangerate.rda create mode 100644 inst/include/commonEFAfunctions.h create mode 100644 man/EFA_Mode_Jumper.Rd create mode 100644 man/IFA_Mode_Jumper.Rd create mode 100644 man/IFA_Mode_Jumper_MixedResponses.Rd create mode 100644 man/bayesefa-package.Rd create mode 100644 man/exchangerate.Rd create mode 100644 man/min2LL_ono.Rd create mode 100644 man/mixedresponse_posterior_prediction.Rd create mode 100644 man/proposal2.Rd create mode 100644 src/EFA_mode_jumper4.cpp create mode 100644 src/IFA_mode_jumper_mixedresponses_v4bounded.cpp create mode 100644 src/IFA_mode_jumper_v5.cpp create mode 100644 src/Makevars create mode 100644 src/Makevars.win create mode 100644 src/RcppExports.cpp create mode 100644 src/commonEFAfunctions.cpp diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..bfa7468 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,24 @@ +Package: bayesefa +Type: Package +Title: Bayesian Exploratory Factor +Version: 0.0.0.4 +Date: 2021-02-05 +Authors@R: c( + person("Albert", "Man", + email = "aman3@illinois.edu", + role = c("aut") + ), + person("Steven Andrew", "Culpepper", + email = "sculpepp@illinois.edu", + role = c("aut", "cre"), + comment = c(ORCID = "0000-0003-4226-6176") + ) + ) +Description: Exploratory Bayesian factor analysis of continuous, mixed-type, and bounded continuous variables using the mode-jumping algorithm of Man and Culpepper (2020) . +License: GPL-3 +Depends: R (>= 3.5.0) +Imports: Rcpp (>= 1.0.1), psych +LinkingTo: Rcpp, RcppArmadillo +Encoding: UTF-8 +LazyData: true +RoxygenNote: 7.1.1 \ No newline at end of file diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..7461c12 --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,10 @@ +# Generated by roxygen2: do not edit by hand + +export(EFA_Mode_Jumper) +export(IFA_Mode_Jumper) +export(IFA_Mode_Jumper_MixedResponses) +export(min2LL_ono) +export(mixedresponse_posterior_prediction) +export(proposal2) +importFrom(Rcpp,sourceCpp) +useDynLib(bayesefa, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..2c946a9 --- /dev/null +++ b/R/RcppExports.R @@ -0,0 +1,476 @@ +# 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) +} + diff --git a/R/bayesefa-internal.R b/R/bayesefa-internal.R new file mode 100644 index 0000000..af77581 --- /dev/null +++ b/R/bayesefa-internal.R @@ -0,0 +1,16 @@ +# Package documentation +#' @aliases bayesefa-package +#' @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. +#' +#' +"_PACKAGE" + + +# The following block is used by usethis to automatically manage +# roxygen namespace tags. Modify with care! +## usethis namespace: start +#' @useDynLib bayesefa, .registration = TRUE +#' @importFrom Rcpp sourceCpp +## usethis namespace: end +NULL \ No newline at end of file diff --git a/R/exchangerate-data.R b/R/exchangerate-data.R new file mode 100644 index 0000000..4226c07 --- /dev/null +++ b/R/exchangerate-data.R @@ -0,0 +1,31 @@ +#' International Exchange Rate Dataset +#' +#' This dataset was considered by West and Harrison (1997) and +#' Lopes and West (2004). The dataset consists of n = 143 monthly +#' first-differences of the exchange rates of 6 international +#' currencies against the British pound, from Jan 1975 to Dec 1986, +#' these currencies are: US dollar (USD), Canadian dollar (CAD), +#' Japanese yen (JPY), French franc (FRF), German (deutsche) mark (DEM), +#' and the Italian lira (ITL). +#' +#' @format A 143 by 7 \code{matrix} of exchange rate time-series. The variables include: +#' \describe{ +#' \item{\code{Month_Year}}{Month and year of exchange rate data.} +#' \item{\code{USD}}{US dollar} +#' \item{\code{CAD}}{Canadian dollar} +#' \item{\code{JPY}}{Japanese yen} +#' \item{\code{FRF}}{French franc} +#' \item{\code{DEM}}{German (deutsche) mark} +#' \item{\code{ITL}}{Italian lira} +#' } +#' @source +#' +#' Lopes, H. F., and West, M. (2004). Bayesian model assessment in factor analysis, Statistica Sinica, 14, 41–67. +#' +#' 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. +#' +#' West, M., and Harrison, J. (1997), Bayesian forecasting and dynamic models (2nd ed.), Berlin, Heidelberg: Springer-Verlag. +#' +#' @author Steven Culpepper +#' +"exchangerate" diff --git a/data/exchangerate.rda b/data/exchangerate.rda new file mode 100644 index 0000000000000000000000000000000000000000..fbff4aac61e0b177b0249684d791c3f4898916cf GIT binary patch literal 5027 zcmbu4DkLIoj6LL~P-_rrbO`#kr4y64mX{Lc3Sc+E!lU&G?_fZch@RD`y#EMAFI{`jG` zodbu#wkCYF`rGlSxH`Qu$|ohQg!LLsyeJYDUsx+A7Wur`$8p}n?{P_q?_;t1POXhR%4ke^!4CRI zTeG}&%hKzHmi5q`vf^BQ2{}y(gY-ar6YkCmY@A0m!U!=i^f?NdM-CJpkoZ?jn4u6{%O6 zvdTvdW(T(<8~ux=u2OjUsBG-uV8-ZsYYaOMJZ6~%HV0N~A3yOiC;0%_jkko9n9<$|WLFhkAD5#3L!I5?0|e3V~!!Edxc>ttKUNT%p5oWnY^WfDPH4^Z3S6j23;JR zZ|Ibdf(p`$mL;;^n=NiZ`$9Y@IR;{6UgSuw4fX3E9;`nI(xetH-jJX&Nj;5~3<-cP zA2s$M%$0^-PDKgS`X=|~&X;9uTmnV#2zHuo3<+T_#IAhvCeHM|o!VQhU^24Pd$O0R zP(r!;cWdtu9g)}eMkSEg^Z+=%+WB{&x-5XOB~YtZXEb+`l?99Z_BxjS!Xc|q*Hsy3 zTSqSH2mfmBYwzy46(oOQ-4MWX`+?nvadmdKOO%F^hJ?SJcErfgP}_HB5hX>U9CZns zck*Emke3rr`~B7BU1e;M(H#)>54Ysz=Do8m?5L!2ywZglA^WHO>;~_4L?moO}W^)bIWyPjnfZO)5?S z8szQyX|UKMmb?V$>ex$c5b`rghd!6eH+Dv=v%j^Ij4;sy+%gSXIf3wYowqH?4_(8K zC-P0J3)4F(H)^eOz{G-ad=BP<60ih@%?g ztW_EAsGzYgLWOm3y`o3S}mF+~RM~idT}%eKeYF?Y5t47WqCBw{h7Q z3;PJ;imHQGQjQ$d4{t5*s5u1Yo{L8fe|8s-I{M_ayYm|?nUhQr+kDO_lVd+VEIH`4 zSttR`LuzKp?nlYy=<~+g&OIX?3+Bz|59JDdA)ZCZ*f1l8&ee0`JOFp=m=N-Hi$%Ur z$+z|;9Pq|m%jnu|^<`;?9z>mAz|`$yTtmeM@C77nLH=X?)C6CYEVI z#V+`TzwwLoKx>tQ)5aWQ14wk}j5X|o-3{_1jVBCHRkwZXNW;!5 z3l;>>3pxMA!Sz7aU97VeFie+Jyu3NMZL)b5{w3&InRp`CLIv=}B)-toF5*E`s(ySG?IP zIWPPSHUPd%Z}UOhxX0^8*cQ|HCQu9ATvs-gR-za-r9twvXT{pWMAU09?@5nkpe6qU zc0!5QV(Ee?QrF;35Y^g|M>LV5zg`4Q#BRFu&T$WiBSS;MU7>_$IED-BZwVCE#*wu- z?scD))+lOHF|+aB^J z*Uwe{$P6hWC^TR{1IgXj#~z%ODRyAX95Z88hqtBH^<4+a zoR~pC4&{*JjWXi$kk<)|KKN&g)3WlE%cIKNpoX!0`i)Y5e0bABZ zXoSFGufHlEQAEa!ET&+7#GGIOi%t8BNaFQqTuEOzqOIOVKpU2Ig;I~W%qVjdJX&hkE*wi@XqI)Gu_KtV z*Jk>YCu{Eyr|pyq_h-LESX37)H^xuUT|zvCxs!R>#>oBhFs8)uEI|f$4Uc85$KMY_ zwThcm&VerMs;}IaJS`2k*+|z#6MMD;7z9mh^B>uVNVu0|A8OQoZBX}I;u=6lnJ`Q$ zbcnM~(ho-%ckGJN8IL=&va%a6nR-CQZs{avS?G+bZ1^)aEB*hOV^&p9B=rr&5iI;nb~v`_nEpU zo$sUh5{Y-S7Y?qa9-6julLbPO**f1a;v2t!uoHb~)ZZ^D7_}Da)-ubS+M>(G} z6=+ydo@)0J!B(MGj8L-5raW^|Wf(4$3j(p54Ele%63uM#buVG%c@rx^@so;D(*2k+ z&cdZnoMIe>Ym`c!`D|@Y9^s5|U-G?ok=whwJg1Vo7uLc}uAPMP`1MEnMg$HNrXt0N*#xgnlNr16XipXtzU&Z5cb5j^N~XY`q}^Y@tjdjP!84EaC` zZw=HV1PpU{<9VNtuf}k1X)T|JV{1HM^&4D@TgVFV)?{ndOT){A21AHuq7xACYZ1hK zs_Hlf=Wesdg&uvrKX=Nv)mvK{&_Hi+f`H}p$D60VP}~5wf?ivn$r<=hPQm!GegRLQ z^VSqmrlx+*NTsb?Mjt=x{gMwId(h0M?Y_KFaAMd94wsv3(9!r?T0t76pX)M&Y9kF za0Q{C(Z@*)n8u0yHC8zdGWatI`!aAw4$BT-S90pYa#~%v*+sO@^oQ$`M%T5AGt}nk zdv}thWeo329 z=En!TsU4<>@i4&9?}d<3mQM~aFI29QR#KoYm`1;pTI1H4Q#@Jy4++M*)fbZ*C^(?K zwaR}uo2gl=*Fg~So$Xiq8wq5vZkF18evv_(I)#NBG<+o-5>rr z)S6;>;5Phmc}DGSw03so8v?zCZuglZo{xw>%|kRB=RY|CYOs_vxv$o;M;hEV@}uya z>7vM_hP7JE{?8R-V-PC!Ht8;-m^Mj|fWVCZopeTghRY1NeY@yVo=WzardJ!loYsvx z2Yopd^b*g0CHt*GX;EcUAs1B#C&i7pQ{aENvHs{=aR16+_=ENSX{mlviOzn(N`%O^ z<9#OdO$MyVH!#=4XG7@|n&iNwr@qV|?hDbfhi+AZUIq|W9XCCJ?bo;0bWQ0Qnz!=* zl=UBL#`IHnJ)sC8G%(&@Hlua6`NX$OFx~KL4acF9@YMk?Qe1iNOk0R2g%NrD)>F}h zx)aZP#b|L{e8y2d=0+>e?^m3!#0GCB+dQUU}`{ z`78-5;Iz_lm4z&z7Yw#26VgsoyFCWr57T2tb6{d(suh-{2hK>H)w&Z!v;7^VI{WA;}LVUf~fKy zf@$XW1$@kX?rMddtXNAFd)tyHO|{$R=xKkHys7uO8IK+-~X6?Jue7}f+ zX_rN`E2*gBZnb0Z$(}llii*E^`kPheu7!EM(N}8KOiS`na`qZXP9VJ_#Lx0JGCFhl c$f?|hR$E_hRdG}4?A)k#D#_@BG#%Z40Crk5m;e9( literal 0 HcmV?d00001 diff --git a/inst/include/commonEFAfunctions.h b/inst/include/commonEFAfunctions.h new file mode 100644 index 0000000..2ba8335 --- /dev/null +++ b/inst/include/commonEFAfunctions.h @@ -0,0 +1,28 @@ +#include + +#ifndef COMMONEFAFUNCTIONS_H +#define COMMONEFAFUNCTIONS_H + +Rcpp::List proposal2(arma::uvec& new_r_idx,arma::mat& lambda_mean, arma::mat& f_mean); + +Rcpp::List mode_jump(const arma::mat& X, arma::mat& lambda_mean, arma::mat& f_mean, + arma::mat& invClam, const arma::vec sigma, arma::uvec r_idx, double my_gamma); + +arma::mat rmvnorm(unsigned int n, const arma::vec& mu, const arma::mat& S); + +double sim_gamma_type(double x,double alph,double g,double d); + +void update_F_matrix(const arma::mat& Y,const arma::mat& I_K,arma::mat& F,arma::mat& Lambda, + arma::vec& psis_inv); + +void update_invClam(const arma::mat& Lambda,arma::mat& invClam); + +arma::uvec my_setdiff(arma::uvec& x, const arma::uvec& y); + +void update_Lambda_loadings_hard_zero(const arma::mat& Y,arma::uvec& r_idx,arma::mat& F, + arma::mat& Lambda,arma::vec& psis_inv, + const arma::mat invClam,double my_gamma); + +arma::mat kappa_initialize(const arma::vec& Ms); + +#endif diff --git a/man/EFA_Mode_Jumper.Rd b/man/EFA_Mode_Jumper.Rd new file mode 100644 index 0000000..6e6eb0b --- /dev/null +++ b/man/EFA_Mode_Jumper.Rd @@ -0,0 +1,67 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{EFA_Mode_Jumper} +\alias{EFA_Mode_Jumper} +\title{Exploratory Factor Analysis of Continuous Response Data} +\usage{ +EFA_Mode_Jumper(Y, M, gamma, burnin, chain_length = 10000L) +} +\arguments{ +\item{Y}{A N by J matrix of mean-centered, continuous variables.} + +\item{M}{An integer specifying the number of factors.} + +\item{gamma}{The value of the mode-jumping tuning parameter. Man and Culpepper (2020) used gamma = 0.5.} + +\item{burnin}{Number of burn-in iterations to discard.} + +\item{chain_length}{The total number of iterations (burn-in + post-burn-in).} +} +\value{ +A list that contains nsamples = chain_length - burnin array draws from the posterior distribution: +\itemize{ +\item \code{LAMBDA}: A J by M by nsamples array of sampled loading matrices. +\item \code{PSIs}: A J by nsamples matrix of vector of variable uniquenesses. +\item \code{ROW_OUT}: A matrix of sampled row indices of founding variables for mode-jumping algorithm. +\item \code{F_OUT}: An array of sampled factor scores. +\item \code{ACCEPTED}: Acceptance rates for mode-jumping Metropolis-Hastings (MH) steps. +} +} +\description{ +Implement the Man and Culpepper (2020) mode-jumping algorithm to factor analyze continuous response data. +} +\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) + +} +\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. +} +\author{ +Steven Andrew Culpepper, Albert Man +} diff --git a/man/IFA_Mode_Jumper.Rd b/man/IFA_Mode_Jumper.Rd new file mode 100644 index 0000000..4b554b1 --- /dev/null +++ b/man/IFA_Mode_Jumper.Rd @@ -0,0 +1,97 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{IFA_Mode_Jumper} +\alias{IFA_Mode_Jumper} +\title{Exploratory Item Factor Analysis for Ordinal Response Data} +\usage{ +IFA_Mode_Jumper(Y, M, gamma, Ms, sdMH, burnin, chain_length = 10000L) +} +\arguments{ +\item{Y}{A N by J matrix of item responses.} + +\item{M}{The number of factors.} + +\item{gamma}{Mode-jumping tuning parameter. Man & Culpepper used 0.5.} + +\item{Ms}{A vector of the number of score categories. Note 2 = two parameter normal ogive, >2 ordinal normal ogive.} + +\item{sdMH}{A vector of tuning parameters of length J for implementing the Cowles (1996) Metropolis-Hastings threshold sampler.} + +\item{burnin}{Number of burn-in iterations to discard.} + +\item{chain_length}{The total number of iterations (burn-in + post-burn-in).} +} +\value{ +A list that contains nsamples = chain_length - burnin array draws from the posterior distribution: +\itemize{ +\item \code{LAMBDA}: A J by M by nsamples array of sampled loading matrices on the standardized metric. +\item \code{PSIs}: A J by nsamples matrix of vector of variable uniquenesses on the standardized metric. +\item \code{ROW_OUT}: A matrix of sampled row indices of founding variables for mode-jumping algorithm. +\item \code{THRESHOLDS}: An array of sampled thresholds. +\item \code{INTERCEPTS}: Sampled variable thresholds on the standardized metric. +\item \code{ACCEPTED}: Acceptance rates for mode-jumping Metropolis-Hastings (MH) steps. +\item \code{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. +} +} +\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. +} +\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) + +} +\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. +} +\author{ +Albert X Man, Steven Andrew Culpepper +} diff --git a/man/IFA_Mode_Jumper_MixedResponses.Rd b/man/IFA_Mode_Jumper_MixedResponses.Rd new file mode 100644 index 0000000..000d594 --- /dev/null +++ b/man/IFA_Mode_Jumper_MixedResponses.Rd @@ -0,0 +1,61 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{IFA_Mode_Jumper_MixedResponses} +\alias{IFA_Mode_Jumper_MixedResponses} +\title{Exploratory Factor Analysis of Mixed-Type Responses} +\usage{ +IFA_Mode_Jumper_MixedResponses( + Y, + M, + gamma, + Ms, + sdMH, + bounds, + burnin, + chain_length = 10000L +) +} +\arguments{ +\item{Y}{A N by J matrix of mixed-type item responses.} + +\item{M}{An interger specifying the number of factors.} + +\item{gamma}{The value of the mode-jumping tuning parameter. Man and Culpepper (2020) used gamma = 0.5.} + +\item{Ms}{model indicator where 0 = "bounded", 1 = "continuous", 2 = "binary", >2 = "ordinal".} + +\item{sdMH}{A J vector of tuning parameters for the Cowles (1996) Metropolis-Hastings sampler for ordinal data latent thresholds.} + +\item{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.} + +\item{burnin}{Number of burn-in iterations to discard.} + +\item{chain_length}{The total number of iterations (burn-in + post-burn-in).} +} +\value{ +A list that contains nsamples = chain_length - burnin array draws from the posterior distribution: +\itemize{ +\item \code{LAMBDA}: A J by M by nsamples array of sampled loading matrices on the standardized metric. +\item \code{PSIs}: A J by nsamples matrix of vector of variable uniquenesses on the standardized metric. +\item \code{ROW_OUT}: A matrix of sampled row indices of founding variables for mode-jumping algorithm. +\item \code{THRESHOLDS}: An array of sampled thresholds. +\item \code{INTERCEPTS}: Sampled variable thresholds on the standardized metric. +\item \code{ACCEPTED}: Acceptance rates for mode-jumping Metropolis-Hastings (MH) steps. +\item \code{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. +\item \code{LAMBDA_unst}: An array of unstandardized loadings. +\item \code{PSIs_inv_unst}: A matrix of unstandardized uniquenesses. +\item \code{THRESHOLDS_unst}: Unstandardized thresholds. +\item \code{INTERCEPTS_unst}: Unstandardized intercepts. +} +} +\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. +} +\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. +} +\author{ +Albert X Man, Steven Andrew Culpepper +} diff --git a/man/bayesefa-package.Rd b/man/bayesefa-package.Rd new file mode 100644 index 0000000..8cae65f --- /dev/null +++ b/man/bayesefa-package.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bayesefa-internal.R +\docType{package} +\name{bayesefa-package} +\alias{bayesefa-package} +\alias{_PACKAGE} +\title{bayesefa: Bayesian Exploratory Factor} +\description{ +Exploratory Bayesian factor analysis of continuous, mixed-type, and bounded continuous variables using the mode-jumping algorithm of Man and Culpepper (2020) . +} +\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. +} +\author{ +\strong{Maintainer}: Steven Andrew Culpepper \email{sculpepp@illinois.edu} (\href{https://orcid.org/0000-0003-4226-6176}{ORCID}) + +Authors: +\itemize{ + \item Albert Man \email{aman3@illinois.edu} +} + +} diff --git a/man/exchangerate.Rd b/man/exchangerate.Rd new file mode 100644 index 0000000..bd35d70 --- /dev/null +++ b/man/exchangerate.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/exchangerate-data.R +\docType{data} +\name{exchangerate} +\alias{exchangerate} +\title{International Exchange Rate Dataset} +\format{ +A 143 by 7 \code{matrix} of exchange rate time-series. The variables include: +\describe{ +\item{\code{Month_Year}}{Month and year of exchange rate data.} +\item{\code{USD}}{US dollar} +\item{\code{CAD}}{Canadian dollar} +\item{\code{JPY}}{Japanese yen} +\item{\code{FRF}}{French franc} +\item{\code{DEM}}{German (deutsche) mark} +\item{\code{ITL}}{Italian lira} +} +} +\source{ +Lopes, H. F., and West, M. (2004). Bayesian model assessment in factor analysis, Statistica Sinica, 14, 41–67. + +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. + +West, M., and Harrison, J. (1997), Bayesian forecasting and dynamic models (2nd ed.), Berlin, Heidelberg: Springer-Verlag. +} +\usage{ +exchangerate +} +\description{ +This dataset was considered by West and Harrison (1997) and +Lopes and West (2004). The dataset consists of n = 143 monthly +first-differences of the exchange rates of 6 international +currencies against the British pound, from Jan 1975 to Dec 1986, +these currencies are: US dollar (USD), Canadian dollar (CAD), +Japanese yen (JPY), French franc (FRF), German (deutsche) mark (DEM), +and the Italian lira (ITL). +} +\author{ +Steven Culpepper +} +\keyword{datasets} diff --git a/man/min2LL_ono.Rd b/man/min2LL_ono.Rd new file mode 100644 index 0000000..02e93ed --- /dev/null +++ b/man/min2LL_ono.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{min2LL_ono} +\alias{min2LL_ono} +\title{Compute Deviance for Ordinal Probit Factor Model} +\usage{ +min2LL_ono(N, J, Y, MISS, as, bs, theta, Ms, Kaps) +} +\arguments{ +\item{N}{The number of observations. (> 0)} + +\item{J}{The number of items. (> 0)} + +\item{Y}{A N by J matrix of item responses.} + +\item{MISS}{A N by J matrix of missing data indicators.} + +\item{as}{A matrix of item loadings.} + +\item{bs}{A vector of threshold parameters.} + +\item{theta}{A matrix of factor scores.} + +\item{Ms}{A vector of the number of score categories. Note 2 = two parameter normal ogive, >2 ordinal normal ogive.} + +\item{Kaps}{A matrix of category thresholds.} +} +\value{ +-2LL. +} +\description{ +Internal function to compute -2LL using unstandardized parameters. +} +\author{ +Steven Andrew Culpepper +} diff --git a/man/mixedresponse_posterior_prediction.Rd b/man/mixedresponse_posterior_prediction.Rd new file mode 100644 index 0000000..31cec36 --- /dev/null +++ b/man/mixedresponse_posterior_prediction.Rd @@ -0,0 +1,37 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{mixedresponse_posterior_prediction} +\alias{mixedresponse_posterior_prediction} +\title{Posterior Predictions of Item Factor Analysis with Mixed Response Types} +\usage{ +mixedresponse_posterior_prediction( + OUTPUT, + Y, + Ms, + variable_predict_flag, + bounds, + n_mcmc_iterations = 10L +) +} +\arguments{ +\item{OUTPUT}{A list of output from IFA_Mode_Jumper_MixedResponses.} + +\item{Y}{A N by J matrix of item responses for predictions. Variables to predict are indicated in Y by NAs.} + +\item{Ms}{model indicator where 0 = "bounded", 1 = "continuous", 2 = "binary", >2 = "ordinal".} + +\item{variable_predict_flag}{A J vector. 0 = do not predict the variable; 1 = predict the variable.} + +\item{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.} + +\item{n_mcmc_iterations}{The number of Gibbs iterations for sampling posterior predictions for factor scores and missing data. The default is 10 iterations.} +} +\value{ +array of predictions for all posterior samples provided in OUTPUT. +} +\description{ +Generate posterior predictions for new variables using posterior samples. +} +\author{ +Steven Andrew Culpepper +} diff --git a/man/proposal2.Rd b/man/proposal2.Rd new file mode 100644 index 0000000..56aecd3 --- /dev/null +++ b/man/proposal2.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{proposal2} +\alias{proposal2} +\title{Rotate loadings and factor scores to a permuted positive lower triangle} +\usage{ +proposal2(new_r_idx, lambda, factors) +} +\arguments{ +\item{new_r_idx}{vector of row indices of length equal to the number of factors.} + +\item{lambda}{Loading matrix.} + +\item{factors}{A n x m matrix of factor scores that is rotated.} +} +\value{ +A list of rotated loadings and factors. +} +\description{ +Rotates loading matrix according to specified founding variable row indices. +} +\author{ +Albert X Man +} diff --git a/src/EFA_mode_jumper4.cpp b/src/EFA_mode_jumper4.cpp new file mode 100644 index 0000000..81b7020 --- /dev/null +++ b/src/EFA_mode_jumper4.cpp @@ -0,0 +1,204 @@ +#include +#include + + + +//unique to EFA code +//sample unique factor variances +//' @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 +// [[Rcpp::export]] +void update_uniquenesses(const arma::mat& Y, + arma::mat& F, + arma::mat& Lambda, + arma::vec& psis_inv) { + // One uniqueness variance per parameter + // Shared between residuals and loadings + + // unsigned int M = Lambda.n_cols; + int J = Lambda.n_rows; + unsigned int N = F.n_rows; + + double b1 = 1; + double b2 = 1; + arma::mat Y_hat = F * Lambda.t(); + arma::mat Y_resids = Y - Y_hat; + arma::colvec d2 = trans( sum( pow(Y_resids, 2), 0) ); + + double u; + int n_obs; + for (int j = 0; j < J; j++) { + u = R::runif(0.0,1.0); + n_obs = N; + // qgamma is parameterized as shape, scale rather than shape and rate + psis_inv(j) = (R::qgamma(u,double(n_obs + 2.0*b1)/2.0,2.0/(2.0*b2 + d2[j]),1,0)); + } +} + + + + + +//' @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) +//' +// [[Rcpp::export]] +Rcpp::List EFA_Mode_Jumper(const arma::mat& Y, + unsigned int M, + double gamma, + unsigned int burnin, + unsigned int chain_length=10000){ + unsigned int N = Y.n_rows; + unsigned int J = Y.n_cols; + double my_gamma = gamma; + unsigned int chain_m_burn = chain_length-burnin; + unsigned int tmburn; + arma::mat I_K = arma::eye(M,M); + + //Saving output + arma::cube LAMBDA(J,M,chain_m_burn); + arma::cube F_OUT(N,M,chain_m_burn); + arma::mat PSIs(J,chain_m_burn); + arma::umat ROW_OUT(M, chain_m_burn); + // arma::mat MH_PROBS(4, chain_m_burn); + arma::vec ACCEPTED(chain_m_burn); + arma::vec LIKELIHOOD(chain_m_burn); + + // Initialize r_idx as PLT + arma::uvec r_idx; + // arma::uvec one_to_J = arma::shuffle(arma::linspace(0, J-1, J)); + arma::uvec one_to_J = arma::linspace(0, J-1, J); + r_idx = sort(one_to_J.subvec(0, M-1)); + + arma::mat F = arma::randn(N,M); + // arma::mat Cstart=inv(I_K+F.t()*F); + + arma::mat Lambda = arma::randn(J,M); + for(unsigned int j=0;j(J); + // arma::vec marginal_like_vec; + + for(unsigned int t = 0; t < chain_length; ++t){ + + accepted = 0; + update_Lambda_loadings_hard_zero(Y, r_idx, F, Lambda, psis_inv, invClam, my_gamma); + update_F_matrix(Y,I_K, F, Lambda, psis_inv); + update_invClam(Lambda, invClam); + update_uniquenesses(Y, F, Lambda, psis_inv); + //perform mode jump step + L = mode_jump(Y, Lambda, F, invClam, pow(psis_inv, -0.5), r_idx, my_gamma); + arma::mat temp = L["lambda"]; + Lambda = temp; + arma::uvec temp2 = L["r_idx"]; + r_idx = temp2; + arma::mat temp3 = L["f_mat"]; + F = temp3; + double temp5 = L["accepted"]; + accepted = temp5; + // double temp6 = L["marginal_like"]; + // marginal_like = temp6; + + // Save output after burn-in + if(t>(burnin-1)){ + tmburn = t-burnin; + ACCEPTED(tmburn) = accepted; + LAMBDA.slice(tmburn) = Lambda; + F_OUT.slice(tmburn) = F; + PSIs.col(tmburn) = 1./psis_inv; + ROW_OUT.col(tmburn) = r_idx; + // LIKELIHOOD(tmburn) = marginal_like; + // MH_PROBS.col(tmburn) = metropolis_probs; + } + } + + return Rcpp::List::create(Rcpp::Named("LAMBDA",LAMBDA), + Rcpp::Named("PSIs",PSIs), + Rcpp::Named("ROW_OUT", ROW_OUT), + Rcpp::Named("F_OUT", F_OUT), + // Rcpp::Named("MH_PROBS", MH_PROBS), + Rcpp::Named("ACCEPTED", ACCEPTED)//,Rcpp::Named("LIKELIHOOD", LIKELIHOOD) + ); +} + diff --git a/src/IFA_mode_jumper_mixedresponses_v4bounded.cpp b/src/IFA_mode_jumper_mixedresponses_v4bounded.cpp new file mode 100644 index 0000000..ffdfe2d --- /dev/null +++ b/src/IFA_mode_jumper_mixedresponses_v4bounded.cpp @@ -0,0 +1,756 @@ +#include +#include + + + + +//unique to IFA and mixed responses +//' @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 +// [[Rcpp::export]] +double rTruncNorm_bounds(double mean,double sd, double upper, double bound){ + double p0 = R::pnorm(bound,mean,sd, 1, 0); + double p1 = 1-p0; + double uZ = R::runif(0,1); + double Z = R::qnorm( (1-upper)*uZ*p0+upper*(p0 + uZ*p1), mean, sd, 1, 0); + return(Z); +} + +//unique to mixed response +//' @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 +// [[Rcpp::export]] +void update_uniquenesses_mixed(const arma::mat& Y, + const arma::vec& Ms, + arma::mat& F, + arma::mat& Lambda, + arma::vec& psis_inv, + const arma::vec& continuous_indicator) { + // One uniqueness variance per parameter + // Shared between residuals and loadings + + int J = Lambda.n_rows; + unsigned int N = F.n_rows; + + double b1 = 1; + double b2 = 1; + + arma::mat Y_hat = F * Lambda.t(); + arma::mat Y_resids = Y - Y_hat; + arma::colvec d2 = trans( sum( pow(Y_resids, 2), 0) ); + + double u; + int n_obs; + for (int j = 0; j < J; j++) { + if (continuous_indicator(j) == 1) { + u = R::runif(0.0,1.0); + n_obs = N; + + // qgamma is parameterized as shape, scale rather than shape and rate + psis_inv(j) = (R::qgamma(u,double(n_obs + 2.0*b1)/2.0,2.0/(2.0*b2 + d2[j]),1,0)); + } + } +} + + +//' @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 +// [[Rcpp::export]] +Rcpp::List update_WKappaZ_NA_mixed(const arma::mat& Y,//binary data + const arma::mat MISS,//missing indicators + arma::mat& Z, + const arma::mat& as,//lambda matrix is called "as" here + const arma::vec& bs,//intercept is called "bs" here + const arma::mat& theta,//F is called "theta" here + const arma::vec& Ms, + const arma::mat& Kaps,//thresholds is called "Kaps" here + const arma::vec& sdMH, + const arma::vec& psis_inv, + const arma::mat& bounds) { + + unsigned int N = Y.n_rows; + unsigned int J = Y.n_cols; + // double Yij; + // arma::vec eta(N); + arma::vec pk_on(4);//Vector that includes probabilities from Cowles (1996) in order: + //P(k_j-1,new|k_j), P(k_j+1|k_j), P(k_j-1|k_j,new), P(k_j+1,new|k_j,new) + arma::vec KapsMH(max(Ms)+1); + arma::vec p2pno(2);//,p3pno(2); //P(0|eta_ij) for 2pno draw of truncated normals + double mpart,ipart; //Ratio for m and i parts of MH step + // arma::vec W(N); + arma::vec oneN = arma::ones(N); + double uMH, uR, uZ;//, u3pno,u4pno; + arma::vec pky_on(4); + arma::vec pky_new(2); + arma::vec MHaccept=arma::zeros(J); + arma::mat Kaps_new = Kaps; + // double sumz_eta,uNC; + + for(unsigned int j=0;j(N); + //Calculate eta_ij for all of the models and individuals with data + //old code that computed eta for only nonmissing values + // for(unsigned int i=0;i2){ + //Cowles MH sampling for items with > 3 categories + // if(Ms(j)>3){ + + //Sample MH Threshold candidates // + //Assuming that Kaps is initialized with -Inf,0,kappas,Inf in a matrix + KapsMH = Kaps.col(j); + + for(unsigned int m=0; m2 = "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 +// [[Rcpp::export]] +Rcpp::List IFA_Mode_Jumper_MixedResponses(const arma::mat& Y, + unsigned int M, + double gamma, + const arma::vec& Ms,const arma::vec& sdMH, + const arma::mat& bounds, + unsigned int burnin, + unsigned int chain_length=10000){ + + //"Ms" : 1 = continuous, 2 = 2pno, >2 = ono + unsigned int N = Y.n_rows; + unsigned int J = Y.n_cols; + double my_gamma = gamma; + // arma::uvec which_j_is_continuous = find( Ms<2 ); + arma::vec continuous_indicator = arma::zeros(J); + // continuous_indicator.elem(which_j_is_continuous) = arma::ones(N,J); + MISS.elem( arma::find_nonfinite(Y) ).fill(0.0); + + //Saving output + arma::cube LAMBDA(J,M,chain_m_burn); + // arma::cube F_OUT(N,M,chain_m_burn); + arma::mat PSIs(J,chain_m_burn); + arma::umat ROW_OUT(M, chain_m_burn); + // arma::mat MH_PROBS(4, chain_m_burn); + arma::vec ACCEPTED(chain_m_burn); + // arma::vec LIKELIHOOD(chain_m_burn); + arma::cube THRESHOLDS_OUT(max(Ms)-1,J, chain_m_burn); + arma::mat INTERCEPT_OUT(J, chain_m_burn); + + //saving unstandardized output + arma::cube LAMBDA_unst(J,M,chain_m_burn); + arma::mat PSIs_inv_unst(J,chain_m_burn); + arma::cube THRESHOLDS_unst(max(Ms)+1,J, chain_m_burn); + arma::mat INTERCEPT_unst(J, chain_m_burn); + + // Initialize r_idx as PLT + arma::uvec r_idx; + arma::uvec one_to_J = arma::linspace(0, J-1, J); + r_idx = sort(one_to_J.subvec(0, M-1)); + + arma::mat thresholds = kappa_initialize(Ms); + arma::mat F = arma::zeros(N,M);//arma::randn(N,M); + // arma::mat Cstart=inv(I_K+F.t()*F); + + // arma::mat Z = arma::randn(N,J); + + arma::mat Lambda = arma::randn(J,M); + for(unsigned int j=0;j(N,J); + + arma::mat Z(N,J); + arma::rowvec meanY=arma::mean(Y); + arma::rowvec sdY=arma::stddev(Y); + for(unsigned int j=0;j2){ + arma::vec pky_new(2); + + for(unsigned int i=0;i(N, J); + arma::vec marginal_like_vec; + + // std::cout << "begin\n"; + + //save MH acceptance for thresholds + arma::vec P_ACCEPT= arma::zeros(J); + arma::vec ACCEPT= arma::zeros(J); + + arma::uvec thresholdindices = arma::linspace(1, max(Ms)-1, max(Ms)-1); + + for(unsigned int t = 0; t < chain_length; ++t){ + + accepted = 0; + + Rcpp::List step1Z = update_WKappaZ_NA_mixed(Y,//binary data + MISS,//missing indicators + Z, + Lambda,//lambda matrix is called "as" here + intercept,//intercept is called "bs" here + F,//F is called "theta" here + Ms, + thresholds,//thresholds is called "Kaps" here + sdMH, + psis_inv, + bounds); + + thresholds = Rcpp::as(step1Z[0]); + ACCEPT= Rcpp::as(step1Z[1]); + + update_intercept_mixed(intercept, Z, F, Lambda, psis_inv, intercept_var); + // std::cout << intercept_var << "\n"; + Z_centered = Z - repmat(intercept, 1, N).t(); + update_F_matrix(Z_centered, I_K, F, Lambda, psis_inv); + update_Lambda_loadings_hard_zero(Z_centered, r_idx, F, Lambda, psis_inv, invClam, my_gamma); + update_invClam(Lambda, invClam); + update_uniquenesses_mixed(Z_centered, Ms, F, Lambda, psis_inv,continuous_indicator);//note changed Z -> Z_centered 10/15/20 + + // Mode Jumping Step + L = mode_jump(Z_centered, Lambda, F, invClam, pow(psis_inv, -0.5), r_idx, my_gamma); + arma::mat temp = L["lambda"]; + Lambda = temp; + arma::uvec temp2 = L["r_idx"]; + r_idx = temp2; + arma::mat temp3 = L["f_mat"]; + F = temp3; + double temp5 = L["accepted"]; + accepted = temp5; + + // Save output after burn-in + if(t>(burnin-1)){ + + // Correct scale invariance, put loadings, intercept, and residual variance + // into standardized representation + arma::mat A1 = diagmat(pow(diagvec(Lambda * Lambda.t()) + pow(psis_inv, -1), -0.5)); + arma::mat new_Lambda = A1 * Lambda; + arma::vec new_intercept = A1 * intercept; + arma::mat new_threshold_transposed = A1 * (thresholds.rows(thresholdindices)).t(); + + arma::mat Lambda_Lambda_t = new_Lambda * new_Lambda.t(); + arma::vec psis_out = 1 - Lambda_Lambda_t.diag(); + + tmburn = t-burnin; + THRESHOLDS_OUT.slice(tmburn) = new_threshold_transposed.t();//thresholds.rows(thresholdindices); + INTERCEPT_OUT.col(tmburn) = new_intercept; + ACCEPTED(tmburn) = accepted; + LAMBDA.slice(tmburn) = A1 * Lambda; + // F_OUT.slice(tmburn) = F; + // PSIs.col(tmburn) = 1./(A*pow(psis_inv, -1)); + PSIs.col(tmburn) = psis_out; + ROW_OUT.col(tmburn) = r_idx; + + THRESHOLDS_unst.slice(tmburn) = thresholds; + INTERCEPT_unst.col(tmburn) = intercept; + LAMBDA_unst.slice(tmburn) = Lambda; + PSIs_inv_unst.col(tmburn) = psis_inv; + + P_ACCEPT = (tmburn*P_ACCEPT + ACCEPT)/(tmburn + 1); + + } + } + + return Rcpp::List::create(Rcpp::Named("LAMBDA",LAMBDA), + Rcpp::Named("PSIs",PSIs), + Rcpp::Named("ROW_OUT", ROW_OUT), + Rcpp::Named("THRESHOLDS", THRESHOLDS_OUT), + Rcpp::Named("INTERCEPTS", INTERCEPT_OUT), + Rcpp::Named("ACCEPTED", ACCEPTED), + Rcpp::Named("MHACCEPT",P_ACCEPT), + Rcpp::Named("LAMBDA_unst",LAMBDA_unst), + Rcpp::Named("PSIs_inv_unst",PSIs_inv_unst), + Rcpp::Named("THRESHOLDS_unst", THRESHOLDS_unst), + Rcpp::Named("INTERCEPTS_unst", INTERCEPT_unst) + ); +} + +//' @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 +// [[Rcpp::export]] +double sampleY_given_Z(const arma::vec& threshold,const double& Msj,const double& Z,const arma::vec& bounds){ + + double Y=0.; + if(Msj>1){ + Y=0; + unsigned int thres_index=1; + while( (Z>threshold(thres_index))&&(thres_indexbounds(1)){ + Y = bounds(1); + } + } + return Y; +} + + +//' @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 +// [[Rcpp::export]] +arma::cube mixedresponse_posterior_prediction(const Rcpp::List& OUTPUT, + const arma::mat& Y, + const arma::vec& Ms, + const arma::vec& variable_predict_flag, + const arma::mat& bounds, + unsigned n_mcmc_iterations=10){ + unsigned int N = Y.n_rows; + unsigned int J = Y.n_cols; + arma::mat MISS = arma::ones(N,J); + MISS.elem( arma::find_nonfinite(Y) ).fill(0.0); + + //extracting samples + arma::cube LAMBDA_unst = OUTPUT["LAMBDA_unst"]; + arma::mat PSIs_inv_unst = OUTPUT["PSIs_inv_unst"]; + arma::cube THRESHOLDS_unst = OUTPUT["THRESHOLDS_unst"]; + arma::mat INTERCEPTS_unst = OUTPUT["INTERCEPTS_unst"]; + + arma::uvec which_vars_to_predict = find(variable_predict_flag==1); + unsigned int n_vars_to_predict = which_vars_to_predict.n_elem; + unsigned int n_mcmc_samples = LAMBDA_unst.n_slices; + unsigned int M = LAMBDA_unst.n_cols; + arma::mat I_K = arma::eye(M,M); + + arma::cube predicted_Y=arma::zeros(N,n_vars_to_predict,n_mcmc_samples); + + //need a function to generate Y|Z, + + for(unsigned int k=0;k(N,M); + arma::mat Z(N,J); + arma::mat Z_centered(N,J); + + arma::vec pky_new(2); + + //for loop to iterate between F samples and Z samples + for(unsigned int Gibbs_its=0;Gibbs_its1){ + for(unsigned int i=0;i1){ + double Zij = Z(i,j); + predicted_Y(i,jpred,k) = sampleY_given_Z(threshold.col(j),Ms(j),Zij,(bounds.row(j)).t() ); + } + if(Ms(j)==1){ + predicted_Y(i,jpred,k)=Z(i,j); + } + if(Ms(j)==0){ + double Zij = Z(i,j); + predicted_Y(i,jpred,k) = sampleY_given_Z(threshold.col(j),Ms(j),Zij,(bounds.row(j)).t() ); + } + } + } + } + return predicted_Y; +} + diff --git a/src/IFA_mode_jumper_v5.cpp b/src/IFA_mode_jumper_v5.cpp new file mode 100644 index 0000000..156c934 --- /dev/null +++ b/src/IFA_mode_jumper_v5.cpp @@ -0,0 +1,615 @@ +#include +#include + + + + + +//function that is unique to IFA +//will need a function to sample for ordered categories too +//' @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 +// [[Rcpp::export]] +double rTruncNorm_lb(double mean,double sd, double b_lb){ + double p0 = R::pnorm(b_lb,mean,sd, 1, 0); + double p1 = 1-p0; + double uZ = R::runif(0,1); + double Z = R::qnorm(p0 + uZ*p1, mean, sd, 1, 0); + return(Z); +} + + + +//' @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 +// [[Rcpp::export]] +Rcpp::List update_WKappaZ_NA(const arma::mat& Y,//binary data + const arma::mat MISS,//missing indicators + arma::mat& Z, + const arma::mat& as,//lambda matrix is called "as" here + const arma::vec& bs,//intercept is called "bs" here + const arma::mat& theta,//F is called "theta" here + const arma::vec& Ms, + const arma::mat& Kaps,//thresholds is called "Kaps" here + const arma::vec& sdMH) { + + unsigned int N = Y.n_rows; + unsigned int J = Y.n_cols; + // double Yij; + // arma::vec eta(N); + arma::vec pk_on(4);//Vector that includes probabilities from Cowles (1996) in order: + //P(k_j-1,new|k_j), P(k_j+1|k_j), P(k_j-1|k_j,new), P(k_j+1,new|k_j,new) + arma::vec KapsMH(max(Ms)+1); + arma::vec p2pno(2);//,p3pno(2); //P(0|eta_ij) for 2pno draw of truncated normals + double mpart,ipart; //Ratio for m and i parts of MH step + // arma::vec W(N); + arma::vec oneN = arma::ones(N); + double uMH, uR, uZ;//, u3pno,u4pno; + // double Phi_eta;//,sj,tj,n0dot,n1dot,n01,n10; + // double us,ug,s_temp,pg_temp,ps_temp,g_temp,pg,ps; + // arma::vec cs_new=arma::zeros(J); + // arma::vec ss_new=arma::zeros(J); + arma::vec pky_on(4); + arma::vec pky_new(2); + arma::vec MHaccept=arma::zeros(J); + arma::mat Kaps_new = Kaps; + // double sumz_eta,uNC; + + for(unsigned int j=0;j(N); + //Calculate eta_ij for all of the models and individuals with data + //old code that computed eta for only nonmissing values + // for(unsigned int i=0;i2){ + // if(model(j)=="ONO"){ + // if(Ms(j)==2){ + // goto Model_is_2PNO; + // } + //Cowles MH sampling for items with > 3 categories + // if(Ms(j)>3){ + + //Sample MH Threshold candidates // + //Assuming that Kaps is initialized with -Inf,0,kappas,Inf in a matrix + KapsMH = Kaps.col(j); + + for(unsigned int m=0; m 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 +// [[Rcpp::export]] +double min2LL_ono(unsigned int N,unsigned int J,const arma::mat& Y,const arma::mat& MISS,const arma::mat& as, + const arma::vec& bs,const arma::mat& theta, + const arma::vec& Ms,const arma::mat& Kaps){ + + double m2ll=0.0; + double p_eta,p_eta0,yij; + + for(unsigned int j=0;j1){ + for(unsigned int i=0;i2 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 +// [[Rcpp::export]] +Rcpp::List IFA_Mode_Jumper(const arma::mat& Y, + unsigned int M,double gamma, + const arma::vec& Ms,const arma::vec& sdMH, + unsigned int burnin,unsigned int chain_length=10000){ + unsigned int N = Y.n_rows; + unsigned int J = Y.n_cols; + double my_gamma = gamma; + unsigned int chain_m_burn = chain_length-burnin; + unsigned int tmburn; + arma::mat I_K = arma::eye(M,M); + + arma::mat MISS = arma::ones(N,J); + MISS.elem( arma::find_nonfinite(Y) ).fill(0.0); + + //Saving output + arma::cube LAMBDA(J,M,chain_m_burn); + // arma::cube F_OUT(N,M,chain_m_burn); + arma::mat PSIs(J,chain_m_burn); + arma::umat ROW_OUT(M, chain_m_burn); + // arma::mat MH_PROBS(4, chain_m_burn); + arma::vec ACCEPTED(chain_m_burn); + // arma::vec LIKELIHOOD(chain_m_burn); + arma::cube THRESHOLDS_OUT(max(Ms)-1,J, chain_m_burn); + arma::mat INTERCEPT_OUT(J, chain_m_burn); + + // Initialize r_idx as PLT + arma::uvec r_idx; + arma::uvec one_to_J = arma::linspace(0, J-1, J); + r_idx = sort(one_to_J.subvec(0, M-1)); + + arma::mat thresholds = kappa_initialize(Ms); + arma::mat F = arma::zeros(N,M);//arma::randn(N,M); + arma::mat Cstart=inv(I_K+F.t()*F); + + arma::mat Lambda = arma::randn(J,M); + for(unsigned int j=0;j(N,J); + + arma::mat Z(N,J); + for(unsigned int j=0;j2){ + arma::vec pky_new(2); + + for(unsigned int i=0;i(N, J); + arma::vec marginal_like_vec; + + // std::cout << "begin\n"; + + //save MH acceptance for thresholds + arma::vec P_ACCEPT= arma::zeros(J); + arma::vec ACCEPT= arma::zeros(J); + + arma::uvec thresholdindices = arma::linspace(1, max(Ms)-1, max(Ms)-1); + + for(unsigned int t = 0; t < chain_length; ++t){ + accepted = 0; + + + //thresholds need to be in the right format for this function + //the first and last thresholds are -inf and inf + //thresholds for items with fewer categories are nan + Rcpp::List step1Z = update_WKappaZ_NA(Y,//binary data + MISS,//missing indicators + Z, + Lambda,//lambda matrix is called "as" here + intercept,//intercept is called "bs" here + F,//F is called "theta" here + Ms, + thresholds,//thresholds is called "Kaps" here + sdMH); + + thresholds = Rcpp::as(step1Z[0]); + ACCEPT= Rcpp::as(step1Z[1]); + update_intercept(intercept, Z, F, Lambda, psis_inv, intercept_var); + + Z_centered = Z - repmat(intercept, 1, N).t(); + update_F_matrix(Z_centered, I_K, F, Lambda, psis_inv); + update_Lambda_loadings_hard_zero(Z_centered, r_idx, F, Lambda, psis_inv, invClam, my_gamma); + + + // std::cout << psis_inv(1) << "\n"; + update_invClam(Lambda, invClam); + + + // Mode Jumping Step + L = mode_jump(Z_centered, Lambda, F, invClam, pow(psis_inv, -0.5), r_idx, my_gamma); + arma::mat temp = L["lambda"]; + Lambda = temp; + arma::uvec temp2 = L["r_idx"]; + r_idx = temp2; + + arma::mat temp3 = L["f_mat"]; + F = temp3; + + double temp5 = L["accepted"]; + accepted = temp5; + + // double temp6 = L["marginal_like"]; + // marginal_like = temp6; + + /* */ + + // Save output after burn-in + if(t>(burnin-1)){ + + arma::mat A1 = diagmat(pow(diagvec(Lambda * Lambda.t()) + pow(psis_inv, -1), -0.5)); + // arma::mat A2 = diagmat(pow(diagvec(Lambda * Lambda.t()) + pow(psis_inv, -1), -1)); + arma::mat new_Lambda = A1 * Lambda; + arma::vec new_intercept = A1 * intercept; + arma::mat new_threshold_transposed = A1 * (thresholds.rows(thresholdindices)).t(); + + arma::mat Lambda_Lambda_t = new_Lambda * new_Lambda.t(); + arma::vec psis_out = 1 - Lambda_Lambda_t.diag(); + + tmburn = t-burnin; + THRESHOLDS_OUT.slice(tmburn) = new_threshold_transposed.t();//thresholds.rows(thresholdindices); + INTERCEPT_OUT.col(tmburn) = new_intercept; + ACCEPTED(tmburn) = accepted; + LAMBDA.slice(tmburn) = A1 * Lambda; + // F_OUT.slice(tmburn) = F; + // PSIs.col(tmburn) = 1./(A*pow(psis_inv, -1)); + PSIs.col(tmburn) = psis_out; + ROW_OUT.col(tmburn) = r_idx; + // LIKELIHOOD(tmburn) = marginal_like; + // MH_PROBS.col(tmburn) = metropolis_probs; + + P_ACCEPT = (tmburn*P_ACCEPT + ACCEPT)/(tmburn + 1); + /* */ + } + } + return Rcpp::List::create(Rcpp::Named("LAMBDA",LAMBDA), + // Rcpp::Named("PSIs",PSIs), + Rcpp::Named("ROW_OUT", ROW_OUT), + // Rcpp::Named("F_OUT", F_OUT), + Rcpp::Named("THRESHOLDS", THRESHOLDS_OUT), + Rcpp::Named("INTERCEPTS", INTERCEPT_OUT), + // Rcpp::Named("MH_PROBS", MH_PROBS), + Rcpp::Named("ACCEPTED", ACCEPTED), + Rcpp::Named("MHACCEPT",P_ACCEPT)//, + // Rcpp::Named("LIKELIHOOD", LIKELIHOOD) + ); +} + + + + + + + + + + diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..186dc73 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,14 @@ + +## With R 3.1.0 or later, you can uncomment the following line to tell R to +## enable compilation with C++11 (where available) +## +## Also, OpenMP support in Armadillo prefers C++11 support. However, for wider +## availability of the package we do not yet enforce this here. It is however +## recommended for client packages to set it. +## +## And with R 3.4.0, and RcppArmadillo 0.7.960.*, we turn C++11 on as OpenMP +## support within Armadillo prefers / requires it +CXX_STD = CXX11 + +PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -I../inst/include/ +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/src/Makevars.win b/src/Makevars.win new file mode 100644 index 0000000..186dc73 --- /dev/null +++ b/src/Makevars.win @@ -0,0 +1,14 @@ + +## With R 3.1.0 or later, you can uncomment the following line to tell R to +## enable compilation with C++11 (where available) +## +## Also, OpenMP support in Armadillo prefers C++11 support. However, for wider +## availability of the package we do not yet enforce this here. It is however +## recommended for client packages to set it. +## +## And with R 3.4.0, and RcppArmadillo 0.7.960.*, we turn C++11 on as OpenMP +## support within Armadillo prefers / requires it +CXX_STD = CXX11 + +PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -I../inst/include/ +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp new file mode 100644 index 0000000..26b1313 --- /dev/null +++ b/src/RcppExports.cpp @@ -0,0 +1,385 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include +#include + +using namespace Rcpp; + +// update_uniquenesses +void update_uniquenesses(const arma::mat& Y, arma::mat& F, arma::mat& Lambda, arma::vec& psis_inv); +RcppExport SEXP _bayesefa_update_uniquenesses(SEXP YSEXP, SEXP FSEXP, SEXP LambdaSEXP, SEXP psis_invSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type F(FSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type Lambda(LambdaSEXP); + Rcpp::traits::input_parameter< arma::vec& >::type psis_inv(psis_invSEXP); + update_uniquenesses(Y, F, Lambda, psis_inv); + return R_NilValue; +END_RCPP +} +// EFA_Mode_Jumper +Rcpp::List EFA_Mode_Jumper(const arma::mat& Y, unsigned int M, double gamma, unsigned int burnin, unsigned int chain_length); +RcppExport SEXP _bayesefa_EFA_Mode_Jumper(SEXP YSEXP, SEXP MSEXP, SEXP gammaSEXP, SEXP burninSEXP, SEXP chain_lengthSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< unsigned int >::type M(MSEXP); + Rcpp::traits::input_parameter< double >::type gamma(gammaSEXP); + Rcpp::traits::input_parameter< unsigned int >::type burnin(burninSEXP); + Rcpp::traits::input_parameter< unsigned int >::type chain_length(chain_lengthSEXP); + rcpp_result_gen = Rcpp::wrap(EFA_Mode_Jumper(Y, M, gamma, burnin, chain_length)); + return rcpp_result_gen; +END_RCPP +} +// rTruncNorm_bounds +double rTruncNorm_bounds(double mean, double sd, double upper, double bound); +RcppExport SEXP _bayesefa_rTruncNorm_bounds(SEXP meanSEXP, SEXP sdSEXP, SEXP upperSEXP, SEXP boundSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< double >::type mean(meanSEXP); + Rcpp::traits::input_parameter< double >::type sd(sdSEXP); + Rcpp::traits::input_parameter< double >::type upper(upperSEXP); + Rcpp::traits::input_parameter< double >::type bound(boundSEXP); + rcpp_result_gen = Rcpp::wrap(rTruncNorm_bounds(mean, sd, upper, bound)); + return rcpp_result_gen; +END_RCPP +} +// update_uniquenesses_mixed +void update_uniquenesses_mixed(const arma::mat& Y, const arma::vec& Ms, arma::mat& F, arma::mat& Lambda, arma::vec& psis_inv, const arma::vec& continuous_indicator); +RcppExport SEXP _bayesefa_update_uniquenesses_mixed(SEXP YSEXP, SEXP MsSEXP, SEXP FSEXP, SEXP LambdaSEXP, SEXP psis_invSEXP, SEXP continuous_indicatorSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type Ms(MsSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type F(FSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type Lambda(LambdaSEXP); + Rcpp::traits::input_parameter< arma::vec& >::type psis_inv(psis_invSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type continuous_indicator(continuous_indicatorSEXP); + update_uniquenesses_mixed(Y, Ms, F, Lambda, psis_inv, continuous_indicator); + return R_NilValue; +END_RCPP +} +// update_WKappaZ_NA_mixed +Rcpp::List update_WKappaZ_NA_mixed(const arma::mat& Y, const arma::mat MISS, arma::mat& Z, const arma::mat& as, const arma::vec& bs, const arma::mat& theta, const arma::vec& Ms, const arma::mat& Kaps, const arma::vec& sdMH, const arma::vec& psis_inv, const arma::mat& bounds); +RcppExport SEXP _bayesefa_update_WKappaZ_NA_mixed(SEXP YSEXP, SEXP MISSSEXP, SEXP ZSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP thetaSEXP, SEXP MsSEXP, SEXP KapsSEXP, SEXP sdMHSEXP, SEXP psis_invSEXP, SEXP boundsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< const arma::mat >::type MISS(MISSSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type Z(ZSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type as(asSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type bs(bsSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type theta(thetaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type Ms(MsSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Kaps(KapsSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sdMH(sdMHSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type psis_inv(psis_invSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type bounds(boundsSEXP); + rcpp_result_gen = Rcpp::wrap(update_WKappaZ_NA_mixed(Y, MISS, Z, as, bs, theta, Ms, Kaps, sdMH, psis_inv, bounds)); + return rcpp_result_gen; +END_RCPP +} +// update_intercept_mixed +void update_intercept_mixed(arma::vec& intercept, arma::mat& Z, arma::mat& F, arma::mat& Lambda, arma::vec& psis_inv, double& intercept_var); +RcppExport SEXP _bayesefa_update_intercept_mixed(SEXP interceptSEXP, SEXP ZSEXP, SEXP FSEXP, SEXP LambdaSEXP, SEXP psis_invSEXP, SEXP intercept_varSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::vec& >::type intercept(interceptSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type Z(ZSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type F(FSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type Lambda(LambdaSEXP); + Rcpp::traits::input_parameter< arma::vec& >::type psis_inv(psis_invSEXP); + Rcpp::traits::input_parameter< double& >::type intercept_var(intercept_varSEXP); + update_intercept_mixed(intercept, Z, F, Lambda, psis_inv, intercept_var); + return R_NilValue; +END_RCPP +} +// IFA_Mode_Jumper_MixedResponses +Rcpp::List IFA_Mode_Jumper_MixedResponses(const arma::mat& Y, unsigned int M, double gamma, const arma::vec& Ms, const arma::vec& sdMH, const arma::mat& bounds, unsigned int burnin, unsigned int chain_length); +RcppExport SEXP _bayesefa_IFA_Mode_Jumper_MixedResponses(SEXP YSEXP, SEXP MSEXP, SEXP gammaSEXP, SEXP MsSEXP, SEXP sdMHSEXP, SEXP boundsSEXP, SEXP burninSEXP, SEXP chain_lengthSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< unsigned int >::type M(MSEXP); + Rcpp::traits::input_parameter< double >::type gamma(gammaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type Ms(MsSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sdMH(sdMHSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type bounds(boundsSEXP); + Rcpp::traits::input_parameter< unsigned int >::type burnin(burninSEXP); + Rcpp::traits::input_parameter< unsigned int >::type chain_length(chain_lengthSEXP); + rcpp_result_gen = Rcpp::wrap(IFA_Mode_Jumper_MixedResponses(Y, M, gamma, Ms, sdMH, bounds, burnin, chain_length)); + return rcpp_result_gen; +END_RCPP +} +// sampleY_given_Z +double sampleY_given_Z(const arma::vec& threshold, const double& Msj, const double& Z, const arma::vec& bounds); +RcppExport SEXP _bayesefa_sampleY_given_Z(SEXP thresholdSEXP, SEXP MsjSEXP, SEXP ZSEXP, SEXP boundsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::vec& >::type threshold(thresholdSEXP); + Rcpp::traits::input_parameter< const double& >::type Msj(MsjSEXP); + Rcpp::traits::input_parameter< const double& >::type Z(ZSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type bounds(boundsSEXP); + rcpp_result_gen = Rcpp::wrap(sampleY_given_Z(threshold, Msj, Z, bounds)); + return rcpp_result_gen; +END_RCPP +} +// mixedresponse_posterior_prediction +arma::cube mixedresponse_posterior_prediction(const Rcpp::List& OUTPUT, const arma::mat& Y, const arma::vec& Ms, const arma::vec& variable_predict_flag, const arma::mat& bounds, unsigned n_mcmc_iterations); +RcppExport SEXP _bayesefa_mixedresponse_posterior_prediction(SEXP OUTPUTSEXP, SEXP YSEXP, SEXP MsSEXP, SEXP variable_predict_flagSEXP, SEXP boundsSEXP, SEXP n_mcmc_iterationsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const Rcpp::List& >::type OUTPUT(OUTPUTSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type Ms(MsSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type variable_predict_flag(variable_predict_flagSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type bounds(boundsSEXP); + Rcpp::traits::input_parameter< unsigned >::type n_mcmc_iterations(n_mcmc_iterationsSEXP); + rcpp_result_gen = Rcpp::wrap(mixedresponse_posterior_prediction(OUTPUT, Y, Ms, variable_predict_flag, bounds, n_mcmc_iterations)); + return rcpp_result_gen; +END_RCPP +} +// rTruncNorm_lb +double rTruncNorm_lb(double mean, double sd, double b_lb); +RcppExport SEXP _bayesefa_rTruncNorm_lb(SEXP meanSEXP, SEXP sdSEXP, SEXP b_lbSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< double >::type mean(meanSEXP); + Rcpp::traits::input_parameter< double >::type sd(sdSEXP); + Rcpp::traits::input_parameter< double >::type b_lb(b_lbSEXP); + rcpp_result_gen = Rcpp::wrap(rTruncNorm_lb(mean, sd, b_lb)); + return rcpp_result_gen; +END_RCPP +} +// update_WKappaZ_NA +Rcpp::List update_WKappaZ_NA(const arma::mat& Y, const arma::mat MISS, arma::mat& Z, const arma::mat& as, const arma::vec& bs, const arma::mat& theta, const arma::vec& Ms, const arma::mat& Kaps, const arma::vec& sdMH); +RcppExport SEXP _bayesefa_update_WKappaZ_NA(SEXP YSEXP, SEXP MISSSEXP, SEXP ZSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP thetaSEXP, SEXP MsSEXP, SEXP KapsSEXP, SEXP sdMHSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< const arma::mat >::type MISS(MISSSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type Z(ZSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type as(asSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type bs(bsSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type theta(thetaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type Ms(MsSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Kaps(KapsSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sdMH(sdMHSEXP); + rcpp_result_gen = Rcpp::wrap(update_WKappaZ_NA(Y, MISS, Z, as, bs, theta, Ms, Kaps, sdMH)); + return rcpp_result_gen; +END_RCPP +} +// update_intercept +void update_intercept(arma::vec& intercept, arma::mat& Z, arma::mat& F, arma::mat& Lambda, arma::vec& psis_inv, double& intercept_var); +RcppExport SEXP _bayesefa_update_intercept(SEXP interceptSEXP, SEXP ZSEXP, SEXP FSEXP, SEXP LambdaSEXP, SEXP psis_invSEXP, SEXP intercept_varSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::vec& >::type intercept(interceptSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type Z(ZSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type F(FSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type Lambda(LambdaSEXP); + Rcpp::traits::input_parameter< arma::vec& >::type psis_inv(psis_invSEXP); + Rcpp::traits::input_parameter< double& >::type intercept_var(intercept_varSEXP); + update_intercept(intercept, Z, F, Lambda, psis_inv, intercept_var); + return R_NilValue; +END_RCPP +} +// min2LL_ono +double min2LL_ono(unsigned int N, unsigned int J, const arma::mat& Y, const arma::mat& MISS, const arma::mat& as, const arma::vec& bs, const arma::mat& theta, const arma::vec& Ms, const arma::mat& Kaps); +RcppExport SEXP _bayesefa_min2LL_ono(SEXP NSEXP, SEXP JSEXP, SEXP YSEXP, SEXP MISSSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP thetaSEXP, SEXP MsSEXP, SEXP KapsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< unsigned int >::type N(NSEXP); + Rcpp::traits::input_parameter< unsigned int >::type J(JSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type MISS(MISSSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type as(asSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type bs(bsSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type theta(thetaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type Ms(MsSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Kaps(KapsSEXP); + rcpp_result_gen = Rcpp::wrap(min2LL_ono(N, J, Y, MISS, as, bs, theta, Ms, Kaps)); + return rcpp_result_gen; +END_RCPP +} +// IFA_Mode_Jumper +Rcpp::List IFA_Mode_Jumper(const arma::mat& Y, unsigned int M, double gamma, const arma::vec& Ms, const arma::vec& sdMH, unsigned int burnin, unsigned int chain_length); +RcppExport SEXP _bayesefa_IFA_Mode_Jumper(SEXP YSEXP, SEXP MSEXP, SEXP gammaSEXP, SEXP MsSEXP, SEXP sdMHSEXP, SEXP burninSEXP, SEXP chain_lengthSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< unsigned int >::type M(MSEXP); + Rcpp::traits::input_parameter< double >::type gamma(gammaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type Ms(MsSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sdMH(sdMHSEXP); + Rcpp::traits::input_parameter< unsigned int >::type burnin(burninSEXP); + Rcpp::traits::input_parameter< unsigned int >::type chain_length(chain_lengthSEXP); + rcpp_result_gen = Rcpp::wrap(IFA_Mode_Jumper(Y, M, gamma, Ms, sdMH, burnin, chain_length)); + return rcpp_result_gen; +END_RCPP +} +// proposal2 +Rcpp::List proposal2(arma::uvec& new_r_idx, arma::mat& lambda, arma::mat& factors); +RcppExport SEXP _bayesefa_proposal2(SEXP new_r_idxSEXP, SEXP lambdaSEXP, SEXP factorsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::uvec& >::type new_r_idx(new_r_idxSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type lambda(lambdaSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type factors(factorsSEXP); + rcpp_result_gen = Rcpp::wrap(proposal2(new_r_idx, lambda, factors)); + return rcpp_result_gen; +END_RCPP +} +// mode_jump +Rcpp::List mode_jump(const arma::mat& X, arma::mat& lambda_mean, arma::mat& f_mean, arma::mat& invClam, const arma::vec sigma, arma::uvec r_idx, double my_gamma); +RcppExport SEXP _bayesefa_mode_jump(SEXP XSEXP, SEXP lambda_meanSEXP, SEXP f_meanSEXP, SEXP invClamSEXP, SEXP sigmaSEXP, SEXP r_idxSEXP, SEXP my_gammaSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type lambda_mean(lambda_meanSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type f_mean(f_meanSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type invClam(invClamSEXP); + Rcpp::traits::input_parameter< const arma::vec >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< arma::uvec >::type r_idx(r_idxSEXP); + Rcpp::traits::input_parameter< double >::type my_gamma(my_gammaSEXP); + rcpp_result_gen = Rcpp::wrap(mode_jump(X, lambda_mean, f_mean, invClam, sigma, r_idx, my_gamma)); + return rcpp_result_gen; +END_RCPP +} +// rmvnorm +arma::mat rmvnorm(unsigned int n, const arma::vec& mu, const arma::mat& S); +RcppExport SEXP _bayesefa_rmvnorm(SEXP nSEXP, SEXP muSEXP, SEXP SSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< unsigned int >::type n(nSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type mu(muSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type S(SSEXP); + rcpp_result_gen = Rcpp::wrap(rmvnorm(n, mu, S)); + return rcpp_result_gen; +END_RCPP +} +// sim_gamma_type +double sim_gamma_type(double x, double alph, double g, double d); +RcppExport SEXP _bayesefa_sim_gamma_type(SEXP xSEXP, SEXP alphSEXP, SEXP gSEXP, SEXP dSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< double >::type x(xSEXP); + Rcpp::traits::input_parameter< double >::type alph(alphSEXP); + Rcpp::traits::input_parameter< double >::type g(gSEXP); + Rcpp::traits::input_parameter< double >::type d(dSEXP); + rcpp_result_gen = Rcpp::wrap(sim_gamma_type(x, alph, g, d)); + return rcpp_result_gen; +END_RCPP +} +// update_F_matrix +void update_F_matrix(const arma::mat& Y, const arma::mat& I_K, arma::mat& F, arma::mat& Lambda, arma::vec& psis_inv); +RcppExport SEXP _bayesefa_update_F_matrix(SEXP YSEXP, SEXP I_KSEXP, SEXP FSEXP, SEXP LambdaSEXP, SEXP psis_invSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type I_K(I_KSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type F(FSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type Lambda(LambdaSEXP); + Rcpp::traits::input_parameter< arma::vec& >::type psis_inv(psis_invSEXP); + update_F_matrix(Y, I_K, F, Lambda, psis_inv); + return R_NilValue; +END_RCPP +} +// update_invClam +void update_invClam(const arma::mat& Lambda, arma::mat& invClam); +RcppExport SEXP _bayesefa_update_invClam(SEXP LambdaSEXP, SEXP invClamSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Lambda(LambdaSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type invClam(invClamSEXP); + update_invClam(Lambda, invClam); + return R_NilValue; +END_RCPP +} +// my_setdiff +arma::uvec my_setdiff(arma::uvec& x, const arma::uvec& y); +RcppExport SEXP _bayesefa_my_setdiff(SEXP xSEXP, SEXP ySEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::uvec& >::type x(xSEXP); + Rcpp::traits::input_parameter< const arma::uvec& >::type y(ySEXP); + rcpp_result_gen = Rcpp::wrap(my_setdiff(x, y)); + return rcpp_result_gen; +END_RCPP +} +// update_Lambda_loadings_hard_zero +void update_Lambda_loadings_hard_zero(const arma::mat& Y, arma::uvec& r_idx, arma::mat& F, arma::mat& Lambda, arma::vec& psis_inv, const arma::mat invClam, double my_gamma); +RcppExport SEXP _bayesefa_update_Lambda_loadings_hard_zero(SEXP YSEXP, SEXP r_idxSEXP, SEXP FSEXP, SEXP LambdaSEXP, SEXP psis_invSEXP, SEXP invClamSEXP, SEXP my_gammaSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< arma::uvec& >::type r_idx(r_idxSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type F(FSEXP); + Rcpp::traits::input_parameter< arma::mat& >::type Lambda(LambdaSEXP); + Rcpp::traits::input_parameter< arma::vec& >::type psis_inv(psis_invSEXP); + Rcpp::traits::input_parameter< const arma::mat >::type invClam(invClamSEXP); + Rcpp::traits::input_parameter< double >::type my_gamma(my_gammaSEXP); + update_Lambda_loadings_hard_zero(Y, r_idx, F, Lambda, psis_inv, invClam, my_gamma); + return R_NilValue; +END_RCPP +} +// kappa_initialize +arma::mat kappa_initialize(const arma::vec& Ms); +RcppExport SEXP _bayesefa_kappa_initialize(SEXP MsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::vec& >::type Ms(MsSEXP); + rcpp_result_gen = Rcpp::wrap(kappa_initialize(Ms)); + return rcpp_result_gen; +END_RCPP +} + +static const R_CallMethodDef CallEntries[] = { + {"_bayesefa_update_uniquenesses", (DL_FUNC) &_bayesefa_update_uniquenesses, 4}, + {"_bayesefa_EFA_Mode_Jumper", (DL_FUNC) &_bayesefa_EFA_Mode_Jumper, 5}, + {"_bayesefa_rTruncNorm_bounds", (DL_FUNC) &_bayesefa_rTruncNorm_bounds, 4}, + {"_bayesefa_update_uniquenesses_mixed", (DL_FUNC) &_bayesefa_update_uniquenesses_mixed, 6}, + {"_bayesefa_update_WKappaZ_NA_mixed", (DL_FUNC) &_bayesefa_update_WKappaZ_NA_mixed, 11}, + {"_bayesefa_update_intercept_mixed", (DL_FUNC) &_bayesefa_update_intercept_mixed, 6}, + {"_bayesefa_IFA_Mode_Jumper_MixedResponses", (DL_FUNC) &_bayesefa_IFA_Mode_Jumper_MixedResponses, 8}, + {"_bayesefa_sampleY_given_Z", (DL_FUNC) &_bayesefa_sampleY_given_Z, 4}, + {"_bayesefa_mixedresponse_posterior_prediction", (DL_FUNC) &_bayesefa_mixedresponse_posterior_prediction, 6}, + {"_bayesefa_rTruncNorm_lb", (DL_FUNC) &_bayesefa_rTruncNorm_lb, 3}, + {"_bayesefa_update_WKappaZ_NA", (DL_FUNC) &_bayesefa_update_WKappaZ_NA, 9}, + {"_bayesefa_update_intercept", (DL_FUNC) &_bayesefa_update_intercept, 6}, + {"_bayesefa_min2LL_ono", (DL_FUNC) &_bayesefa_min2LL_ono, 9}, + {"_bayesefa_IFA_Mode_Jumper", (DL_FUNC) &_bayesefa_IFA_Mode_Jumper, 7}, + {"_bayesefa_proposal2", (DL_FUNC) &_bayesefa_proposal2, 3}, + {"_bayesefa_mode_jump", (DL_FUNC) &_bayesefa_mode_jump, 7}, + {"_bayesefa_rmvnorm", (DL_FUNC) &_bayesefa_rmvnorm, 3}, + {"_bayesefa_sim_gamma_type", (DL_FUNC) &_bayesefa_sim_gamma_type, 4}, + {"_bayesefa_update_F_matrix", (DL_FUNC) &_bayesefa_update_F_matrix, 5}, + {"_bayesefa_update_invClam", (DL_FUNC) &_bayesefa_update_invClam, 2}, + {"_bayesefa_my_setdiff", (DL_FUNC) &_bayesefa_my_setdiff, 2}, + {"_bayesefa_update_Lambda_loadings_hard_zero", (DL_FUNC) &_bayesefa_update_Lambda_loadings_hard_zero, 7}, + {"_bayesefa_kappa_initialize", (DL_FUNC) &_bayesefa_kappa_initialize, 1}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_bayesefa(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/commonEFAfunctions.cpp b/src/commonEFAfunctions.cpp new file mode 100644 index 0000000..0bf0c6d --- /dev/null +++ b/src/commonEFAfunctions.cpp @@ -0,0 +1,361 @@ +#include +#include "commonEFAfunctions.h" +using namespace Rcpp; + +// [[Rcpp::depends(RcppArmadillo)]] + + +Rcpp::List proposal2(arma::uvec& new_r_idx,arma::mat& lambda_mean, arma::mat& f_mean); + +Rcpp::List mode_jump(const arma::mat& X, arma::mat& lambda_mean, arma::mat& f_mean, + arma::mat& invClam, const arma::vec sigma, arma::uvec r_idx, double my_gamma); + +arma::mat rmvnorm(unsigned int n, const arma::vec& mu, const arma::mat& S); + +double sim_gamma_type(double x,double alph,double g,double d); + +void update_F_matrix(const arma::mat& Y,const arma::mat& I_K,arma::mat& F,arma::mat& Lambda, + arma::vec& psis_inv); + +void update_invClam(const arma::mat& Lambda,arma::mat& invClam); + +arma::uvec my_setdiff(arma::uvec& x, const arma::uvec& y); + +void update_Lambda_loadings_hard_zero(const arma::mat& Y,arma::uvec& r_idx,arma::mat& F, + arma::mat& Lambda,arma::vec& psis_inv, + const arma::mat invClam,double my_gamma); + +arma::mat kappa_initialize(const arma::vec& Ms); + +//' @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 +// [[Rcpp::export]] +Rcpp::List proposal2(arma::uvec& new_r_idx, // Input in R format, aka counting starts from 1 not 0 + arma::mat& lambda, + arma::mat& factors) { + arma::mat Q, R; + //extract lambda_r + arma::mat lambda_sub = lambda.rows(new_r_idx - 1); + + //compute QR decomposition + qr(Q, R, lambda_sub.t()); + arma::vec sgn = sign(R.diag()); + //compute Q for rotating to new row index configuration + Q = Q * diagmat(sgn); + arma::mat lambda_new = lambda * Q; + arma::mat f_new = factors * Q; + + return Rcpp::List::create(Rcpp::Named("lambda_new",lambda_new), + Rcpp::Named("f_new",f_new)); +} + + +//' @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 +// [[Rcpp::export]] +Rcpp::List mode_jump(const arma::mat& X, + arma::mat& lambda_mean, + arma::mat& f_mean, + arma::mat& invClam, + const arma::vec sigma, + arma::uvec r_idx, // Input in Rcpp format; index starts at 0 not 1 + double my_gamma) { + int J = lambda_mean.n_rows; + int M = lambda_mean.n_cols; + bool accepted = false; + + //generate and sample proposed row indices for PLT configuration + arma::uvec new_r_idx; + arma::uvec one_to_J = arma::shuffle(arma::linspace(1, J, J)); + new_r_idx = sort(one_to_J.subvec(0, M-1)); + + //rotate current lambda and F to proposed configuration + Rcpp::List forward_prop = proposal2(new_r_idx, lambda_mean, f_mean); + + //computations for current state + arma::mat lambda_sub = lambda_mean.rows(r_idx); + double log_like_old = 0; + // Lambda prior variance + //compute the log determinant for current state + log_like_old += accu(log(pow(lambda_sub.diag(), my_gamma))); + + //computations for proposed rotation + //extract rotated lambda and F + arma::mat lambda_proposed = forward_prop["lambda_new"]; + arma::mat F_proposed = forward_prop["f_new"]; + + //compute logged determinant + arma::mat lambda_sub_proposed = lambda_proposed.rows(new_r_idx-1); + double log_like_prop = 0; + log_like_prop += accu(log(pow(lambda_sub_proposed.diag(), my_gamma))); + + // double marginal_like = NULL; + // double marginal_like_prop = NULL; + + //compute acceptance probability + // double A = std::min(1.0, exp(log_like_prop - log_like_old)); + // if (!std::isfinite(A)) { + // std::cout << "NA acceptance probability\n"; + // return Rcpp::List::create(Rcpp::Named("lambda",lambda_mean), + // Rcpp::Named("r_idx",r_idx), + // Rcpp::Named("f_mat",f_mean), + // Rcpp::Named("accepted",accepted), + // Rcpp::Named("marginal_like",log_like_old)); + // } + + double u = arma::randu(); + + if (log(u) < log_like_prop - log_like_old) { + // if (u < A) { + arma::mat temp = forward_prop["f_new"]; + f_mean = temp; + lambda_mean = lambda_proposed; + r_idx = new_r_idx - 1; + accepted = true; + } + return Rcpp::List::create(Rcpp::Named("lambda",lambda_mean), + Rcpp::Named("r_idx",r_idx), + Rcpp::Named("f_mat",f_mean), + Rcpp::Named("accepted",accepted)); +} + +//' @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)) +//' +// [[Rcpp::export]] +arma::mat rmvnorm(unsigned int n, const arma::vec& mu, const arma::mat& S) { + unsigned int ncols = S.n_cols; + arma::mat Y(n, ncols); + Y.imbue( norm_rand ) ; + return arma::repmat(mu, 1, n).t() + Y * arma::chol(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 +//' +// [[Rcpp::export]] +double sim_gamma_type(double x,double alph,double g,double d){ + //update u given x + double v1 = R::runif(0.0,1.0); + double u = v1*exp(-pow(x-g,2)/d); + //update x given u + double v2 = R::runif(0.0,1.0); + arma::vec acand(2); + acand(0) = 0.; + acand(1) = g - sqrt(-1.*d*log(u)); + double a = arma::max(acand); + double b = g + sqrt(-1.*d*log(u)); + double x1 = pow(v2*pow(b,alph)+(1-v2)*pow(a,alph),1./alph); + return x1; +} + +//update F +//this could be updated to sample F within an i loop, which may be advantageous for ordinal samples +//' @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 +// [[Rcpp::export]] +void update_F_matrix(const arma::mat& Y, + const arma::mat& I_K, + arma::mat& F, + arma::mat& Lambda, + arma::vec& psis_inv){ + unsigned int N = F.n_rows; + unsigned int M = F.n_cols; + arma::mat Psi_inv = arma::diagmat(psis_inv); + arma::mat Psi_inv_Lam = Psi_inv*Lambda; + arma::mat VF = inv(I_K + Lambda.t()*Psi_inv_Lam); + arma::mat mF = Y*Psi_inv_Lam*VF; + F = mF + arma::randn(N,M)*arma::chol(VF); + +} + + +//update precision for elements of lambda +//note there is just single element for the prior +//' @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 +// [[Rcpp::export]] +void update_invClam(const arma::mat& Lambda, + arma::mat& invClam) { + // No spike-and-slab, only column-wise variance + unsigned int J = Lambda.n_rows; + int M = Lambda.n_cols; + double u = R::runif(0.0,1.0); + double d = accu(pow(Lambda, 2)); + invClam.fill(R::qgamma(u,double(J*M-((M-1)*M/2) + 2)/2.0,2.0/(2.0 + d),1,0)); +} + + +//' @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 +//' +// [[Rcpp::export]] +arma::uvec my_setdiff(arma::uvec& x, const arma::uvec& y){ + for (size_t j = 0; j < y.n_elem; j++) { + arma::uword q1 = arma::conv_to::from(arma::find(x == y(j))); + x.shed_row(q1); + } + return x; +} + +//update loadings given Y, row indices, F, and uniquenesses +//' @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 +// [[Rcpp::export]] +void update_Lambda_loadings_hard_zero(const arma::mat& Y, + arma::uvec& r_idx, + arma::mat& F, + arma::mat& Lambda, + arma::vec& psis_inv, + const arma::mat invClam, + double my_gamma){ + unsigned int J = Lambda.n_rows; + unsigned int M = F.n_cols; + //update Lambda + // Lambda.fill(0); + arma::mat FpF = F.t()*F; + arma::mat FpY = F.t()*Y; + arma::rowvec lambdaj(M); + + // Constrained case + unsigned int j; + for(unsigned int m=0;m 0) { + arma::vec m1 = m_1toj(arma::span(0,m-1)); + arma::mat S1 = C_1toj(arma::span(0,m-1),arma::span(0,m-1)); + arma::vec S12= C_1toj(arma::span(0,m-1),m); + arma::vec mb = m1+S12*(lambdajj-m_1toj(m))/C_1toj(m,m); + arma::mat Sb = S1-S12*S12.t()/C_1toj(m,m); + arma::rowvec lambda_1tjm1 = rmvnorm(1,mb,Sb); + Lambda(j,arma::span(0,m-1))=lambda_1tjm1; + } + Lambda(j,m) = lambdajj; + + //fill remaining elements in PLT rows with zeros + if ((m+1) < M) { + arma::rowvec fill_zeros = arma::zeros(M-(m+1)); + Lambda(j,arma::span(m+1,M-1)) = fill_zeros; + } + } + + // Unconstrained case + arma::uvec one_to_J = arma::linspace(0, J-1, J); + arma::uvec not_r_idx = my_setdiff(one_to_J, r_idx);//find indices not in r + unsigned int I = not_r_idx.n_elem; + for(unsigned int i=0; iMs(j)){ + KAP0(m,j)=arma::datum::nan; + } + } + } + return KAP0; +} +