R/update_probit.R

In bomarkussen/probit: Item response analysis with random effects.

#' update method for probit-class
#'
#' @description
#' update method for \code{probit}-object as generated by \code{\link{probit}}.
#'
#'
#' @param object Object of class \code{probit}.
#' @param fixed Update of right hand side of model formula for fixed effects. Defaults to \code{fixed=NULL}, which corresponds to no update.
#' @param random List of updates of right hand side of formulas for the random effects. Defaults to \code{random=NULL}, which corresponds to no update.
#' @param A Transformation matrix of predicted random effects. Defaults to \code{A=NULL}, which corresponds to identity transformation of residing random effects.
#' @param dependence Text string (\code{"marginal" or "joint"}) deciding whether random effects are assumed independent or with a common joint normal distribution. Defaults to \code{dependence=NULL}, which corresponds to no update.
#' @param data Date frame with data on the wide format. Defaults to \code{data=NULL}, which corresponds to \code{data} taken from the call object.
#' @param B Number of simulations in minimization step. Defaults to \code{B=NULL}, which corresponds to \code{B} taken from the call object.
#' @param BB Number of simulations per subject in \code{anova} or \code{update}. Defaults to \code{BB=NULL}, which corresponds to \code{BB} taken from the call object.
#' @param maxit Maximal number of minimization-maximization steps. Defaults to \code{maxit=20}.
#' @param sig.level Significance level at which the iterative stochastic optimization will be stopped. Defaults to \code{sig.level=0.60}.
#' @param verbose Numeric controlling amount of convergence diagnostics. Default: \code{verbose=0} corresponding to no output.
#' @param estimate.models Boolean deciding if model parameters are estimated in \code{update}.
#'
#' @return \code{\link{probit-class}} object.
#'
#' @export
update.probit <- function(object,fixed=NULL,random=NULL,A=NULL,dependence=NULL,data=NULL,B=NULL,BB=NULL,maxit=20,sig.level=0.6,verbose=0,estimate.models=TRUE) {
  # update fixed effects?
  if (is.null(fixed)) {
    fixed <- formula(object$fixed)
    m.fixed <- object$m.fixed
  } else {
    fixed <- update(object$fixed,fixed)
    m.fixed <- NULL
    if (!estimate.models) warning("Fixed effects model had been updated, but it is not reestimated")
  }

  # update random effect models
  if (is.null(random)) {
    random <- object$random
    m.random <- object$m.random
  } else {
    new_random <- object$random
    random_ii  <- match(sapply(random,function(x){all.vars(x[[2]])}),
                        sapply(object$random,function(x){all.vars(x[[2]])}))
    for (i in 1:length(random_ii)) {
      if (is.na(random_ii[i])) {
        new_random <- c(new_random,random[[i]])
      } else {
        new_random[[random_ii[i]]] <- update(new_random[[random_ii[i]]],random[[i]])
      }
    }
    random <- new_random
    m.random <- NULL
    if (!estimate.models) warning("Random effects models have been updated, but they are not reestimated")
  }

  # remove non-used random effects
  random <- random[is.element(sapply(random,function(x){all.vars(x[[2]])}),
                              all.vars(delete.response(terms(fixed))))]
  q <- length(random)
  if (q<1) stop("Present implementation assumes at least one random effect")

  # update values of dependence, B, BB, data?
  if (is.null(dependence)) {dependence <- object$dependence}
  if (is.null(B))          {B <- object$B}
  if (is.null(BB))         {BB <- object$BB}
  if (is.null(data)) {data <- object$data} else {data <- as_tibble(data)}

  # fix dependence
  if (dependence!="marginal") dependence <- "joint"

  # recode ordinal variables as numeric
  for (i in object$items.ordinal) {
    data[[i]] <- as.numeric(data[[i]])
  }

  # find random effects
  new_random.eff <- sapply(random,function(x){all.vars(x[[2]])})
  old_random.eff <- sapply(object$random,function(x){all.vars(x[[2]])})

  # remove observations with missing explanatory variables
  i <- complete.cases(select(data,any_of(setdiff(
    c(all.vars(fixed),unlist(sapply(random,function(x){all.vars(x)}))),
    c(object$items.interval,object$items.ordinal,new_random.eff)))))
  if (sum(!i)>0) {
    data <- data[i,]
    warning("Remove ",sum(!i)," observations with non-complete explanatory variables. NOTE: This may interact badly with update().")
  }

  # transformation matrix
  if (is.null(A)) {
    A <- matrix(0,length(new_random.eff),length(old_random.eff))
    A[outer(new_random.eff,old_random.eff,"==")] <- 1
  } else {
    # check dimensions of A matrix
    if (nrow(A)!=length(new_random.eff)) stop("Number of rows in A must match number of new random effects")
    if (ncol(A)!=length(old_random.eff)) stop("Number of columns in A must match number of old random effects")
  }

  # find subjects
  new_subjects <- unique(data[[object$subject.name]])
  old_subjects <- unique(object$data[[object$subject.name]])

  # initiate mu at model predictions
  old_mu <- matrix(0,length(new_subjects),ncol(object$mu))
  data.short <- data %>% group_by(!!as.name(object$subject.name)) %>% slice_head(n=1)
  for (i in 1:ncol(object$mu)) {
    old_mu[,i] <- predict(object$m.random[[i]],newdata=data.short)
  }
#  # initiate psi at the mean (doesn't necessarily make any sense)
#  old_psi <- matrix(colMeans(object$psi),length(new_subjects),ncol(object$psi),byrow=TRUE)
  # initiate psi at Gamma
  old_psi <- matrix(object$Gamma[upper.tri(object$Gamma,diag=TRUE)],length(new_subjects),ncol(object$psi),byrow=TRUE)
  # reset (mu,psi) for known subjects at their present value
  ii <- is.element(new_subjects,old_subjects)
  old_mu[ii,]  <- object$mu[match(new_subjects,old_subjects)[ii]]
  old_psi[ii,] <- object$psi[match(new_subjects,old_subjects)[ii]]

  # update mu and psi
  mu  <- old_mu%*%t(A)
  psi <- matrix(0,length(new_subjects),q*(q+1)/2)
  for (s in 1:length(new_subjects)) {
    tmp <- matrix(0,length(old_random.eff),length(old_random.eff))
    tmp[upper.tri(tmp,diag = TRUE)] <- old_psi[s,]
    tmp <- chol(solve(A%*%solve(t(tmp)%*%tmp)%*%t(A) + diag(1e-4,nrow=nrow(A))))
    psi[s,] <- tmp[upper.tri(tmp,diag = TRUE)]
  }

  # update Gamma
  Gamma <- chol(solve(A%*%solve(t(object$Gamma)%*%object$Gamma)%*%t(A) +
                        diag(1e-4,nrow=nrow(A))))

#  # estimate random effects models
#  # TO DO: really update m.random like this??
#  if (estimate.models) {
#    m.random <- vector("list",length(new_random.eff))
#    names(m.random) <- new_random.eff
#    data.short <- data %>% group_by(!!as.name(object$subject.name)) %>% slice_head(n=1)
#    U <- as.data.frame(cbind(new_subjects,mu))
#    names(U) <- c(object$subject.name,new_random.eff)
#    data.short <- full_join(data.short,U,by=object$subject.name)
#    for (i in 1:q) m.random[[i]] <- lm(random[[i]],data=data.short)
#  }

  # refit and return
  return(probit:::MM_probit(maxit=maxit,sig.level=sig.level,verbose=verbose,
                   fixed=fixed,response.name=object$response.name,weight.name=object$weight.name,
                   item.name=object$item.name,items.interval=object$items.interval,items.ordinal=object$items.ordinal,
                   ordinal.levels=object$ordinal.levels,
                   subject.name=object$subject.name,random=random,dependence=dependence,
                   m.fixed=m.fixed,sigma2=object$sigma2,eta=object$eta,m.random=m.random,Gamma=Gamma,
                   mu=mu,psi=psi,
                   B=B,BB=BB,
                   data=data,
                   estimate.models=estimate.models))
}