R/sirtr.plot.icc.R

In eatRest: eatRest

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Author:    Alexander Robitzsch
#			 a.robitzsch@bifie.at
#			 https://www.bifie.at/user/robitzsch-alexander
#			 
# Quelle:  	 Kopie aus sirtr_0.18-15
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



#----------------------------------------------------------------------------------------------------
np.dich <- function( dat , theta , thetagrid , progress = F , bwscale = 1.1 , method = "binomial"){
    #.....................................................................
    # 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          
    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 == T){ cat( "Item")}
    for (ii in 1:I){
        dfr <- data.frame( dat[,ii ] , theta )
        missing.resp[ , ii] <- 1* ( rowSums( is.na(dfr) ) > 0 )
        dfr <- na.omit(dfr)
        y <- dfr[,1]
        x <- dfr[,2]
        if (method == "normal" ){
            h1 <- ksmooth( x , y , bandwidth = bwscale * sd( x ) * length(x)^(-1/5)  , x.points = thetagrid , kernel="normal")
            estimate[ii,] <- h1$y
                    }
        if (method == "binomial"){      
                hsel <- h.select(x = x , y = y )
                h1 <- 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 == T){ if ( ii %% 20 == 0) { cat("" , ii , "\n" ) } else { cat( " " , ii ) } 
        flush.console() }
        }
    if (method == "normal"){ lower <- upper <- NULL }
    if (progress == T){ cat( "\n")}
    return( list( "dat" = dat , "thetagrid" = thetagrid , "theta" = theta , 
            "estimate" = estimate , "lower" = lower , "upper" = upper , "missing.resp" = missing.resp ) )
    }
#----------------------------------------------------------------------------------------------------


#-------------------------------------------------------------------------------------------#
# plotting Rasch ICCs and their corresponding nonparametric estimates                       #
plot.rasch.np <- function( np.dich.obj , b , infit = NULL , outfit = NULL , nsize = 100 , askplot = T , progress = T , 
                            bands = F , plot.b = F , shade = F , shadecol="burlywood1" , icc.pdfs = NULL ){
    # 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                                      #

    # 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){
	
		if ( !is.null ( icc.pdfs ) ) pdf ( icc.pdfs[ii] )
		
		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") ) }
        plot( thetagrid , estimate[ii,] , type="n"  , ylim = c( - .02 , 1.02 ) ,
                    xlab = expression( theta ) , ylab = expression( P( theta )) , 
                    main = plotname )
        abline( h = 0 , lty=3 )
        abline( h = 1 , lty=3 )
        rug( theta[ missing.resp[,ii] == 0 ]  , col= "cyan" )
        # shade area
        if (shade == T){ .shade.area( xval = thetagrid  , yval1 = p[,ii] , yval2 = estimate[ii,] , 
                            eps=0 , shadecol= shadecol )
                      }
        if (bands ==T & !is.null(upper) ){ 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 ) )
        lines( thetagrid[ind] , estimate[  ii, ind ] , lwd= 2 , col = 1 ,  lty= 5 )
        lines( thetagrid , p[ , ii ] , lwd= 2, col = 2 )
        if (plot.b == T){   lines( c( b[ii] , b[ii] ) , c( 0 , .5) , col=3 , lty=4 , lwd=2)
                            lines( c( -90 , b[ii] ) , c( .5 , .5) , col=3 , lty=4 ,  lwd=2 )        
                            th80 <- qlogis(.8) + b[ii]
                            lines( c( th80 , th80 ) , c( 0 , .8 ) , col=3 , lty=3 ,  lwd=2 )                                    
                            lines( c(  th80 , 90 ) , c( .8 , .8 ) , col=3 , lty=3 ,  lwd=2 ) 
                            lines( c( b[ii] , th80) , c(0,0) , col=3 , lwd=2 )             
                         }
        par(ask = askplot)
        if ( !is.null ( icc.pdfs ) ) dev.off()
		}
    par(ask=F)
    }
#-------------------------------------------------------------------------------------------#