## File Name: tetrachoric2.R
## File Version: 1.343
tetrachoric2 <- function( dat, method="Ol", delta=.007, maxit=1000000,
cor.smooth=TRUE, progress=TRUE)
{
# process data
dat <- as.matrix(dat)
if (method=="Ol"){
res <- polychoric2( dat=dat, cor.smooth=cor.smooth, maxiter=maxit)
}
if ( method !="Ol" ){
# calculate tau
tau <- - stats::qnorm( colMeans(dat,na.rm=TRUE ) )
dat.resp <- 1 - is.na(dat)
dat[ dat.resp==0] <- 9
I <- ncol(dat)
# calculate frequencies
dfr <- data.frame( item1=rep(1:I,I), item2=rep(1:I, each=I ) )
h1 <- t( ( dat==1 ) ) %*% ( dat==0 )
dfr$f11 <- matrix( t( ( dat==1 ) ) %*% ( dat==1 ), ncol=1, byrow=TRUE ) + .5
dfr$f10 <- matrix( t( ( dat==1 ) ) %*% ( dat==0 ), ncol=1, byrow=TRUE ) + .5
dfr$f01 <- matrix( t( ( dat==0 ) ) %*% ( dat==1 ), ncol=1, byrow=TRUE ) + .5
dfr$f00 <- matrix( t( ( dat==0 ) ) %*% ( dat==0 ), ncol=1, byrow=TRUE ) + .5
dfr$ftot <- dfr$f11 + dfr$f10 + dfr$f01 + dfr$f00
dfr$p11 <- dfr$f11 / dfr$ftot
dfr$pi1 <- ( dfr$f11 + dfr$f10 - 1 ) / dfr$ftot
dfr$pi2 <- ( dfr$f11 + dfr$f01 - 1) / dfr$ftot
# subdata of dfr
dfr <- dfr[ dfr$item1 > dfr$item2, ]
dfr <- dfr[ dfr$ftot > 0, ]
dfr$qi1 <- stats::qnorm( dfr$pi1)
dfr$qi2 <- stats::qnorm( dfr$pi2)
#******************
# method of Bonett
if ( method %in% c("Bo","Di") ){
dfr$pmin <- ifelse( dfr$pi1 < dfr$pi2, dfr$pi1, dfr$pi2 )
dfr$c <- ( 1 - abs( dfr$pi1 - dfr$pi2 ) / 5 - ( 0.5 - dfr$pmin)^2 ) / 2
dfr$omega <- ( dfr$f00 * dfr$f11 ) / ( dfr$f01 * dfr$f10)
dfr$r0 <- cos( pi / ( 1 + dfr$omega^( dfr$c ) ) )
}
#*****
# method of Divgi
if ( method=="Di"){
dfr2 <- as.matrix(dfr)
numdiffparm <- 1e-5
dfr$r0 <- sirt_rcpp_tetrachoric2_rcpp( dfr=dfr2, h=numdiffparm, maxiter=maxit )
}
#******************
# method of Tucker
if ( method=="Tu"){
# functions defined by Cengiz Zopluoglu
L <- function(r,h,k) {(1/(2*pi*sqrt(1-r^2)))*exp(-((h^2-2*h*k*r+k^2)/(2*(1-r^2))))}
S <- function(x,r,h,k) { x-(L(r,h,k)*delta) }
A0 <- dfr$A0 <- dfr$p11 - dfr$pi1 * dfr$pi2
dfr$r0 <- delta / 2
dfr$iter <- 0
dfr$iter <- 0
dfr$conv <- 0
ii <- 0
vars <- c("A0","r0","iter","conv")
if (progress){ cat(" 0 : " ) }
while( ( mean( dfr$conv) < 1 ) & ( ii < maxit ) ){
# iterations
dfr0 <- dfr
#***
ii <- ii+1
dfr$A0 <- dfr$A0 - delta * L( dfr$r0, dfr$qi1, dfr$qi2 )
dfr$r0 <- dfr$r0 + delta
dfr$iter <- ii
ind <- which( dfr$A0 < 0 )
if ( length(ind) > 0 ){
h1 <- dfr$r0 - delta/2 + dfr$A0 / L( dfr0$r0, dfr$qi1, dfr$qi2 )
dfr$r0[ind] <- h1[ind]
dfr$conv[ind] <- 1
}
i2 <- which( dfr$conv==1 & dfr0$conv==1 )
if (length(i2) > 0){ dfr[ i2, vars] <- dfr0[ i2, vars] }
if (progress){
if (ii %% 100==0 ){ cat(".") ; utils::flush.console() }
if (ii %% 1000==0 ){ cat("\n",ii, ": ") }
}
}
#******************
cat("\n")
}
TC <- matrix(NA, I, I )
diag(TC) <- 1
TC[ as.matrix(dfr[, c("item1","item2") ] ) ] <- dfr$r0
TC[ as.matrix(dfr[, c("item2","item1") ] ) ] <- dfr$r0
if (cor.smooth){
TC <- sirt_import_psych_cor.smooth(x=TC)
}
rownames(TC) <- colnames(TC) <- colnames(dat)
res <- list(tau=tau, rho=TC )
} # method !="Ol"
return(res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.