11# # File Name: tam.mml.wle2.R
2- # # File Version: 0.853
2+ # # File Version: 0.869
33
44# ###############################################################
55tam.mml.wle2 <- function ( tamobj , score.resp = NULL , WLE = TRUE , adj = .3 , Msteps = 20 ,
6- convM = .0001 , progress = TRUE , output.prob = FALSE , pid = NULL )
6+ convM = .0001 , progress = TRUE , output.prob = FALSE , pid = NULL ,
7+ theta_init = NULL )
78{
89 CALL <- match.call()
10+ iweights <- NULL
911 # --- process input data
1012 res <- tam_mml_wle_proc_input_data( tamobj = tamobj , score.resp = score.resp )
1113 AXsi <- res $ AXsi
@@ -37,23 +39,33 @@ tam.mml.wle2 <- function( tamobj, score.resp=NULL, WLE=TRUE, adj=.3, Msteps=20,
3739 col.index <- rep( 1 : nitems , each = maxK )
3840 cResp <- resp [, col.index ]* resp.ind [, col.index ]
3941 cResp <- 1 * t( t(cResp )== rep(0 : (maxK - 1 ), nitems ) )
40- cB <- t( matrix ( aperm( B , c(2 ,1 ,3 ) ), nrow = dim(B )[3 ], byrow = TRUE ) )
42+ if (! is.null(iweights )){
43+ B1 <- B * iweights
44+ } else {
45+ B1 <- B
46+ iweights <- rep(1 , nitems )
47+ }
48+
49+ cB <- t( matrix ( aperm( B1 , c(2 ,1 ,3 ) ), nrow = dim(B )[3 ], byrow = TRUE ) )
4150 cB [is.na(cB )] <- 0
4251 # Compute person sufficient statistics (total score on each dimension)
4352 PersonScores <- cResp %*% cB
4453
4554 # Compute possible maximum score for each item on each dimension
46- maxBi <- apply(B , 3 , tam_rowMaxs , na.rm = TRUE )
55+ maxBi <- apply(B1 , 3 , tam_rowMaxs , na.rm = TRUE )
56+
4757
4858 # Compute possible maximum score for each person on each dimension
4959 PersonMax <- resp.ind %*% maxBi
5060 PersonMax [ PersonMax == 0 ] <- 2 * adj
5161
5262 # Adjust perfect scores for each person on each dimension
53- PersonScores [PersonScores == PersonMax ] <- PersonScores [PersonScores == PersonMax ] - adj
63+ ind_max <- which(PersonScores == PersonMax )
64+ PersonScores [ind_max ] <- PersonScores [ind_max ] - adj
5465
5566 # Adjust zero scores for each person on each dimension
56- PersonScores [PersonScores == 0 ] <- PersonScores [PersonScores == 0 ] + adj
67+ ind0 <- which(PersonScores == 0 )
68+ PersonScores [ind0 ] <- PersonScores [ind0 ] + adj
5769
5870 # Calculate Axsi. Only need to do this once.
5971 # for (i in 1:nitems) {
@@ -63,7 +75,12 @@ tam.mml.wle2 <- function( tamobj, score.resp=NULL, WLE=TRUE, adj=.3, Msteps=20,
6375 # }
6476
6577 # Initialise theta (WLE) values for all students
66- theta <- log((PersonScores + .5 )/ (PersonMax - PersonScores + 1 )) # log of odds ratio of raw score
78+ if (is.null(theta_init )){
79+ theta <- log((PersonScores + .5 )/ (PersonMax - PersonScores + 1 ))
80+ # log of odds ratio of raw score
81+ } else {
82+ theta <- as.matrix( theta_init )
83+ }
6784
6885 # #####################################
6986 # Compute WLE
@@ -92,12 +109,13 @@ tam.mml.wle2 <- function( tamobj, score.resp=NULL, WLE=TRUE, adj=.3, Msteps=20,
92109 xsi = xsi , theta = theta , nnodes = nstud , maxK = maxK , recalc = FALSE ,
93110 use_rcpp = TRUE , maxcat = max(maxK ), avoid_outer = TRUE )
94111 rprobsWLE <- resWLE $ rprobs
112+
95113 rprobsWLEL <- matrix (rprobsWLE , nrow = nitems * maxK , ncol = nstud )
96114 rprobsWLEL [is.na(rprobsWLEL )] <- 0
97115
98116 resB <- tam_rcpp_wle_suffstat( RPROBS = rprobsWLEL , CBL = BL , CBB = BBL ,
99- CBBB = BBBL , cndim = ndim , cnitems = nitems , cmaxK = maxK , cnstud = nstud ,
100- resp_ind = resp.ind )
117+ CBBB = BBBL , cndim = ndim , cnitems = nitems , cmaxK = maxK ,
118+ cnstud = nstud , resp_ind = resp.ind )
101119 B_bari <- array (resB $ B_bari , dim = c(nstud , nitems ,ndim ))
102120 BB_bari <- array (resB $ BB_bari , dim = c(nstud , nitems , ndim , ndim ))
103121 BBB_bari <- array (resB $ BBB_bari , dim = c(nstud , nitems , ndim ))
@@ -145,7 +163,7 @@ tam.mml.wle2 <- function( tamobj, score.resp=NULL, WLE=TRUE, adj=.3, Msteps=20,
145163 # dampening the increment
146164 for ( d1 in 1 : ndim ){
147165 # increment[,d1] <- ifelse( abs(increment[,d1]) > 3, sign( increment[,d1] )*3, increment[,d1] )
148- ci <- ceiling( abs(increment [,d1 ]) / ( abs( old_increment [,d1 ]) + 10 ^ ( - 10 ) ) )
166+ ci <- ceiling( abs(increment [,d1 ]) / ( abs( old_increment [,d1 ]) + 1e -10 ) )
149167 increment [,d1 ] <- ifelse( abs( increment [,d1 ]) > abs(old_increment [,d1 ]),
150168 increment [,d1 ]/ (2 * ci ),
151169 increment [,d1 ] )
@@ -165,15 +183,16 @@ tam.mml.wle2 <- function( tamobj, score.resp=NULL, WLE=TRUE, adj=.3, Msteps=20,
165183 Miter <- Miter + 1
166184
167185 if (progress ){
168- cat( paste( " Iteration in WLE/MLE estimation " , Miter ,
169- " | Maximal change " , round( max(abs(increment )), 4 ), " \n " ) )
170- utils :: flush.console()
171- }
186+ cat( paste( " Iteration in WLE/MLE estimation " , Miter ,
187+ " | Maximal change " , round( max(abs(increment )), 4 ), " \n " ))
188+ utils :: flush.console()
189+ }
172190 } # end of Newton-Raphson
173191
174192 res <- tam_mml_wle_postproc( ndim = ndim , err_inv = err_inv , theta = theta , pid = pid ,
175193 resp.ind = resp.ind , PersonScores = PersonScores , PersonMax = PersonMax ,
176- adj = adj , WLE = WLE , rprobsWLE = rprobsWLE , output.prob = output.prob , progress = progress ,
177- pweights = pweights , CALL = CALL , B = B , score.resp = score.resp )
194+ adj = adj , WLE = WLE , rprobsWLE = rprobsWLE , output.prob = output.prob ,
195+ progress = progress , pweights = pweights , CALL = CALL , B = B ,
196+ score.resp = score.resp )
178197 return (res )
179198}
0 commit comments