R/np.dich.R

In alexanderrobitzsch/sirt: Supplementary Item Response Theory Models

## 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 )
    }
#-------------------------------------------------------------------------------------------#