R/multord_fit.R

In MultOrd: Multivariate Ordinal Regression Models

multord.fit <- function(rho){
  if(!all(sapply(rho$y, is.ordered))) stop("Responses need to be ordered factors", call.=FALSE)
  # split y according to missingness pattern
  rho$y.NA.ind <- split.NA.pattern(rho$y)
  if(is.null(rho$PL.lag)) rho$PL.lag <- rho$ndim

  ## number of thresholds per outcome
  rho$ntheta <- sapply(1:rho$ndim,function(j) nlevels(rho$y[, j]) -1)#rho$ncat - 1

  rho$threshold.values <- if (is.null(rho$threshold.values)) {
                              if ((rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")) && (rho$intercept.type == "fixed")) {
                              lapply(1:rho$ndim, function(j) rep(NA,rho$ntheta[j]))
                          } else if((rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")) && (rho$intercept.type == "flexible")){
                              lapply(1:rho$ndim, function(j) if(rho$ntheta[j] == 1) 0 else c(0,rep(NA,max(rho$ntheta[j]-1,0))))
                          } else if ((rho$error.structure$type %in% c("covGeneral")) && (rho$intercept.type == "flexible")) {
                              lapply(1:rho$ndim, function(j) if(rho$ntheta[j] == 1) 0 else c(0,1,rep(NA,max(rho$ntheta[j]-2,0))))
                          } else if ((rho$error.structure$type %in% c("covGeneral")) && (rho$intercept.type == "fixed")) {
                              lapply(1:rho$ndim, function(j) if(rho$ntheta[j] == 1) 0 else c(0,rep(NA,max(rho$ntheta[j]-1,0))))
                          }#else if(rho$model == "CMORcov2"){
                              #lapply(1:rho$ndim, function(j) if(rho$ntheta[j] == 1) 0 else c(0,rep(NA,max(rho$ntheta[j]-2,0),1)))
                            #}
                        } else rho$threshold.values

  #values of fixed (non-NA) threshold parameters
  rho$threshold.values.fixed <- lapply(1:rho$ndim, function(j){ind <- !is.na(rho$threshold.values[[j]])
                                                                  rho$threshold.values[[j]][ind]})
  ## number of non-fixed thresholds per rater
  rho$npar.theta <- sapply(1:rho$ndim, function(j) sum(is.na(rho$threshold.values[[j]])))

  #checks if binary outcome is present
  rho$binary <- any(sapply(1:rho$ndim, function(j) length(rho$threshold.values[[j]]) == 1))

  rho$threshold <- set.threshold.type(rho)

  #number of columns in the covariate matrix
  rho$p <- ncol(rho$x[[1]])

  # if(is.null(rho$coef.constraints)) rho$coef.constraints <- matrix(1:rho$ndim,ncol = rho$p, nrow = rho$ndim)
  # if(is.null(rho$coef.values)) rho$coef.values <- matrix(NA,ncol = rho$p, nrow = rho$ndim)

  rho$ind.coef <- getInd.coef(rho$coef.constraints, rho$coef.values)

  if(is.null(rho$threshold.constraints)) rho$threshold.constraints <- 1:rho$ndim

  rho$n <- nrow(rho$y)

  #number of total coefficients
  rho$npar.betas <- max(rho$ind.coef, na.rm = TRUE)#sum(unique(c(rho$ind.coef)) != 0)#sum(rho$p)

  ##INCLUDE CHECKS here
  checkArgs(rho)
##############################################################################################
  rho$getInd.thresholds <- switch(rho$threshold,
                                  flexible = getInd.thresholds.flexible,
                                  fix1first = getInd.thresholds.fix1,
                                  fix2first = getInd.thresholds.fix2,
                                  fix2firstlast = getInd.thresholds.fix2,
                                  fixall = getInd.thresholds.fixall)
  rho$ind.thresholds <- rho$getInd.thresholds(rho$threshold.constraints,rho)
  rho$npar.theta.opt <- rho$npar.theta
  rho$npar.theta.opt[duplicated(rho$threshold.constraints)] <- 0
  #number of flexible threshold parameters (in optimizer)
  rho$npar.thetas <- sum(rho$npar.theta.opt)#max(unlist(rho$ind.thresholds))

  #first indices of coefficients in parameter vector for each rater
  rho$first.ind.beta <- rho$npar.thetas + getInd.coef(rho$coef.constraints, rho$coef.values)[, 1]
  #rho$first.ind.theta <- rep(0, rho$ndim)
  #rho$first.ind.theta[rho$npar.theta > 0] <- sapply(rho$ind.thresholds, "[[", 1)
  rho$first.ind.theta <- sapply(rho$ind.thresholds, "[", 1)


  # degrees of freedom for the t-distribution.
  rho$df.t <- switch(rho$link,
                     probit = Inf,
                     logit  = as.integer(8))
  # variance parameter (set to 1)
  rho$sd.y <- switch(rho$link,
                     probit = 1,
                     logit =  pi/sqrt(3) * sqrt((rho$df.t - 2)/rho$df.t))
  rho$inf.value <- 1000

  ###############################
  # error.structure
  ###############################
  #number of correlation parameters for a matrix
  rho$npar.cor <- switch(rho$error.structure$type,
                              corGeneral = (rho$ndim) * (rho$ndim - 1)/2,
                              covGeneral = (rho$ndim) * (rho$ndim - 1)/2,
                              corEqui = ncol(rho$error.structure$x),
                              corAR1 = ncol(rho$error.structure$x))# = 1)
  #number of sigmas
  rho$ncor.levels <- switch(rho$error.structure$type,
                              corGeneral = if(rho$error.structure$formula == ~1)
                                              1 else nlevels(rho$error.structure$x),
                              covGeneral = if(rho$error.structure$formula == ~1)
                                              1 else nlevels(rho$error.structure$x),
                              corEqui    = 1,
                              corAR1     = 1)

  if (rho$error.structure$type == "covGeneral") {
    rho$npar.cor.sd <- rho$ndim
    rho$start.lvar <- rep(0, rho$npar.cor.sd)
  } else {
    rho$npar.cor.sd <- 0
    rho$start.lvar <- NULL
  }
  #number of parameters for all sigmas
  rho$npar.sigmas <- rho$ncor.levels * (rho$npar.cor + rho$npar.cor.sd)

  rho$error.structure$levels <- switch(rho$error.structure$type,
                                       corGeneral = levels(rho$error.structure$x),
                                       covGeneral = levels(rho$error.structure$x),
                                       corEqui = NULL,
                                       corAR1 = NULL)#levels(rho$error.structure$x))

  if (is.null(rho$start.values)) {
    rho$start <- c(getStart.values(rho),
                   rep(0,rho$npar.betas),
                   rep(0,rho$npar.sigmas))
  }

  rho$transf.thresholds <- switch(rho$threshold,
                                  flexible      = transf.thresholds.flexible,
                                  fix1first     = transf.thresholds.fix1.first,
                                  fix2first     = transf.thresholds.fix2.first,
                                  fix2firstlast = transf.thresholds.fix2.firstlast,
                                  fixall        = transf.thresholds.fixall)

  rho$transf.sigmas <- switch(rho$error.structure$type,
                              corGeneral = transf.sigmas.spheric,
                              covGeneral = transf.sigmas.spheric,
                              corEqui    = transf.sigmas.corEqui,
                              corAR1     = transf.sigmas.corAR1)

  rho$transf.par <- switch(rho$error.structure$type,
                           corGeneral = transf.par.cor,
                           corEqui    = transf.par.cor,
                           covGeneral = transf.par.cov,
                           corAR1     = transf.par.cor)

  rho$ObjFun <- switch(rho$link,
                       probit =  PL.probit,
                       logit  =  PL.logit)
  rho$optRes <- suppressWarnings(optimx(rho$start, function(x)  rho$ObjFun(x, rho),
                                method = rho$solver,
                                hessian = FALSE,
                                #itnmax = 10 * length(rho$start)^2,
                                #itnmax = 5,
                                control = list(maxit=200000, trace = 1, kkt = FALSE)))
  if (rho$optRes["convcode"] != 0){
    stop("NO/FALSE CONVERGENCE - choose a different optimizer")
  }
  if (rho$optRes["fevals"] >= 200000){
    warning("reached function evalution limit")
  }
  rho$optpar <- unlist(rho$optRes[1:length(rho$start)])
  rho$objective <- unlist(rho$optRes["value"])
  if (rho$se) {
    rho <- PL_se(rho)
  }

  res <- list()
  res <- multord.finalize(rho)
  res$rho <- rho

  class(res) <- "multord"

  return(res)
}
##########################################################
######                      PL             ###############
##########################################################
# k<-1
# h<-1
# par <- rho$start
PL.probit <- function(par, rho){
  tmp <- rho$transf.par(par, rho)
  pred.upper <- tmp$U
  pred.lower <- tmp$L
  sigmas <- tmp$sigmas
  vecPL <- sapply(1:length(rho$y.NA.ind), function(k) {
    q <- as.numeric(strsplit(names(rho$y.NA.ind[k]), "")[[1]])# turn "101"to 1 0 1)
    ind <- rho$y.NA.ind[[k]]
    if (rho$error.structure$type %in% c("corGeneral", "covGeneral")) { #,"corAR1"
      lev <- match(rho$error.structure$x[ind], rho$error.structure$levels)
    }
    Uq <-  matrix(pred.upper[ind, ], nrow = length(ind))
    Lq <-  matrix(pred.lower[ind, ], nrow = length(ind))
    if (sum(q) == 1){
      pr <- pnorm(Uq[, q == 1]) - pnorm(Lq[, q == 1])
      pr[pr < .Machine$double.eps] <- .Machine$double.eps
      sum(rho$weights[ind] * log(pr))
    } else {
      combis <- combn(which(q == 1), 2)
      #combis <- combis[,which((combis[2,] - combis[1,])  <= 2* rho$PL.lag)]
      combis <- combis[,which((combis[2,] - combis[1,])  <= 2* rho$PL.lag), drop = FALSE]
      #?if rho$PL.lag = 0 ?not allowed
      lprk <- 0
      for (h in seq_len(ncol(combis))){#1:ncol(combis)){
        if(rho$error.structure$type %in% c("corGeneral", "covGeneral")){ #, "corAR1"
          r <- sapply(sigmas[lev],'[', combis[1, h], combis[2, h])
        } else if(rho$error.structure$type %in% c("corAR1")){
          r <- sapply(sigmas[ind],'[', combis[1, h], combis[2, h])
        } else r <- sigmas[ind] ###?ind
        prh <- rectbiv.norm.prob(Uq[, combis[, h]], Lq[, combis[, h]], r)
        prh[prh < .Machine$double.eps] <- .Machine$double.eps
        lprk <- lprk + sum(rho$weights[ind] * log(prh))
      }
      lprk
    }
  })
  -sum(vecPL)
}

PL.logit <- function(par, rho){
  tmp <- rho$transf.par(par, rho)
  pred.upper <- tmp$U
  pred.lower <- tmp$L
  sigmas <- tmp$sigmas
  vecPL <- sapply(1:length(rho$y.NA.ind), function(k){
    q <- as.numeric(strsplit(names(rho$y.NA.ind[k]), "")[[1]])# turn "1_0_1"to 1 0 1)
    ind <- rho$y.NA.ind[[k]]
    if (rho$error.structure$type %in% c("corGeneral", "covGeneral")) { #,"corAR1"
      lev <- match(rho$error.structure$x[ind], rho$error.structure$levels)
    }
    U <-  matrix(pred.upper[ind, ], nrow = length(ind))
    L <-  matrix(pred.lower[ind, ], nrow = length(ind))
    if (is.null(dim(U)[2]) & is.null(dim(L)[2])) {
      dim(U) <- dim(L) <- c(1, length(U))
    }
    if (sum(q) == 1){
      #combis <- which(q == 1)
      pr <- pt(U[, q == 1], df =  rho$df.t) - pt(L[, q == 1], df = rho$df.t)
      pr[pr < .Machine$double.eps] <- .Machine$double.eps
      sum(rho$weights[ind] * log(pr))
    } else {
      combis <- combn(which(q == 1), 2)
      #combis <- combis[,which((combis[2,] - combis[1,])  <= 2* rho$PL.lag)]
      combis <- combis[,which((combis[2,] - combis[1,])  <= 2 * rho$PL.lag), drop = FALSE]
      #?if rho$PL.lag = 0 ?not allowed
      lprk <- 0
      for (h in seq_len(ncol(combis))){#1:ncol(combis)){
        prh <- NULL
        if(rho$error.structure$type %in% c("corGeneral", "covGeneral")){ #, "corAR1"
          r <- sapply(sigmas[lev],'[', combis[1, h], combis[2, h])
        } else if(rho$error.structure$type %in% c("corAR1")){
          r <- sapply(sigmas[ind],'[', combis[1, h], combis[2, h])
        } else r <- sigmas[ind] ###?ind
        for (i in 1:length(ind)) {
        prh[i] <- biv.nt.prob2(df=rho$df.t,
                               lower = L[i, combis[,h]],
                               upper = U[i, combis[,h]],
                               r = r[i])
        }
        prh[prh < .Machine$double.eps] <- .Machine$double.eps
        lprk <- lprk + sum(rho$weights[ind] * log(prh))
      }
      lprk
    }
  })
  -sum(vecPL)
}

.onLoad <- function(library, pkg)
{
  library.dynam("MultOrd", pkg, library)
  invisible()
}