R/RcppExports.R

Defines functions Generate_Choice Total_Tabulate riwishart rwishart rmvnorm cIRT probitHLM TwoPLChoicemcmc center_matrix direct_sum

Documented in center_matrix cIRT direct_sum Generate_Choice probitHLM riwishart rmvnorm rwishart Total_Tabulate TwoPLChoicemcmc

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

#' Direct Sum of Matrices
#'
#' Computes the direct sum of all matrices passed in via the list.
#' 
#' @param x A `field<matrix>` or `list` containing matrices
#' 
#' @return 
#' Matrix containing the direct sum of all matrices in the list.
#' 
#' @author 
#' James Joseph Balamuta
#' 
#' @details
#' Consider matrix \eqn{A} (\eqn{M \times N}{M x N}) and
#' \eqn{B} (\eqn{K \times P}{k x p}). A direct sum is a diagonal matrix 
#' \eqn{A (+) B} with dimensions \eqn{(m + k) x (n + p)}.
#' 
#' @export
#' @examples
#' 
#' x = list(matrix(0, nrow = 5, ncol = 3),
#'          matrix(1, nrow = 5, ncol = 3))
#' direct_sum(x)
#'
#' x = list(matrix(rnorm(15), nrow = 5, ncol = 3),
#'          matrix(rnorm(30), nrow = 5, ncol = 6),
#'          matrix(rnorm(18), nrow = 2, ncol = 9))
#' direct_sum(x)
direct_sum <- function(x) {
    .Call(`_cIRT_direct_sum`, x)
}

#' Center a Matrix
#'
#' Obtains the mean of each column of the matrix and subtracts it from the
#' given matrix in a centering operation.
#' 
#' @param x A `matrix` with any dimensions
#' 
#' @return
#' A `matrix` with the same dimensions of X that has been centered. 
#' 
#' @author 
#' James Joseph Balamuta
#' 
#' @details 
#' The application of this function to a matrix mimics the use of a 
#' centering matrix given by:
#' 
#' \deqn{{C_n} = {I_n} - \frac{1}{n}{11^T}} 
#' 
#' @seealso
#' [cIRT()]
#' 
#' @export
#' @examples
#' nobs = 500
#' nvars = 20
#' x = matrix(rnorm(nobs * nvars), nrow = nobs, ncol = nvars) 
#' r_centered = scale(x) 
#' arma_centered1 = center_matrix(x)
center_matrix <- function(x) {
    .Call(`_cIRT_center_matrix`, x)
}

#' Two Parameter Choice IRT Model MCMC
#'
#' Performs an MCMC routine for a two parameter IRT Model using Choice Data
#'
#' @param unique_subject_ids A `vector` with length \eqn{N \times 1}{N x 1}
#'                           containing unique subject IDs.
#' @param subject_ids        A `vector` with length \eqn{NK \times 1}{N*K x 1}
#'                           containing subject IDs.
#' @param choices_nk         A `vector` with length \eqn{NK \times 1}{N*K x 1}
#'                           containing subject choices.
#' @param fixed_effects      A `matrix` with dimensions
#'                           \eqn{NK \times P_1}{N*K x P_1} containing
#'                           fixed effect design matrix without theta.
#' @param B                  A \eqn{V} dimensional column `vector`
#'                           relating \eqn{\theta_i}{theta_i} and
#'                           \eqn{\zeta_i}{zeta_i}.
#' @param rv_effects_design  A `matrix` with dimensions
#'                           \eqn{NK \times V}{N*K x V} containing
#'                           random effect variables.
#' @param gamma              A `vector` with dimensions \eqn{P \times 1}{P x 1}
#'                           containing fixed parameter estimates,
#'                           where \eqn{P = P_1 + P_2}
#' @param beta               A `vector` with dimensions \eqn{P_2}
#'                           containing random parameter estimates.
#' @param zeta_rv            A `matrix` with dimensions \eqn{N \times V}{N x V}
#'                           containing random parameter estimates.
#' @param Sigma_zeta_inv     A `matrix` with dimensions
#'                           \eqn{P_2 \times P_2}{P_2 x P_2}.
#' @param Y                  A `matrix` of dimensions \eqn{N \times J}{N x J}
#'                           for Dichotomous item responses
#' @param theta0             A `vector` of length \eqn{N \times 1}{N x 1}
#'                           for latent theta.
#' @param a0                 A `vector` of length \eqn{J}
#'                           for item discriminations.
#' @param b0                 A `vector` of length \eqn{J}
#'                           for item locations.
#' @param mu_xi0             A `vector` of dimension 2 (i.e. c(0,1)) that is a
#'                           prior for item parameter means.
#' @param Sig_xi0            A `matrix` of dimension 2x2 (i.e. diag(2)) that
#'                           is a prior for item parameter vc matrix.
#'
#' @return
#' A `list` that contains:
#' \describe{
#'   \item{\code{ai1}}{A `vector` of length J}
#'   \item{\code{bi1}}{A `vector` of length J}
#'   \item{\code{theta1}}{A `vector` of length N}
#'   \item{\code{Z_c}}{A `matrix` of length NK}
#'   \item{\code{Wzeta_0}}{A `matrix` of length NK}
#' }
#'
#' @author
#' Steven Andrew Culpepper and James Joseph Balamuta
#'
#' @seealso [cIRT()], [rmvnorm()], and [riwishart()]
#' @export
#' @examples \dontrun{
#' # Call with the following data:
#' TwoPLChoicemcmc(cogDAT, theta0, a0, b0, mu_xi0, Sig_xi0)
#' }
TwoPLChoicemcmc <- function(unique_subject_ids, subject_ids, choices_nk, fixed_effects, B, rv_effects_design, gamma, beta, zeta_rv, Sigma_zeta_inv, Y, theta0, a0, b0, mu_xi0, Sig_xi0) {
    .Call(`_cIRT_TwoPLChoicemcmc`, unique_subject_ids, subject_ids, choices_nk, fixed_effects, B, rv_effects_design, gamma, beta, zeta_rv, Sigma_zeta_inv, Y, theta0, a0, b0, mu_xi0, Sig_xi0)
}

#' Probit Hierarchical Level Model
#'
#' Performs modeling procedure for a Probit Hierarchical Level Model.
#'
#' @param unique_subject_ids   A `vector` with length N x 1 containing
#'                             unique subject IDs.
#' @param subject_ids          A `vector` with length N*K x 1 containing
#'                             subject IDs.
#' @param choices_nk           A `vector` with length N*K x 1 containing
#'                             subject choices.
#' @param fixed_effects_design A `matrix` with dimensions N*K x P containing
#'                             fixed effect variables.
#' @param rv_effects_design    A `matrix` with dimensions N*K x V containing
#'                             random effect variables.
#' @param B_elem_plus1         A `V[[1]]` dimensional column `vector`
#'                             indicating which zeta_i relate to theta_i.
#' @param gamma                A `vector` with dimensions P_1 x 1 containing
#'                             fixed parameter estimates.
#' @param beta                 A `vector` with dimensions P_2 x 1 containing
#'                             random parameter estimates.
#' @param theta                A `vector` with dimensions N x 1 containing
#'                             subject understanding estimates.
#' @param zeta_rv              A `matrix` with dimensions N x V containing
#'                             random parameter estimates.
#' @param WtW                  A `field<matrix>` P x P x N contains the
#'                             caching for direct sum.
#' @param Z_c                  A `vector` with dimensions N*K x 1
#' @param Wzeta_0              A `vector` with dimensions N*K x 1
#' @param inv_Sigma_gamma      A `matrix` with dimensions P x P that is the
#'                             prior inverse sigma matrix for gamma.
#' @param mu_gamma             A `vector` with length P x 1 that is the prior
#'                             mean vector for gamma.
#' @param Sigma_zeta_inv       A `matrix` with dimensions V x V that is the
#'                             prior inverse sigma matrix for zeta.
#' @param S0                   A `matrix` with dimensions V x V that is the
#'                             prior sigma matrix for zeta.
#' @param mu_beta              A `vector` with dimensions P_2 x 1, that is
#'                             the mean of beta.
#' @param sigma_beta_inv       A `matrix` with dimensions P_2 x P_2, that is
#'                             the inverse sigma matrix of beta.
#' @return
#' A `list` that contains:
#' \describe{
#'   \item{\code{zeta_1}}{A `vector` of length N}
#'   \item{\code{sigma_zeta_inv_1}}{A `matrix` of dimensions V x V}
#'   \item{\code{gamma_1}}{A `vector` of length P}
#'   \item{\code{beta_1}}{A `vector` of length V}
#'   \item{\code{B}}{A `matrix` of length V}
#' }
#' 
#' @author
#' Steven Andrew Culpepper and James Joseph Balamuta
#'
#' @details
#' The function is implemented to decrease the amount of vectorizations
#' necessary.
#'
#' @seealso
#' [rwishart()] and [TwoPLChoicemcmc()]
#'
#' @export
probitHLM <- function(unique_subject_ids, subject_ids, choices_nk, fixed_effects_design, rv_effects_design, B_elem_plus1, gamma, beta, theta, zeta_rv, WtW, Z_c, Wzeta_0, inv_Sigma_gamma, mu_gamma, Sigma_zeta_inv, S0, mu_beta, sigma_beta_inv) {
    .Call(`_cIRT_probitHLM`, unique_subject_ids, subject_ids, choices_nk, fixed_effects_design, rv_effects_design, B_elem_plus1, gamma, beta, theta, zeta_rv, WtW, Z_c, Wzeta_0, inv_Sigma_gamma, mu_gamma, Sigma_zeta_inv, S0, mu_beta, sigma_beta_inv)
}

#' Generic Implementation of Choice IRT MCMC
#'
#' Builds a model using MCMC
#'
#' @param subject_ids    A `vector` that contains subject IDs for each line of
#'                       data in the choice vector (e.g. For 1 subject that
#'                       made 5 choices, we would have the number 1 appear
#'                       five times consecutively.)
#' @param fixed_effects  A `matrix` with NK x P1 dimensions that acts as the
#'                       design matrix for terms WITHOUT theta.
#' @param B_elem_plus1   A `V[[1]]` dimensional column `vector` indicating
#'                       which zeta_i relate to theta_i.
#' @param rv_effects     A `matrix` with NK x V dimensions for random effects
#'                       design matrix.
#' @param trial_matrix   A `matrix` with N x J dimensions, where J denotes the
#'                       number of items presented. The matrix MUST contain
#'                       only 1's and 0's.
#' @param choices_nk     A `vector` with NK length that contains the choice
#'                       value e.g. 0 or 1.
#' @param chain_length   An `int` that controls how many MCMC draws there are.
#'                       (> 0)
#' @param burnit         An `int` that describes how many MCMC draws should be
#'                       discarded.
#'
#' @return
#' A `list` that contains:
#'
#' \describe{
#'   \item{\code{as}}{A `matrix` of dimension chain_length x J}
#'   \item{\code{bs}}{A `matrix` of dimension chain_length x J}
#'   \item{\code{gs}}{A `matrix` of dimension chain_length x P_1}
#'   \item{\code{Sigma_zeta_inv}}{An `array` of dimension V x V x chain_length}
#'   \item{\code{betas}}{A `matrix` of dimension chain_length x P_2}
#' }
#'
#' @author
#' Steven Andrew Culpepper and James Joseph Balamuta
#'
#' @seealso
#' [TwoPLChoicemcmc()], [probitHLM()], [center_matrix()],
#' [rmvnorm()], [rwishart()], and [riwishart()]
#'
#' @export
#' @examples
#' \dontrun{
#' # Variables
#' # Y = trial matix
#' # C = KN vector of binary choices
#' # N = #of subjects
#' # J = # of items
#' # K= # of choices
#' # atrue = true item discriminations
#' # btrue = true item locations
#' # thetatrue = true thetas/latent performance
#' # gamma = fixed effects coefficients
#' # Sig = random-effects variance-covariance
#' # subid = id variable for subjects
#' 
#' # Load the Package
#' library(cIRT)
#' 
#' # Load the Data
#' data(trial_matrix)
#' data(choice_matrix)
#'
#' # Thurstone design matrices
#' all_nopractice = subset(all_data_trials, experiment_loop.thisN > -1)
#' hard_items = choice_matrix$hard_q_id
#' easy_items = choice_matrix$easy_q_id
#'
#' D_easy = model.matrix( ~ -1 + factor(easy_items))
#' D_hard = -1 * model.matrix( ~ -1 + factor(hard_items))[, -c(5, 10, 15)]
#'
#' # Defining effect-coded contrasts
#' high_contrasts = rbind(-1, diag(4))
#' rownames(high_contrasts) = 12:16
#' low_contrasts = rbind(-1, diag(2))
#' rownames(low_contrasts) = 4:6
#'
#' # Creating high & low factors
#' high = factor(choice_matrix[, 'high_value'])
#' low = factor(choice_matrix[, 'low_value'])
#' contrasts(high) = high_contrasts
#' contrasts(low) = low_contrasts
#'
#' fixed_effects = model.matrix( ~ high + low)
#' fixed_effects_base = fixed_effects[, 1]
#' fixed_effects_int = model.matrix( ~ high * low)
#'
#'
#' # Model with Thurstone D Matrix
#' system.time({
#'  out_model_thurstone = cIRT(
#'    choice_matrix[, 'subject_id'],
#'    cbind(fixed_effects[, -1], D_easy, D_hard),
#'    c(1:ncol(fixed_effects)),
#'    as.matrix(fixed_effects),
#'    as.matrix(trial_matrix),
#'    choice_matrix[, 'choose_hard_q'],
#'    20000,
#'    25000
#'  )
#' })
#'
#'
#' vlabels_thurstone = colnames(cbind(fixed_effects[, -1], D_easy, D_hard))
#' G_thurstone = t(apply(
#'  out_model_thurstone$gs0,
#'  2,
#'  FUN = quantile,
#'  probs = c(.5, .025, .975)
#' ))
#'
#' rownames(G_thurstone) = vlabels_thurstone
#' B_thurstone = t(apply(
#'  out_model_thurstone$beta,
#'  2,
#'  FUN = quantile,
#'  probs = c(.5, 0.025, .975)
#' ))
#'
#' rownames(B_thurstone) = colnames(fixed_effects)
#'
#' S_thurstone = solve(
#'   apply(out_model_thurstone$Sigma_zeta_inv, c(1, 2), FUN = mean)
#' )
#'
#' inv_sd = diag(1 / sqrt(diag(solve(
#'  apply(out_model_thurstone$Sigma_zeta_inv, c(1, 2),
#'        FUN = mean)
#' ))))
#'
#' inv_sd %*% S_thurstone %*% inv_sd
#' apply(out_model_thurstone$as, 2, FUN = mean)
#' apply(out_model_thurstone$bs, 2, FUN = mean)
#' }
#' @export
cIRT <- function(subject_ids, fixed_effects, B_elem_plus1, rv_effects, trial_matrix, choices_nk, burnit, chain_length = 10000L) {
    .Call(`_cIRT_cIRT`, subject_ids, fixed_effects, B_elem_plus1, rv_effects, trial_matrix, choices_nk, burnit, chain_length)
}

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

#' Generate Random Wishart Distribution
#'
#' Creates a random wishart distribution when given degrees of freedom
#' and a sigma matrix. 
#' 
#' @param df An `integer`, which gives the degrees of freedom of the Wishart.
#'           (> 0) 
#' @param S  A `matrix` with dimensions m x m that provides Sigma, the
#'           covariance matrix. 
#'           
#' @return 
#' A `matrix` that is a Wishart distribution, aka the sample covariance
#' matrix of a Multivariate Normal Distribution 
#' 
#' @author 
#' James Joseph Balamuta
#'  
#' @seealso 
#' [riwishart()] and [probitHLM()]
#' 
#' @export 
#' @examples 
#' # Call with the following data:
#' rwishart(3, diag(2))
#'
#' # Validation
#' set.seed(1337)
#' S = toeplitz((10:1)/10)
#' n = 10000
#' o = array(dim = c(10,10,n))
#' for(i in 1:n){
#' o[,,i] = rwishart(20, S)
#' }
#' mR = apply(o, 1:2, mean)
#' Va = 20*(S^2 + tcrossprod(diag(S)))
#' vR = apply(o, 1:2, var)
#' stopifnot(all.equal(vR, Va, tolerance = 1/16))
rwishart <- function(df, S) {
    .Call(`_cIRT_rwishart`, df, S)
}

#' Generate Random Inverse Wishart Distribution
#'
#' Creates a random inverse wishart distribution when given
#' degrees of freedom and a sigma matrix. 
#' 
#' @param df An `integer` that represents the degrees of freedom. (> 0) 
#' @param S  A `matrix` with dimensions m x m that provides Sigma, 
#'           the covariance matrix. 
#' 
#' @return 
#' A `matrix` that is an inverse wishart distribution. 
#' 
#' @author
#' James Joseph Balamuta 
#'  
#' @seealso 
#' [rwishart()] and [TwoPLChoicemcmc()]
#' 
#' @export 
#' @examples
#' #Call with the following data:
#' riwishart(3, diag(2))
riwishart <- function(df, S) {
    .Call(`_cIRT_riwishart`, df, S)
}

#' Calculate Tabulated Total Scores
#'
#' Internal function to -2LL
#' 
#' @param N An `integer`, which gives the number of observations.  (> 0)
#' @param J An `integer`, which gives the number of items.  (> 0)
#' @param Y A N by J `matrix` of item responses.
#' 
#' @return 
#' A vector of tabulated total scores.
#' 
#' @author 
#' Steven Andrew Culpepper
#' 
#' @export
Total_Tabulate <- function(N, J, Y) {
    .Call(`_cIRT_Total_Tabulate`, N, J, Y)
}

#' Generate Observed Data from choice model
#'
#' Generates observed cognitive and choice data from the IRT-Thurstone model.
#' 
#' @param N                  An `integer` for the number of observations.
#' @param J                  An `integer` for the number of items.
#' @param K                  An `integer` for the number of paired comparisons.
#' @param theta              A `vector` of latent cognitive variables.
#' @param as                 A `vector` of length J with item discriminations.
#' @param bs                 A `vector` of length J with item locations.
#' @param zeta               A `matrix` with dimensions N x V containing
#'                           random parameter estimates. 
#' @param gamma              A `vector` with dimensions P x 1 containing
#'                           fixed parameter estimates, where \eqn{P = P_1 + P_2} 
#' @param X                  A `matrix` with dimensions N*K x P_1 containing 
#'                           fixed effect design matrix without theta. 
#' @param W                  A `matrix` with dimensions N*K x V containing
#'                           random effect variables. 
#' @param subject_ids        A `vector` with length NK x 1 containing
#'                           subject-choice IDs. 
#' @param unique_subject_ids A `vector` with length N x 1 containing
#'                           unique subject IDs. 
#' 
#' @return 
#' A \code{list} that contains: 
#' 
#' \describe{
#'    \item{\code{Y}}{A `matrix` of dimension N by J}
#'    \item{\code{C}}{A `vector` of length NK}
#' } 
#' 
#' @author 
#' Steven Andrew Culpepper and James Joseph Balamuta 
#' 
#' @export
Generate_Choice <- function(N, J, K, theta, as, bs, zeta, gamma, X, W, subject_ids, unique_subject_ids) {
    .Call(`_cIRT_Generate_Choice`, N, J, K, theta, as, bs, zeta, gamma, X, W, subject_ids, unique_subject_ids)
}

Try the cIRT package in your browser

Any scripts or data that you put into this service are public.

cIRT documentation built on March 18, 2022, 7:38 p.m.