Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

robustify estimation a little more

  • Loading branch information...
commit 989b54275bde5cd0ed34ce6e542fbddd94af326c 1 parent 7efa84b
@philchalmers authored
View
4 NEWS
@@ -1,3 +1,7 @@
+* item probability functions now only permit permissible values, and models may converge even when
+ the log-likelihood decreases during estimation. In the EM if the model does not have a strictly
+ increasing log-likelihood then a warning message will be printed
+
* removed infit and outfit stats since they only are applicable to Rasch models (and are available in
the eRm package)
View
45 R/02-item_methods.R
@@ -56,6 +56,14 @@ setMethod(
)
setMethod(
+ f = "print",
+ signature = signature(x = 'GroupPars'),
+ definition = function(x, ...){
+ cat('Object of class:', class(x))
+ }
+)
+
+setMethod(
f = "show",
signature = signature(object = 'dich'),
definition = function(object){
@@ -111,6 +119,14 @@ setMethod(
}
)
+setMethod(
+ f = "show",
+ signature = signature(object = 'GroupPars'),
+ definition = function(object){
+ print(object)
+ }
+)
+
#----------------------------------------------------------------------------
#LogLik
@@ -118,8 +134,7 @@ setMethod(
f = "LogLik",
signature = signature(x = 'dich', Theta = 'matrix'),
definition = function(x, Theta, EM=FALSE, prior=NULL){
- itemtrace <- ProbTrace(x=x, Theta=Theta)
- itemtrace[itemtrace < 1e-8] <- 1e-8
+ itemtrace <- ProbTrace(x=x, Theta=Theta)
Prior <- rep(1, nrow(itemtrace))
if(EM) Prior <- prior
LL <- (-1) * sum(x@rs * log(itemtrace) * Prior)
@@ -132,8 +147,7 @@ setMethod(
f = "LogLik",
signature = signature(x = 'graded', Theta = 'matrix'),
definition = function(x, Theta, EM = FALSE, prior = NULL){
- itemtrace <- ProbTrace(x=x, Theta=Theta)
- itemtrace[itemtrace < 1e-8] <- 1e-8
+ itemtrace <- ProbTrace(x=x, Theta=Theta)
Prior <- rep(1, nrow(itemtrace))
if(EM) Prior <- prior
LL <- (-1) * sum(x@rs * log(itemtrace) * Prior)
@@ -146,8 +160,7 @@ setMethod(
f = "LogLik",
signature = signature(x = 'rating', Theta = 'matrix'),
definition = function(x, Theta, EM = FALSE, prior = NULL){
- itemtrace <- ProbTrace(x=x, Theta=Theta)
- itemtrace[itemtrace < 1e-8] <- 1e-8
+ itemtrace <- ProbTrace(x=x, Theta=Theta)
Prior <- rep(1, nrow(itemtrace))
if(EM) Prior <- prior
LL <- (-1) * sum(x@rs * log(itemtrace) * Prior)
@@ -161,7 +174,6 @@ setMethod(
signature = signature(x = 'gpcm', Theta = 'matrix'),
definition = function(x, Theta, EM = FALSE, prior = NULL){
itemtrace <- ProbTrace(x=x, Theta=Theta)
- itemtrace[itemtrace < 1e-8] <- 1e-8
Prior <- rep(1, nrow(itemtrace))
if(EM) Prior <- prior
LL <- (-1) * sum(x@rs * log(itemtrace) * Prior)
@@ -175,7 +187,6 @@ setMethod(
signature = signature(x = 'nominal', Theta = 'matrix'),
definition = function(x, Theta, EM = FALSE, prior = NULL){
itemtrace <- ProbTrace(x=x, Theta=Theta)
- itemtrace[itemtrace < 1e-8] <- 1e-8
Prior <- rep(1, nrow(itemtrace))
if(EM) Prior <- prior
LL <- (-1) * sum(x@rs * log(itemtrace) * Prior)
@@ -189,7 +200,6 @@ setMethod(
signature = signature(x = 'partcomp', Theta = 'matrix'),
definition = function(x, Theta, EM = FALSE, prior = NULL){
itemtrace <- ProbTrace(x=x, Theta=Theta)
- itemtrace[itemtrace < 1e-8] <- 1e-8
Prior <- rep(1, nrow(itemtrace))
if(EM) Prior <- prior
LL <- (-1) * sum(x@rs * log(itemtrace) * Prior)
@@ -203,7 +213,6 @@ setMethod(
signature = signature(x = 'mcm', Theta = 'matrix'),
definition = function(x, Theta, EM = FALSE, prior = NULL){
itemtrace <- ProbTrace(x=x, Theta=Theta)
- itemtrace[itemtrace < 1e-8] <- 1e-8
Prior <- rep(1, nrow(itemtrace))
if(EM) Prior <- prior
LL <- (-1) * sum(x@rs * log(itemtrace) * Prior)
@@ -475,15 +484,17 @@ P.poly <- function(a, d, Theta, itemexp = FALSE, D)
P <- matrix(0,nrow(Theta),ncat)
for(i in ncat:1)
P[ ,i] <- Pk[ ,i] - Pk[ ,i+1]
- Pk <- P
- }
+ P[P < 1e-8] <- 1e-8
+ P[(1 - P) < 1e-8] <- 1 - 1e-8
+ Pk <- P
+ }
return(Pk)
}
# Trace lines for mirt models
P.mirt <- function(a, d, Theta, g = 0, u = 1, D)
{
- traces <- .Call("traceLinePts", a, d, g, u, Theta, D) ###FIX RCPP CODE
+ traces <- .Call("traceLinePts", a, d, g, u, Theta, D)
return(traces)
}
@@ -495,6 +506,8 @@ P.comp <- function(a, d, Theta, g, u = 1, D)
for(i in 1:nfact)
P <- P * P.mirt(a[i], d[i], Theta[ ,i, drop=FALSE], g=0, u=1, D=D)
P <- g + (u - g) * P
+ P[P < 1e-8] <- 1e-8
+ P[(1 - P) < 1e-8] <- 1 - 1e-8
P
}
@@ -507,6 +520,8 @@ P.nominal <- function(a, ak, d, Theta, D, returnNum = FALSE){
for(i in 1:ncat)
numerator[ ,i] <- exp(D * ak[i] * (Theta %*% a) + D * d[i])
P <- numerator/rowSums(numerator)
+ P[P < 1e-8] <- 1e-8
+ P[(1 - P) < 1e-8] <- 1 - 1e-8
if(returnNum) return(numerator)
return(P)
}
@@ -521,6 +536,8 @@ P.gpcm <- function(a, d, Theta, D){
for(i in 1:ncat)
numerator[ ,i] <- exp(D * k[i] * (Theta %*% a) + D * d[i])
P <- numerator/rowSums(numerator)
+ P[P < 1e-8] <- 1e-8
+ P[(1 - P) < 1e-8] <- 1 - 1e-8
return(P)
}
@@ -541,6 +558,8 @@ P.mcm <- function(a, ak, d, t, Theta, D){
t[1] <- 1 - sum(t[2:length(t)])
T <- matrix(t, nrow(P), ncat, byrow = TRUE)
P <- C0 * T + (1 - C0) * numerator/denominator
+ P[P < 1e-8] <- 1e-8
+ P[(1 - P) < 1e-8] <- 1 - 1e-8
return(P)
}
View
2  R/04-derivatives.R
@@ -24,7 +24,6 @@ setMethod(
hess <- matrix(0,nfact+3, nfact+3)
if(x@par[parlength] == 1){ #'3PL'
Pstar <- P.mirt(a,d,Theta,0,1, D=x@D)
- Pstar[Pstar < 1e-8] <- 1e-8
Qstar <- 1 - Pstar
da <- rep(0,nfact)
dd <- sum((1-g)*Pstar*Qstar*(dat/P - (f-dat)/Q)*Prior)
@@ -603,7 +602,6 @@ setMethod(
ind <- x@nfact + 1
stop('Derivatives for mcm items not yet written')
grad <- hess <- NULL
-
return(list(grad=grad, hess=hess))
}
)
View
13 R/EMstep.group.R
@@ -45,6 +45,7 @@ EM.group <- function(pars, constrain, PrepList, list, Theta, debug)
}
stagecycle <- 1
converge <- 1
+ LLwarn <- FALSE
inverse_fail_count <- 1
##
L <- c()
@@ -73,7 +74,7 @@ EM.group <- function(pars, constrain, PrepList, list, Theta, debug)
for(g in 1:ngroups)
r[[g]] <- PrepList[[g]]$tabdata[, ncol(PrepList[[g]]$tabdata)]
#EM
- for (cycles in 1:NCYCLES){
+ for (cycles in 1:NCYCLES){
#priors
for(g in 1:ngroups){
gstructgrouppars[[g]] <- ExtractGroupPars(pars[[g]][[J+1]])
@@ -104,6 +105,7 @@ EM.group <- function(pars, constrain, PrepList, list, Theta, debug)
}
LL <- LL + sum(r[[g]]*log(rlist[[g]]$expected))
}
+ if(LL < lastLL && cycles > 1) LLwarn <- TRUE
for(g in 1:ngroups){
for(i in 1:J){
tmp <- c(itemloc[i]:(itemloc[i+1] - 1))
@@ -170,8 +172,10 @@ EM.group <- function(pars, constrain, PrepList, list, Theta, debug)
ind1 <- ind2 + 1
}
grad <- g %*% L
- hess <- L %*% h %*% L
- grad <- grad[1, estpars & !redun_constr]
+ hess <- L %*% h %*% L
+ grad <- grad[1, estpars & !redun_constr]
+ if(any(is.na(grad)))
+ stop('Model did not converge (unacceptable gradient caused by extreme parameter values)')
Hess <- Matrix(hess[estpars & !redun_constr, estpars & !redun_constr], sparse = TRUE)
inv.Hess <- try(solve(Hess), silent = TRUE)
if(class(inv.Hess) == 'try-error'){
@@ -212,6 +216,9 @@ EM.group <- function(pars, constrain, PrepList, list, Theta, debug)
listpars[[g]][[i]] <- pars[[g]][[i]]@par
} #END EM
+ if(LLwarn)
+ warning('Log-likelihood did not strictly decrease during estimation.
+ Solution may not be a maximum.')
ret <- list(pars=pars, cycles = cycles, info=matrix(0), longpars=longpars, converge=converge,
logLik=LL, rlist=rlist)
ret
View
7 R/MHRM.group.R
@@ -175,8 +175,11 @@ MHRM.group <- function(pars, constrain, PrepList, list, debug)
grad <- g %*% L
ave.h <- (-1)*L %*% h %*% L
grad <- grad[1, estpars & !redun_constr]
- ave.h <- ave.h[estpars & !redun_constr, estpars & !redun_constr]
- if(is.na(attr(gtheta0[[1]],"log.lik"))) stop('Estimation halted. Model did not converge.')
+ ave.h <- ave.h[estpars & !redun_constr, estpars & !redun_constr]
+ if(any(is.na(grad)))
+ stop('Model did not converge (unacceptable gradient caused by extreme parameter values)')
+ if(is.na(attr(gtheta0[[1]],"log.lik")))
+ stop('Estimation halted. Model did not converge.')
if(verbose){
if((cycles + 1) %% 10 == 0){
if(cycles < BURNIN)
View
12 R/confmirt.R
@@ -26,6 +26,15 @@
#' \code{anova} function, where a Chi-squared difference test and AIC/BIC
#' difference values are displayed.
#'
+#' @section Congergence:
+#'
+#' The MHRM algorithm often is more stable than the EM counterpart in \code{mirt} and
+#' congergence of the algorithm should be taken with caution. When the number of iterations reaches
+#' a large number (e.g., greater than 1500) or when \code{Max Change = .2500} values are repeatedly printed
+#' to the console too often (indicating that the parameters were being constrained since they are naturally
+#' moving in steps greater than 0.25, and hence unstable) then the model may either be ill defined or have a
+#' very flat likelihood surface, and genuine maximum likelihood estimates may be difficult to find.
+#'
#' @section Confirmatory IRT:
#'
#' Specification of the confirmatory item factor analysis model follows many of
@@ -38,8 +47,7 @@
#'
#' Specifying a number as the second input to confmirt an exploratory IRT model is estimated and
#' can be viewed as a stochastic analogue of \code{mirt}, with much of the same behaviour and
-#' specifications.
-#' Rotation and target matrix options will be used in this subroutine and will be
+#' specifications. Rotation and target matrix options will be used in this subroutine and will be
#' passed to the returned object for use in generic functions such as \code{summary()} and
#' \code{fscores}. Again, factor means and variances are fixed to ensure proper identification. See
#' \code{\link{mirt}} for more details.
View
14 man/confmirt.Rd
@@ -227,6 +227,20 @@
the \code{anova} function, where a Chi-squared difference
test and AIC/BIC difference values are displayed.
}
+\section{Congergence}{
+ The MHRM algorithm often is more stable than the EM
+ counterpart in \code{mirt} and congergence of the
+ algorithm should be taken with caution. When the number
+ of iterations reaches a large number (e.g., greater than
+ 1500) or when \code{Max Change = .2500} values are
+ repeatedly printed to the console too often (indicating
+ that the parameters were being constrained since they are
+ naturally moving in steps greater than 0.25, and hence
+ unstable) then the model may either be ill defined or
+ have a very flat likelihood surface, and genuine maximum
+ likelihood estimates may be difficult to find.
+}
+
\section{Confirmatory IRT}{
Specification of the confirmatory item factor analysis
model follows many of the rules in the SEM framework for
View
5 src/traceLinePts.cpp
@@ -34,8 +34,11 @@ RcppExport SEXP traceLinePts(SEXP Ra, SEXP Rd, SEXP Rg, SEXP Ru, SEXP RTheta, SE
}
z(j) += d(0) * D(0);
}
- for (i = 0; i < nquad; i++)
+ for (i = 0; i < nquad; i++){
P(i) = g(0) + (u(0) - g(0)) * (exp(z(i))/(1 + exp(z(i))));
+ if(P(i) < 1e-8) P(i) = 1e-8;
+ if((1.0 - P(i)) < 1e-8) P(i) = 1.0 - 1e-8;
+ }
return(P);
Please sign in to comment.
Something went wrong with that request. Please try again.