From f5984931e96efe46826ea648e64c1864638ef2ff Mon Sep 17 00:00:00 2001 From: James J Balamuta Date: Tue, 22 Jan 2019 15:24:54 -0600 Subject: [PATCH] Update fourPNO (#1) fourPNO: 1.0.5 modernizing R code --- .Rbuildignore | 1 + .travis.yml | 25 +- DESCRIPTION | 20 +- NAMESPACE | 8 +- NEWS.md | 17 +- R/RcppExports.R | 548 +++++++----- R/fourPNO-internal.R | 130 +-- README.Rmd | 61 +- README.md | 71 +- cran-comments.md | 11 +- inst/CITATION | 18 +- man/Gibbs_2PNO.Rd | 12 +- man/Gibbs_4PNO.Rd | 45 +- man/Total_Tabulate.Rd | 6 +- man/Y_4pno_simulate.Rd | 6 +- man/fourPNO-package.Rd | 8 +- man/kappa_initialize.Rd | 23 - man/min2LL_4pno.Rd | 2 +- man/rmvnorm.Rd | 11 +- man/update_2pno.Rd | 46 - man/update_WKappaZ_NA.Rd | 50 -- man/update_ab_NA.Rd | 37 - man/update_ab_norestriction.Rd | 37 - man/update_theta.Rd | 35 - src/.clang-format | 26 + src/Makevars | 9 +- src/Makevars.win | 8 +- src/RcppExports.cpp | 116 --- src/fourpno_101315.cpp | 1479 ++++++++++++++++++-------------- 29 files changed, 1428 insertions(+), 1438 deletions(-) delete mode 100644 man/kappa_initialize.Rd delete mode 100644 man/update_2pno.Rd delete mode 100644 man/update_WKappaZ_NA.Rd delete mode 100644 man/update_ab_NA.Rd delete mode 100644 man/update_ab_norestriction.Rd delete mode 100644 man/update_theta.Rd create mode 100644 src/.clang-format diff --git a/.Rbuildignore b/.Rbuildignore index b72160a..ffc6dd1 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -4,3 +4,4 @@ ^README-.*\.png$ ^\.travis\.yml$ ^cran-comments\.md$ +^src/\.clang-format$ diff --git a/.travis.yml b/.travis.yml index 40f26e0..fa7c088 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,3 +1,26 @@ -language: R +# R for travis: see documentation at https://docs.travis-ci.com/user/languages/r + +###### Default build via devtools::use_travis() + +language: r sudo: false cache: packages + +###### Custom Options + +# Containers have 2 CPUs by default. Speed up the build by using both. +# c.f. https://docs.travis-ci.com/user/reference/overview/#virtualisation-environment-vs-operating-system +env: + global: + - MAKEFLAGS="-j 2" + - _R_CHECK_FORCE_SUGGESTS_=false + +# Run R package on both release and developer versions +jobs: + include: + - r: release + - r: devel + +# If one fails, avoid running the other. +matrix: + fast_finish: true diff --git a/DESCRIPTION b/DESCRIPTION index cb68381..8dd6d42 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,17 +1,21 @@ Package: fourPNO Type: Package Title: Bayesian 4 Parameter Item Response Model -Version: 1.0.4.1000 -Authors@R: c(person("Steven Andrew", "Culpepper", role = c("aut","cre"), email = - "sculpepp@illinois.edu")) +Version: 1.0.5 +Authors@R: c(person("Steven Andrew", "Culpepper", + email = "sculpepp@illinois.edu", + role = c("aut", "cre", "cph"), + comment = c(ORCID = "0000-0003-4226-6176") + )) Description: Estimate Barton & Lord's (1981) four parameter IRT model with lower and upper asymptotes using Bayesian formulation described by Culpepper (2016) . URL: https://github.com/tmsalab/fourPNO BugReports: https://github.com/tmsalab/fourPNO/issues License: GPL (>= 2) -Imports: Rcpp (>= 0.12.10) -LinkingTo: Rcpp (>= 0.12.10), RcppArmadillo (>= 0.7.800) -Depends: R (>= 3.0.2) -RoxygenNote: 6.0.1 - +Depends: R (>= 3.5.0) +Imports: Rcpp (>= 1.0.0) +LinkingTo: Rcpp, RcppArmadillo (>= 0.9.200) +RoxygenNote: 6.1.1 +Roxygen: list(markdown = TRUE) +Encoding: UTF-8 diff --git a/NAMESPACE b/NAMESPACE index 27ca53d..30ee174 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,13 +4,7 @@ export(Gibbs_2PNO) export(Gibbs_4PNO) export(Total_Tabulate) export(Y_4pno_simulate) -export(kappa_initialize) export(min2LL_4pno) export(rmvnorm) -export(update_2pno) -export(update_WKappaZ_NA) -export(update_ab_NA) -export(update_ab_norestriction) -export(update_theta) importFrom(Rcpp,evalCpp) -useDynLib("fourPNO", .registration=TRUE) +useDynLib(fourPNO, .registration=TRUE) diff --git a/NEWS.md b/NEWS.md index 980ac0a..2598d64 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,19 @@ -# fourPNO 1.0.4.1000 +# fourPNO 1.0.5 -- Switched generation of native registration to Rcpp. +## Changes + +- Increased _R_ dependency to 3.5.0. +- Autogenerate native registration with Rcpp 1.0.0. +- Suppressed internal functions from being exported. + +## Documentation + +- Enabled markdown in documentation. +- Cleaned up documentation entries. + +## Testing + +- Enabled TMSA Lab's travis-ci configuration. # fourPNO 1.0.4 diff --git a/R/RcppExports.R b/R/RcppExports.R index b1be89e..6055880 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,205 +1,187 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -#' Generate Random Multivariate Normal Distribution -#' -#' Creates a random Multivariate Normal when given number of obs, mean, and sigma. -#' @usage rmvnorm(n, mu, 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 sigma A \code{matrix} with dimensions m x m that provides the covariance matrix. -#' @return A \code{matrix} that is a Multivariate Normal distribution -#' @author James J Balamuta -#' @examples -#' #Call with the following data: -#' rmvnorm(2, c(0,0), diag(2)) -#' @export -rmvnorm <- function(n, mu, sigma) { - .Call(`_fourPNO_rmvnorm`, n, mu, sigma) -} - #' Initialize Thresholds -#' +#' #' Internal function for initializing item thresholds. -#' @param Ms A \code{vector} with the number of scale values. -#' @return A \code{matrix} that is a Multivariate Normal distribution -#' @seealso \code{\link{Gibbs_4PNO}} -#' @author Steven Andrew Culpepper -#' @export -kappa_initialize <- function(Ms) { - .Call(`_fourPNO_kappa_initialize`, Ms) -} +#' +#' @param Ms A `vector` with the number of scale values. +#' +#' @return +#' A `matrix` that is a Multivariate Normal distribution +#' +#' @author +#' Steven Andrew Culpepper +#' +#' @seealso +#' [Gibbs_4PNO()] +#' +#' @noRd +NULL #' Internal Function for Updating Theta in Gibbs Sampler -#' +#' #' Update theta in Gibbs sampler -#' @param N An \code{int}, which gives the number of observations. (> 0) -#' @param Z A \code{matrix} N by J of continuous augmented data. -#' @param as A \code{vector} of item discrimination parameters. -#' @param bs A \code{vector} of item threshold parameters. -#' @param theta A \code{vector} of prior thetas. -#' @param mu_theta The prior mean for theta. +#' +#' @param N An `int`, which gives the number of observations. +#' (> 0) +#' @param Z A `matrix` N by J of continuous augmented data. +#' @param as A `vector` of item discrimination parameters. +#' @param bs A `vector` of item threshold parameters. +#' @param theta A `vector` of prior thetas. +#' @param mu_theta The prior mean for theta. #' @param Sigma_theta_inv The prior inverse variance for theta. -#' @return A \code{vector} of thetas. -#' @seealso \code{\link{Gibbs_4PNO}} -#' @author Steven Andrew Culpepper -#' @export -update_theta <- function(N, Z, as, bs, theta, mu_theta, Sigma_theta_inv) { - .Call(`_fourPNO_update_theta`, N, Z, as, bs, theta, mu_theta, Sigma_theta_inv) -} +#' +#' @return +#' A `vector` of thetas. +#' +#' @author +#' Steven Andrew Culpepper +#' +#' @seealso +#' [Gibbs_4PNO()] +#' +#' @noRd +NULL #' Update a and b Parameters of 2PNO, 3PNO, 4PNO -#' +#' #' Update item slope and threshold -#' @param N An \code{int}, which gives the number of observations. (> 0) -#' @param J An \code{int}, which gives the number of items. (> 0) -#' @param Z A \code{matrix} N by J of continuous augmented data. -#' @param as A \code{vector} of item discrimination parameters. -#' @param bs A \code{vector} of item threshold parameters. -#' @param theta A \code{vector} of prior thetas. -#' @param mu_xi A two dimensional \code{vector} of prior item parameter means. -#' @param Sigma_xi_inv A two dimensional identity \code{matrix} of prior item parameter VC matrix. -#' @return A list of item parameters. -#' @seealso \code{\link{Gibbs_4PNO}} -#' @author Steven Andrew Culpepper -#' @export -update_ab_NA <- function(N, J, Z, as, bs, theta, mu_xi, Sigma_xi_inv) { - .Call(`_fourPNO_update_ab_NA`, N, J, Z, as, bs, theta, mu_xi, Sigma_xi_inv) -} +#' +#' @param N An `int`, which gives the number of observations. (> 0) +#' @param J An `int`, which gives the number of items. (> 0) +#' @param Z A `matrix` N by J of continuous augmented data. +#' @param as A `vector` of item discrimination parameters. +#' @param bs A `vector` of item threshold parameters. +#' @param theta A `vector` of prior thetas. +#' @param mu_xi A two dimensional `vector` of prior item parameter +NULL + +#' parameter VC matrix. +#' +#' @return +#' A `list` of item parameters. +#' +#' @author +#' Steven Andrew Culpepper +#' +#' @seealso +#' [Gibbs_4PNO()] +#' +#' @noRd +NULL #' Update a and b Parameters of 4pno without alpha > 0 Restriction -#' +#' #' Update item slope and threshold -#' @param N An \code{int}, which gives the number of observations. (> 0) -#' @param J An \code{int}, which gives the number of items. (> 0) -#' @param Z A \code{matrix} N by J of continuous augmented data. -#' @param as A \code{vector} of item discrimination parameters. -#' @param bs A \code{vector} of item threshold parameters. -#' @param theta A \code{vector} of prior thetas. -#' @param mu_xi A two dimensional \code{vector} of prior item parameter means. -#' @param Sigma_xi_inv A two dimensional identity \code{matrix} of prior item parameter VC matrix. -#' @return A list of item parameters. -#' @seealso \code{\link{Gibbs_4PNO}} -#' @author Steven Andrew Culpepper -#' @export -update_ab_norestriction <- function(N, J, Z, as, bs, theta, mu_xi, Sigma_xi_inv) { - .Call(`_fourPNO_update_ab_norestriction`, N, J, Z, as, bs, theta, mu_xi, Sigma_xi_inv) -} +#' +#' @param N An `int`, which gives the number of observations. (> 0) +#' @param J An `int`, which gives the number of items. (> 0) +#' @param Z A `matrix` N by J of continuous augmented data. +#' @param as A `vector` of item discrimination parameters. +#' @param bs A `vector` of item threshold parameters. +#' @param theta A `vector` of prior thetas. +#' @param mu_xi A two dimensional `vector` of prior item parameter +NULL + +#' +#' @return +#' A `list` of item parameters. +#' +#' @author +#' Steven Andrew Culpepper +#' +#' @seealso +#' [Gibbs_4PNO()] +#' +#' @noRd +NULL #' Update Lower and Upper Asymptote Parameters of 4PNO -#' +#' #' Internal function to update item lower and upper asymptote -#' @param Y A N by J \code{matrix} of item responses. -#' @param Ysum A \code{vector} of item total scores. -#' @param Z A \code{matrix} N by J of continuous augmented data. -#' @param as A \code{vector} of item discrimination parameters. -#' @param bs A \code{vector} of item threshold parameters. -#' @param gs A \code{vector} of item lower asymptote parameters. -#' @param ss A \code{vector} of item upper asymptote parameters. -#' @param theta A \code{vector} of prior thetas. -#' @param Kaps A \code{matrix} for item thresholds (used for internal computations). -#' @param alpha_c The lower asymptote prior 'a' parameter. -#' @param beta_c The lower asymptote prior 'b' parameter. -#' @param alpha_s The upper asymptote prior 'a' parameter. -#' @param beta_s The upper asymptote prior 'b' parameter. -#' @param gwg_reps The number of Gibbs within Gibbs MCMC samples for marginal distribution of gamma. -#' @return A list of item threshold parameters. -#' @seealso \code{\link{Gibbs_4PNO}} -#' @author Steven Andrew Culpepper -#' @export -update_WKappaZ_NA <- function(Y, Ysum, Z, as, bs, gs, ss, theta, Kaps, alpha_c, beta_c, alpha_s, beta_s, gwg_reps) { - .Call(`_fourPNO_update_WKappaZ_NA`, Y, Ysum, Z, as, bs, gs, ss, theta, Kaps, alpha_c, beta_c, alpha_s, beta_s, gwg_reps) -} - -#' Compute 4PNO Deviance -#' -#' Internal function to -2LL -#' @param N An \code{int}, which gives the number of observations. (> 0) -#' @param J An \code{int}, which gives the number of items. (> 0) -#' @param Y A N by J \code{matrix} of item responses. -#' @param as A \code{vector} of item discrimination parameters. -#' @param bs A \code{vector} of item threshold parameters. -#' @param gs A \code{vector} of item lower asymptote parameters. -#' @param ss A \code{vector} of item upper asymptote parameters. -#' @param theta A \code{vector} of prior thetas. -#' @return -2LL. -#' @seealso \code{\link{Gibbs_4PNO}} -#' @author Steven Andrew Culpepper -#' @export -min2LL_4pno <- function(N, J, Y, as, bs, gs, ss, theta) { - .Call(`_fourPNO_min2LL_4pno`, N, J, Y, as, bs, gs, ss, theta) -} - -#' Simulate from 4PNO Model -#' -#' Generate item responses under the 4PNO -#' @param N An \code{int}, which gives the number of observations. (> 0) -#' @param J An \code{int}, which gives the number of items. (> 0) -#' @param as A \code{vector} of item discrimination parameters. -#' @param bs A \code{vector} of item threshold parameters. -#' @param gs A \code{vector} of item lower asymptote parameters. -#' @param ss A \code{vector} of item upper asymptote parameters. -#' @param theta A \code{vector} of prior thetas. -#' @return A N by J \code{matrix} of dichotomous item responses. -#' @seealso \code{\link{Gibbs_4PNO}} -#' @author Steven Andrew Culpepper -#' @export -Y_4pno_simulate <- function(N, J, as, bs, gs, ss, theta) { - .Call(`_fourPNO_Y_4pno_simulate`, N, J, as, bs, gs, ss, theta) -} - -#' Calculate Tabulated Total Scores -#' -#' Internal function to -2LL -#' @param N An \code{int}, which gives the number of observations. (> 0) -#' @param J An \code{int}, which gives the number of items. (> 0) -#' @param Y A N by J \code{matrix} of item responses. -#' @return A vector of tabulated total scores. -#' @seealso \code{\link{Gibbs_4PNO}} -#' @author Steven Andrew Culpepper -#' @export -Total_Tabulate <- function(N, J, Y) { - .Call(`_fourPNO_Total_Tabulate`, N, J, Y) -} +#' +#' @param Y A N by J `matrix` of item responses. +#' @param Ysum A `vector` of item total scores. +#' @param Z A `matrix` N by J of continuous augmented data. +#' @param as A `vector` of item discrimination parameters. +#' @param bs A `vector` of item threshold parameters. +#' @param gs A `vector` of item lower asymptote parameters. +#' @param ss A `vector` of item upper asymptote parameters. +#' @param theta A `vector` of prior thetas. +#' @param Kaps A `matrix` for item thresholds +#' (used for internal computations). +#' @param alpha_c The lower asymptote prior 'a' parameter. +#' @param beta_c The lower asymptote prior 'b' parameter. +#' @param alpha_s The upper asymptote prior 'a' parameter. +#' @param beta_s The upper asymptote prior 'b' parameter. +#' @param gwg_reps The number of Gibbs within Gibbs MCMC samples for +#' marginal distribution of gamma. +#' +#' @return +#' A `list` of item threshold parameters. +#' +#' @author +#' Steven Andrew Culpepper +#' +#' @seealso +#' [Gibbs_4PNO()] +#' +#' @noRd +NULL #' Gibbs Implementation of 4PNO -#' +#' #' Internal function to -2LL -#' @param Y A N by J \code{matrix} of item responses. -#' @param mu_xi A two dimensional \code{vector} of prior item parameter means. -#' @param Sigma_xi_inv A two dimensional identity \code{matrix} of prior item parameter VC matrix. +#' +#' @param Y A N by J `matrix` of item responses. +#' @param mu_xi A two dimensional `vector` of prior item parameter +#' means. +#' @param Sigma_xi_inv A two dimensional identity `matrix` of prior item +#' parameter VC matrix. #' @param mu_theta The prior mean for theta. #' @param Sigma_theta_inv The prior inverse variance for theta. -#' @param alpha_c The lower asymptote prior 'a' parameter. -#' @param beta_c The lower asymptote prior 'b' parameter. -#' @param alpha_s The upper asymptote prior 'a' parameter. -#' @param beta_s The upper asymptote prior 'b' parameter. -#' @param burnin The number of MCMC samples to discard. -#' @param cTF A J dimensional \code{vector} indicating which lower asymptotes to estimate. 0 = exclude lower asymptote and 1 = include lower asymptote. -#' @param sTF A J dimensional \code{vector} indicating which upper asymptotes to estimate. 0 = exclude upper asymptote and 1 = include upper asymptote. -#' @param gwg_reps The number of Gibbs within Gibbs MCMC samples for marginal distribution of gamma. Values between 5 to 10 are adequate. -#' @param chain_length The number of MCMC samples. -#' @return Samples from posterior. -#' @author Steven Andrew Culpepper +#' @param alpha_c The lower asymptote prior 'a' parameter. +#' @param beta_c The lower asymptote prior 'b' parameter. +#' @param alpha_s The upper asymptote prior 'a' parameter. +#' @param beta_s The upper asymptote prior 'b' parameter. +#' @param burnin The number of MCMC samples to discard. +#' @param cTF A J dimensional `vector` indicating which +#' lower asymptotes to estimate. +#' 0 = exclude lower asymptote and +#' 1 = include lower asymptote. +#' @param sTF A J dimensional `vector` indicating which +#' upper asymptotes to estimate. +#' 0 = exclude upper asymptote and +#' 1 = include upper asymptote. +#' @param gwg_reps The number of Gibbs within Gibbs MCMC samples for +#' marginal distribution of gamma. Values between +#' 5 to 10 are adequate. +#' @param chain_length The number of MCMC samples. +#' +#' @return +#' Samples from posterior. +#' +#' @author +#' Steven Andrew Culpepper +#' #' @export #' @examples -#' require(fourPNO) -#' #' # Simulate small 4PNO dataset to demonstrate function #' J = 5 #' N = 100 -#' +#' #' # Population item parameters #' as_t = rnorm(J,mean=2,sd=.5) #' bs_t = rnorm(J,mean=0,sd=.5) -#' +#' #' # Sampling gs and ss with truncation #' gs_t = rbeta(J,1,8) #' ps_g = pbeta(1-gs_t,1,8) #' ss_t = qbeta(runif(J)*ps_g,1,8) #' theta_t <- rnorm(N) #' Y_t = Y_4pno_simulate(N,J,as=as_t,bs=bs_t,gs=gs_t,ss=ss_t,theta=theta_t) -#' +#' #' # Setting prior parameters #' mu_theta=0 #' Sigma_theta_inv=1 @@ -207,62 +189,80 @@ Total_Tabulate <- function(N, J, Y) { #' alpha_c=alpha_s=beta_c=beta_s=1 #' Sigma_xi_inv = solve(2*matrix(c(1,0,0,1),2,2)) #' burnin = 1000 -#' +#' #' # Execute Gibbs sampler -#' out_t = Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta,Sigma_theta_inv,alpha_c,beta_c,alpha_s, -#' beta_s,burnin,rep(1,J),rep(1,J),gwg_reps=5,chain_length=burnin*2) -#' +#' out_t = Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta, +#' Sigma_theta_inv,alpha_c,beta_c,alpha_s, +#' beta_s,burnin,rep(1,J),rep(1,J), +#' gwg_reps=5,chain_length=burnin*2) +#' #' # Summarizing posterior distribution -#' OUT = cbind(apply(out_t$AS[,-c(1:burnin)],1,mean),apply(out_t$BS[,-c(1:burnin)],1,mean), -#' apply(out_t$GS[,-c(1:burnin)],1,mean),apply(out_t$SS[,-c(1:burnin)],1,mean), -#' apply(out_t$AS[,-c(1:burnin)],1,sd),apply(out_t$BS[,-c(1:burnin)],1,sd), -#' apply(out_t$GS[,-c(1:burnin)],1,sd),apply(out_t$SS[,-c(1:burnin)],1,sd) ) -#' +#' OUT = cbind(apply(out_t$AS[,-c(1:burnin)],1,mean), +#' apply(out_t$BS[,-c(1:burnin)],1,mean), +#' apply(out_t$GS[,-c(1:burnin)],1,mean), +#' apply(out_t$SS[,-c(1:burnin)],1,mean), +#' apply(out_t$AS[,-c(1:burnin)],1,sd), +#' apply(out_t$BS[,-c(1:burnin)],1,sd), +#' apply(out_t$GS[,-c(1:burnin)],1,sd), +#' apply(out_t$SS[,-c(1:burnin)],1,sd) ) +#' #' OUT = cbind(1:J,OUT) -#' colnames(OUT) = c('Item','as','bs','gs','ss','as_sd','bs_sd','gs_sd','ss_sd') -#' print(OUT,digits=3) -Gibbs_4PNO <- function(Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, alpha_c, beta_c, alpha_s, beta_s, burnin, cTF, sTF, gwg_reps, chain_length = 10000L) { - .Call(`_fourPNO_Gibbs_4PNO`, Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, alpha_c, beta_c, alpha_s, beta_s, burnin, cTF, sTF, gwg_reps, chain_length) -} +#' colnames(OUT) = +NULL #' Update 2PNO Model Parameters -#' +#' #' Internal function to update 2PNO parameters -#' @param N The number of observations. -#' @param J The number of items. -#' @param Y A N by J \code{matrix} of item responses. -#' @param Z A \code{matrix} N by J of continuous augmented data. -#' @param as A \code{vector} of item discrimination parameters. -#' @param bs A \code{vector} of item threshold parameters. -#' @param theta A \code{vector} of prior thetas. -#' @param Kaps A \code{matrix} for item thresholds (used for internal computations). -#' @param mu_xi Prior mean for item parameters. -#' @param Sigma_xi_inv Prior item parameter inverse variance-covariance matrix. -#' @param mu_theta Prior mean for theta. +#' +#' @param N The number of observations. +#' @param J The number of items. +#' @param Y A N by J `matrix` of item responses. +#' @param Z A `matrix` N by J of continuous augmented data. +#' @param as A `vector` of item discrimination parameters. +#' @param bs A `vector` of item threshold parameters. +#' @param theta A `vector` of prior thetas. +#' @param Kaps A `matrix` for item thresholds +#' (used for internal computations). +#' @param mu_xi Prior mean for item parameters. +#' @param Sigma_xi_inv Prior item parameter inverse variance-covariance +#' matrix. +#' @param mu_theta Prior mean for theta. #' @param Sigma_theta_inv Prior inverse variance for theta. -#' @return A list of item parameters. -#' @seealso \code{\link{Gibbs_4PNO}} -#' @author Steven Andrew Culpepper -#' @export -update_2pno <- function(N, J, Y, Z, as, bs, theta, Kaps, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv) { - .Call(`_fourPNO_update_2pno`, N, J, Y, Z, as, bs, theta, Kaps, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv) -} +#' +#' @return +#' A `list` of item parameters. +#' +#' @author +#' Steven Andrew Culpepper +#' +#' @seealso +#' [Gibbs_2PNO()] +#' +#' @noRd +NULL #' Gibbs Implementation of 2PNO -#' +#' #' Implement Gibbs 2PNO Sampler -#' @param Y A N by J \code{matrix} of item responses. -#' @param mu_xi A two dimensional \code{vector} of prior item parameter means. -#' @param Sigma_xi_inv A two dimensional identity \code{matrix} of prior item parameter VC matrix. -#' @param mu_theta The prior mean for theta. -#' @param Sigma_theta_inv The prior inverse variance for theta. -#' @param burnin The number of MCMC samples to discard. -#' @param chain_length The number of MCMC samples. -#' @return Samples from posterior. -#' @author Steven Andrew Culpepper +#' +#' @param Y A N by J `matrix` of item responses. +#' @param mu_xi A two dimensional `vector` of prior item parameter +#' means. +#' @param Sigma_xi_inv A two dimensional identity `matrix` of prior item +#' parameter VC matrix. +#' @param mu_theta The prior mean for theta. +#' @param Sigma_theta_inv The prior inverse variance for theta. +#' @param burnin The number of MCMC samples to discard. +#' @param chain_length The number of MCMC samples. +#' +#' @return +#' Samples from posterior. +#' +#' @author +#' Steven Andrew Culpepper +#' #' @export #' @examples -#' require(fourPNO) #' # simulate small 2PNO dataset to demonstrate function #' J = 5 #' N = 100 @@ -270,14 +270,14 @@ update_2pno <- function(N, J, Y, Z, as, bs, theta, Kaps, mu_xi, Sigma_xi_inv, mu #' # Population item parameters #' as_t = rnorm(J,mean=2,sd=.5) #' bs_t = rnorm(J,mean=0,sd=.5) -#' +#' #' # Sampling gs and ss with truncation #' gs_t = rbeta(J,1,8) #' ps_g = pbeta(1-gs_t,1,8) #' ss_t = qbeta(runif(J)*ps_g,1,8) #' theta_t = rnorm(N) #' Y_t = Y_4pno_simulate(N,J,as=as_t,bs=bs_t,gs=gs_t,ss=ss_t,theta=theta_t) -#' +#' #' # Setting prior parameters #' mu_theta = 0 #' Sigma_theta_inv = 1 @@ -285,19 +285,123 @@ update_2pno <- function(N, J, Y, Z, as, bs, theta, Kaps, mu_xi, Sigma_xi_inv, mu #' alpha_c = alpha_s = beta_c = beta_s = 1 #' Sigma_xi_inv = solve(2*matrix(c(1,0,0,1), 2, 2)) #' burnin = 1000 -#' +#' #' # Execute Gibbs sampler. This should take about 15.5 minutes -#' out_t <- Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta,Sigma_theta_inv,alpha_c,beta_c,alpha_s, -#' beta_s,burnin,rep(1,J),rep(1,J),gwg_reps=5,chain_length=burnin*2) -#' +#' out_t = Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta,Sigma_theta_inv, +#' alpha_c,beta_c,alpha_s, beta_s,burnin, +#' rep(1,J),rep(1,J),gwg_reps=5,chain_length=burnin*2) +#' #' # Summarizing posterior distribution -#' OUT = cbind(apply(out_t$AS[,-c(1:burnin)],1,mean),apply(out_t$BS[,-c(1:burnin)],1,mean), -#' apply(out_t$GS[,-c(1:burnin)],1,mean),apply(out_t$SS[,-c(1:burnin)],1,mean), -#' apply(out_t$AS[,-c(1:burnin)],1,sd),apply(out_t$BS[,-c(1:burnin)],1,sd), -#' apply(out_t$GS[,-c(1:burnin)],1,sd),apply(out_t$SS[,-c(1:burnin)],1,sd) ) -#' OUT = cbind(1:J, OUT) -#' colnames(OUT) = c('Item','as','bs','gs','ss','as_sd','bs_sd','gs_sd','ss_sd') -#' print(OUT, digits=3) +#' OUT = +NULL + +#' apply(out_t$GS[,-c(1:burnin)],1,mean),apply(out_t$SS[,-c(1:burnin)],1,mean), +#' apply(out_t$AS[,-c(1:burnin)],1,sd),apply(out_t$BS[,-c(1:burnin)],1,sd), ' +NULL + +#' Generate Random Multivariate Normal Distribution +#' +#' Creates a random Multivariate Normal when given number of +#' obs, mean, and sigma. +#' +#' @param n An `int`, which gives the number of observations. (> 0) +#' @param mu A `vector` length m that represents the means of +#' the normals. +#' @param sigma A `matrix` with dimensions m x m that provides the +#' covariance matrix. +#' +#' @return +#' A `matrix` that is a Multivariate Normal distribution +#' +#' @author +#' James J Balamuta +#' +#' @export +#' @examples +#' # Call with the following data: +#' rmvnorm(2, c(0,0), diag(2)) +rmvnorm <- function(n, mu, sigma) { + .Call(`_fourPNO_rmvnorm`, n, mu, sigma) +} + +#' Compute 4PNO Deviance +#' +#' Internal function to -2LL +#' +#' @param N An `int`, which gives the number of observations. (> 0) +#' @param J An `int`, which gives the number of items. (> 0) +#' @param Y A N by J `matrix` of item responses. +#' @param as A `vector` of item discrimination parameters. +#' @param bs A `vector` of item threshold parameters. +#' @param gs A `vector` of item lower asymptote parameters. +#' @param ss A `vector` of item upper asymptote parameters. +#' @param theta A `vector` of prior thetas. +#' +#' @return +#' -2LL. +#' @author +#' Steven Andrew Culpepper +#' +#' @seealso +#' [Gibbs_4PNO()] +#' +#' @export +min2LL_4pno <- function(N, J, Y, as, bs, gs, ss, theta) { + .Call(`_fourPNO_min2LL_4pno`, N, J, Y, as, bs, gs, ss, theta) +} + +#' Simulate from 4PNO Model +#' +#' Generate item responses under the 4PNO +#' +#' @param N An `int`, which gives the number of observations. (> 0) +#' @param J An `int`, which gives the number of items. (> 0) +#' @param as A `vector` of item discrimination parameters. +#' @param bs A `vector` of item threshold parameters. +#' @param gs A `vector` of item lower asymptote parameters. +#' @param ss A `vector` of item upper asymptote parameters. +#' @param theta A `vector` of prior thetas. +#' +#' @return +#' A N by J `matrix` of dichotomous item responses. +#' +#' @author +#' Steven Andrew Culpepper +#' +#' @seealso +#' [Gibbs_4PNO()] +#' +#' @export +Y_4pno_simulate <- function(N, J, as, bs, gs, ss, theta) { + .Call(`_fourPNO_Y_4pno_simulate`, N, J, as, bs, gs, ss, theta) +} + +#' Calculate Tabulated Total Scores +#' +#' Internal function to -2LL +#' +#' @param N An `int`, which gives the number of observations. (> 0) +#' @param J An `int`, 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 +#' +#' @seealso +#' [Gibbs_4PNO()] +#' +#' @export +Total_Tabulate <- function(N, J, Y) { + .Call(`_fourPNO_Total_Tabulate`, N, J, Y) +} + +Gibbs_4PNO <- function(Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, alpha_c, beta_c, alpha_s, beta_s, burnin, cTF, sTF, gwg_reps, chain_length = 10000L) { + .Call(`_fourPNO_Gibbs_4PNO`, Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, alpha_c, beta_c, alpha_s, beta_s, burnin, cTF, sTF, gwg_reps, chain_length) +} + Gibbs_2PNO <- function(Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, burnin, chain_length = 10000L) { .Call(`_fourPNO_Gibbs_2PNO`, Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, burnin, chain_length) } diff --git a/R/fourPNO-internal.R b/R/fourPNO-internal.R index ac0cdcc..2edc622 100644 --- a/R/fourPNO-internal.R +++ b/R/fourPNO-internal.R @@ -1,131 +1,3 @@ -.Random.seed <- -c(403L, 10L, -916551914L, 1672741312L, 1551819139L, 1470003105L, --242550896L, 1313623150L, -395847143L, -441564501L, 731944602L, --1666989476L, 1679130463L, 413471029L, -934399796L, -795239902L, -26596477L, -1884263049L, 1508164254L, 1963004248L, -1674377029L, -270505561L, 729794024L, 230146070L, -1626690671L, -23194557L, --1993364654L, 1711206372L, 1982475815L, -2016762243L, 389117812L, -755690330L, 1913044389L, 1218607551L, -1217777082L, 1080434896L, -1389132531L, 598080529L, 1450773056L, 937643550L, 1762691305L, --1627207141L, 988903978L, -721654452L, -1329320177L, 1027797573L, -1641239548L, -1832220078L, -1027934259L, -1623858809L, -446745490L, --1351062456L, -1107945909L, -232969303L, 1781375736L, -1602337242L, --501080255L, -2126207789L, 1217341570L, -683227468L, 659644983L, -1872390445L, -2028729788L, 426391402L, 1267147733L, 2071672175L, -1473540214L, 1870257056L, 2105228387L, 1996704001L, -578698128L, --275095538L, 1105150073L, 1514569739L, -2022305478L, -1971482564L, -1836029823L, -773149483L, -376794900L, 1117850178L, -1029452387L, --949159913L, 1630764414L, -1411374280L, -1275642853L, 954828153L, --239193144L, -609927178L, 654115377L, 346305635L, -348160974L, --284927612L, -232808505L, -1033088355L, 2132775444L, -1560502086L, -1760513285L, -1070450977L, 43531110L, -2118804752L, -1325422061L, --1332853839L, 1983278688L, -1731486722L, -1944018039L, 1764508859L, --907042678L, 1934238572L, -1748473617L, -1300995291L, 959393244L, -639941554L, -416472403L, -822044057L, -1359631986L, -1481245080L, --1230632341L, 2083546121L, 48511640L, 56163782L, -1886397151L, --1413952077L, 1509116194L, 1142514452L, -316736617L, 1179476877L, -899686052L, 160305290L, -1246144779L, 143518031L, 1370998358L, -1890974976L, 1963984963L, 371567585L, -1365757104L, -1131371090L, --1117541543L, -131514005L, -1163958054L, -1386757732L, -1542984289L, -1539325941L, 1243508364L, -1414739102L, -743644355L, -958229577L, -906483806L, 1488197016L, -967483141L, 371219609L, -1688421080L, -508726230L, -884600751L, -575986941L, 1618514578L, 1998495780L, --494432537L, 410837821L, -1219535948L, -1757525990L, 1506966373L, --410245633L, 1757735046L, 901682192L, -348698957L, 603288913L, --1343263488L, 1553102814L, 1925723177L, 2038141915L, 601098474L, -1855178380L, 1493187023L, -1653101179L, -698853828L, 1638841106L, --1785325555L, 335981383L, 952680750L, -1773623160L, -1249792629L, --204892823L, -993085640L, -1195820570L, -2006385279L, 1460318483L, --734982462L, -992728716L, 1971409271L, -202325395L, 981585156L, -599484458L, 55584149L, -604354001L, 1009063734L, -872544416L, -53317027L, 384919233L, 559013552L, -823469618L, -1033212615L, --1095909045L, -1848851718L, -2077076484L, 1181930559L, 1162964245L, --1460447700L, -1343997182L, 800679901L, 1286853335L, -174789186L, --342997000L, 1452983259L, 804339257L, -299432056L, -2011588810L, --1892465295L, 496566179L, 38627058L, 842157636L, -149286649L, -1530783965L, -2059600428L, -202600198L, 1612474437L, 132591775L, -292272678L, -1960317264L, 1210094419L, 875153120L, 157857916L, --2070079792L, 1409849762L, 527956008L, 1202080524L, 1799784892L, --1841519822L, -1526445376L, -4658124L, -1852129592L, 1658582202L, -2140378768L, -103633172L, -819406252L, -108887406L, -332281568L, -511021516L, 845005536L, 466796738L, -2062005208L, 245438380L, -2147166908L, -124331022L, 1198282288L, 570658724L, -886321736L, --1606887590L, -1338814320L, 577653628L, -1529732844L, 90032450L, -1289228544L, 1558601148L, -1653989040L, 337746210L, -1674962360L, -2139189260L, -1040406948L, -242611758L, -624478720L, -742202028L, -735620328L, 515739898L, -1767954736L, -1009128660L, 1380200692L, --1452353134L, 481197408L, -1126802420L, 801815776L, -1969704350L, --1766628184L, -730858420L, 674826108L, -1368352462L, -1604492688L, --1837440188L, -614579240L, 1963283322L, 1070831664L, 589746108L, --1013730732L, 2045521986L, -1920749920L, 1433584316L, -2143941680L, -771308898L, 38948584L, -162485044L, -1738923524L, 1825193458L, -1277977728L, -271701516L, -1085945208L, -1883120838L, 626462032L, -2127031404L, 636474772L, -989306798L, -424296992L, -1880845108L, -669859488L, 1703276930L, 779348200L, -2059941140L, -1731254980L, -1876422770L, 1921913456L, -1933141468L, 1637017656L, -1120810726L, --59423536L, -1528166084L, -1046449260L, -639154302L, -476443712L, --1461905860L, 1834746256L, -129279966L, -1930537720L, -939251572L, -550072092L, 2080806034L, 1625609600L, 447537172L, 133121832L, --1757723078L, -364523056L, 2102753644L, 1039778740L, -1772385774L, -1237977056L, -1697519156L, 978385120L, 880643106L, -1732996312L, -52328972L, 1135514428L, -1192513550L, -654834960L, 1843681860L, -15318616L, 590673722L, -1358376976L, 1684251324L, 657167636L, -491461826L, -924271520L, 1982932348L, -473751600L, 1469444386L, -1577522856L, -359901300L, 2079360188L, -1379419342L, 251274816L, --1743493964L, 1186534472L, 7755194L, -841025008L, -537707284L, -1297949012L, -2035797230L, -912298592L, 772561612L, 598678112L, -987082690L, 1946055464L, 847968684L, 198409532L, -118211726L, --1599771344L, 1665738532L, 964502840L, 440649434L, 1227030416L, --708814084L, -544597100L, -1711557566L, -668003072L, 385181756L, -703389008L, 1333609634L, 1936427848L, 1959591692L, -2101821732L, -857648978L, -915286400L, 832936020L, -1606786072L, -723850502L, -1555191120L, 49495852L, 1497801588L, 1446634386L, -302689696L, --1968834676L, 979720032L, -688323614L, 1272880040L, -179467060L, -1580366332L, 1898469298L, 2070963696L, -773225020L, -61744808L, --2096460422L, -570032720L, 1338549180L, -1610731436L, -1173811262L, --177954784L, -866651972L, -425158064L, 1089062882L, -1420575128L, --760566196L, 1179750652L, -189507598L, 823813376L, -250043532L, --1386487800L, 405721530L, 695606736L, -177056788L, -1610232940L, -1022811090L, 1566451552L, -647954100L, -109792992L, 523825026L, -366069480L, -346105620L, 960034364L, 229782386L, -264735120L, -657056804L, 389869368L, 1774909594L, 2009525968L, -728347460L, --809293164L, 329327234L, -716705600L, -1163057860L, -718325616L, -84176563L, -1961948716L, -655120158L, -1349979785L, 1844441281L, --608891834L, 298570628L, 566887733L, -830635233L, 1863689848L, --206824618L, -1441078813L, 1593079717L, 1192105538L, -1098720208L, -223608921L, -641962549L, -1077223876L, 807501930L, -843695457L, --518355159L, -769311810L, 1302768300L, 1978546173L, 1760378695L, --42450896L, -1944811186L, 1880150395L, -1622628643L, -1284120822L, --267889112L, -563151215L, 1514886787L, -1812330620L, -1258324558L, --1306660089L, -1153005231L, -838151978L, -382386764L, 1791367621L, --148874993L, 952952712L, 867341990L, -1411283373L, -1540556875L, -399807762L, 441298656L, 2057892681L, 1676045691L, 1082107980L, -955593178L, 920259535L, -24803943L, 334863790L, 693512892L, 1315064621L, --1816643369L, 1167671584L, -1347205506L, -592375733L, 1819503181L, --327026566L, -669562696L, -440243615L, -929640429L, -1541210572L, --631875518L, -1576871849L, 1349614497L, 712473702L, -149842844L, --659312747L, 154716351L, 2085271640L, -837227402L, 2136503427L, -576967877L, -1438871070L, 1058396240L, 112404857L, -646926677L, -1895458332L, 1154486602L, -2817729L, 119596937L, 1116948382L, --1804364212L, -847831011L, 31519143L, 1627092176L, 237172334L, -1540477211L, -1699462531L, 1635031146L, -455727864L, 1479193009L, --2112207325L, 1413497636L, 908567506L, 180592679L, -1364608399L, -1529901238L, 330745940L, -84821787L, 396973295L, -105506776L, -1892953222L, 1686693555L, -1380767723L, -1050162574L, 2023659200L, -668560169L, 119187867L, 665562732L, 74568314L, -148384849L, -514341191L, --505272114L, 387744284L, 259344781L, -544077833L, -1547531520L, -207434078L, -589065685L, -1142883027L, -2046172774L, 1873878872L, --1018087743L, -25753229L, -535731564L, -867150046L, -771034569L, -325179137L, 1195593222L, 1797387332L, 1067648373L, 1042794719L, -342490168L, 210246678L, -1234133981L, -371667995L, 495975938L, -1432376432L, 1659707929L, -475427829L, -1911058692L, 1983407914L, --617936545L, 1099809001L, -1861552130L, 1147681772L, 60608317L, --1886924025L, -84702480L, 1314615950L, -258667333L, 1493207197L, -779686986L, -1890087320L, 419029585L, -1391419069L, -611313596L, -1186962162L, -1316986169L, 2086674577L, -500175338L, -636113804L, --1255905531L, -849337009L, 353033183L) - - -#' @useDynLib "fourPNO", .registration=TRUE +#' @useDynLib fourPNO, .registration=TRUE #' @importFrom Rcpp evalCpp "_PACKAGE" diff --git a/README.Rmd b/README.Rmd index 8589f5a..a8e2951 100644 --- a/README.Rmd +++ b/README.Rmd @@ -12,11 +12,62 @@ knitr::opts_chunk$set( ) ``` -[![Travis-CI Build Status](https://travis-ci.org/tmsalab/fourPNO.svg?branch=master)](https://travis-ci.org/tmsalab/fourPNO) -[![CRAN RStudio mirror downloads](http://cranlogs.r-pkg.org/badges/fourPNO)](http://www.r-pkg.org/pkg/fourPNO) -[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/fourPNO)](https://cran.r-project.org/package=fourPNO) +[![Build Status](https://travis-ci.org/tmsalab/fourPNO.svg)](https://travis-ci.org/tmsalab/fourPNO) +[![Package-License](http://img.shields.io/badge/license-GPL%20(%3E=2)-brightgreen.svg?style=flat)](http://www.gnu.org/licenses/gpl-2.0.html) +[![CRAN Version Badge](http://www.r-pkg.org/badges/version/fourPNO)](https://cran.r-project.org/package=fourPNO) +[![CRAN Status](https://cranchecks.info/badges/worst/fourPNO)](https://cran.r-project.org/web/checks/check_results_fourPNO.html) +[![RStudio CRAN Mirror's Monthly Downloads](http://cranlogs.r-pkg.org/badges/fourPNO?color=brightgreen)](http://www.r-pkg.org/pkg/fourPNO) +[![RStudio CRAN Mirror's Total Downloads](http://cranlogs.r-pkg.org/badges/grand-total/fourPNO?color=brightgreen)](http://www.r-pkg.org/pkg/fourPNO) +[![Coverage status](https://codecov.io/gh/tmsalab/fourPNO/branch/master/graph/badge.svg)](https://codecov.io/github/tmsalab/fourPNO?branch=master) # `fourPNO` R package -Estimate Barton & Lord's four parameter IRT model with lower and upper -asymptotes using Bayesian formulation described by Culpepper (2015). +Estimate Barton & Lord's (1981) +four parameter IRT model with lower and upper asymptotes using Bayesian +formulation described by Culpepper (2016) . + +## Installation + +You can install `fourPNO` from CRAN using: + +```{r cran-installation, eval = FALSE} +install.packages("fourPNO") +``` + +Or, you can be on the cutting-edge development version on GitHub using: + +```{r gh-installation, eval = FALSE} +if(!requireNamespace("devtools")) install.packages("devtools") +devtools::install_github("tmsalab/fourPNO") +``` + +## Usage + +To use the `fourPNO` package, load it into _R_ using: + +```{r example, message = FALSE} +library("fourPNO") +``` + +Inside the package, the estimation routines can be viewed as: + +- `Gibbs_2PNO()` +- `Gibbs_4PNO()` + +## Author + +Steven Andrew Culpepper + +## Citing the `fourPNO` package + +To ensure future development of the package, please cite `fourPNO` +package if used during an analysis or simulation study. Citation information +for the package may be acquired by using in *R*: + +```{r citation-info, eval = FALSE} +citation("fourPNO") +``` + +## License + +GPL (>= 2) diff --git a/README.md b/README.md index b7f7ac1..cdca617 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,70 @@ -[![Travis-CI Build -Status](https://travis-ci.org/tmsalab/fourPNO.svg?branch=master)](https://travis-ci.org/tmsalab/fourPNO) -[![CRAN RStudio mirror -downloads](http://cranlogs.r-pkg.org/badges/fourPNO)](http://www.r-pkg.org/pkg/fourPNO) -[![CRAN\_Status\_Badge](http://www.r-pkg.org/badges/version/fourPNO)](https://cran.r-project.org/package=fourPNO) +[![Build +Status](https://travis-ci.org/tmsalab/fourPNO.svg)](https://travis-ci.org/tmsalab/fourPNO) +[![Package-License](http://img.shields.io/badge/license-GPL%20\(%3E=2\)-brightgreen.svg?style=flat)](http://www.gnu.org/licenses/gpl-2.0.html) +[![CRAN Version +Badge](http://www.r-pkg.org/badges/version/fourPNO)](https://cran.r-project.org/package=fourPNO) +[![CRAN +Status](https://cranchecks.info/badges/worst/fourPNO)](https://cran.r-project.org/web/checks/check_results_fourPNO.html) +[![RStudio CRAN Mirror’s Monthly +Downloads](http://cranlogs.r-pkg.org/badges/fourPNO?color=brightgreen)](http://www.r-pkg.org/pkg/fourPNO) +[![RStudio CRAN Mirror’s Total +Downloads](http://cranlogs.r-pkg.org/badges/grand-total/fourPNO?color=brightgreen)](http://www.r-pkg.org/pkg/fourPNO) +[![Coverage +status](https://codecov.io/gh/tmsalab/fourPNO/branch/master/graph/badge.svg)](https://codecov.io/github/tmsalab/fourPNO?branch=master) # `fourPNO` R package -Estimate Barton & Lord’s four parameter IRT model with lower and upper -asymptotes using Bayesian formulation described by Culpepper (2015). +Estimate Barton & Lord’s (1981) + four parameter IRT +model with lower and upper asymptotes using Bayesian formulation +described by Culpepper (2016) +. + +## Installation + +You can install `fourPNO` from CRAN using: + +``` r +install.packages("fourPNO") +``` + +Or, you can be on the cutting-edge development version on GitHub using: + +``` r +if(!requireNamespace("devtools")) install.packages("devtools") +devtools::install_github("tmsalab/fourPNO") +``` + +## Usage + +To use the `fourPNO` package, load it into *R* using: + +``` r +library("fourPNO") +``` + +Inside the package, the estimation routines can be viewed as: + + - `Gibbs_2PNO()` + - `Gibbs_4PNO()` + +## Author + +Steven Andrew Culpepper + +## Citing the `fourPNO` package + +To ensure future development of the package, please cite `fourPNO` +package if used during an analysis or simulation study. Citation +information for the package may be acquired by using in *R*: + +``` r +citation("fourPNO") +``` + +## License + +GPL (\>= 2) diff --git a/cran-comments.md b/cran-comments.md index 4746a12..4736a42 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,15 +1,8 @@ ## Test environments -* local OS X install, R 3.4.0 -* ubuntu 12.04 (on travis-ci), R 3.4.0 +* local OS X install, R 3.5.2 +* ubuntu 14.04 (on travis-ci), R 3.5.2 * win-builder (devel and release) -## Feedback - -The previously submitted version of this package had a CRAN maintainer request -the addition paper reference material including the year and doi of publications -referenced within the `DESCRIPTION` file's description section. -This was addressed. - ## R CMD check results 0 errors | 0 warnings | 1 note diff --git a/inst/CITATION b/inst/CITATION index f843ea6..43457a3 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -1,4 +1,14 @@ -citHeader("To cite package fourPNO in publications use:") +citHeader("To cite the package 'fourPNO' in publications use:") + +citEntry(entry = "Manual", + title = "{fourPNO: Bayesian 4 Parameter Item Response Model}", + author = personList(as.person("Steven Andrew Culpepper") + ), + year = 2019, + textVersion = paste("Culpepper, S.A. (2019)", + "fourPNO: Bayesian 4 Parameter Item Response Model.", + "URL https://cran.r-project.org/package=fourPNO.") +) citEntry(entry = "Article", title = "Revisiting the 4-Parameter Item Response Model: Bayesian Estimation and Application", @@ -11,5 +21,9 @@ citEntry(entry = "Article", issn = "1860-0980", doi = "10.1007/s11336-015-9477-6", url = "http://dx.doi.org/10.1007/s11336-015-9477-6", - textVersion = "Culpepper, S. A. (2016). Revisiting the 4-parameter item response model: Bayesian estimation and application, Psychometrika, 81, 1142-1163. DOI: 10.1007/s11336-015-9477-6.") + textVersion = paste("Culpepper, S. A. (2016).", + "Revisiting the 4-parameter item response model: Bayesian estimation and application,", + "Psychometrika, 81, 1142-1163.", + "DOI: 10.1007/s11336-015-9477-6.") +) diff --git a/man/Gibbs_2PNO.Rd b/man/Gibbs_2PNO.Rd index 9a7ddec..27388bc 100644 --- a/man/Gibbs_2PNO.Rd +++ b/man/Gibbs_2PNO.Rd @@ -10,9 +10,11 @@ Gibbs_2PNO(Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, burnin, \arguments{ \item{Y}{A N by J \code{matrix} of item responses.} -\item{mu_xi}{A two dimensional \code{vector} of prior item parameter means.} +\item{mu_xi}{A two dimensional \code{vector} of prior item parameter +means.} -\item{Sigma_xi_inv}{A two dimensional identity \code{matrix} of prior item parameter VC matrix.} +\item{Sigma_xi_inv}{A two dimensional identity \code{matrix} of prior item +parameter VC matrix.} \item{mu_theta}{The prior mean for theta.} @@ -29,7 +31,6 @@ Samples from posterior. Implement Gibbs 2PNO Sampler } \examples{ -require(fourPNO) # simulate small 2PNO dataset to demonstrate function J = 5 N = 100 @@ -54,8 +55,9 @@ Sigma_xi_inv = solve(2*matrix(c(1,0,0,1), 2, 2)) burnin = 1000 # Execute Gibbs sampler. This should take about 15.5 minutes -out_t <- Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta,Sigma_theta_inv,alpha_c,beta_c,alpha_s, - beta_s,burnin,rep(1,J),rep(1,J),gwg_reps=5,chain_length=burnin*2) +out_t = Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta,Sigma_theta_inv, + alpha_c,beta_c,alpha_s, beta_s,burnin, + rep(1,J),rep(1,J),gwg_reps=5,chain_length=burnin*2) # Summarizing posterior distribution OUT = cbind(apply(out_t$AS[,-c(1:burnin)],1,mean),apply(out_t$BS[,-c(1:burnin)],1,mean), diff --git a/man/Gibbs_4PNO.Rd b/man/Gibbs_4PNO.Rd index 0c1fe1a..2218d28 100644 --- a/man/Gibbs_4PNO.Rd +++ b/man/Gibbs_4PNO.Rd @@ -4,15 +4,18 @@ \alias{Gibbs_4PNO} \title{Gibbs Implementation of 4PNO} \usage{ -Gibbs_4PNO(Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, alpha_c, beta_c, - alpha_s, beta_s, burnin, cTF, sTF, gwg_reps, chain_length = 10000L) +Gibbs_4PNO(Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, alpha_c, + beta_c, alpha_s, beta_s, burnin, cTF, sTF, gwg_reps, + chain_length = 10000L) } \arguments{ \item{Y}{A N by J \code{matrix} of item responses.} -\item{mu_xi}{A two dimensional \code{vector} of prior item parameter means.} +\item{mu_xi}{A two dimensional \code{vector} of prior item parameter +means.} -\item{Sigma_xi_inv}{A two dimensional identity \code{matrix} of prior item parameter VC matrix.} +\item{Sigma_xi_inv}{A two dimensional identity \code{matrix} of prior item +parameter VC matrix.} \item{mu_theta}{The prior mean for theta.} @@ -28,11 +31,19 @@ Gibbs_4PNO(Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, alpha_c, beta_c, \item{burnin}{The number of MCMC samples to discard.} -\item{cTF}{A J dimensional \code{vector} indicating which lower asymptotes to estimate. 0 = exclude lower asymptote and 1 = include lower asymptote.} +\item{cTF}{A J dimensional \code{vector} indicating which +lower asymptotes to estimate. +0 = exclude lower asymptote and +1 = include lower asymptote.} -\item{sTF}{A J dimensional \code{vector} indicating which upper asymptotes to estimate. 0 = exclude upper asymptote and 1 = include upper asymptote.} +\item{sTF}{A J dimensional \code{vector} indicating which +upper asymptotes to estimate. +0 = exclude upper asymptote and +1 = include upper asymptote.} -\item{gwg_reps}{The number of Gibbs within Gibbs MCMC samples for marginal distribution of gamma. Values between 5 to 10 are adequate.} +\item{gwg_reps}{The number of Gibbs within Gibbs MCMC samples for +marginal distribution of gamma. Values between +5 to 10 are adequate.} \item{chain_length}{The number of MCMC samples.} } @@ -43,8 +54,6 @@ Samples from posterior. Internal function to -2LL } \examples{ -require(fourPNO) - # Simulate small 4PNO dataset to demonstrate function J = 5 N = 100 @@ -69,14 +78,20 @@ Sigma_xi_inv = solve(2*matrix(c(1,0,0,1),2,2)) burnin = 1000 # Execute Gibbs sampler -out_t = Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta,Sigma_theta_inv,alpha_c,beta_c,alpha_s, - beta_s,burnin,rep(1,J),rep(1,J),gwg_reps=5,chain_length=burnin*2) +out_t = Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta, + Sigma_theta_inv,alpha_c,beta_c,alpha_s, + beta_s,burnin,rep(1,J),rep(1,J), + gwg_reps=5,chain_length=burnin*2) # Summarizing posterior distribution -OUT = cbind(apply(out_t$AS[,-c(1:burnin)],1,mean),apply(out_t$BS[,-c(1:burnin)],1,mean), - apply(out_t$GS[,-c(1:burnin)],1,mean),apply(out_t$SS[,-c(1:burnin)],1,mean), - apply(out_t$AS[,-c(1:burnin)],1,sd),apply(out_t$BS[,-c(1:burnin)],1,sd), - apply(out_t$GS[,-c(1:burnin)],1,sd),apply(out_t$SS[,-c(1:burnin)],1,sd) ) +OUT = cbind(apply(out_t$AS[,-c(1:burnin)],1,mean), + apply(out_t$BS[,-c(1:burnin)],1,mean), + apply(out_t$GS[,-c(1:burnin)],1,mean), + apply(out_t$SS[,-c(1:burnin)],1,mean), + apply(out_t$AS[,-c(1:burnin)],1,sd), + apply(out_t$BS[,-c(1:burnin)],1,sd), + apply(out_t$GS[,-c(1:burnin)],1,sd), + apply(out_t$SS[,-c(1:burnin)],1,sd) ) OUT = cbind(1:J,OUT) colnames(OUT) = c('Item','as','bs','gs','ss','as_sd','bs_sd','gs_sd','ss_sd') diff --git a/man/Total_Tabulate.Rd b/man/Total_Tabulate.Rd index 92b9c4f..be90ba4 100644 --- a/man/Total_Tabulate.Rd +++ b/man/Total_Tabulate.Rd @@ -7,9 +7,9 @@ Total_Tabulate(N, J, Y) } \arguments{ -\item{N}{An \code{int}, which gives the number of observations. (> 0)} +\item{N}{An \code{int}, which gives the number of observations. (> 0)} -\item{J}{An \code{int}, which gives the number of items. (> 0)} +\item{J}{An \code{int}, which gives the number of items. (> 0)} \item{Y}{A N by J \code{matrix} of item responses.} } @@ -20,7 +20,7 @@ A vector of tabulated total scores. Internal function to -2LL } \seealso{ -\code{\link{Gibbs_4PNO}} +\code{\link[=Gibbs_4PNO]{Gibbs_4PNO()}} } \author{ Steven Andrew Culpepper diff --git a/man/Y_4pno_simulate.Rd b/man/Y_4pno_simulate.Rd index 23caa40..2be44db 100644 --- a/man/Y_4pno_simulate.Rd +++ b/man/Y_4pno_simulate.Rd @@ -7,9 +7,9 @@ Y_4pno_simulate(N, J, as, bs, gs, ss, theta) } \arguments{ -\item{N}{An \code{int}, which gives the number of observations. (> 0)} +\item{N}{An \code{int}, which gives the number of observations. (> 0)} -\item{J}{An \code{int}, which gives the number of items. (> 0)} +\item{J}{An \code{int}, which gives the number of items. (> 0)} \item{as}{A \code{vector} of item discrimination parameters.} @@ -28,7 +28,7 @@ A N by J \code{matrix} of dichotomous item responses. Generate item responses under the 4PNO } \seealso{ -\code{\link{Gibbs_4PNO}} +\code{\link[=Gibbs_4PNO]{Gibbs_4PNO()}} } \author{ Steven Andrew Culpepper diff --git a/man/fourPNO-package.Rd b/man/fourPNO-package.Rd index d5f11c3..64d38f7 100644 --- a/man/fourPNO-package.Rd +++ b/man/fourPNO-package.Rd @@ -6,9 +6,9 @@ \alias{fourPNO-package} \title{fourPNO: Bayesian 4 Parameter Item Response Model} \description{ -Estimate Barton & Lord's (1981) -four parameter IRT model with lower and upper asymptotes using Bayesian -formulation described by Culpepper (2016) . +Estimate Barton & Lord's (1981) + four parameter IRT model with lower and upper asymptotes using Bayesian + formulation described by Culpepper (2016) . } \seealso{ Useful links: @@ -19,6 +19,6 @@ Useful links: } \author{ -\strong{Maintainer}: Steven Andrew Culpepper \email{sculpepp@illinois.edu} +\strong{Maintainer}: Steven Andrew Culpepper \email{sculpepp@illinois.edu} (0000-0003-4226-6176) [copyright holder] } diff --git a/man/kappa_initialize.Rd b/man/kappa_initialize.Rd deleted file mode 100644 index fe5fb0f..0000000 --- a/man/kappa_initialize.Rd +++ /dev/null @@ -1,23 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{kappa_initialize} -\alias{kappa_initialize} -\title{Initialize Thresholds} -\usage{ -kappa_initialize(Ms) -} -\arguments{ -\item{Ms}{A \code{vector} with the number of scale values.} -} -\value{ -A \code{matrix} that is a Multivariate Normal distribution -} -\description{ -Internal function for initializing item thresholds. -} -\seealso{ -\code{\link{Gibbs_4PNO}} -} -\author{ -Steven Andrew Culpepper -} diff --git a/man/min2LL_4pno.Rd b/man/min2LL_4pno.Rd index be30caf..28c2cbe 100644 --- a/man/min2LL_4pno.Rd +++ b/man/min2LL_4pno.Rd @@ -30,7 +30,7 @@ min2LL_4pno(N, J, Y, as, bs, gs, ss, theta) Internal function to -2LL } \seealso{ -\code{\link{Gibbs_4PNO}} +\code{\link[=Gibbs_4PNO]{Gibbs_4PNO()}} } \author{ Steven Andrew Culpepper diff --git a/man/rmvnorm.Rd b/man/rmvnorm.Rd index 2049c84..655b5d7 100644 --- a/man/rmvnorm.Rd +++ b/man/rmvnorm.Rd @@ -9,18 +9,21 @@ rmvnorm(n, mu, sigma) \arguments{ \item{n}{An \code{int}, which gives the number of observations. (> 0)} -\item{mu}{A \code{vector} length m that represents the means of the normals.} +\item{mu}{A \code{vector} length m that represents the means of +the normals.} -\item{sigma}{A \code{matrix} with dimensions m x m that provides the covariance matrix.} +\item{sigma}{A \code{matrix} with dimensions m x m that provides the +covariance matrix.} } \value{ A \code{matrix} that is a Multivariate Normal distribution } \description{ -Creates a random Multivariate Normal when given number of obs, mean, and sigma. +Creates a random Multivariate Normal when given number of +obs, mean, and sigma. } \examples{ -#Call with the following data: +# Call with the following data: rmvnorm(2, c(0,0), diag(2)) } \author{ diff --git a/man/update_2pno.Rd b/man/update_2pno.Rd deleted file mode 100644 index 9cc8e7c..0000000 --- a/man/update_2pno.Rd +++ /dev/null @@ -1,46 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{update_2pno} -\alias{update_2pno} -\title{Update 2PNO Model Parameters} -\usage{ -update_2pno(N, J, Y, Z, as, bs, theta, Kaps, mu_xi, Sigma_xi_inv, mu_theta, - Sigma_theta_inv) -} -\arguments{ -\item{N}{The number of observations.} - -\item{J}{The number of items.} - -\item{Y}{A N by J \code{matrix} of item responses.} - -\item{Z}{A \code{matrix} N by J of continuous augmented data.} - -\item{as}{A \code{vector} of item discrimination parameters.} - -\item{bs}{A \code{vector} of item threshold parameters.} - -\item{theta}{A \code{vector} of prior thetas.} - -\item{Kaps}{A \code{matrix} for item thresholds (used for internal computations).} - -\item{mu_xi}{Prior mean for item parameters.} - -\item{Sigma_xi_inv}{Prior item parameter inverse variance-covariance matrix.} - -\item{mu_theta}{Prior mean for theta.} - -\item{Sigma_theta_inv}{Prior inverse variance for theta.} -} -\value{ -A list of item parameters. -} -\description{ -Internal function to update 2PNO parameters -} -\seealso{ -\code{\link{Gibbs_4PNO}} -} -\author{ -Steven Andrew Culpepper -} diff --git a/man/update_WKappaZ_NA.Rd b/man/update_WKappaZ_NA.Rd deleted file mode 100644 index 61c85da..0000000 --- a/man/update_WKappaZ_NA.Rd +++ /dev/null @@ -1,50 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{update_WKappaZ_NA} -\alias{update_WKappaZ_NA} -\title{Update Lower and Upper Asymptote Parameters of 4PNO} -\usage{ -update_WKappaZ_NA(Y, Ysum, Z, as, bs, gs, ss, theta, Kaps, alpha_c, beta_c, - alpha_s, beta_s, gwg_reps) -} -\arguments{ -\item{Y}{A N by J \code{matrix} of item responses.} - -\item{Ysum}{A \code{vector} of item total scores.} - -\item{Z}{A \code{matrix} N by J of continuous augmented data.} - -\item{as}{A \code{vector} of item discrimination parameters.} - -\item{bs}{A \code{vector} of item threshold parameters.} - -\item{gs}{A \code{vector} of item lower asymptote parameters.} - -\item{ss}{A \code{vector} of item upper asymptote parameters.} - -\item{theta}{A \code{vector} of prior thetas.} - -\item{Kaps}{A \code{matrix} for item thresholds (used for internal computations).} - -\item{alpha_c}{The lower asymptote prior 'a' parameter.} - -\item{beta_c}{The lower asymptote prior 'b' parameter.} - -\item{alpha_s}{The upper asymptote prior 'a' parameter.} - -\item{beta_s}{The upper asymptote prior 'b' parameter.} - -\item{gwg_reps}{The number of Gibbs within Gibbs MCMC samples for marginal distribution of gamma.} -} -\value{ -A list of item threshold parameters. -} -\description{ -Internal function to update item lower and upper asymptote -} -\seealso{ -\code{\link{Gibbs_4PNO}} -} -\author{ -Steven Andrew Culpepper -} diff --git a/man/update_ab_NA.Rd b/man/update_ab_NA.Rd deleted file mode 100644 index 03514d2..0000000 --- a/man/update_ab_NA.Rd +++ /dev/null @@ -1,37 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{update_ab_NA} -\alias{update_ab_NA} -\title{Update a and b Parameters of 2PNO, 3PNO, 4PNO} -\usage{ -update_ab_NA(N, J, Z, as, bs, theta, mu_xi, Sigma_xi_inv) -} -\arguments{ -\item{N}{An \code{int}, which gives the number of observations. (> 0)} - -\item{J}{An \code{int}, which gives the number of items. (> 0)} - -\item{Z}{A \code{matrix} N by J of continuous augmented data.} - -\item{as}{A \code{vector} of item discrimination parameters.} - -\item{bs}{A \code{vector} of item threshold parameters.} - -\item{theta}{A \code{vector} of prior thetas.} - -\item{mu_xi}{A two dimensional \code{vector} of prior item parameter means.} - -\item{Sigma_xi_inv}{A two dimensional identity \code{matrix} of prior item parameter VC matrix.} -} -\value{ -A list of item parameters. -} -\description{ -Update item slope and threshold -} -\seealso{ -\code{\link{Gibbs_4PNO}} -} -\author{ -Steven Andrew Culpepper -} diff --git a/man/update_ab_norestriction.Rd b/man/update_ab_norestriction.Rd deleted file mode 100644 index 4fe4dda..0000000 --- a/man/update_ab_norestriction.Rd +++ /dev/null @@ -1,37 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{update_ab_norestriction} -\alias{update_ab_norestriction} -\title{Update a and b Parameters of 4pno without alpha > 0 Restriction} -\usage{ -update_ab_norestriction(N, J, Z, as, bs, theta, mu_xi, Sigma_xi_inv) -} -\arguments{ -\item{N}{An \code{int}, which gives the number of observations. (> 0)} - -\item{J}{An \code{int}, which gives the number of items. (> 0)} - -\item{Z}{A \code{matrix} N by J of continuous augmented data.} - -\item{as}{A \code{vector} of item discrimination parameters.} - -\item{bs}{A \code{vector} of item threshold parameters.} - -\item{theta}{A \code{vector} of prior thetas.} - -\item{mu_xi}{A two dimensional \code{vector} of prior item parameter means.} - -\item{Sigma_xi_inv}{A two dimensional identity \code{matrix} of prior item parameter VC matrix.} -} -\value{ -A list of item parameters. -} -\description{ -Update item slope and threshold -} -\seealso{ -\code{\link{Gibbs_4PNO}} -} -\author{ -Steven Andrew Culpepper -} diff --git a/man/update_theta.Rd b/man/update_theta.Rd deleted file mode 100644 index cb5124c..0000000 --- a/man/update_theta.Rd +++ /dev/null @@ -1,35 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{update_theta} -\alias{update_theta} -\title{Internal Function for Updating Theta in Gibbs Sampler} -\usage{ -update_theta(N, Z, as, bs, theta, mu_theta, Sigma_theta_inv) -} -\arguments{ -\item{N}{An \code{int}, which gives the number of observations. (> 0)} - -\item{Z}{A \code{matrix} N by J of continuous augmented data.} - -\item{as}{A \code{vector} of item discrimination parameters.} - -\item{bs}{A \code{vector} of item threshold parameters.} - -\item{theta}{A \code{vector} of prior thetas.} - -\item{mu_theta}{The prior mean for theta.} - -\item{Sigma_theta_inv}{The prior inverse variance for theta.} -} -\value{ -A \code{vector} of thetas. -} -\description{ -Update theta in Gibbs sampler -} -\seealso{ -\code{\link{Gibbs_4PNO}} -} -\author{ -Steven Andrew Culpepper -} diff --git a/src/.clang-format b/src/.clang-format new file mode 100644 index 0000000..4c7303a --- /dev/null +++ b/src/.clang-format @@ -0,0 +1,26 @@ + +################################################################################ +# James Joseph Balamuta +# balamut2@illinois.edu +################################################################################ +# All cpp files should be formatted according to the style rules +# enforced by this file using clang-format. +# +# Official style options documentation: +# https://clang.llvm.org/docs/ClangFormatStyleOptions.html +# Style options with interactive examples: +# https://clangformat.com/ +################################################################################ +### Example use +# shopt -s extglob +# clang-format -i -style=file src/!(RcppExports).@(c|h|cpp) +# +################################################################################ +BasedOnStyle: LLVM +IndentWidth: 4 +AlignTrailingComments: true +UseTab: Never +BreakBeforeBraces: Linux +AllowShortIfStatementsOnASingleLine: false +IndentCaseLabels: false + diff --git a/src/Makevars b/src/Makevars index 22ebc63..332db97 100644 --- a/src/Makevars +++ b/src/Makevars @@ -1 +1,8 @@ -PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) +# Enable C++11 +CXX_STD=CXX11 + +# Enable the header file and OpenMP +PKG_CXXFLAGS=-I../inst/include $(SHLIB_OPENMP_CXXFLAGS) + +# Specify the required linking setup +PKG_LIBS=$(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/src/Makevars.win b/src/Makevars.win index 9542baf..332db97 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -1,2 +1,8 @@ +# Enable C++11 +CXX_STD=CXX11 -PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) +# Enable the header file and OpenMP +PKG_CXXFLAGS=-I../inst/include $(SHLIB_OPENMP_CXXFLAGS) + +# Specify the required linking setup +PKG_LIBS=$(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 1c78a85..7aeca70 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -19,94 +19,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// kappa_initialize -arma::mat kappa_initialize(const arma::vec& Ms); -RcppExport SEXP _fourPNO_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 -} -// update_theta -Rcpp::List update_theta(unsigned int N, const arma::mat& Z, const arma::vec& as, const arma::vec& bs, arma::vec& theta, const double& mu_theta, const double& Sigma_theta_inv); -RcppExport SEXP _fourPNO_update_theta(SEXP NSEXP, SEXP ZSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP thetaSEXP, SEXP mu_thetaSEXP, SEXP Sigma_theta_invSEXP) { -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::mat& >::type Z(ZSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type as(asSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type bs(bsSEXP); - Rcpp::traits::input_parameter< arma::vec& >::type theta(thetaSEXP); - Rcpp::traits::input_parameter< const double& >::type mu_theta(mu_thetaSEXP); - Rcpp::traits::input_parameter< const double& >::type Sigma_theta_inv(Sigma_theta_invSEXP); - rcpp_result_gen = Rcpp::wrap(update_theta(N, Z, as, bs, theta, mu_theta, Sigma_theta_inv)); - return rcpp_result_gen; -END_RCPP -} -// update_ab_NA -Rcpp::List update_ab_NA(unsigned int N, unsigned int J, const arma::mat& Z, arma::vec& as, arma::vec& bs, const arma::vec& theta, const arma::vec& mu_xi, const arma::mat& Sigma_xi_inv); -RcppExport SEXP _fourPNO_update_ab_NA(SEXP NSEXP, SEXP JSEXP, SEXP ZSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP thetaSEXP, SEXP mu_xiSEXP, SEXP Sigma_xi_invSEXP) { -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 Z(ZSEXP); - Rcpp::traits::input_parameter< arma::vec& >::type as(asSEXP); - Rcpp::traits::input_parameter< arma::vec& >::type bs(bsSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type theta(thetaSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type mu_xi(mu_xiSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type Sigma_xi_inv(Sigma_xi_invSEXP); - rcpp_result_gen = Rcpp::wrap(update_ab_NA(N, J, Z, as, bs, theta, mu_xi, Sigma_xi_inv)); - return rcpp_result_gen; -END_RCPP -} -// update_ab_norestriction -Rcpp::List update_ab_norestriction(unsigned int N, unsigned int J, const arma::mat& Z, arma::vec& as, arma::vec& bs, const arma::vec& theta, const arma::vec& mu_xi, const arma::mat& Sigma_xi_inv); -RcppExport SEXP _fourPNO_update_ab_norestriction(SEXP NSEXP, SEXP JSEXP, SEXP ZSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP thetaSEXP, SEXP mu_xiSEXP, SEXP Sigma_xi_invSEXP) { -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 Z(ZSEXP); - Rcpp::traits::input_parameter< arma::vec& >::type as(asSEXP); - Rcpp::traits::input_parameter< arma::vec& >::type bs(bsSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type theta(thetaSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type mu_xi(mu_xiSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type Sigma_xi_inv(Sigma_xi_invSEXP); - rcpp_result_gen = Rcpp::wrap(update_ab_norestriction(N, J, Z, as, bs, theta, mu_xi, Sigma_xi_inv)); - return rcpp_result_gen; -END_RCPP -} -// update_WKappaZ_NA -Rcpp::List update_WKappaZ_NA(const arma::mat& Y, const arma::vec& Ysum, arma::mat& Z, const arma::vec& as, const arma::vec& bs, const arma::vec& gs, const arma::vec& ss, const arma::vec& theta, const arma::mat& Kaps, double alpha_c, double beta_c, double alpha_s, double beta_s, unsigned int gwg_reps); -RcppExport SEXP _fourPNO_update_WKappaZ_NA(SEXP YSEXP, SEXP YsumSEXP, SEXP ZSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP gsSEXP, SEXP ssSEXP, SEXP thetaSEXP, SEXP KapsSEXP, SEXP alpha_cSEXP, SEXP beta_cSEXP, SEXP alpha_sSEXP, SEXP beta_sSEXP, SEXP gwg_repsSEXP) { -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::vec& >::type Ysum(YsumSEXP); - Rcpp::traits::input_parameter< arma::mat& >::type Z(ZSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type as(asSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type bs(bsSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type gs(gsSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type ss(ssSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type theta(thetaSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type Kaps(KapsSEXP); - Rcpp::traits::input_parameter< double >::type alpha_c(alpha_cSEXP); - Rcpp::traits::input_parameter< double >::type beta_c(beta_cSEXP); - Rcpp::traits::input_parameter< double >::type alpha_s(alpha_sSEXP); - Rcpp::traits::input_parameter< double >::type beta_s(beta_sSEXP); - Rcpp::traits::input_parameter< unsigned int >::type gwg_reps(gwg_repsSEXP); - rcpp_result_gen = Rcpp::wrap(update_WKappaZ_NA(Y, Ysum, Z, as, bs, gs, ss, theta, Kaps, alpha_c, beta_c, alpha_s, beta_s, gwg_reps)); - return rcpp_result_gen; -END_RCPP -} // min2LL_4pno double min2LL_4pno(unsigned int N, unsigned int J, const arma::mat& Y, const arma::vec& as, const arma::vec& bs, const arma::vec& gs, const arma::vec& ss, const arma::vec& theta); RcppExport SEXP _fourPNO_min2LL_4pno(SEXP NSEXP, SEXP JSEXP, SEXP YSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP gsSEXP, SEXP ssSEXP, SEXP thetaSEXP) { @@ -179,28 +91,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// update_2pno -Rcpp::List update_2pno(unsigned int N, unsigned int J, const arma::mat& Y, arma::mat& Z, const arma::vec& as, const arma::vec& bs, const arma::vec& theta, const arma::mat& Kaps, const arma::vec& mu_xi, const arma::mat& Sigma_xi_inv, const double& mu_theta, const double& Sigma_theta_inv); -RcppExport SEXP _fourPNO_update_2pno(SEXP NSEXP, SEXP JSEXP, SEXP YSEXP, SEXP ZSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP thetaSEXP, SEXP KapsSEXP, SEXP mu_xiSEXP, SEXP Sigma_xi_invSEXP, SEXP mu_thetaSEXP, SEXP Sigma_theta_invSEXP) { -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< arma::mat& >::type Z(ZSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type as(asSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type bs(bsSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type theta(thetaSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type Kaps(KapsSEXP); - Rcpp::traits::input_parameter< const arma::vec& >::type mu_xi(mu_xiSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type Sigma_xi_inv(Sigma_xi_invSEXP); - Rcpp::traits::input_parameter< const double& >::type mu_theta(mu_thetaSEXP); - Rcpp::traits::input_parameter< const double& >::type Sigma_theta_inv(Sigma_theta_invSEXP); - rcpp_result_gen = Rcpp::wrap(update_2pno(N, J, Y, Z, as, bs, theta, Kaps, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv)); - return rcpp_result_gen; -END_RCPP -} // Gibbs_2PNO Rcpp::List Gibbs_2PNO(const arma::mat& Y, const arma::vec& mu_xi, const arma::mat& Sigma_xi_inv, const double& mu_theta, const double& Sigma_theta_inv, unsigned int burnin, unsigned int chain_length); RcppExport SEXP _fourPNO_Gibbs_2PNO(SEXP YSEXP, SEXP mu_xiSEXP, SEXP Sigma_xi_invSEXP, SEXP mu_thetaSEXP, SEXP Sigma_theta_invSEXP, SEXP burninSEXP, SEXP chain_lengthSEXP) { @@ -221,16 +111,10 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_fourPNO_rmvnorm", (DL_FUNC) &_fourPNO_rmvnorm, 3}, - {"_fourPNO_kappa_initialize", (DL_FUNC) &_fourPNO_kappa_initialize, 1}, - {"_fourPNO_update_theta", (DL_FUNC) &_fourPNO_update_theta, 7}, - {"_fourPNO_update_ab_NA", (DL_FUNC) &_fourPNO_update_ab_NA, 8}, - {"_fourPNO_update_ab_norestriction", (DL_FUNC) &_fourPNO_update_ab_norestriction, 8}, - {"_fourPNO_update_WKappaZ_NA", (DL_FUNC) &_fourPNO_update_WKappaZ_NA, 14}, {"_fourPNO_min2LL_4pno", (DL_FUNC) &_fourPNO_min2LL_4pno, 8}, {"_fourPNO_Y_4pno_simulate", (DL_FUNC) &_fourPNO_Y_4pno_simulate, 7}, {"_fourPNO_Total_Tabulate", (DL_FUNC) &_fourPNO_Total_Tabulate, 3}, {"_fourPNO_Gibbs_4PNO", (DL_FUNC) &_fourPNO_Gibbs_4PNO, 14}, - {"_fourPNO_update_2pno", (DL_FUNC) &_fourPNO_update_2pno, 12}, {"_fourPNO_Gibbs_2PNO", (DL_FUNC) &_fourPNO_Gibbs_2PNO, 7}, {NULL, NULL, 0} }; diff --git a/src/fourpno_101315.cpp b/src/fourpno_101315.cpp index 9907379..d3b6460 100644 --- a/src/fourpno_101315.cpp +++ b/src/fourpno_101315.cpp @@ -1,459 +1,574 @@ #include //' Generate Random Multivariate Normal Distribution -//' -//' Creates a random Multivariate Normal when given number of obs, mean, and sigma. -//' @usage rmvnorm(n, mu, 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 sigma A \code{matrix} with dimensions m x m that provides the covariance matrix. -//' @return A \code{matrix} that is a Multivariate Normal distribution -//' @author James J Balamuta -//' @examples -//' #Call with the following data: -//' rmvnorm(2, c(0,0), diag(2)) +//' +//' Creates a random Multivariate Normal when given number of +//' obs, mean, and sigma. +//' +//' @param n An `int`, which gives the number of observations. (> 0) +//' @param mu A `vector` length m that represents the means of +//' the normals. +//' @param sigma A `matrix` with dimensions m x m that provides the +//' covariance matrix. +//' +//' @return +//' A `matrix` that is a Multivariate Normal distribution +//' +//' @author +//' James J Balamuta +//' //' @export +//' @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& sigma) { - unsigned int ncols = sigma.n_cols; - Rcpp::RNGScope scope; - arma::mat Y(n, ncols); - Y.imbue( norm_rand ) ; - return arma::repmat(mu, 1, n).t() + Y * arma::chol(sigma); +arma::mat rmvnorm(unsigned int n, const arma::vec &mu, const arma::mat &sigma) +{ + unsigned int ncols = sigma.n_cols; + Rcpp::RNGScope scope; + arma::mat Y(n, ncols); + Y.imbue(norm_rand); + return arma::repmat(mu, 1, n).t() + Y * arma::chol(sigma); } //' Initialize Thresholds -//' +//' //' Internal function for initializing item thresholds. -//' @param Ms A \code{vector} with the number of scale values. -//' @return A \code{matrix} that is a Multivariate Normal distribution -//' @seealso \code{\link{Gibbs_4PNO}} -//' @author Steven Andrew Culpepper -//' @export -// [[Rcpp::export]] -arma::mat kappa_initialize(const arma::vec& Ms){ - unsigned int J = Ms.n_elem; - unsigned int M = max(Ms); - arma::mat KAP0(M+1,J); - - (KAP0.row(0)).fill(-arma::datum::inf); - KAP0.row(1).fill(0.0); - - for(unsigned int j=0;jMs(j)){ - KAP0(m,j)=arma::datum::nan; - } - - } - - } - - return KAP0; +//' +//' @param Ms A `vector` with the number of scale values. +//' +//' @return +//' A `matrix` that is a Multivariate Normal distribution +//' +//' @author +//' Steven Andrew Culpepper +//' +//' @seealso +//' [Gibbs_4PNO()] +//' +//' @noRd +arma::mat kappa_initialize(const arma::vec &Ms) +{ + unsigned int J = Ms.n_elem; + unsigned int M = max(Ms); + arma::mat KAP0(M + 1, J); + + (KAP0.row(0)).fill(-arma::datum::inf); + KAP0.row(1).fill(0.0); + + for (unsigned int j = 0; j < J; j++) { + for (unsigned int m = 2; m < M + 1; m++) { + + if (m < Ms(j)) { + KAP0(m, j) = (m - 1) * 2; + } + + if (m == Ms(j)) { + KAP0(m, j) = arma::datum::inf; + } + + if (m > Ms(j)) { + KAP0(m, j) = arma::datum::nan; + } + } + } + + return KAP0; } //' Internal Function for Updating Theta in Gibbs Sampler -//' +//' //' Update theta in Gibbs sampler -//' @param N An \code{int}, which gives the number of observations. (> 0) -//' @param Z A \code{matrix} N by J of continuous augmented data. -//' @param as A \code{vector} of item discrimination parameters. -//' @param bs A \code{vector} of item threshold parameters. -//' @param theta A \code{vector} of prior thetas. -//' @param mu_theta The prior mean for theta. +//' +//' @param N An `int`, which gives the number of observations. +//' (> 0) +//' @param Z A `matrix` N by J of continuous augmented data. +//' @param as A `vector` of item discrimination parameters. +//' @param bs A `vector` of item threshold parameters. +//' @param theta A `vector` of prior thetas. +//' @param mu_theta The prior mean for theta. //' @param Sigma_theta_inv The prior inverse variance for theta. -//' @return A \code{vector} of thetas. -//' @seealso \code{\link{Gibbs_4PNO}} -//' @author Steven Andrew Culpepper -//' @export -// [[Rcpp::export]] -Rcpp::List update_theta(unsigned int N,const arma::mat& Z,const arma::vec& as,const arma::vec& bs, - arma::vec& theta,const double& mu_theta,const double& Sigma_theta_inv) { - - double apb=dot(as,bs); - arma::vec oneN = arma::ones(N); - double vartheta = 1.0/(dot(as,as)+Sigma_theta_inv); -// apb = dot(as,bs); - arma::vec zetas(N,arma::fill::randn); -// zetas.imbue( norm_rand ); - - theta = zetas*sqrt(vartheta) + vartheta*(Z*as + oneN*(apb+mu_theta*Sigma_theta_inv)); - -return Rcpp::List::create(//Rcpp::Named("theta",theta), - Rcpp::Named("apb",apb) -// Rcpp::Named("ajIakIna",ajIakIna), -// Rcpp::Named("AZ",AZ), -// Rcpp::Named("Sigma_theta_star",Sigma_theta_star), -// Rcpp::Named("mu_theta_star",mu_theta_star) - ); +//' +//' @return +//' A `vector` of thetas. +//' +//' @author +//' Steven Andrew Culpepper +//' +//' @seealso +//' [Gibbs_4PNO()] +//' +//' @noRd +Rcpp::List update_theta(unsigned int N, const arma::mat &Z, const arma::vec &as, + const arma::vec &bs, arma::vec &theta, + const double &mu_theta, const double &Sigma_theta_inv) +{ + + double apb = dot(as, bs); + arma::vec oneN = arma::ones(N); + double vartheta = 1.0 / (dot(as, as) + Sigma_theta_inv); + // apb = dot(as,bs); + arma::vec zetas(N, arma::fill::randn); + // zetas.imbue( norm_rand ); + + theta = zetas * sqrt(vartheta) + + vartheta * (Z * as + oneN * (apb + mu_theta * Sigma_theta_inv)); + + return Rcpp::List::create( // Rcpp::Named("theta",theta), + Rcpp::Named("apb", apb) + // Rcpp::Named("ajIakIna",ajIakIna), + // Rcpp::Named("AZ",AZ), + // Rcpp::Named("Sigma_theta_star",Sigma_theta_star), + // Rcpp::Named("mu_theta_star",mu_theta_star) + ); } - - //' Update a and b Parameters of 2PNO, 3PNO, 4PNO -//' +//' //' Update item slope and threshold -//' @param N An \code{int}, which gives the number of observations. (> 0) -//' @param J An \code{int}, which gives the number of items. (> 0) -//' @param Z A \code{matrix} N by J of continuous augmented data. -//' @param as A \code{vector} of item discrimination parameters. -//' @param bs A \code{vector} of item threshold parameters. -//' @param theta A \code{vector} of prior thetas. -//' @param mu_xi A two dimensional \code{vector} of prior item parameter means. -//' @param Sigma_xi_inv A two dimensional identity \code{matrix} of prior item parameter VC matrix. -//' @return A list of item parameters. -//' @seealso \code{\link{Gibbs_4PNO}} -//' @author Steven Andrew Culpepper -//' @export -// [[Rcpp::export]] -Rcpp::List update_ab_NA(unsigned int N,unsigned int J,const arma::mat& Z,arma::vec& as,arma::vec& bs, - const arma::vec& theta,const arma::vec& mu_xi,const arma::mat& Sigma_xi_inv) { - - arma::vec m1(N); - m1.fill(-1.0); - arma::mat X_theta = join_rows(theta,m1); - arma::mat XX_xi = X_theta.t()*X_theta; - arma::mat XZ = X_theta.t()*Z; - - arma::mat Sig_xi_star = inv(Sigma_xi_inv + XX_xi); - arma::vec mu_xi_star(2); - double pa, mb_a, vb_a,ua,ub; - - for(unsigned j=0;j 0) +//' @param J An `int`, which gives the number of items. (> 0) +//' @param Z A `matrix` N by J of continuous augmented data. +//' @param as A `vector` of item discrimination parameters. +//' @param bs A `vector` of item threshold parameters. +//' @param theta A `vector` of prior thetas. +//' @param mu_xi A two dimensional `vector` of prior item parameter +//means. ' @param Sigma_xi_inv A two dimensional identity `matrix` of prior item +//' parameter VC matrix. +//' +//' @return +//' A `list` of item parameters. +//' +//' @author +//' Steven Andrew Culpepper +//' +//' @seealso +//' [Gibbs_4PNO()] +//' +//' @noRd +Rcpp::List update_ab_NA(unsigned int N, unsigned int J, const arma::mat &Z, + arma::vec &as, arma::vec &bs, const arma::vec &theta, + const arma::vec &mu_xi, const arma::mat &Sigma_xi_inv) +{ + + arma::vec m1(N); + m1.fill(-1.0); + arma::mat X_theta = join_rows(theta, m1); + arma::mat XX_xi = X_theta.t() * X_theta; + arma::mat XZ = X_theta.t() * Z; + + arma::mat Sig_xi_star = inv(Sigma_xi_inv + XX_xi); + arma::vec mu_xi_star(2); + double pa, mb_a, vb_a, ua, ub; + + for (unsigned j = 0; j < J; j++) { + + mu_xi_star = Sig_xi_star * (Sigma_xi_inv * mu_xi + XZ.col(j)); + + ua = R::runif(0.0, 1.0); + pa = R::pnorm(0.0, mu_xi_star(0), sqrt(Sig_xi_star(0, 0)), 1, 0); + as(j) = R::qnorm(pa + ua * (1.0 - pa), mu_xi_star(0), + sqrt(Sig_xi_star(0, 0)), 1, 0); + + ub = R::runif(0.0, 1.0); + mb_a = mu_xi_star(1) + + Sig_xi_star(0, 1) / Sig_xi_star(1, 1) * (as(j) - mu_xi_star(0)); + vb_a = Sig_xi_star(1, 1) - + Sig_xi_star(0, 1) * Sig_xi_star(0, 1) / Sig_xi_star(0, 0); + bs(j) = R::qnorm(ub, mb_a, sqrt(vb_a), 1, 0); + } + + return Rcpp::List::create( // Rcpp::Named("XX_xi")=XX_xi, + // Rcpp::Named("XZ")=XZ, + // Rcpp::Named("as")=as, + // Rcpp::Named("bs")=bs, + Rcpp::Named("pa") = pa //, + // Rcpp::Named("ua")=ua, + // Rcpp::Named("ub")=ub, + // Rcpp::Named("mb_a")=mb_a, + // Rcpp::Named("vb_a")=vb_a, + // Rcpp::Named("Sig_xi_star")=Sig_xi_star, + // Rcpp::Named("mu_xi_star")=mu_xi_star + ); } - - - - //' Update a and b Parameters of 4pno without alpha > 0 Restriction -//' +//' //' Update item slope and threshold -//' @param N An \code{int}, which gives the number of observations. (> 0) -//' @param J An \code{int}, which gives the number of items. (> 0) -//' @param Z A \code{matrix} N by J of continuous augmented data. -//' @param as A \code{vector} of item discrimination parameters. -//' @param bs A \code{vector} of item threshold parameters. -//' @param theta A \code{vector} of prior thetas. -//' @param mu_xi A two dimensional \code{vector} of prior item parameter means. -//' @param Sigma_xi_inv A two dimensional identity \code{matrix} of prior item parameter VC matrix. -//' @return A list of item parameters. -//' @seealso \code{\link{Gibbs_4PNO}} -//' @author Steven Andrew Culpepper -//' @export -// [[Rcpp::export]] -Rcpp::List update_ab_norestriction(unsigned int N,unsigned int J,const arma::mat& Z,arma::vec& as,arma::vec& bs, - const arma::vec& theta,const arma::vec& mu_xi,const arma::mat& Sigma_xi_inv) { - - arma::vec m1(N); - m1.fill(-1.0); - arma::mat X_theta = join_rows(theta,m1); - arma::mat XX_xi = X_theta.t()*X_theta; - arma::mat XZ = X_theta.t()*Z; - - arma::mat Sig_xi_star = inv(Sigma_xi_inv + XX_xi); - arma::vec mu_xi_star(2); - double pa, mb_a, vb_a,ua,ub; - - for(unsigned j=0;j 0) +//' @param J An `int`, which gives the number of items. (> 0) +//' @param Z A `matrix` N by J of continuous augmented data. +//' @param as A `vector` of item discrimination parameters. +//' @param bs A `vector` of item threshold parameters. +//' @param theta A `vector` of prior thetas. +//' @param mu_xi A two dimensional `vector` of prior item parameter +//means. ' @param Sigma_xi_inv A two dimensional identity `matrix` of prior ' +//item parameter VC matrix. +//' +//' @return +//' A `list` of item parameters. +//' +//' @author +//' Steven Andrew Culpepper +//' +//' @seealso +//' [Gibbs_4PNO()] +//' +//' @noRd +Rcpp::List update_ab_norestriction(unsigned int N, unsigned int J, + const arma::mat &Z, arma::vec &as, + arma::vec &bs, const arma::vec &theta, + const arma::vec &mu_xi, + const arma::mat &Sigma_xi_inv) +{ + + arma::vec m1(N); + m1.fill(-1.0); + arma::mat X_theta = join_rows(theta, m1); + arma::mat XX_xi = X_theta.t() * X_theta; + arma::mat XZ = X_theta.t() * Z; + + arma::mat Sig_xi_star = inv(Sigma_xi_inv + XX_xi); + arma::vec mu_xi_star(2); + double pa, mb_a, vb_a, ua, ub; + + for (unsigned j = 0; j < J; j++) { + mu_xi_star = Sig_xi_star * (Sigma_xi_inv * mu_xi + XZ.col(j)); + + ua = R::runif(0.0, 1.0); + pa = 0; + as(j) = R::qnorm(pa + ua * (1.0 - pa), mu_xi_star(0), + sqrt(Sig_xi_star(0, 0)), 1, 0); + + ub = R::runif(0.0, 1.0); + mb_a = mu_xi_star(1) + + Sig_xi_star(0, 1) / Sig_xi_star(1, 1) * (as(j) - mu_xi_star(0)); + vb_a = Sig_xi_star(1, 1) - + Sig_xi_star(0, 1) * Sig_xi_star(0, 1) / Sig_xi_star(0, 0); + bs(j) = R::qnorm(ub, mb_a, sqrt(vb_a), 1, 0); + } + + return Rcpp::List::create( // Rcpp::Named("XX_xi")=XX_xi, + Rcpp::Named("ua") = ua); } //' Update Lower and Upper Asymptote Parameters of 4PNO -//' +//' //' Internal function to update item lower and upper asymptote -//' @param Y A N by J \code{matrix} of item responses. -//' @param Ysum A \code{vector} of item total scores. -//' @param Z A \code{matrix} N by J of continuous augmented data. -//' @param as A \code{vector} of item discrimination parameters. -//' @param bs A \code{vector} of item threshold parameters. -//' @param gs A \code{vector} of item lower asymptote parameters. -//' @param ss A \code{vector} of item upper asymptote parameters. -//' @param theta A \code{vector} of prior thetas. -//' @param Kaps A \code{matrix} for item thresholds (used for internal computations). -//' @param alpha_c The lower asymptote prior 'a' parameter. -//' @param beta_c The lower asymptote prior 'b' parameter. -//' @param alpha_s The upper asymptote prior 'a' parameter. -//' @param beta_s The upper asymptote prior 'b' parameter. -//' @param gwg_reps The number of Gibbs within Gibbs MCMC samples for marginal distribution of gamma. -//' @return A list of item threshold parameters. -//' @seealso \code{\link{Gibbs_4PNO}} -//' @author Steven Andrew Culpepper -//' @export -// [[Rcpp::export]] -Rcpp::List update_WKappaZ_NA(const arma::mat& Y,const arma::vec& Ysum,arma::mat& Z,const arma::vec& as, - const arma::vec& bs,const arma::vec& gs,const arma::vec& ss,const arma::vec& theta,const arma::mat& Kaps, - double alpha_c,double beta_c,double alpha_s,double beta_s,unsigned int gwg_reps){ - - unsigned int N = Y.n_rows; - unsigned int J = Y.n_cols; - - arma::vec eta(N); - - arma::vec p4pno(2); - arma::mat W=arma::zeros(N,J); - arma::vec oneN = arma::ones(N); - double us,ug, uZ, u4pno; - double Phi_eta; - double pg_temp,ps_temp,s_temp,g_temp; - arma::vec gs_new=arma::zeros(J); - arma::vec ss_new=arma::zeros(J); - double n1dot,n0dot,n01,n10; - double ps,pg; - - for(unsigned int j=0;j u4pno and Y(i,j)==1.0){ - W(i,j) = 1.0; +//' +//' @param Y A N by J `matrix` of item responses. +//' @param Ysum A `vector` of item total scores. +//' @param Z A `matrix` N by J of continuous augmented data. +//' @param as A `vector` of item discrimination parameters. +//' @param bs A `vector` of item threshold parameters. +//' @param gs A `vector` of item lower asymptote parameters. +//' @param ss A `vector` of item upper asymptote parameters. +//' @param theta A `vector` of prior thetas. +//' @param Kaps A `matrix` for item thresholds +//' (used for internal computations). +//' @param alpha_c The lower asymptote prior 'a' parameter. +//' @param beta_c The lower asymptote prior 'b' parameter. +//' @param alpha_s The upper asymptote prior 'a' parameter. +//' @param beta_s The upper asymptote prior 'b' parameter. +//' @param gwg_reps The number of Gibbs within Gibbs MCMC samples for +//' marginal distribution of gamma. +//' +//' @return +//' A `list` of item threshold parameters. +//' +//' @author +//' Steven Andrew Culpepper +//' +//' @seealso +//' [Gibbs_4PNO()] +//' +//' @noRd +Rcpp::List update_WKappaZ_NA(const arma::mat &Y, const arma::vec &Ysum, + arma::mat &Z, const arma::vec &as, + const arma::vec &bs, const arma::vec &gs, + const arma::vec &ss, const arma::vec &theta, + const arma::mat &Kaps, double alpha_c, + double beta_c, double alpha_s, double beta_s, + unsigned int gwg_reps) +{ + + unsigned int N = Y.n_rows; + unsigned int J = Y.n_cols; + + arma::vec eta(N); + + arma::vec p4pno(2); + arma::mat W = arma::zeros(N, J); + arma::vec oneN = arma::ones(N); + double us, ug, uZ, u4pno; + double Phi_eta; + double pg_temp, ps_temp, s_temp, g_temp; + arma::vec gs_new = arma::zeros(J); + arma::vec ss_new = arma::zeros(J); + double n1dot, n0dot, n01, n10; + double ps, pg; + + for (unsigned int j = 0; j < J; j++) { + + // Calculate eta_ij for all of the models and individuals with data + eta = as(j) * theta - oneN * bs(j); + + //////////////////////////////////////////////// + // 4PNO model + //////////////////////////////////////////////// + // if(model(j)=="4PNO"){ + + for (unsigned int i = 0; i < N; i++) { + + // if(MISS(i,j)==1.0){ + uZ = R::runif(0.0, 1.0); + u4pno = R::runif(0.0, 1.0); + Phi_eta = R::pnorm(eta(i), 0.0, 1.0, 1, 0); + + if ((1.0 - ss(j)) * Phi_eta / + (gs(j) + (1.0 - ss(j) - gs(j)) * Phi_eta) > + u4pno and + Y(i, j) == 1.0) { + W(i, j) = 1.0; } - - if(ss(j)*Phi_eta/( 1.0-gs(j) - (1.0-ss(j)-gs(j))*Phi_eta) > u4pno and Y(i,j)==0.0){ - W(i,j)=1.0; + + if (ss(j) * Phi_eta / + (1.0 - gs(j) - (1.0 - ss(j) - gs(j)) * Phi_eta) > + u4pno and + Y(i, j) == 0.0) { + W(i, j) = 1.0; } -// sj = (1-W(i))*Y(i,j) + sj; - - p4pno(0) = R::pnorm(Kaps(W(i,j) ,j) ,eta(i),1.0, 1, 0); - p4pno(1) = R::pnorm(Kaps(W(i,j)+1.0,j) ,eta(i),1.0, 1, 0); - Z(i,j) = R::qnorm(p4pno(0) + uZ*( p4pno(1)-p4pno(0) ),eta(i),1.0,1,0) ; -// } - } - - //update gs and ss - us=R::runif(0,1); - ug=R::runif(0,1); - - n1dot = sum(W.col(j)); - n0dot = N-n1dot; - n01 = Ysum(j) - dot(W.col(j),Y.col(j) ); - n10 = n1dot - dot(W.col(j),Y.col(j)); - - //draw g from marginal distribution w Gibbs within Gibbs - s_temp = R::rbeta(n10+alpha_s,n1dot-n10+beta_s); - for(unsigned int r =0;r 0) -//' @param J An \code{int}, which gives the number of items. (> 0) -//' @param Y A N by J \code{matrix} of item responses. -//' @param as A \code{vector} of item discrimination parameters. -//' @param bs A \code{vector} of item threshold parameters. -//' @param gs A \code{vector} of item lower asymptote parameters. -//' @param ss A \code{vector} of item upper asymptote parameters. -//' @param theta A \code{vector} of prior thetas. -//' @return -2LL. -//' @seealso \code{\link{Gibbs_4PNO}} -//' @author Steven Andrew Culpepper +//' +//' @param N An `int`, which gives the number of observations. (> 0) +//' @param J An `int`, which gives the number of items. (> 0) +//' @param Y A N by J `matrix` of item responses. +//' @param as A `vector` of item discrimination parameters. +//' @param bs A `vector` of item threshold parameters. +//' @param gs A `vector` of item lower asymptote parameters. +//' @param ss A `vector` of item upper asymptote parameters. +//' @param theta A `vector` of prior thetas. +//' +//' @return +//' -2LL. +//' @author +//' Steven Andrew Culpepper +//' +//' @seealso +//' [Gibbs_4PNO()] +//' //' @export // [[Rcpp::export]] -double min2LL_4pno(unsigned int N,unsigned int J,const arma::mat& Y,const arma::vec& as, - const arma::vec& bs,const arma::vec& gs,const arma::vec& ss,const arma::vec& theta){ - - double m2ll=0.0; - double eta,p_eta; - - for(unsigned int j=0;j 0) -//' @param J An \code{int}, which gives the number of items. (> 0) -//' @param as A \code{vector} of item discrimination parameters. -//' @param bs A \code{vector} of item threshold parameters. -//' @param gs A \code{vector} of item lower asymptote parameters. -//' @param ss A \code{vector} of item upper asymptote parameters. -//' @param theta A \code{vector} of prior thetas. -//' @return A N by J \code{matrix} of dichotomous item responses. -//' @seealso \code{\link{Gibbs_4PNO}} -//' @author Steven Andrew Culpepper +//' +//' @param N An `int`, which gives the number of observations. (> 0) +//' @param J An `int`, which gives the number of items. (> 0) +//' @param as A `vector` of item discrimination parameters. +//' @param bs A `vector` of item threshold parameters. +//' @param gs A `vector` of item lower asymptote parameters. +//' @param ss A `vector` of item upper asymptote parameters. +//' @param theta A `vector` of prior thetas. +//' +//' @return +//' A N by J `matrix` of dichotomous item responses. +//' +//' @author +//' Steven Andrew Culpepper +//' +//' @seealso +//' [Gibbs_4PNO()] +//' //' @export // [[Rcpp::export]] -arma::mat Y_4pno_simulate(unsigned int N,unsigned int J,const arma::vec& as,const arma::vec& bs, - const arma::vec& gs,const arma::vec& ss,const arma::vec& theta){ - - double u,eta,p_eta; - arma::mat Ysim(N,J); - - for(unsigned int j=0;j u ){ - Ysim(i,j) = 1.0; +arma::mat Y_4pno_simulate(unsigned int N, unsigned int J, const arma::vec &as, + const arma::vec &bs, const arma::vec &gs, + const arma::vec &ss, const arma::vec &theta) +{ + + double u, eta, p_eta; + arma::mat Ysim(N, J); + + for (unsigned int j = 0; j < J; j++) { + for (unsigned int i = 0; i < N; i++) { + + u = R::runif(0.0, 1.0); + eta = as(j) * theta(i) - bs(j); + p_eta = R::pnorm(eta, 0.0, 1.0, 1, 0); + + if (gs(j) + (1.0 - ss(j) - gs(j)) * p_eta > u) { + Ysim(i, j) = 1.0; } - - else{ - Ysim(i,j)=.0; + + else { + Ysim(i, j) = .0; } - } - } - return Ysim; + } + } + return Ysim; } -//' Calculate Tabulated Total Scores -//' +//' Calculate Tabulated Total Scores +//' //' Internal function to -2LL -//' @param N An \code{int}, which gives the number of observations. (> 0) -//' @param J An \code{int}, which gives the number of items. (> 0) -//' @param Y A N by J \code{matrix} of item responses. -//' @return A vector of tabulated total scores. -//' @seealso \code{\link{Gibbs_4PNO}} -//' @author Steven Andrew Culpepper +//' +//' @param N An `int`, which gives the number of observations. (> 0) +//' @param J An `int`, 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 +//' +//' @seealso +//' [Gibbs_4PNO()] +//' //' @export // [[Rcpp::export]] -arma::uvec Total_Tabulate(unsigned int N,unsigned int J,const arma::mat Y){ - - arma::vec T(N); -// arma::vec H(J); - T = sum(Y,1); - arma::uvec H = arma::hist(T,arma::linspace(0,J,J+1)); - - return H; +arma::uvec Total_Tabulate(unsigned int N, unsigned int J, const arma::mat Y) +{ + + arma::vec T(N); + // arma::vec H(J); + T = sum(Y, 1); + arma::uvec H = arma::hist(T, arma::linspace(0, J, J + 1)); + + return H; } //' Gibbs Implementation of 4PNO -//' +//' //' Internal function to -2LL -//' @param Y A N by J \code{matrix} of item responses. -//' @param mu_xi A two dimensional \code{vector} of prior item parameter means. -//' @param Sigma_xi_inv A two dimensional identity \code{matrix} of prior item parameter VC matrix. +//' +//' @param Y A N by J `matrix` of item responses. +//' @param mu_xi A two dimensional `vector` of prior item parameter +//' means. +//' @param Sigma_xi_inv A two dimensional identity `matrix` of prior item +//' parameter VC matrix. //' @param mu_theta The prior mean for theta. //' @param Sigma_theta_inv The prior inverse variance for theta. -//' @param alpha_c The lower asymptote prior 'a' parameter. -//' @param beta_c The lower asymptote prior 'b' parameter. -//' @param alpha_s The upper asymptote prior 'a' parameter. -//' @param beta_s The upper asymptote prior 'b' parameter. -//' @param burnin The number of MCMC samples to discard. -//' @param cTF A J dimensional \code{vector} indicating which lower asymptotes to estimate. 0 = exclude lower asymptote and 1 = include lower asymptote. -//' @param sTF A J dimensional \code{vector} indicating which upper asymptotes to estimate. 0 = exclude upper asymptote and 1 = include upper asymptote. -//' @param gwg_reps The number of Gibbs within Gibbs MCMC samples for marginal distribution of gamma. Values between 5 to 10 are adequate. -//' @param chain_length The number of MCMC samples. -//' @return Samples from posterior. -//' @author Steven Andrew Culpepper +//' @param alpha_c The lower asymptote prior 'a' parameter. +//' @param beta_c The lower asymptote prior 'b' parameter. +//' @param alpha_s The upper asymptote prior 'a' parameter. +//' @param beta_s The upper asymptote prior 'b' parameter. +//' @param burnin The number of MCMC samples to discard. +//' @param cTF A J dimensional `vector` indicating which +//' lower asymptotes to estimate. +//' 0 = exclude lower asymptote and +//' 1 = include lower asymptote. +//' @param sTF A J dimensional `vector` indicating which +//' upper asymptotes to estimate. +//' 0 = exclude upper asymptote and +//' 1 = include upper asymptote. +//' @param gwg_reps The number of Gibbs within Gibbs MCMC samples for +//' marginal distribution of gamma. Values between +//' 5 to 10 are adequate. +//' @param chain_length The number of MCMC samples. +//' +//' @return +//' Samples from posterior. +//' +//' @author +//' Steven Andrew Culpepper +//' //' @export //' @examples -//' require(fourPNO) -//' //' # Simulate small 4PNO dataset to demonstrate function //' J = 5 //' N = 100 -//' +//' //' # Population item parameters //' as_t = rnorm(J,mean=2,sd=.5) //' bs_t = rnorm(J,mean=0,sd=.5) -//' +//' //' # Sampling gs and ss with truncation //' gs_t = rbeta(J,1,8) //' ps_g = pbeta(1-gs_t,1,8) //' ss_t = qbeta(runif(J)*ps_g,1,8) //' theta_t <- rnorm(N) //' Y_t = Y_4pno_simulate(N,J,as=as_t,bs=bs_t,gs=gs_t,ss=ss_t,theta=theta_t) -//' +//' //' # Setting prior parameters //' mu_theta=0 //' Sigma_theta_inv=1 @@ -461,213 +576,245 @@ arma::uvec Total_Tabulate(unsigned int N,unsigned int J,const arma::mat Y){ //' alpha_c=alpha_s=beta_c=beta_s=1 //' Sigma_xi_inv = solve(2*matrix(c(1,0,0,1),2,2)) //' burnin = 1000 -//' +//' //' # Execute Gibbs sampler -//' out_t = Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta,Sigma_theta_inv,alpha_c,beta_c,alpha_s, -//' beta_s,burnin,rep(1,J),rep(1,J),gwg_reps=5,chain_length=burnin*2) -//' +//' out_t = Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta, +//' Sigma_theta_inv,alpha_c,beta_c,alpha_s, +//' beta_s,burnin,rep(1,J),rep(1,J), +//' gwg_reps=5,chain_length=burnin*2) +//' //' # Summarizing posterior distribution -//' OUT = cbind(apply(out_t$AS[,-c(1:burnin)],1,mean),apply(out_t$BS[,-c(1:burnin)],1,mean), -//' apply(out_t$GS[,-c(1:burnin)],1,mean),apply(out_t$SS[,-c(1:burnin)],1,mean), -//' apply(out_t$AS[,-c(1:burnin)],1,sd),apply(out_t$BS[,-c(1:burnin)],1,sd), -//' apply(out_t$GS[,-c(1:burnin)],1,sd),apply(out_t$SS[,-c(1:burnin)],1,sd) ) -//' +//' OUT = cbind(apply(out_t$AS[,-c(1:burnin)],1,mean), +//' apply(out_t$BS[,-c(1:burnin)],1,mean), +//' apply(out_t$GS[,-c(1:burnin)],1,mean), +//' apply(out_t$SS[,-c(1:burnin)],1,mean), +//' apply(out_t$AS[,-c(1:burnin)],1,sd), +//' apply(out_t$BS[,-c(1:burnin)],1,sd), +//' apply(out_t$GS[,-c(1:burnin)],1,sd), +//' apply(out_t$SS[,-c(1:burnin)],1,sd) ) +//' //' OUT = cbind(1:J,OUT) -//' colnames(OUT) = c('Item','as','bs','gs','ss','as_sd','bs_sd','gs_sd','ss_sd') -//' print(OUT,digits=3) +//' colnames(OUT) = +//c('Item','as','bs','gs','ss','as_sd','bs_sd','gs_sd','ss_sd') ' +//print(OUT,digits=3) // [[Rcpp::export]] -Rcpp::List Gibbs_4PNO(const arma::mat& Y,const arma::vec& mu_xi,const arma::mat& Sigma_xi_inv, - const double& mu_theta,const double& Sigma_theta_inv,double alpha_c,double beta_c,double alpha_s, - double beta_s,unsigned int burnin,const arma::vec& cTF,const arma::vec& sTF, - unsigned int gwg_reps,unsigned int chain_length=10000){ - - //arma::vec theta,arma::vec as,arma::vec bs, - unsigned int N = Y.n_rows; - unsigned int J = Y.n_cols; +Rcpp::List Gibbs_4PNO(const arma::mat &Y, const arma::vec &mu_xi, + const arma::mat &Sigma_xi_inv, const double &mu_theta, + const double &Sigma_theta_inv, double alpha_c, + double beta_c, double alpha_s, double beta_s, + unsigned int burnin, const arma::vec &cTF, + const arma::vec &sTF, unsigned int gwg_reps, + unsigned int chain_length = 10000) +{ + + // arma::vec theta,arma::vec as,arma::vec bs, + unsigned int N = Y.n_rows; + unsigned int J = Y.n_cols; // unsigned int K = max(area);//assuming area is coded from 1 to # of areas - - //save model deviances over iteractions - arma::vec deviance(chain_length-burnin); - - //sum of Y's columns & vc matrix - arma::vec oneN = arma::ones(N); - arma::vec Ysum = Y.t()*oneN; - - double tmburn; - arma::mat Ysimt(N,J); - arma::umat HistS(J+1,chain_length-burnin); - - //compute average kappas, thetas, as, and bs - arma::vec EAPtheta = arma::zeros(N); - - //Setting 2 categories per item - arma::vec Ms(J); - Ms.fill(2.0); - - //Savinging as,bs,gs,kappas, means and vcs of theta - arma::mat AS(J,chain_length); - arma::mat BS(J,chain_length); - arma::mat GSS(J,chain_length); - arma::mat SSS(J,chain_length); - arma::vec ms_thetas(chain_length); - arma::vec SD_thetas(chain_length); - - -//need to initialize, theta, as, bs, kappas, Z; eventually gs - arma::vec theta = arma::randn(N); - arma::vec oneJ = arma::ones(J); - arma::vec as = arma::ones(J) + 0.5*arma::randu(J); - arma::vec bs = arma::randn(J); - arma::vec uc = arma::randu(J); - arma::vec gs = (oneJ-arma::sqrt(uc) ) % cTF; - arma::vec ss = (arma::randu(J)%(oneJ-gs) ) % sTF; - arma::mat KAPS = kappa_initialize(Ms); - arma::mat Z = arma::zeros(N,J); - - - //Start chain - for(unsigned int t = 0; t < chain_length; t++){ - //Generate Z matrix. Note that Z will be overwritten to disk here. - Rcpp::List step1Z = update_WKappaZ_NA(Y,Ysum,Z,as,bs,gs,ss,theta,KAPS,alpha_c,beta_c,alpha_s,beta_s,gwg_reps); - - //update value for ss and gs - ss = Rcpp::as(step1Z[0]) % sTF; - gs = Rcpp::as(step1Z[1]) % cTF; - - //Update a, b; as and bs should be automatically updated here by writing to disk - Rcpp::List step2ab = update_ab_NA(N,J,Z,as,bs,theta,mu_xi,Sigma_xi_inv); - // Rcpp::List ste2ab = update_ab_norestriction(N,J,Z,as,bs,theta,mu_xi,Sigma_xi_inv);//I(alpha>0) removed - - //Update theta. Theta is stored in memory and function should update each iteration too. - Rcpp::List step3theta = update_theta(N,Z,as,bs,theta,mu_theta,Sigma_theta_inv); - //saving means and vcs of thetas + + // save model deviances over iteractions + arma::vec deviance(chain_length - burnin); + + // sum of Y's columns & vc matrix + arma::vec oneN = arma::ones(N); + arma::vec Ysum = Y.t() * oneN; + + double tmburn; + arma::mat Ysimt(N, J); + arma::umat HistS(J + 1, chain_length - burnin); + + // compute average kappas, thetas, as, and bs + arma::vec EAPtheta = arma::zeros(N); + + // Setting 2 categories per item + arma::vec Ms(J); + Ms.fill(2.0); + + // Savinging as,bs,gs,kappas, means and vcs of theta + arma::mat AS(J, chain_length); + arma::mat BS(J, chain_length); + arma::mat GSS(J, chain_length); + arma::mat SSS(J, chain_length); + arma::vec ms_thetas(chain_length); + arma::vec SD_thetas(chain_length); + + // need to initialize, theta, as, bs, kappas, Z; eventually gs + arma::vec theta = arma::randn(N); + arma::vec oneJ = arma::ones(J); + arma::vec as = arma::ones(J) + 0.5 * arma::randu(J); + arma::vec bs = arma::randn(J); + arma::vec uc = arma::randu(J); + arma::vec gs = (oneJ - arma::sqrt(uc)) % cTF; + arma::vec ss = (arma::randu(J) % (oneJ - gs)) % sTF; + arma::mat KAPS = kappa_initialize(Ms); + arma::mat Z = arma::zeros(N, J); + + // Start chain + for (unsigned int t = 0; t < chain_length; t++) { + // Generate Z matrix. Note that Z will be overwritten to disk here. + Rcpp::List step1Z = + update_WKappaZ_NA(Y, Ysum, Z, as, bs, gs, ss, theta, KAPS, alpha_c, + beta_c, alpha_s, beta_s, gwg_reps); + + // update value for ss and gs + ss = Rcpp::as(step1Z[0]) % sTF; + gs = Rcpp::as(step1Z[1]) % cTF; + + // Update a, b; as and bs should be automatically updated here by + // writing to disk + Rcpp::List step2ab = + update_ab_NA(N, J, Z, as, bs, theta, mu_xi, Sigma_xi_inv); + // Rcpp::List ste2ab = + // update_ab_norestriction(N,J,Z,as,bs,theta,mu_xi,Sigma_xi_inv);//I(alpha>0) + // removed + + // Update theta. Theta is stored in memory and function should update + // each iteration too. + Rcpp::List step3theta = + update_theta(N, Z, as, bs, theta, mu_theta, Sigma_theta_inv); + // saving means and vcs of thetas SD_thetas(t) = stddev(theta); ms_thetas(t) = mean(theta); - - //Storing output for as, bs, and kappas - AS.col(t) = as; - BS.col(t) = bs; - GSS.col(t) = gs; - SSS.col(t) = ss; - - // /* - if(t>burnin-1){ - tmburn = t-burnin; - //Compute EAPs here - EAPtheta = (tmburn*EAPtheta + theta)/(tmburn + 1.0); - - //Compute Simulated VC matrix - Ysimt = Y_4pno_simulate(N,J,as,bs,gs,ss,theta); - HistS.col(tmburn)=Total_Tabulate(N,J,Ysimt); - //Compute -2LL for DIC or Bayes Factors - deviance(tmburn) = min2LL_4pno(N,J,Y,as,bs,gs,ss,theta); + + // Storing output for as, bs, and kappas + AS.col(t) = as; + BS.col(t) = bs; + GSS.col(t) = gs; + SSS.col(t) = ss; + + // /* + if (t > burnin - 1) { + tmburn = t - burnin; + // Compute EAPs here + EAPtheta = (tmburn * EAPtheta + theta) / (tmburn + 1.0); + + // Compute Simulated VC matrix + Ysimt = Y_4pno_simulate(N, J, as, bs, gs, ss, theta); + HistS.col(tmburn) = Total_Tabulate(N, J, Ysimt); + // Compute -2LL for DIC or Bayes Factors + deviance(tmburn) = min2LL_4pno(N, J, Y, as, bs, gs, ss, theta); } - - - } - - return Rcpp::List::create(Rcpp::Named("EAPtheta",EAPtheta), - Rcpp::Named("HistS",HistS), - Rcpp::Named("AS",AS), - Rcpp::Named("BS",BS), - Rcpp::Named("GS",GSS), - Rcpp::Named("SS",SSS), - Rcpp::Named("ms_thetas",ms_thetas), - Rcpp::Named("SD_thetas",SD_thetas), - Rcpp::Named("Ds",deviance) - ); + } + + return Rcpp::List::create( + Rcpp::Named("EAPtheta", EAPtheta), Rcpp::Named("HistS", HistS), + Rcpp::Named("AS", AS), Rcpp::Named("BS", BS), Rcpp::Named("GS", GSS), + Rcpp::Named("SS", SSS), Rcpp::Named("ms_thetas", ms_thetas), + Rcpp::Named("SD_thetas", SD_thetas), Rcpp::Named("Ds", deviance)); } - + //' Update 2PNO Model Parameters -//' +//' //' Internal function to update 2PNO parameters -//' @param N The number of observations. -//' @param J The number of items. -//' @param Y A N by J \code{matrix} of item responses. -//' @param Z A \code{matrix} N by J of continuous augmented data. -//' @param as A \code{vector} of item discrimination parameters. -//' @param bs A \code{vector} of item threshold parameters. -//' @param theta A \code{vector} of prior thetas. -//' @param Kaps A \code{matrix} for item thresholds (used for internal computations). -//' @param mu_xi Prior mean for item parameters. -//' @param Sigma_xi_inv Prior item parameter inverse variance-covariance matrix. -//' @param mu_theta Prior mean for theta. +//' +//' @param N The number of observations. +//' @param J The number of items. +//' @param Y A N by J `matrix` of item responses. +//' @param Z A `matrix` N by J of continuous augmented data. +//' @param as A `vector` of item discrimination parameters. +//' @param bs A `vector` of item threshold parameters. +//' @param theta A `vector` of prior thetas. +//' @param Kaps A `matrix` for item thresholds +//' (used for internal computations). +//' @param mu_xi Prior mean for item parameters. +//' @param Sigma_xi_inv Prior item parameter inverse variance-covariance +//' matrix. +//' @param mu_theta Prior mean for theta. //' @param Sigma_theta_inv Prior inverse variance for theta. -//' @return A list of item parameters. -//' @seealso \code{\link{Gibbs_4PNO}} -//' @author Steven Andrew Culpepper -//' @export -// [[Rcpp::export]] -Rcpp::List update_2pno(unsigned int N,unsigned int J,const arma::mat& Y,arma::mat& Z,const arma::vec& as, - const arma::vec& bs,const arma::vec& theta, const arma::mat& Kaps,const arma::vec& mu_xi, - const arma::mat& Sigma_xi_inv,const double& mu_theta,const double& Sigma_theta_inv){ - - arma::vec eta(N); - arma::vec p4pno(2); - arma::vec oneN = arma::ones(N); - double uZ; - - //parms for alpha and beta - arma::vec m1(N); - m1.fill(-1.0); - arma::mat X_theta = join_rows(theta,m1); - arma::mat XX_xi = X_theta.t()*X_theta; - arma::mat Sig_xi_star = inv(Sigma_xi_inv + XX_xi); - arma::vec mu_xi_star(2); - double pa,mb_a,vb_a,ua,ub; - arma::vec as1(J),bs1(J); - - - for(unsigned int j=0;j(N); + double uZ; + + // parms for alpha and beta + arma::vec m1(N); + m1.fill(-1.0); + arma::mat X_theta = join_rows(theta, m1); + arma::mat XX_xi = X_theta.t() * X_theta; + arma::mat Sig_xi_star = inv(Sigma_xi_inv + XX_xi); + arma::vec mu_xi_star(2); + double pa, mb_a, vb_a, ua, ub; + arma::vec as1(J), bs1(J); + + for (unsigned int j = 0; j < J; j++) { + // Calculate eta_ij for all of the models and individuals with data + eta = as(j) * theta - oneN * bs(j); + // Update Z, W + for (unsigned int i = 0; i < N; i++) { + uZ = R::runif(0.0, 1.0); + p4pno(0) = R::pnorm(Kaps(Y(i, j), j), eta(i), 1.0, 1, 0); + p4pno(1) = R::pnorm(Kaps(Y(i, j) + 1.0, j), eta(i), 1.0, 1, 0); + Z(i, j) = R::qnorm(p4pno(0) + uZ * (p4pno(1) - p4pno(0)), eta(i), + 1.0, 1, 0); + } + + // update alpha,beta + arma::mat XZ_col = X_theta.t() * Z.col(j); + mu_xi_star = Sig_xi_star * (Sigma_xi_inv * mu_xi + XZ_col); + ua = R::runif(0.0, 1.0); + pa = R::pnorm(0.0, mu_xi_star(0), sqrt(Sig_xi_star(0, 0)), 1, 0); + as1(j) = R::qnorm(pa + ua * (1.0 - pa), mu_xi_star(0), + sqrt(Sig_xi_star(0, 0)), 1, 0); + ub = R::runif(0.0, 1.0); + mb_a = mu_xi_star(1) + + Sig_xi_star(0, 1) / Sig_xi_star(1, 1) * (as(j) - mu_xi_star(0)); + vb_a = Sig_xi_star(1, 1) - + Sig_xi_star(0, 1) * Sig_xi_star(0, 1) / Sig_xi_star(0, 0); + bs1(j) = R::qnorm(ub, mb_a, sqrt(vb_a), 1, 0); + } + + // update theta + double apb = dot(as1, bs1); + double vartheta = 1.0 / (dot(as1, as1) + Sigma_theta_inv); + arma::vec zetas(N, arma::fill::randn); + arma::vec theta1 = + zetas * sqrt(vartheta) + + vartheta * (Z * as1 + oneN * (apb + mu_theta * Sigma_theta_inv)); + + return Rcpp::List::create(Rcpp::Named("as1", as1), Rcpp::Named("bs1", bs1), + Rcpp::Named("theta1", theta1)); } //' Gibbs Implementation of 2PNO -//' +//' //' Implement Gibbs 2PNO Sampler -//' @param Y A N by J \code{matrix} of item responses. -//' @param mu_xi A two dimensional \code{vector} of prior item parameter means. -//' @param Sigma_xi_inv A two dimensional identity \code{matrix} of prior item parameter VC matrix. -//' @param mu_theta The prior mean for theta. -//' @param Sigma_theta_inv The prior inverse variance for theta. -//' @param burnin The number of MCMC samples to discard. -//' @param chain_length The number of MCMC samples. -//' @return Samples from posterior. -//' @author Steven Andrew Culpepper +//' +//' @param Y A N by J `matrix` of item responses. +//' @param mu_xi A two dimensional `vector` of prior item parameter +//' means. +//' @param Sigma_xi_inv A two dimensional identity `matrix` of prior item +//' parameter VC matrix. +//' @param mu_theta The prior mean for theta. +//' @param Sigma_theta_inv The prior inverse variance for theta. +//' @param burnin The number of MCMC samples to discard. +//' @param chain_length The number of MCMC samples. +//' +//' @return +//' Samples from posterior. +//' +//' @author +//' Steven Andrew Culpepper +//' //' @export //' @examples -//' require(fourPNO) //' # simulate small 2PNO dataset to demonstrate function //' J = 5 //' N = 100 @@ -675,14 +822,14 @@ Rcpp::List update_2pno(unsigned int N,unsigned int J,const arma::mat& Y,arma::ma //' # Population item parameters //' as_t = rnorm(J,mean=2,sd=.5) //' bs_t = rnorm(J,mean=0,sd=.5) -//' +//' //' # Sampling gs and ss with truncation //' gs_t = rbeta(J,1,8) //' ps_g = pbeta(1-gs_t,1,8) //' ss_t = qbeta(runif(J)*ps_g,1,8) //' theta_t = rnorm(N) //' Y_t = Y_4pno_simulate(N,J,as=as_t,bs=bs_t,gs=gs_t,ss=ss_t,theta=theta_t) -//' +//' //' # Setting prior parameters //' mu_theta = 0 //' Sigma_theta_inv = 1 @@ -690,99 +837,101 @@ Rcpp::List update_2pno(unsigned int N,unsigned int J,const arma::mat& Y,arma::ma //' alpha_c = alpha_s = beta_c = beta_s = 1 //' Sigma_xi_inv = solve(2*matrix(c(1,0,0,1), 2, 2)) //' burnin = 1000 -//' +//' //' # Execute Gibbs sampler. This should take about 15.5 minutes -//' out_t <- Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta,Sigma_theta_inv,alpha_c,beta_c,alpha_s, -//' beta_s,burnin,rep(1,J),rep(1,J),gwg_reps=5,chain_length=burnin*2) -//' +//' out_t = Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta,Sigma_theta_inv, +//' alpha_c,beta_c,alpha_s, beta_s,burnin, +//' rep(1,J),rep(1,J),gwg_reps=5,chain_length=burnin*2) +//' //' # Summarizing posterior distribution -//' OUT = cbind(apply(out_t$AS[,-c(1:burnin)],1,mean),apply(out_t$BS[,-c(1:burnin)],1,mean), -//' apply(out_t$GS[,-c(1:burnin)],1,mean),apply(out_t$SS[,-c(1:burnin)],1,mean), -//' apply(out_t$AS[,-c(1:burnin)],1,sd),apply(out_t$BS[,-c(1:burnin)],1,sd), -//' apply(out_t$GS[,-c(1:burnin)],1,sd),apply(out_t$SS[,-c(1:burnin)],1,sd) ) -//' OUT = cbind(1:J, OUT) -//' colnames(OUT) = c('Item','as','bs','gs','ss','as_sd','bs_sd','gs_sd','ss_sd') -//' print(OUT, digits=3) +//' OUT = +//cbind(apply(out_t$AS[,-c(1:burnin)],1,mean),apply(out_t$BS[,-c(1:burnin)],1,mean), +//' apply(out_t$GS[,-c(1:burnin)],1,mean),apply(out_t$SS[,-c(1:burnin)],1,mean), +//' apply(out_t$AS[,-c(1:burnin)],1,sd),apply(out_t$BS[,-c(1:burnin)],1,sd), ' +//apply(out_t$GS[,-c(1:burnin)],1,sd),apply(out_t$SS[,-c(1:burnin)],1,sd) ) ' +//OUT = cbind(1:J, OUT) ' colnames(OUT) = +//c('Item','as','bs','gs','ss','as_sd','bs_sd','gs_sd','ss_sd') ' print(OUT, +//digits=3) // [[Rcpp::export]] -Rcpp::List Gibbs_2PNO(const arma::mat& Y,const arma::vec& mu_xi,const arma::mat& Sigma_xi_inv, - const double& mu_theta,const double& Sigma_theta_inv,unsigned int burnin,unsigned int chain_length=10000){ - - //arma::vec theta,arma::vec as,arma::vec bs, - unsigned int N = Y.n_rows; - unsigned int J = Y.n_cols; +Rcpp::List Gibbs_2PNO(const arma::mat &Y, const arma::vec &mu_xi, + const arma::mat &Sigma_xi_inv, const double &mu_theta, + const double &Sigma_theta_inv, unsigned int burnin, + unsigned int chain_length = 10000) +{ + + // arma::vec theta,arma::vec as,arma::vec bs, + unsigned int N = Y.n_rows; + unsigned int J = Y.n_cols; // unsigned int K = max(area);//assuming area is coded from 1 to # of areas - - //save model deviances over iteractions - arma::vec deviance(chain_length-burnin); - - //sum of Y's columns & vc matrix - arma::vec oneN = arma::ones(N); - arma::vec Ysum = Y.t()*oneN; - - double tmburn; - arma::mat Ysimt(N,J); - arma::umat HistS(J+1,chain_length-burnin); - - //compute average kappas, thetas, as, and bs - arma::vec EAPtheta = arma::zeros(N); - - //Setting 2 categories per item - arma::vec Ms(J); - Ms.fill(2.0); - - //Savinging as,bs,gs,kappas, means and vcs of theta - arma::mat AS(J,chain_length); - arma::mat BS(J,chain_length); - arma::vec ms_thetas(chain_length); - arma::vec SD_thetas(chain_length); - -//need to initialize, theta, as, bs, kappas, Z; eventually gs - arma::vec theta = arma::randn(N);//arma::zeros(N); - arma::vec oneJ = arma::ones(J); - arma::vec as = 0.5*arma::randu(J) + arma::ones(J); - arma::vec bs = arma::randn(J); - arma::vec gs = arma::zeros(J); - arma::vec ss = arma::zeros(J); - arma::mat KAPS = kappa_initialize(Ms); - arma::mat Z = arma::zeros(N,J); - - //Start chain - for(unsigned int t = 0; t < chain_length; t++){ - //Generate Z matrix. Note that Z will be overwritten to disk here. - Rcpp::List step1Z = update_2pno(N,J,Y,Z,as,bs,theta,KAPS,mu_xi,Sigma_xi_inv,mu_theta,Sigma_theta_inv); - //update value for as, bs, theta - as = Rcpp::as(step1Z[0]) ; - bs = Rcpp::as(step1Z[1]) ; - theta = Rcpp::as(step1Z[2]) ; + + // save model deviances over iteractions + arma::vec deviance(chain_length - burnin); + + // sum of Y's columns & vc matrix + arma::vec oneN = arma::ones(N); + arma::vec Ysum = Y.t() * oneN; + + double tmburn; + arma::mat Ysimt(N, J); + arma::umat HistS(J + 1, chain_length - burnin); + + // compute average kappas, thetas, as, and bs + arma::vec EAPtheta = arma::zeros(N); + + // Setting 2 categories per item + arma::vec Ms(J); + Ms.fill(2.0); + + // Savinging as,bs,gs,kappas, means and vcs of theta + arma::mat AS(J, chain_length); + arma::mat BS(J, chain_length); + arma::vec ms_thetas(chain_length); + arma::vec SD_thetas(chain_length); + + // need to initialize, theta, as, bs, kappas, Z; eventually gs + arma::vec theta = arma::randn(N); // arma::zeros(N); + arma::vec oneJ = arma::ones(J); + arma::vec as = 0.5 * arma::randu(J) + arma::ones(J); + arma::vec bs = arma::randn(J); + arma::vec gs = arma::zeros(J); + arma::vec ss = arma::zeros(J); + arma::mat KAPS = kappa_initialize(Ms); + arma::mat Z = arma::zeros(N, J); + + // Start chain + for (unsigned int t = 0; t < chain_length; t++) { + // Generate Z matrix. Note that Z will be overwritten to disk here. + Rcpp::List step1Z = + update_2pno(N, J, Y, Z, as, bs, theta, KAPS, mu_xi, Sigma_xi_inv, + mu_theta, Sigma_theta_inv); + // update value for as, bs, theta + as = Rcpp::as(step1Z[0]); + bs = Rcpp::as(step1Z[1]); + theta = Rcpp::as(step1Z[2]); SD_thetas(t) = stddev(theta); ms_thetas(t) = mean(theta); - - //Storing output for as, bs, and kappas - AS.col(t) = as; - BS.col(t) = bs; - - if(t>burnin-1){ - tmburn = t-burnin; - //Compute EAPs here - EAPtheta = (tmburn*EAPtheta + theta)/(tmburn + 1.0); - - //Compute Simulated VC matrix - Ysimt = Y_4pno_simulate(N,J,as,bs,gs,ss,theta); - HistS.col(tmburn)=Total_Tabulate(N,J,Ysimt); - //Compute -2LL for DIC or Bayes Factors - deviance(tmburn) = min2LL_4pno(N,J,Y,as,bs,gs,ss,theta); + + // Storing output for as, bs, and kappas + AS.col(t) = as; + BS.col(t) = bs; + + if (t > burnin - 1) { + tmburn = t - burnin; + // Compute EAPs here + EAPtheta = (tmburn * EAPtheta + theta) / (tmburn + 1.0); + + // Compute Simulated VC matrix + Ysimt = Y_4pno_simulate(N, J, as, bs, gs, ss, theta); + HistS.col(tmburn) = Total_Tabulate(N, J, Ysimt); + // Compute -2LL for DIC or Bayes Factors + deviance(tmburn) = min2LL_4pno(N, J, Y, as, bs, gs, ss, theta); } - - - } - - return Rcpp::List::create(Rcpp::Named("EAPtheta",EAPtheta), - Rcpp::Named("HistS",HistS), - Rcpp::Named("AS",AS), - Rcpp::Named("BS",BS), - Rcpp::Named("ms_thetas",ms_thetas), - Rcpp::Named("SD_thetas",SD_thetas), - Rcpp::Named("Ds",deviance) - ); + } + + return Rcpp::List::create( + Rcpp::Named("EAPtheta", EAPtheta), Rcpp::Named("HistS", HistS), + Rcpp::Named("AS", AS), Rcpp::Named("BS", BS), + Rcpp::Named("ms_thetas", ms_thetas), + Rcpp::Named("SD_thetas", SD_thetas), Rcpp::Named("Ds", deviance)); }