### Discrete Graded Response Models. ### These are similar to the DiBello-Samejima models, except that the ### two Samejima curves have different discriminations. Also, the ### discriminations go with the output states, not the input ### variables. ### skillLevels --- list of labels for the skill level variable. ### obsLevels --- labels for obserable values. ### skillLevels and obsLevels are assumed to be ordered highest to ### lowest. ### Unlike the DiBello-Samejima rule, lnAlphas, betas, and ### lnAlphas --- set of log slopes, one for each obs level ### beta --- difficulty (-intercept) parameter ### rule --- Function for computing effective theta. This should have the ### signature function(thetas,alphas,beta) ### The first one in the list representes the difference between the ### highest and next lowest obs state, and so forth. calcDPCTable <- function (skillLevels, obsLevels, lnAlphas, betas, rules="Compensatory", link="partialCredit", linkScale=NULL,Q=TRUE, tvals=lapply(skillLevels, function (sl) effectiveThetas(length(sl)))) { k <- length(obsLevels) if (k < 2) { stop("There must be at least two obsLevels to caluclate a CPT, got", obsLevels) } if (!is.list(lnAlphas)) lnAlphas <- list(lnAlphas) if (length(lnAlphas) != k-1) lnAlphas <- rep(lnAlphas,k-1) if (!is.list(betas)) betas <- list(betas) if (length(betas) != k-1) betas <- rep(betas,k-1) if (!is.list(rules)) rules <- list(rules) if (length(rules) != k-1) rules <- rep(rules,k-1) if (!is.list(tvals)) tvals <- list(tvals) p <- length(skillLevels) if (length(Q)==1) { Q <- matrix(TRUE,k-1,min(p,1)) } else if (!is.matrix(Q) || nrow(Q) != k-1 || ncol(Q) != p) { stop("Q must be a",k-1,"by",p,"matrix.") } thetas <- do.call("expand.grid",tvals) if (nrow(thetas) == 0L) { thetas <- data.frame(X0=0) } et <- matrix(0,nrow(thetas),k-1) #Take care of no parent case for (kk in 1:(k-1)) { et[,kk] <- do.call(rules[[kk]], list(thetas[,Q[kk,],drop=FALSE], exp(lnAlphas[[kk]]),betas[[kk]])) } do.call(link,list(et,linkScale,obsLevels)) } calcDPCFrame <- function (skillLevels, obsLevels, lnAlphas, betas, rules="Compensatory", link="partialCredit", linkScale=NULL,Q=TRUE, tvals=lapply(skillLevels, function (sl) effectiveThetas(length(sl)))) { markers <- expand.grid(skillLevels) tab <- calcDPCTable(skillLevels,paste(obsLevels), lnAlphas,betas,rules,link, linkScale,Q,tvals) if (length(skillLevels)>0L) { result <- data.frame(markers,tab) } else { result <- data.frame(tab) } if (is.numeric(obsLevels) || names(result)[length(skillLevels)+1]!=paste(obsLevels[1])) { ## R is "helpfully" fixing our numeric labels. Need to insist. names(result) <- c(names(skillLevels),paste(obsLevels)) } result } ### Link Functions ### These are a function of three arguments ### et -- a table of effective thetas, one for row for each ### configuration of parent variables, and one column for ### each child state except for the last. ### k -- the number of states of the child variable. ### linkScale --- a scale parameter used by some link functions. ### obsLevel --- a list of names of the observables, assume to be ### sorted from highest to lowest. ### It returns a conditional probility table. partialCredit <- function (et,linkScale=NULL,obsLevels=NULL) { m <- ncol(et)+1 zt <- apply(et,1,function(x) rev(cumsum(rev(x)))) if (m==2) { zt <- matrix(zt,length(et),1) } else { zt <- t(zt) } pt <- cbind(exp(1.7*zt),1) pt <- sweep(pt,1,apply(pt,1,sum),"/") probs <- pt[,1:m,drop=FALSE] probs <- ifelse(probs<0,0,probs) if (!is.null(obsLevels)) { dimnames(probs) <- list(NULL,obsLevels) } probs } gradedResponse <- function (et,linkScale=NULL,obsLevels=NULL) { m <- ncol(et)+1 zt <- 1/(1+exp(-1.7*et)) ## The cummax corrects for problems where the ## Pr(X> m+1) is greater than Pr(X>m) zt <- apply(zt,1,cummax) if (m==2) { zt <- matrix(zt,length(et),1) } else { zt <- t(zt) } pt <- cbind(0,zt,1) probs <- pt[,2:(m+1),drop=FALSE]-pt[,1:m,drop=FALSE] # probs <- ifelse(probs<0,0,probs) if (!is.null(obsLevels)) { dimnames(probs) <- list(NULL,obsLevels) } probs } normalLink <- function (et,linkScale=NULL,obsLevels=NULL) { if (is.null(linkScale)) stop("Link Scale parameter required for normal link.") m <- ncol(et)+1 cuts <- qnorm( ((m-1):1)/m) ## Only play attention to the first column. pt <- pnorm(outer(-et[,1],cuts,"+")/linkScale) pt <- cbind(1,pt,0) probs <- pt[,1:m,drop=FALSE]-pt[,1+(1:m),drop=FALSE] if (!is.null(obsLevels)) { dimnames(probs) <- list(NULL,obsLevels) } probs }