## File Name: np.dich.R
## File Version: 0.30
#----------------------------------------------------------------------------------------------------
np.dich <- function( dat, theta, thetagrid, progress=FALSE,
bwscale=1.1, method="normal"){
#.....................................................................
# INPUT:
# dat ... data frame
# theta ... estimated theta values
# thetagrid ... grid on which ability is nonparametrically evaluated #
# method ... binomial and normal (nonparametric binomial or normal nonparametric regression)
# bwscale ... bandwidth scale h=bwscale * ^N^{-1/5) #
#......................................................................
# recode values of infinity to missing
TAM::require_namespace_msg("sm")
theta[ theta %in% c( -Inf, Inf ) ] <- NA
I <- ncol(dat)
# display
cat("Nonparametric ICC Estimation for Dichotomous Data \n" )
cat(I, "Items \n" )
# matrix of lower and upper confidence limits and estimates
lower <- upper <- estimate <- matrix( 0, nrow=I, ncol=length(thetagrid) )
missing.resp <- dat
if (progress){ cat( "Item")}
for (ii in 1:I){
dfr <- data.frame( dat[,ii ], theta )
missing.resp[, ii] <- 1* ( rowSums( is.na(dfr) ) > 0 )
dfr <- stats::na.omit(dfr)
y <- dfr[,1]
x <- dfr[,2]
if (method=="normal" ){
h1 <- stats::ksmooth( x, y, bandwidth=bwscale * sd( x ) * length(x)^(-1/5),
x.points=thetagrid, kernel="normal")
estimate[ii,] <- h1$y
}
if (method=="binomial"){
hsel <- sm::h.select(x=x, y=y )
h1 <- sm::sm.binomial( x, y, N=rep(1, length(x) ), h=hsel,
eval.points=thetagrid, display="none" )
estimate[ ii,] <- h1$estimate
lower[ii,] <- h1$lower
upper[ii,] <- h1$upper
}
if (progress ){ if ( ii %% 20==0) { cat("", ii, "\n" ) } else { cat( "", ii ) }
utils::flush.console() }
}
if (method=="normal"){ lower <- upper <- NULL }
if (progress){ cat( "\n")}
res <- list( "dat"=dat, "thetagrid"=thetagrid, "theta"=theta,
"estimate"=estimate, "lower"=lower, "upper"=upper,
"missing.resp"=missing.resp )
class(res) <- "np.dich"
return(res)
}
#----------------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------
# plotting function for shades between 2 curves
.shade.area <- function( xval, yval1, yval2, shadecol="gray",eps=0){
# INPUT:
# xval .. x values
# yval1/yval2 ... vectors of y curves
# shadecol ... color for shading the area
yvaldiff <- 1*( yval1 - yval2 )
xx <- c( xval, rev(xval) )
yy <- c( yval1 - eps*yvaldiff, rev( yval2 ) - eps * yvaldiff )
graphics::polygon(xx, yy, col=shadecol, border=shadecol )
}
#--------------------------------------------------------------------------------------------------
##NS # S3method(plot, rasch.np)
#-------------------------------------------------------------------------------------------#
# plotting Rasch ICCs and their corresponding nonparametric estimates #
plot.np.dich <- function( x, b, infit=NULL, outfit=NULL, nsize=100,
askplot=TRUE, progress=TRUE,
bands=FALSE, plot.b=FALSE, shade=FALSE, shadecol="burlywood1", ...){
# INPUT: #
# np.dich.obj... object generated by np.dich #
# theta ... ability estimates obtained from Rasch analysis #
# b ... Rasch item difficulties #
# thetalim ... limits of theta values #
# nsize ... number of grid points #
# askplot ... should before every plot asked to go to the next one? #
# progress ... display iteration process of nonparametric estimation #
# bands ... should confidence bands be drawn? #
# plot.b ... plot item difficulties #
P <- NULL
np.dich.obj <- object <- x
# extract informations from object np.dich
thetagrid <- np.dich.obj$thetagrid
theta <- np.dich.obj$theta
nsize <- length(thetagrid)
dat <- np.dich.obj$dat
estimate <- np.dich.obj$estimate
lower <- np.dich.obj$lower
upper <- np.dich.obj$upper
missing.resp <- np.dich.obj$missing.resp
I <- ncol(dat)
n <- nrow(dat)
p <- .prob.rasch( thetagrid, b )
for (ii in 1:I){
plotname <- paste( colnames(dat)[ii], ": N=",
sum((1-is.na(dat))[,ii]),
"; b=", formatC( b[ii], 2, format="f"), sep="")
if ( !is.null(infit) ){ plotname <- paste( plotname, "\n Infit=", formatC( infit[ii],2, format="f") ) }
if ( !is.null(infit) ){ plotname <- paste( plotname, " Outfit=", formatC( outfit[ii],2, format="f") ) }
graphics::plot( thetagrid, estimate[ii,], type="n", ylim=c( - .02, 1.02 ),
xlab=expression( theta ), ylab=expression( P( theta )),
main=plotname )
graphics::abline( h=0, lty=3 )
graphics::abline( h=1, lty=3 )
# graphics::rug( theta[ missing.resp[,ii]==0 ], col="cyan" )
# shade area
if (shade ){
.shade.area( xval=thetagrid, yval1=p[,ii], yval2=estimate[ii,],
eps=0, shadecol=shadecol )
}
if (bands & !is.null(upper) ){
graphics::segments( thetagrid, lower[ii,], thetagrid,
upper[ii,], col="gray")
}
# select display
theta.ii <- theta[ missing.resp[,ii]==0 ]
ind <- which( thetagrid > min(theta.ii) & thetagrid < max(theta.ii ) )
graphics::lines( thetagrid[ind], estimate[ ii, ind ], lwd=2, col=1, lty=5 )
graphics::lines( thetagrid, p[, ii ], lwd=2, col=2 )
if (plot.b ){
graphics::lines( c( b[ii], b[ii] ), c( 0, .5), col=3, lty=4, lwd=2)
graphics::lines( c( -90, b[ii] ), c( .5, .5), col=3, lty=4, lwd=2 )
th80 <- stats::qlogis(.8) + b[ii]
graphics::lines( c( th80, th80 ), c( 0, .8 ), col=3, lty=3, lwd=2 )
graphics::lines( c( th80, 90 ), c( .8, .8 ), col=3, lty=3, lwd=2 )
graphics::lines( c( b[ii], th80), c(0,0), col=3, lwd=2 )
}
graphics::par(ask=askplot)
}
graphics::par(! askplot )
}
#-------------------------------------------------------------------------------------------#
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.