Nothing
#' Draw values from multinomial distribution
#'
#' Draws values from a mutinomial distribution and matches them with a row in the pmatrix (attribute profile matrix)
#' to obtain the alpha values
#' @param pmatrix a matrix indicating all possible attribute profiles
#' @param y number of random vectors to draw
#' @param nus0 a numeric vector of length \code{nclasses} for nu parameter estimates
#' @return a numeric vector of alpha values
DrawAlphas <- function(pmatrix, y, nus0){
pmatrix[which(rmultinom(y, 1, nus0) ==1),]
}
#' Draw latent variable values
#'
#' Produces values of alphas(s) (i.e., latent variable(s)) by using a rejection/acceptance Metropolis-Hastings procedure.
#' @param nattributes a numeric value indicating number of alphas to be drawn
#' @param class0 a numeric vector of previous classes
#' @param nus0 a numeric vector of length \code{nclasses} for nu parameter estimates
#' @param observations a data frame or matrix of dichotomous responses
#' @param nobservations a numeric value indicating number of rows of the observation data frame or matrix
#' @param threshold.values an nclasses by nitems numeric matrix with appropriate item threshold values
#' @param pmatrix a numeric nclasses by nattributes matrix of all possible attribute profiles
#' @return a numeric vector of classes indicating one of the possible attribute profile for each student
DrawClasses <- function(nattributes, class0, nus0, observations, nobservations, threshold.values, pmatrix) {
alphas <- t(sapply(rep(1, nobservations), function(y) DrawAlphas(pmatrix, y, nus0)))
class1 <- apply(alphas, 1, function(y) which(apply(pmatrix, 1, function(x) all(x == y))))
probs0 <- InvLogit(threshold.values[class0,])
probs0 <- ifelse(observations == 1, probs0, 1 - probs0)
dcm0 <- apply(probs0, 1, function(x) prod(x, na.rm=T))
probs1 <- InvLogit(threshold.values[class1,])
probs1 <- ifelse(observations == 1, probs1, 1 - probs1)
dcm1 <- apply(probs1, 1, function(x) prod(x, na.rm=T))
accept <- ifelse(dcm0==0,1,dcm1/dcm0)
accept <- ifelse(accept>1,1,accept)
accept <- ifelse(runif(nobservations)<accept,1,0)
class1 <- ifelse(accept,class1,class0)
return(class1)
}
#' Perform one iteration of MCMC procedure
#'
#' If applicable, randomly samples new set of parameter estimates, obtains applicable estimates
#' and uses those to calculate threshold values for both items and latent variables, draws
#' new set of alpha values.
#' @param nattributes numberic value for number of attributes
#' @param class0 The previous value of attribute profile for each respondent
#' @param estimates0 a numeric vector of parameter estimates
#' @param threshold.labels an nclasses by nitems character matrix with appropriate threshold labels
#' @param lambda.equations equations for lambda parameters
#' @param is.pi.r If \code{FALSE} (the default), parameter values are the type of \code{taus} and \code{nus} or
#' \code{lambdas} and \code{gammas} else they are the type \code{pis} and \code{rs} as used in NC-RUM parameterization
#' @param parameter.means a numerical vector of calibrated item and structural parameters
#' @param parameter.acov a numerical matrix of covariances of item and structural parameters
#' @param observations a data frame or matrix of dichotomous responses
#' @param nobservations a numeric value of number of observations
#' @param is.parameter.randomized if true parameter estimates are randomized using acov matrix
#' @param qmatrix a data frame or matrix of 1s and 0s indicating relation between items and attributes.
#' This matrix specifies which items are required for mastery of each attribute (i.e., latent variable).
#' A matrix must be a size of \code{nItems X nAttributes}
#' @param pmatrix a numeric nclasses by nattributes matrix of all possible attribute profiles
#' @return a list of newly sampled classes and parameter estimates
iterate <-function(nattributes, class0, estimates0, threshold.labels, lambda.equations
, is.pi.r, parameter.means, parameter.acov, observations, nobservations
, is.parameter.randomized, qmatrix, pmatrix) {
if (is.parameter.randomized)
estimates0 <- SampleParameterEstimates(parameter.means, parameter.acov)
if (!is.null(lambda.equations)){
threshold.info <- GetThresholdValuesKernel(lambda.equations, estimates0, threshold.labels)
taus <- -threshold.info$taus
threshold.values <- -threshold.info$threshold.values
mus <- GetMusFromGammas(estimates0, pmatrix)
estimates0 <- c(taus, mus)
}
if (is.pi.r){
taus <- GetThresholdValuesKernelPiR(estimates0, qmatrix, pmatrix, threshold.labels)
mus <- GetMusFromGammas(estimates0, pmatrix)
estimates0 <- c(taus, mus)
}
parms <- InitializeParameters(estimates0, nrow(pmatrix))
taus0 <- parms[['taus']]
nus0 <- parms[['nus']]
if (is.null(lambda.equations)){
threshold.values <- GetThresholdValues(threshold.labels, taus0)
}
class1 <- DrawClasses(nattributes, class0, nus0, observations, nobservations, threshold.values, pmatrix)
return(list(new.class = class1, estimates = estimates0))
}
#' Randomly sample parameter estimates
#'
#' Given parameter estimates and asymptotic covaraince matrix of parameter estimates, randomly samples a new
#' set of estimates
#' @param means a numerical vector of calibrated item and structural parameters
#' @param acov a numerical matrix of covariances of item and structural parameters
#' @return A vector of randomly sampled parameter estimate values
SampleParameterEstimates <- function(means, acov){
#TODO: check for non-singularity of acov matrix
if (det(acov) != 0)
L= chol(acov)
else
L = t(ldl(acov)) #chol decomposition of acov (L*L` = acov solve for L, square root of acov matrix)
Z = rnorm(matrix(0, 1, length(means))) #random standard normal variates (theta ~ N(means, acov)
# where theta is vector of all parameters- taus and mus)
Z = matrix(Z, 1, length(means)) #placing in matrix form
simulated.values = c(Z%*%L) + means#getting simulated values (theta = N(0,1)*sqrt(acov) + means)
return(simulated.values)
}
#' Initialize parameter estimates
#'
#' Obtains appropriate parameter estimates for a non-kernel model (e.g., taus and nus)
#' @param estimates a vector of parameter estimates
#' @param nclasses a numeric value representing the number of unique attribute profiles
InitializeParameters <- function(estimates, nclasses){
taus <- unlist(estimates[grep('^tau', names(estimates))])
mus <- estimates[grep('^mu', names(estimates))]
if (length(mus) < nclasses)
mus <- c(mus, rep(0, nclasses -length(mus)))
nus <- GetClassProbsFromMus(mus)
list("taus" = taus, "nus" = nus)
}
#' Performs MCMC routine for DCM
#'
#' @param observations a data frame or matrix of dichotomous responses
#' @param nattributes numeric value of number of attributes
#' @param qmatrix a data frame or matrix of 1s and 0s indicating relation between items and attributes.
#' This matrix specifies which items are required for mastery of each attribute (i.e., latent variable).
#' A matrix must be a size of \code{nItems X nAttributes}
#' @param pmatrix a numeric nclasses by nattributes matrix of all possible attribute profiles
#' @param parameter.means a numerical vector of calibrated item and structural parameters
#' @param parameter.acov a numerical matrix of covariances of item and structural parameters
#' @param nobservations a numeric value indicating number of rows of the observation data frame or matrix
#' @param nreps The number of iterations in MCMC per chain
#' @param initial.class The initial value of attribute profile for each respondent
#' @param nchains The number of chains in MCMC
#' @param threshold.labels an nclasses by nitems character matrix with appropriate item threshold labels
#' @param lambda.equations lambda parameter equations
#' @param is.pi.r If \code{FALSE} (the default), parameter values are the type of \code{taus} and \code{nus} or
#' \code{lambdas} and \code{gammas} else they are the type \code{pis} and \code{rs} as used in NC-RUM parameterization
#' @param is.parameter.randomized if true parameter estimates are randomized using acov matrix
#' @param parameterization.method optional character string of parameterization method used to calibrate parameters
#' @param percent.reps.to.discard The percent of iterations to be discarded
#' @return a list of class and parameter data frame containing all accepted iteraction of MCMC
mcmc <- function(observations, nattributes, qmatrix, pmatrix, parameter.means
, parameter.acov, nobservations, nreps, initial.class, nchains
, threshold.labels, lambda.equations, is.pi.r
, is.parameter.randomized, parameterization.method, percent.reps.to.discard){
chain.class.results <- list()
chain.parameter.results <- list()
cols.to.discard = ceiling(nreps*percent.reps.to.discard*0.01)
# begin chains
for (chain in 1:nchains){
# sample from the parameter estimates (mean vector) and acov matrix (cov matrix)
estimates <- parameter.means
class.out <- matrix(NA, nobservations, nreps)
class.out[, 1] <- initial.class
class0 <- initial.class
estimates.out <- matrix(NA, length(estimates), nreps)
rownames(estimates.out) <- names(estimates)
estimates.out[, 1] <- estimates
# begin replications
for(i in 2:nreps){
# perform one iteration
out <- iterate(nattributes, class0, estimates, threshold.labels, lambda.equations, is.pi.r
, parameter.means, parameter.acov, observations, nobservations
,is.parameter.randomized, qmatrix, pmatrix)
class.out[, i] <- class0 <- out[["new.class"]]
estimates.out[, i] <- out[["estimates"]]
}
chain.class.results[[chain]] <- class.out[, -(1:cols.to.discard)]
chain.parameter.results[[chain]] <- estimates.out[, -(1:cols.to.discard)]
}
names(chain.class.results) <- names(chain.parameter.results) <- 1:nchains
return(list(class.result = chain.class.results, parameter.result = chain.parameter.results))
}
#' Calculate attribute probabilities
#'
#' Calculates the attribute (i.e., latent variable) probabilities for all students
#' @param class.result a list of classes containing all interactions and chains of MCMC
#' @param pmatrix a numeric nclasses by nattributes matrix of all possible attribute profiles
GetAttributesProbabilitiesMCMC <- function(class.result, pmatrix){
#theta_values <- ddply(narrow.results, .(rep, id, group), print(group))
narrow.results <- with(class.result, LongFormatResults(class.result, 'klass'))
alpha.values <- with(narrow.results, ddply(narrow.results, .(chain, rep, id, klass), function(x) pmatrix[x$klass,]))
narrow.alpha.values <- with(alpha.values, melt(alpha.values, id.vars = c("chain", "rep", "id", "klass"), variable.name = "attr.name", value.name = "attr.value"))
narrow.alpha.results <- with(narrow.alpha.values, ddply(narrow.alpha.values
, .(id, attr.name), summarize, means = mean(attr.value), sds = sd(attr.value)))
narrow.alpha.results.melted <- with(narrow.alpha.results, melt(narrow.alpha.results
, id.vars = c("id", "attr.name"), variable.name = "moment", value.name = "moment.value"))
alpha.results <- with(narrow.alpha.results.melted, dcast(narrow.alpha.results.melted, id ~ attr.name + moment, value.var = 'moment.value'))
return(alpha.results)
}
#' Calculate attribute profile probabilities
#'
#' Calculates the attribute profile probabilities for all students
#' @param class.result a list of classes containing all interactions and chains of MCMC
#' @param nclasses a numeric value representing the number of unique attribute profiles
GetClassProbabilitiesMCMC <- function(class.result, nclasses){
narrow.results <- with(class.result, LongFormatResults(class.result, 'klass'))
narrow.class.probs <- with(narrow.results, ddply(narrow.results, .(chain, id, klass), summarize, prob = length(klass)/max(rep)))
narrow.class.probs <- with(narrow.class.probs, ddply(narrow.class.probs, .(id, klass), summarize, prob = mean(prob)))
narrow.class.probs <- with(narrow.class.probs, rbind(narrow.class.probs, data.frame(id = rep(0, nclasses), klass = 1:nclasses, prob = rep(0, nclasses))))
max_ <- with(narrow.class.probs, ddply(narrow.class.probs, .(id), summarize, max.class = klass[which(max(prob, na.rm=T) == prob)]))
class.probs <- with(narrow.class.probs, dcast(narrow.class.probs, id ~ klass, value.var = 'prob', fill = 0))
class.probs.max <- merge(class.probs, max_, by = "id")
return(class.probs.max[!class.probs.max$id == 0,])
}
#' Calculate results for attributes and attribute profiles
#'
#' Calculates results for attributes and attribute profiles for all students based on MCMC output
#' @param parameter.result a list parameter estimates for each iterations and chains
GetParameterResultsMCMC <- function(parameter.result){
parameter.result <- with(parameter.result, LongFormatResults(parameter.result, 'estimate'))
parameter.result <- with(parameter.result, ddply(parameter.result, .(id), summarize, parm.means = mean(estimate), parm.sds = sd(estimate)))
return(parameter.result)
}
#' Tranform dataframe to long format
#'
#' Given a dataframe in wide format, transform to long format
#' @param d a list of MCMC output either classes or parameter estimates
#' @param result.type optional character strings containing result.type: \code{class} or \code{parameter}
LongFormatResults <- function(d, result.type = "class"){
long.d <- ldply(d, function(x) melt(x, value.name = result.type))
long.d <- rename(long.d, c(".id" = "chain", "Var1" = "id", "Var2" = "rep"))
return(long.d)
}
#' Score dichotomous response data using DCM and MCMC
#'
#' Given dichotomous response data, a Q-matrix (relation between items and attributes) and calibrated item and structural parameters
#' \code{ScoreDCM} estimates posterior probabilities of attribute profiles of respondents using a Diagnostic Classification Model (DCM)
#' and Markov Chain Monte Carlo (MCMC) method. The estimation procedure uses the loglinear cognitive diagnostic modeling (LDCM)
#' framework that allows for the estimation of a host of DCMs such as DCM, DINA, C-RUM, NIDO, NIDA, NC-RUM etc.
#' @param observations a data frame or matrix of dichotomous responses in wide format (\code{nobservations X nitems})
#' @param qmatrix a data frame or matrix of 1s and 0s indicating relation between items and attributes. This matrix specifies which items are required
#' for mastery of each attribute (i.e., latent variable). The matrix must be a size of \code{nitems X nattributes}
#' @param parameter.means a numerical vector of calibrated item and structural parameters. Values must be in the order of
#' \code{\link{GetParameterNames}} if parametrization method is Mplus and non-kernel parameters are used. If kernel parameters values
#' are used must be in order of \code{\link{GetKernelParameterNames}}
#' @param parameter.acov optional matrix of covariances of all model parameters. If \code{NULL} (the default) model parameters are not randomized
#' for each iteration of MCMC
#' @param parameterization.method optional character string of parameterization method used to calibrate parameters.
#' If not specified then the default will be set to \code{Mplus}
#' @param is.kernel.parameters If \code{FALSE} (the default), parameter values are of type \code{taus} and \code{nus} else they are of type
#' kernel parameters, i.e., \code{lambdas} and \code{gammas}
#' @param model.type If \code{is.kernel.parameter} is \code{TRUE}, model type must be one of \code{DCM}, \code{DINA}, \code{CRUM}, \code{DINO},
#' \code{NIDO}, \code{NCRUM}. Kernel parameters are different for each model type
#' @param nreps The number of iterations in MCMC per chain. If not specified the default is 1000
#' @param nchains The number of chains in MCMC. If not specified the default is 1
#' @param initial.class The initial value of the attribute profile for each respondent. If not specified the default is 1
#' @param percent.reps.to.discard The percent of iterations to be discarded. If not specifed then the default is 5
#'
#' @return ScoreDCM returns an object of class \code{\link{dcm.scorer.class}}; a list consisting of
#' \describe{
#' \item{inputs}{a list of all input arguments to the funciton like observations, Q-matrix, parameter.means, etc.}
#' \item{mcmc.inputs}{a list of all mcmc inputs}
#' \item{results}{an object of \code{\link{all.results.class}} class; a list consisting of
#' \describe{
#' \item{attribute.profile.result}{an object of \code{\link{attribute.profile.class}}}
#' \item{attribute.result}{an object of \code{\link{attribute.class}}}
#' \item{parameter.result}{an object of \code{\link{parameter.class}}}
#' }
#' }
#' \item{mcmc.outputs}{a list of all mcmc output}
#' }
#' The function \code{\link{summary}} is used to obtain and print a summary of
#' MCMC runs in the form of probabilities of mastering each attribute and attribute profile probabilities.
#' The function \code{\link{plot}} is used to plot the aggregated mean of both attribute mastery (type = "attr.means")
#' and attribute profile probabily (type = "attr.profile.means") across all respondents. Other plot options include
#' attribute mastery (type = "attr.profiles") and attribute profile probabilities (type = "attr.profile.profiles")
#' of individual respondents.
#' @author Margi Dubal \email{margidubal@@gmail.com} & Diane Losardo \email{dlosardo@@amplify.com}
#' @usage ScoreDCM(observations, qmatrix, parameter.means, parameter.acov= NULL,
#' parameterization.method = "Mplus", is.kernel.parameters = FALSE, model.type = NULL,
#' nreps = 1000 , nchains = 1, initial.class = 1, percent.reps.to.discard = 5)
#' @examples
#'
#' \dontrun{
#'
#' data(obervations.test)
#' data(qmatrix.test)
#' parameter.names <- GetParameterNames(qmatrix.test, ncol(qmatrix.test))
#' parameter.names
#' data(parameter.means.DCM.Mplus.test)
#' obj <- ScoreDCM(observations = observations.test, qmatrix = qmatrix.test
#' , parameter.means = parameter.means.DCM.Mplus.test)
#' summary(obj)
#' plot(obj)
#' }
#' @export
ScoreDCM <- function(observations, qmatrix, parameter.means, parameter.acov= NULL,
parameterization.method = "Mplus", is.kernel.parameters = FALSE, model.type = NULL,
nreps = 1000 , nchains = 1, initial.class = 1, percent.reps.to.discard = 5){
kAllModelTypes <- c('DCM', 'DINA', 'CRUM', 'DINO', 'NIDO', 'NCRUM')
is.pi.r <- FALSE
if (!hasArg(observations))
stop("argument \"observations\" is missing, with no default")
if (!hasArg(qmatrix))
stop("argument \"qmatrix\" is missing, with no default")
if (!hasArg(parameter.means))
stop("argument \"parameter.means\" is missing, with no default")
if (is.null(observations) | (!(is.data.frame(observations) | is.matrix(observations))))
stop("argument \"observations\" must be data frame or a matrix")
if (nrow(observations) == 0)
stop("argument \"observations\" can not be empty")
if (!all(sapply(observations, is.numeric)))
stop("argument \"observations\" must be a matrix of all numeric values")
if (all(sapply(observations, function(x) x!=1 & x!=0)))
stop("argument \"observations\" must be a matrix of 0s and 1s")
if (!is.matrix(qmatrix))
stop("argument \"qmatrix\" must be a matrix")
if (nrow(qmatrix) == 0 | ncol(qmatrix) == 0)
stop("argument \"qmatrix\" should not be empty")
if (!all(sapply(qmatrix, is.numeric)))
stop("argument \"qmatrix\" must be a matrix of all numeric values")
if (all(sapply(qmatrix, function(x) x!=1 & x!=0)))
stop("argument \"qmatrix\" must be a matrix of 0s and 1s")
if (nrow(qmatrix) != ncol(observations))
stop(sprintf("argument \"qmatrix\" must be a matrix %d rows", ncol(observations)))
if (!all(is.numeric(parameter.means)))
stop("parameter.means must be vector of numeric values")
if (!all(sapply(qmatrix, is.numeric)))
stop("parameter.acov must be matrix of numeric values")
if (is.kernel.parameters){
if (is.null(model.type))
stop(sprintf("model.type must be one of: %s", paste(kAllModelTypes, collapse = ', ')))
if(!(model.type %in% kAllModelTypes))
stop(sprintf("model.type must be one of: %s", paste(kAllModelTypes, collapse = ', ')))
}
nobservations <- nrow(observations)
nitems <- nrow(qmatrix)
nattributes <- ncol(qmatrix)
nclasses <- 2^nattributes
if (!is.null(rownames(qmatrix))){
item.names <- rownames(qmatrix)
} else{
item.names <- sprintf('item%d', 1:nitems)
}
if (!is.null(colnames(qmatrix))){
attribute.names <- colnames(qmatrix)
} else{
attribute.names <- sprintf('attribute%d', 1:nattributes)
}
dimnames(qmatrix) <- list(item.names, attribute.names)
pmatrix <- GetAttributeProfiles(nattributes, attribute.names)
#required parameters for a given model type
if (!is.kernel.parameters){
required.parms <- GetParameterNames(qmatrix, nattributes, parameterization.method)
names(parameter.means) <- required.parms
if (parameterization.method == "Mplus"){
parameter.means[grep('^tau', names(parameter.means))] <- -parameter.means[grep('^tau', names(parameter.means))]
}
lambda.equations <- NULL
} else {
kernel.parms <- GetKernelParameterNames(qmatrix, nattributes, model.type)
if (model.type == "NCRUM"){
names(parameter.means) <- kernel.parms$pis.rias
is.pi.r <- TRUE
} else {
names(parameter.means) <- c(kernel.parms$lambdas$lambda.name, kernel.parms$gammas)
lambda.equations <- kernel.parms$lambdas$lambda.equations
}
}
if(length(parameter.means) != length(required.parms))
stop(sprintf("parameter means should be a length of %d", length(required.parms)))
if(!is.null(parameter.acov)){
is.parameter.randomized = TRUE
if ((nrow(parameter.acov) != ncol(parameter.acov)) | (ncol(parameter.acov) != length(required.parms)))
stop(sprintf("parameter acov should be a matrix of [%d, %d]", length(required.parms), length(required.parms)))
} else{
is.parameter.randomized = FALSE
}
initial.class <- rep(initial.class, nobservations)
threshold.labels <- GetThresholdLabels(qmatrix, pmatrix)
mcmc.out <- mcmc(observations, nattributes, qmatrix, pmatrix, parameter.means
, parameter.acov, nobservations, nreps, initial.class, nchains
, threshold.labels, lambda.equations, is.pi.r, is.parameter.randomized
, parameterization.method, percent.reps.to.discard)
attribute.profile.out <- attribute.profile.class(results = GetClassProbabilitiesMCMC(mcmc.out$class.result, nclasses), attribute.profile.matrix = pmatrix)
attribute.out <- attribute.class(results = GetAttributesProbabilitiesMCMC(mcmc.out$class.result, pmatrix))
parameter.out <- parameter.class(results = GetParameterResultsMCMC(mcmc.out$parameter.result))
inputs <- list(observations = observations, qmatrix = qmatrix, parameter.means = parameter.means,
nobservations = nobservations, nitems = nitems, nattributes = nattributes, nclasses = nclasses,
attribute.profiles = pmatrix, parameterization.method = parameterization.method,
is.parameter.randomized = is.parameter.randomized, parameter.acov = parameter.acov)
mcmc.inputs <- list(nreps = nreps, nchains = nchains, burn = percent.reps.to.discard,
initial.class = initial.class)
outputs <- all.results.class(attribute.profile.result = attribute.profile.out , attribute.result = attribute.out,
parameter.result= parameter.out)
score.dcm <- dcm.scorer.class(inputs = inputs, mcmc.inputs = mcmc.inputs, results = outputs, mcmc.output = mcmc.out)
return(score.dcm)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.