result <- list("pairs" =
c(row.names(cheaters)), "Ncheat" = nrow(cheaters),
"TotalCompare" = totalCompare, "alpha" = alpha,
"ExactMatch" = cheaters$Nexact, "Zobs" =
cheaters$Zobs, "Zcrit" = Zcrit,
"Mean" = cheaters$Mean, "Variance" =
cheaters$Var, "Probs" = stuProbs)
result
and the only place Var1 appears in the plot function is here
prop.correct <- subset(data.frame(prop.table(table(tmp[,
i+1], tmp$Estimate), margin=2)), Var1 == 1)[, 2:3]
Many thanks,
Harold
> sessionInfo()
R version 2.10.0 (2009-10-26)
i386-pc-mingw32
locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
cheat.fit <- function(dat, key, wrongChoice, alpha = .01, rfa
= c('nr', 'uni', 'bsct'), bonf = c('yes','no'), con = 1e-12,
lower = 0, upper = 50){
bonf <- tolower(bonf)
bonf <- match.arg(bonf)
rfa <- match.arg(rfa)
rfa <- tolower(rfa)
dat <- t(dat)
correctStuMat <- numeric(ncol(dat))
for(i in 1:ncol(dat)){
correctStuMat[i] <-
mean(key==dat[,i], na.rm= TRUE)
}
correctClsMat <- numeric(length(key))
for(i in 1:length(key)){
correctClsMat[i] <-
mean(key[i]==dat[i,], na.rm= TRUE)
}
### this is here for cases if all students in a class
### did not answer the item
correctClsMat[is.na(correctClsMat)] <- 0
pCorr <- function(R,c,q){
numer <- function(R,a,c,q){
result <-
sum((1-(1-R)^a)^(1/a), na.rm= TRUE)-c*q
result
}
denom <- function(R,a,c,q){
result <- sum(na.rm= TRUE,
-((1 - (1 - R)^a)^(1/a) * (log((1 - (1 - R)^a)) * (1/a^2)) +
(1 - (1 - R)^a)^((1/a) - 1) *
((1/a) * ((1 - R)^a * log((1 - R))))))
result
}
aConst <- function(R, c, q, con){
a <- .5 # starting value for a
change <- 1
while(abs(change) > con) {
r1 <- numer(R,a,c,q)
r2 <- denom(R,a,c,q)
change <- r1/r2
a <- a - change
}
a
}
bisect <- function(R, c, q, lower, upper, con){
f <- function(a) sum((1 -
(1-R)^a)^(1/a)) - c * q
if(f(lower) * f(upper) > 0)
stop("endpoints must have opposite signs")
while(abs(lower-upper) > con){
x = .5 * (lower+upper)
if(f(x) *
f(lower) >=0) lower = x
else upper = x
}
.5 * (lower+upper)
}
if(rfa == 'nr'){
if(any(correctClsMat==1))
correctClsMat[correctClsMat==1]<-.9999 else correctClsMat
if(any(correctClsMat==0))
correctClsMat[correctClsMat==0]<-.0001 else correctClsMat
a <- aConst(R,c,q, con)
} else if(rfa == 'uni'){
f <- function(R, a, c, q)
sum((1 - (1-R)^a)^(1/a)) - c * q
a <- uniroot(f,
c(lower,upper), R = R, c = c, q = q)$root
} else if(rfa == 'bsct'){
a <- bisect(R, c, q, lower =
lower, upper = upper, con)
}
result <- (1-(1-R)^a)^(1/a)
result
} # end pCorr function
stuProbs <- matrix(0, ncol=ncol(dat), nrow=nrow(dat))
for(i in 1:ncol(dat)){
if(correctStuMat[i]==1){
stuProbs[,i] <- 1
} else if(correctStuMat[i]==0){
stuProbs[,i] <- 0
} else {
stuProbs[,i] <-
pCorr(correctClsMat, correctStuMat[i], q = length(key))
}
}
stuProbs[which(is.na(dat))] <- NA # this is a
bit kludgy. But, it places NAs in the right place
matchProb <- function(StuMat, wrongMat, numItems){
ind <- combn(c(1:ncol(StuMat)),2) # These are
all the possible combinations
result <- matrix(0, ncol=ncol(ind), nrow=numItems)
for(j in 1:ncol(ind)){
for(i in 1:numItems){
if(is.na(StuMat[i,ind[1,j]]) | is.na(StuMat[i,ind[2,j]]) ) {
result[i,j] <- NA
} else { result[i,j] <- StuMat[i,ind[1,j]] * StuMat[i,ind[2,j]] +
(1-StuMat[i,ind[1,j]]) * (1-StuMat[i,ind[2,j]]) *
sum(wrongChoice[[i]]^2)
}
}
}
result <- data.frame(result)
names(result) <- paste('S',paste(ind[1,],
ind[2,], sep=''), sep='')
result
}
aa <- matchProb(stuProbs, wrongChoice,
numItems = nrow(dat))
matchTotal <- function(dat){
ind <- combn(c(1:ncol(dat)),2) # These are
all the possible combinations
Match <- numeric(ncol(ind))
for(j in 1:ncol(ind)){
Match[j] <- sum(dat[,
ind[1,j]] == dat[, ind[2,j]], na.rm=TRUE) ### added this just now
}
Match <- data.frame(Match)
row.names(Match) <- paste('S',paste(ind[1,],
ind[2,], sep=':'), sep='')
Match
}
exactMatch <- matchTotal(dat)
means <- colSums(aa, na.rm=TRUE)
vars <- colSums(aa*(1-aa), na.rm= TRUE)
cheaters <- (exactMatch - .5 - means)/sqrt(vars)
totalCompare <- nrow(cheaters)
if(bonf == 'yes'){
alpha <- alpha/nrow(cheaters)
Zcrit <- qnorm(alpha, mean =
0, sd = 1, lower.tail = FALSE)
} else {
Zcrit <- qnorm(alpha, mean = 0, sd = 1, lower.tail = FALSE)
}
cheaters <- cbind(data.frame(cheaters), exactMatch)
names(cheaters)[1] <- 'Zobs'
names(cheaters)[2] <- 'Nexact'
cheaters$Zcrit <- Zcrit
cheaters$Mean <- means
cheaters$Var <- vars
cheaters <- subset(cheaters, Zobs >= Zcrit)
result <- list("pairs" =
c(row.names(cheaters)), "Ncheat" = nrow(cheaters),
"TotalCompare" = totalCompare, "alpha" = alpha,
"ExactMatch" = cheaters$Nexact, "Zobs" =
cheaters$Zobs, "Zcrit" = Zcrit,
"Mean" = cheaters$Mean, "Variance" =
cheaters$Var, "Probs" = stuProbs)
result
}
plot.jml <- function(x, ask = TRUE, all = TRUE, item, ...){
par(ask = ask)
xvals <- seq(from=-5, to =5, by=.01)
params <- coef(x)
L <- length(params)
tmp <- data.frame(x$model.frame)
tmp$Raw.Score <- rowSums(tmp)
mle <- summary(scoreCon(coef(x)))$coef[c(1,3)]
tmp <- merge(tmp, mle, by='Raw.Score', sort = FALSE)
if(all){
for(i in 1:L ){
exp <- 1/(1 +
exp(params[i] - xvals))
prop.correct
<- subset(data.frame(prop.table(table(tmp[, i+1],
tmp$Estimate), margin=2)),
Var1 == 1)[, 2:3]
names(prop.correct)[1:2] <- c('Estimate', 'prop')
prop.correct$Estimate <- as.numeric(levels(prop.correct$Estimate))
plot(xvals,
exp, type='l', ylim = c(0,1), xlab = 'Ability', ylab =
'Proportion Correct', ...)
points(prop.correct$Estimate, prop.correct$prop, ...)
title(main =
paste("Plot of Item", i))
}
} else {
par(ask = FALSE)
i <- item
exp <- 1/(1 + exp(params[i] - xvals))
prop.correct <-
subset(data.frame(prop.table(table(tmp[, i+1], tmp$Estimate),
margin=2)),
Var1 == 1)[, 2:3]
names(prop.correct)[1:2] <-
c('Estimate', 'prop')
prop.correct$Estimate <-
as.numeric(levels(prop.correct$Estimate))
plot(xvals, exp, type='l',
ylim = c(0,1), xlab = 'Ability', ylab = 'Proportion Correct', ...)
points(prop.correct$Estimate,
prop.correct$prop, ...)
title(main = paste("Plot of Item", i))
}
par(ask = FALSE)
}
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
