Parameter Estimation

library(CPTtools)
skill1 <- c("High","Med","Low")
skill2 <- c("High","Med","Low")
isCorrect <- c("Correct","Incorrect")
troph <- c("Gold","Silver","None")
tol <- .001

Compensatory Model

True Model

pLevels <- list(skill1)
obsLevels <- troph
rules <- "Compensatory"
link <- "partialCredit"
testname <- paste("Testing:  rules=",paste(rules,collapse=", "),"; link=",link)
trueLnAlphas <- list(gold=log(1),silver=log(.25))
trueBetas <- list(gold=2,silver=-.5)
truedist <- calcDPCTable(pLevels,obsLevels,trueLnAlphas,trueBetas,
                           rules=rules,link=link)
round(truedist,3)
#>       Gold Silver  None
#> [1,] 0.119  0.687 0.195
#> [2,] 0.023  0.685 0.293
#> [3,] 0.004  0.606 0.390

Data Generation

weights <- 1000
counts <- round(sweep(truedist,1,weights,"*"))
round(counts)
#>      Gold Silver None
#> [1,]  119    687  195
#> [2,]   23    685  293
#> [3,]    4    606  390

Prior Model

priorLnAlphas <- list(log(.5),log(.5))
priorBetas <- list(1,-1)
prior <- calcDPCTable(pLevels,obsLevels,priorLnAlphas,priorBetas,
                      rules=rules,link=link)
round(prior,3)
#>       Gold Silver  None
#> [1,] 0.278  0.668 0.054
#> [2,] 0.134  0.732 0.134
#> [3,] 0.054  0.668 0.278

Model fitting

map1 <- mapDPC(counts,pLevels,obsLevels,
               priorLnAlphas,priorBetas,rules,link,gamma=.001)
if (map1$convergence != 0) {
  warning(paste("Optimization did not converge:", map1$message))
}

Parameter Recovery

postLnAlphas <- map1$lnAlphas
names(postLnAlphas) <- names(trueLnAlphas)
cat("True LnAlphas:",sapply(trueLnAlphas,round,3),".\n")
#> True LnAlphas: 0 -1.386 .
cat("Est LnAlphas:",sapply(postLnAlphas,round,3),".\n")
#> Est LnAlphas: -0.004 -1.39 .
all.equal(trueLnAlphas,postLnAlphas,tolerance=tol)
#> [1] "Component \"gold\": Mean absolute difference: 0.00371074"   
#> [2] "Component \"silver\": Mean relative difference: 0.002950396"
postBetas <- map1$betas
names(postBetas) <- names(trueBetas)
cat("True Betas:",sapply(trueBetas,round,3),".\n")
#> True Betas: 2 -0.5 .
cat("Est Betas:",sapply(postBetas,round,3),".\n")
#> Est Betas: 1.995 -0.5 .
all.equal(trueBetas,postBetas,tolerance=tol)
#> [1] "Component \"gold\": Mean relative difference: 0.002453082"

CPT Recovery

fitdist <- calcDPCTable(pLevels,obsLevels,map1$lnAlphas,
                        map1$betas,rules,link)
round(fitdist,3)
#>       Gold Silver  None
#> [1,] 0.119  0.686 0.195
#> [2,] 0.023  0.684 0.293
#> [3,] 0.004  0.606 0.390
cat("Maximum difference is ",max(abs(fitdist-truedist)),".\n")
#> Maximum difference is  0.0004223954 .
cat("KL divergence",cptKL(fitdist,truedist),".\n")
#> KL divergence 1.934476e-06 .