Bugfix: Remove Bitwise Shift and use Xi (#2)

* 1. Switch from using a bitwise shift to using the correct comparison
2. Use the subset variant of X for the subject.

* Ensure cIRT is exported.

* Code formatting in the vignette

* Restrict text to 80 characters in vignette.

* Add a bullet point to NEWS.md

* CITATION spacing [skip-ci]

NEWS.md

@@ -6,6 +6,10 @@
 - Enabled C++11 and OpenMP.
 - Switched to allowing Rcpp to handle native registration
 
+## Bugfix
+
+- Addressed issues in the choice generation procedure.
+
 ## Documentation
 
 - Improved in-line documentation.

R/RcppExports.R

@@ -331,6 +331,7 @@ probitHLM <- function(unique_subject_ids, subject_ids, choices_nk, fixed_effects
 #' 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)
 }

inst/CITATION

@@ -6,7 +6,7 @@ citEntry(entry = "Manual",
                                    as.person("James Balamuta")
                                   ),
          year         = 2015,
-         textVersion  = paste("Culpepper, S.A. and Balamuta, J.J (2015)",
+         textVersion  = paste("Culpepper, S. A. and Balamuta, J. J. (2015)",
                               "cIRT: Choice Item Response Theory.",
                               "URL https://cran.r-project.org/package=cIRT.")
 )

src/project_source.cpp

@@ -629,6 +629,7 @@ probitHLM(const arma::vec &unique_subject_ids, const arma::vec &subject_ids,
 //' apply(out_model_thurstone$as, 2, FUN = mean)
 //' apply(out_model_thurstone$bs, 2, FUN = mean)
 //' }
+//' @export
 // [[Rcpp::export]]
 Rcpp::List cIRT(arma::vec subject_ids, arma::mat fixed_effects,
                 arma::uvec B_elem_plus1, arma::mat rv_effects,

src/total_cpp.cpp

@@ -90,13 +90,14 @@ Rcpp::List Generate_Choice(unsigned int N, unsigned int J, unsigned int K,
     arma::vec Wzeta_0(nk);
     arma::mat etaC_mat(N, K);
 
-    for (unsigned int i = 0; i << N; i++) {
+    for (unsigned int i = 0; i < N; i++) {
         arma::uvec current_subject = find(subject_ids == unique_subject_ids(i));
         arma::mat Wi = W.rows(current_subject);
         arma::mat Xi = X.rows(current_subject);
         Wzeta_0.elem(current_subject) = Wi * trans(zeta.row(i)); // nk x 1
-        etaC_mat.row(i) = gamma.t() * X.t() + zeta.row(i) * Wi.t();
+        etaC_mat.row(i) = gamma.t() * Xi.t() + zeta.row(i) * Wi.t();
     }
+
     arma::vec etaC = X * gamma + Wzeta_0;
     Cs.elem(find(etaC_mat > ZC)).ones();
 

vignettes/Estimating-the-Model-in-the-Paper.Rmd

@@ -10,17 +10,22 @@ vignette: >
 
 # Intro
 
-The objective of this vignette is to provide clarity as to the estimation procedure used to obtain the results in the paper. This vignette breaks down each section in the example displayed in `cIRT()` function. 
+The objective of this vignette is to provide clarity as to the estimation 
+procedure used to obtain the results in the paper. This vignette breaks down 
+each section in the example displayed in `cIRT()` function. 
 
 # Data
 
-Within this vignette, we used two different data sets to generate the model that are located within `cIRT`.
+Within this vignette, we used two different data sets to generate a Choice
+Item Response Theory Model with the rountine located in `cIRT`.
 
 ```{r load_pkg}
 library(cIRT)
 ```
 
-The first data set, `trial_matrix`, contains whether or not the subject correctly identified the spatial rotation. The second dataset, `choice_matrix`, provides information regarding the choice decision subjects were asked to make. 
+The first data set, `trial_matrix`, contains whether or not the subject
+correctly identified the spatial rotation. The second dataset, `choice_matrix`,
+provides information regarding the choice decision subjects were asked to make. 
 
 ```{r load_data}
 data(trial_matrix)
@@ -35,27 +40,28 @@ Here we construct a thurstone design matrix by obtaining the IDs of hard and eas
 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)]
+D_easy = model.matrix( ~ -1 + factor(easy_items))
+D_hard = -1 * model.matrix( ~ -1 + factor(hard_items))[, -c(5, 10, 15)]
 ```
 
 Within this setting, we setup the effect-codes for different constraints. 
+
 ```{r effect_coding}
 # Defining effect-coded contrasts
-high_contrasts <- rbind(-1,diag(4))
+high_contrasts = rbind(-1, diag(4))
 rownames(high_contrasts) = 12:16
-low_contrasts <- rbind(-1,diag(2))
+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'])
+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)
+fixed_effects = model.matrix( ~ high + low)
+fixed_effects_base = fixed_effects[, 1]
+fixed_effects_int = model.matrix( ~ high * low)
 ```
 
 # Modeling the Data
@@ -64,19 +70,22 @@ Generate the cIRT model using a Thurstone Design Matrix generated above.
 
 ```{r model_data}
 # 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)
+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
+  )
 })
 ```
 
-We recommend saving the model object as a `.rda` file even though the total computational time is less than 2.5 minutes.  
+We recommend saving the model object as a `.rda` file even though the total
+computational time is less than 2.5 minutes.  
 ```
 ## Save model output to an rda file.
 # save(out_model_thurstone, file='choiceMCMCoutput.rda')
@@ -90,15 +99,33 @@ We recommend saving the model object as a `.rda` file even though the total comp
 Next up, we obtain the parameter estimates of the model by averaging over the different estimates obtained via the Gibbs sampling technique employed. 
 
 ```{r param_ests}
-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)))))
-corrmat = inv_sd%*%S_thurstone%*%inv_sd
+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)
+))))
+
+corrmat = inv_sd %*% S_thurstone %*% inv_sd
 as = apply(out_model_thurstone$as, 2, FUN = mean)
 bs = apply(out_model_thurstone$bs, 2, FUN = mean)
 ```
@@ -108,18 +135,15 @@ Thus, we have the following results:
 ```{r param_results}
 # gs0
 G_thurstone
-
 # betas
 B_thurstone
-
 # Sigma Thurstone
 S_thurstone
 
-# item parameters
+## Item parameters ----
 
 # a
 as
-
 # b
 bs
 ```