## File Name: rasch.jml.biascorr.R
## File Version: 0.12
#########################################################################################
# analytical bias correction of JML estimated Rasch model
# rasch.jml.biascorr <- function( data, data.resp, theta, b, adjmin=.5 ){
rasch.jml.biascorr <- function( jmlobj, itemfac=NULL #, adjmin=.5
){
mod <- jmlobj
m1 <- mod$data.proc
data <- m1$data
data.resp <- m1$data.resp
theta <- m1$theta
b <- mod$item$itemdiff
I <- length(b)
I1 <- mean( rowSums( data.resp ) )
if (! is.null( itemfac) ){ I1 <- itemfac }
N <- nrow(data)
p_it <- stats::plogis( outer( theta, b, "-" ) )
q_it <- 1 - p_it
u_it <- - data + p_it # d/d(theta) f
v_it <- - u_it # d/d(alpha) f
# auxiliary calculations
u_it_theta <- - p_it * q_it
u_it_alpha <- p_it * q_it
v_it_alpha <- - p_it * q_it
u_it_alpha2 <- p_it * q_it^2 - p_it^2 * q_it
v_it_alpha2 <- - p_it * q_it^2 + p_it^2 * q_it
# terms
s2 <- rowWgtMeans( v_it_alpha, data.resp )
t00 <- rowWgtMeans( v_it^2, data.resp )
t1z <- - 1 / t00
t0 <- - t00 / s2
t1b <- rowWgtMeans( v_it * v_it_alpha, data.resp )
t2b <- rowWgtMeans( v_it_alpha2, data.resp )
t2c2 <- t1c2 <- t2c2 <- s2
t2z <- 1/2 / t1c2
b.bias2 <- b.bias1 <- rep(NA, I )
M1 <- matrix( 1, nrow=N,ncol=I)
for (ii in 1:I){
# ii <- 1
M1.ii <- M1
M1.ii[,ii] <- 0
M1.ii <- M1 - M1.ii
data.resp.ii <- data.resp[,ii]
#******
# bias correction method 1
s1 <- rowWgtMeans( u_it_alpha * M1.ii, data.resp )
s0 <- rowWgtMeans( u_it_theta * M1.ii, data.resp )
I_i <- sum( - ( s0 - s1^2 / s2 ) * data.resp.ii / sum(data.resp.ii) )
t1a <- rowWgtMeans( v_it * u_it_alpha * M1.ii, data.resp )
t2c1 <- t1c1 <- rowWgtMeans( u_it_alpha* M1.ii, data.resp )
t2a <- rowWgtMeans( u_it_alpha2* M1.ii, data.resp )
# term
btheta <- t0 * ( - t1z * ( t1a - t1b * t1c1 / t1c2 )
- t2z * ( t2a - t2b * t2c1 / t2c2 ) )
btheta <- sum( btheta * data.resp.ii / sum(data.resp.ii ) )
b.bias1[ii] <- btheta / I_i
#******
# bias correction method 2
a1 <- rowWgtMeans( u_it^2 * M1.ii, data.resp )
a4 <- a2 <- rowWgtMeans( u_it * v_it * M1.ii, data.resp )
a3 <- rowWgtMeans( v_it^2, data.resp )
I_i <- sum( - ( a1 - a2 * a3 * a4 ) )
U_it <- u_it * M1.ii - a2 / a3 * v_it
V_2it <- v_it^2 + v_it_alpha
t1 <- rowSums( U_it * V_2it * data.resp )
t2 <- 2*rowSums( v_it_alpha * data.resp )
b.bias2[ii] <- - sum( t1/t2 ) / I_i
}
b.analytcorr1 <- b - b.bias1 / I1
b.analytcorr2 <- b - b.bias2 / I1
# ind <- which( abs(b) > adjmin )
# b.analytcorr1b <- mean( ( b.analytcorr1a / b )[ind] )*b
# b.analytcorr2b <- mean( ( b.analytcorr2a / b )[ind] )*b
dfr <- data.frame( "b.JML"=b, "b.JMLcorr"=(I1-1)/I1*b,
b.analytcorr1, # b.analytcorr1b,
b.analytcorr2 #, b.analytcorr2b
)
res <- list( "b.biascorr"=dfr, "b.bias1"=b.bias1, "b.bias2"=b.bias2,
"itemfac"=I1 )
return(res)
}
#########################################################################################
# auxiliary function for bias correction
rowWgtMeans <- function( M, W ){
rowSums( M*W ) / rowSums(W)
}
#........................................................................................
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.