R/utilities.R

## TODO if no intercept in the covGeneral model set default to fix1first
set.threshold.type <- function(rho){
#fixall
  if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))) {#all thresholds are fixed in all dimensions
      if (rho$intercept == FALSE) cat("We suggest to include an intercept in the model (formula = y ~ 1 + ...)")
      type <- "fixall"
#fix2first
  } else if (all(sapply(1:rho$ndim, function(j){
      (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,2)) #all first two thresholds are not NA
      ) || ((length(rho$threshold.values[[j]]) == 1) &&  (which(!is.na(rho$threshold.values[[j]])) == 1))
      }))){
      if (rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")) cat("We suggest to fix only one threshold or the intercept in a correlation model.\n")
      if ((rho$error.structure$type %in% c("covGeneral"))&& (rho$intercept.type == "fixed")) cat("We suggest to fix either two thresholds or one threshold and the intercept in a covGeneral model.\n")
      type <- "fix2first"
#fix2firstlast
  } else if (all(sapply(1:rho$ndim, function(j){
      (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j]))#all first and last two thresholds are not NA
      ) || ((length(rho$threshold.values[[j]]) == 1) &&  (which(!is.na(rho$threshold.values[[j]])) == 1))
      }))){
      if (rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")) cat("We suggest to fix only one threshold or the intercept in a correlation model.\n")
      if ((rho$error.structure$type %in% c("covGeneral"))&& (rho$intercept.type == "fixed")) cat("We suggest to fix either two thresholds or one threshold and the intercept in a covGeneral model.\n")
      type <- "fix2firstlast"
#fix1first
  } else if (all(sapply(1:rho$ndim, function(j) (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && all(which(!is.na(rho$threshold.values[[j]])) == 1))))){ #all first thresholds are not NA (and no additional)
      if ((rho$error.structure$type == "covGeneral") && (rho$intercept.type == "flexible")) stop("Model with covGeneral is not identifiable.
                          Please either fix two thresholds or one threshold and the intercept.\n", call. = FALSE)
      if ((rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui"))&& (rho$intercept.type == "fixed")) cat("We suggest to fix only one threshold or the intercept in a correlation model.\n")
      type <- "fix1first"
#flexible
  } else if (all(sapply(1:rho$ndim, function(j) all(is.na(rho$threshold.values[[j]]))))){#all thresholds NA
      if (rho$error.structure$type == "covGeneral") stop("Model with covGeneral is not identifiable.
                                                        Please either fix two thresholds or one threshold and the intercept.\n", call. = FALSE)
      if ((rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")) && (rho$intercept.type == "flexible")) stop("Model is not identifiable.
                                                                                                                                  Please either fix one threshold or the intercept.", call. = FALSE)
      type <- "flexible"
#ERRORS
  } else stop("Either fix all thresholds in one or more outcome dimensions,
              or consistently in all other outcome dimensions, all first thresholds or none.\n", call. = FALSE)
  if((rho$error.structure$type == "covGeneral") && (rho$binary == TRUE)){
    if(type != "fix1first") stop("In the presence of binary outcomes intercept and one threshold
                                  have to be fixed to some value.\n", call. = FALSE)
  }
  type
}


# #CORRELATION
#   } else if (rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")){
#    #NO INTERCEPT
#         if (rho$intercept.type == "fixed"){
#             if (all(sapply(1:rho$ndim, function(j) all(is.na(rho$threshold.values[[j]]))))){#all thresholds NA
#                 type <- "flexible"
#           } else stop("If you want fix a threshold value please use intercept in formula (remove ~ 0 + ...)", call. = FALSE)
#   #TODO? other parameterization are not allowed for intercept == FALSE?
#
#   #INTERCEPT
#     } else if (rho$intercept == TRUE){
#         if (all(sapply(1:rho$ndim, function(j) all(is.na(rho$threshold.values[[j]]))))){ #all thresholds NA
#             if (all(!is.na(rho$coef.values[,1]))){ #intercepts are fixed to constants
#                 type <- "flexible"
#           } else stop("Remove intercept in formula (~ 0 + ...) or fix the first threshold in threshold.values", call. = FALSE)
#         }
#         if (all(sapply(1:rho$ndim, function(j) (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && all(which(!is.na(rho$threshold.values[[j]])) == 1))))){ #all first thresholds are not NA (and no additional)
#             type = "fix1first"
#       } else if (any(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in at least one dimension
#             if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))  ||  all(
#               (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && (which(!is.na(rho$threshold.values[[j]])) == 1)))))){ #all not NA or only first fixed for each j
#                   type <- "fix1first"
#           } else if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))  ||  all(is.na(rho$threshold.values[[j]]))))){ #all not NA or all NA for each j
#                   type <- "flexible"
#           } else stop("Either fix all thresholds in one or more outcome dimensions,
#                            or consistently in all other outcome dimensions, all first thresholds or none.
#                             In addition, intercept must be included in the model.", call. = FALSE)
#     } else{
#       stop("Model is not identifiable. One of the following could help:
#            1. Remove intercept in formula (~ 0 + ...) or
#            2. Fix all first thresholds in threshold.values or
#            3. Fix all thresholds in threshold.values or
#            4. Either fix all thresholds in one or more outcome dimensions,
#            or consistently in all other outcome dimensions, all first thresholds or none", call. = FALSE)
#     }
#     }
#
# #covGeneral
#   } else if (rho$error.structure$type == "covGeneral"){
#     if (rho$intercept == FALSE) stop("formula needs intercept if you want to fix a threshold in model CMORcov", call. = FALSE)#TODO?
#     if (all(sapply(1:rho$ndim, function(j){
#       (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,2)) #all first two thresholds are not NA
#       ) || ((length(rho$threshold.values[[j]]) == 1) &&  (which(!is.na(rho$threshold.values[[j]])) == 1))
#     }))){
#       type == "fix2first"
#     } else if (all(sapply(1:rho$ndim, function(j){
#       (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j]))#all first and last two thresholds are not NA
#       ) || ((length(rho$threshold.values[[j]]) == 1) &&  (which(!is.na(rho$threshold.values[[j]])) == 1))
#     }))){
#       type == "fix2firstlast"
#     }else if (any(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in at least one dimensions
#       ##CHECK somewhere if all values are fixed in one dimensions + ?warning(sense)
#       if (all(sapply(1:rho$ndim, function(j) #
#         #(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1))) &&  all((which(!is.na(rho$threshold.values[[j]])) == 1))) ||
#         (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) &&  all((which(!is.na(rho$threshold.values[[j]])) == c(1,2)))) ||
#         (all(!is.na(rho$threshold.values[[j]])))))){
#         type <- "fix2first"
#       } else if (all(sapply(1:rho$ndim, function(j)
#         #(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1))) &&  all((which(!is.na(rho$threshold.values[[j]])) == 1))) ||
#         (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) &&  all((which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j])))) ||
#         (all(!is.na(rho$threshold.values[[j]])))))){
#         type <- "fix2firstlast"
#       } else stop("Either fix all thresholds in one or more outcome dimensions,
#                     or consistently in all other outcome dimensions, all first two thresholds or all first and last thresholds", call. = FALSE)
#     } else stop("Either fix first two or first and last threshold in threshold.values if you use a model with covGeneral", call. = FALSE)
#   }
#   type
# }

#data.x <- data.multord$x
#ndim <- rho$ndim
set.error.structure <- function(error.structure, data.x, ndim){
  if (error.structure$formula == ~1 ) {
    if (error.structure$type %in% c("corEqui", "corAR1")) {
     error.structure$x <- model.matrix(error.structure$formula, data.x[[1]])
   } else {
     error.structure$x <- as.factor(model.matrix(error.structure$formula, data.x[[1]]))
   }
  } else if (error.structure$type %in% c("corEqui", "corAR1")){

    tmp <- lapply(1:ndim, function(j) {
      tmp1 <- model.matrix(error.structure$formula, data.x[[j]])
      attribute <- attr(tmp1, "assign")
      tmp2 <- tmp1[match(rownames(data.x[[j]]),rownames(tmp1)),]
      rownames(tmp2) <- rownames(data.x[[j]])
      attr(tmp2, "assign") <- attribute
      tmp2
    }
    )
    error.x <- sapply(1:ncol(tmp[[1]]), function(k){
      xtcol <- do.call(cbind,lapply(tmp, `[`,,k))
      xtcol_final <- apply(xtcol,1,function(i) unique(i[!is.na(i)]))
      if(is.list(xtcol_final)) stop("Covariate dependent variables need to be constant across multiple observations", call.=FALSE)
      xtcol_final
    })

    colnames(error.x) <- colnames(tmp[[1]])
    attr(error.x, "assign") <- attr(tmp[[1]], "assign")

    error.structure$x <- error.x
  # } else if (error.structure$type == "corAR1"){
  #   #depends on one factor (factor of first year)
  #   if (!is.factor(data.x[[1]][,all.vars(error.structure$formula)])) stop("variable must be of type factor in corAR1", call. = FALSE)
  #   tmp <- lapply(1:ndim,function(j) data.x[[j]][,all.vars(error.structure$formula)])
  #   if (!all(sapply(1:ndim, function(j) all(tmp[[1]]==tmp[[j]], na.rm=T)))) stop("Covariate dependent variables need to be constant across multiple observations", call. = FALSE)
  #   error.structure$x <- tmp[[1]]#data.x[[1]][,all.vars(error.structure$formula)]
  } else if (error.structure$type %in%  c("corGeneral", "covGeneral")){
    #depends on one factor
    if (length(all.vars(error.structure$formula[[2]])) > 1) stop("only one factor is allowed in covGeneral and corGeneral", call. = FALSE)
    if (!is.factor(data.x[[1]][,all.vars(error.structure$formula)])) stop("variable must be of type factor in covGeneral and corGeneral", call. = FALSE)
    tmp <- lapply(1:ndim,function(j) data.x[[j]][,all.vars(error.structure$formula)])

    if (!all(sapply(1:ndim, function(j) all(tmp[[1]]==tmp[[j]], na.rm=T)))) stop("Covariate dependent variables need to be constant across multiple observations", call. = FALSE)
    error.structure$x <- as.factor(apply(do.call("cbind.data.frame", tmp), 1, function(x) unique(x[!is.na(x)])))
     } else {
    #depends on one factor
    if (length(all.vars(error.structure$formula[[2]])) > 1) stop("only one factor is allowed in covGeneral and corGeneral", call. = FALSE)
    if (!is.factor(data[,all.vars(error.structure$formula)])) stop("variable must be of type factor in covGeneral and corGeneral", call. = FALSE)
    error.structure$x <- data.x[,all.vars(error.structure$formula)]
  }
  error.structure
}


set.error.structure.CMOR <- function(error.structure, data){
  if (error.structure$type == "corEqui"){
    error.structure$x <- model.matrix(error.structure$formula, data)
  } else if (error.structure$type == "corAR1"){
    stop("corAR1 is not applicable in CMOR", call. = FALSE)
    # if (error.structure$formula == ~1){
    #   error.structure$x <- as.factor(model.matrix(error.structure$formula, data[[1]]))
    # } else{ #depends on one factor (factor of first year)
    #   if (!is.factor(data[[1]][,all.vars(error.structure$formula)])) stop("variable must be of type factor in covGeneral and corGeneral", call. = FALSE)
    #   error.structure$x <- data[[1]][,all.vars(error.structure$formula)]
    # }
  } else if (error.structure$type %in%  c("corGeneral", "covGeneral")){
    if (error.structure$formula == ~ 1){
      error.structure$x <- as.factor(model.matrix(error.structure$formula, data))
    } else{ #depends on one factor
      if (length(all.vars(error.structure$formula[[2]])) > 1) stop("only one factor is allowed in covGeneral and corGeneral", call. = FALSE)
      if (!is.factor(data[,all.vars(error.structure$formula)])) stop("variable must be of type factor in covGeneral and corGeneral", call. = FALSE)
      error.structure$x <- data[,all.vars(error.structure$formula)]
    }
  }
  error.structure
}

checkArgs <- function(rho){
  #CHECK if treshold.values is in line with threshold.constraints
  if (length(rho$threshold.constraints) != rho$ndim) stop("dimensions of threshold.values and number of thresholds do not match", call. = FALSE)
    if (any(sapply(1:rho$ndim, function(j) length(rho$threshold.values[[j]]) != rho$ntheta[j])))
      stop("dimensions of threshold.values and number of thresholds do not match", call. = FALSE)
  for (j in unique(rho$threshold.constraints)){
    ind <- which(rho$threshold.constraints == j)
    if (length(unique(rho$threshold.values[ind]))!=1) stop("If constraints are set on thresholds (by threshold.constraints), threshold.values need to be specified accordingly for these outcome dimensions.
                                                           Maybe dimension do not have the same number of threshold parameters.", call. = FALSE)
  }
  if (nrow(rho$coef.constraints) != rho$ndim) stop("row dimension of coef.constraints and outcome dimension do not match (?factor)", call. = FALSE)

  for (j in 1:ncol(rho$coef.constraints)){
    indj <- unique(rho$coef.constraints[,j])
    indj <- indj[!is.na(indj)]
    lapply(seq_len(length(indj)), function(k) {
      tmpind <- which(rho$coef.constraints[,j] == indj[k])
      tmp <- rho$coef.values[tmpind,j]
      if(length(unique(tmp)) != 1) stop("If constraints are set on the coefficients (by coef.constraints),
                                        coef.values need to be specified accordingly for these outcome dimensions.", call. = FALSE)
    })
    }
}

##########################################
###### AUXILIARY FUNCTIONS ##############
###########################################################################
# rectbiv.norm.prob_TEST <- function(U, L, r) {
#   # computes the rectangle probabilities for biv.normal-distribution
#   if (is.null(dim(U)[2]) & is.null(dim(L)[2])) {
#     dim(U) <- dim(L) <- c(1, length(U))
#   }

#   n <- nrow(U)
#   x <- c(U[, 1],L[, 1],U[, 1],L[, 1])
#   y <- c(U[, 2],U[, 2],L[, 2],L[, 2])
#   r

#    p <- pbivnorm(x,y,r)

#    p[is.nan(p)] <- 0

#    p[seq_len(n)]-p[n + seq_len(n)]-p[2*n + seq_len(n)] + p[3*n + seq_len(n)]
# }

rectbiv.norm.prob <- function(U, L, r) {
  # computes the rectangle probabilities for biv.normal-distribution
  if (is.null(dim(U)[2]) & is.null(dim(L)[2])) {
   dim(U) <- dim(L) <- c(1, length(U))
  }
  p1 <- pbivnorm(U[, 1], U[, 2], r)
  p2 <- pbivnorm(L[, 1], U[, 2], r)
  p3 <- pbivnorm(U[, 1], L[, 2], r)
  p4 <- pbivnorm(L[, 1], L[, 2], r)

  p1[is.nan(p1)] <- 0
  p2[is.nan(p2)] <- 0
  p3[is.nan(p3)] <- 0
  p4[is.nan(p4)] <- 0

  p1-p2-p3+p4
}

biv.nt.prob2 <-function (df, lower, upper,
                         mean=c(0,0), r) {
  if (df == Inf) {
    nu <- 0
  } else {
    nu <- df
  }
  rho <- r
  if (any(lower > upper)) stop("lower > upper integration limits")
  if (any(lower == upper)) return(0)
  infin <- as.integer(c(2, 2))
  prob <- as.double(0)
  a <- .Fortran("smvbvt", prob, nu, lower, upper, infin, rho,
                #PACKAGE = "mnormt")
                PACKAGE = "MultOrd")
  return(a[[1]])
}

casesNA<- function(y) {
  # returns a column vector v, assigning a different number/label to each unique combination of NA pattern
  contrast <- as.matrix((!is.na(y)) + 0)
  contrast.char <- apply(contrast, 1, paste, collapse = "") # make each contrast row as character i.e., 111 has all three response variables
  factor(contrast.char, levels = unique(contrast.char))
}

split.NA.pattern <- function(y) {
  # y should be a data.frame
  # return a list, whose elements are groups of the matrix y according to each unique combination of response variables
  if (!is.data.frame(y)) y <- as.data.frame(y)
  yg <- split(y, casesNA(y)) # grouped y by cases of NA pattern
  lapply(yg, function(x)  which(rownames(y) %in% rownames(x)))
}



getStart.values <- function(rho){
  #   # set starting value for beta and the gammas thresholds
  #   #thetas <-  qlogis((1:rho$ntheta) / (rho$ntheta + 1))
  #   #if (rho$link != "logistic")
  #   #    thetas <- thetas/1.7
  #   mclm <- lapply(1:rho$ndim, function(j)
  #     ordinal:::clm(factor(rho$y[,j]) ~   rho$x[,which(rho$coef.constraints[j,] != 0)], link = "logit"))
#   #thetas <-  lapply(mclm, "[[", "alpha")
#   betas <-  lapply(mclm, "[[", "beta") #Reduce("+", lapply(mclm, "[[", "beta"))
#   #gammas <- lapply(thetas, function(x) c(x[1L], log(diff(x))))
gammas <- sapply(1:rho$ndim, function(j) {
  if (rho$npar.theta.opt[j] != 0){
    theta <- if (rho$ntheta[j] >= 2) polr(rho$y[, j] ~1)$zeta else 0#qlogis((1:rho$ntheta[j])/(rho$ntheta[j] + 1))
    if (!grepl("logit", rho$link)) theta <- theta/1.7
    c(theta[1L], log(diff(theta)))[1:rho$npar.theta.opt[j]]
  } else NULL
})
c(unlist(gammas))#, unlist(betas))# rep(, sum(rho$npar.betas))) }
}


transf.par.cor <- function(par, rho) {
  sigmas <- rho$transf.sigmas(par[rho$npar.thetas + rho$npar.betas + seq_len(rho$npar.sigmas)],rho)
  theta <- rho$transf.thresholds(par[seq_len(rho$npar.thetas)], rho)
  par_beta <- par[rho$npar.thetas + seq(rho$npar.betas)]
  beta <- sapply(1:ncol(rho$coef.constraints), function(j){
    sapply(1:nrow(rho$coef.constraints), function(i,j) ifelse(is.na(rho$ind.coef[i,j]), rho$coef.values[i,j], par_beta[rho$ind.coef[i, j]]), j)
  })
  pred.fixed <- sapply(1:rho$ndim, function(j) rho$x[[j]] %*% beta[j, ])
  theta.lower <- sapply(1:rho$ndim, function(j) c(-rho$inf.value, theta[[j]])[rho$y[, j]])
  theta.upper <- sapply(1:rho$ndim, function(j) c(theta[[j]], rho$inf.value)[rho$y[, j]])
  pred.lower <- (theta.lower - pred.fixed)/rho$sd.y
  pred.upper <- (theta.upper - pred.fixed)/rho$sd.y
  list(U = pred.upper,
       L = pred.lower,
       sigmas = sigmas)
}

transf.par.cov <- function(par, rho) {
  exp.par.sd <- exp(par[rho$npar.thetas + rho$npar.betas +
    rho$npar.cor * rho$ncor.levels +
    seq_len(rho$npar.cor.sd * rho$ncor.levels)])
  sigmas <- rho$transf.sigmas(par[rho$npar.thetas + rho$npar.betas + seq_len(rho$npar.sigmas)],rho,
                                exp.par.sd)
  exp.par.sd <- lapply(1:rho$ncor.levels, function(l) exp.par.sd[ (l-1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd) ])
  theta <- rho$transf.thresholds(par[seq_len(rho$npar.thetas)], rho)
  par_beta <- par[rho$npar.thetas + seq_len(rho$npar.betas)]
  beta <- sapply(1:ncol(rho$coef.constraints), function(j){
    sapply(1:nrow(rho$coef.constraints), function(i,j)
      ifelse(is.na(rho$ind.coef[i,j]), rho$coef.values[i,j], par_beta[rho$ind.coef[i, j]]), j)
  })
  pred.fixed <- sapply(1:rho$ndim, function(j) rho$x[[j]] %*% beta[j, ])
  theta.lower <- sapply(1:rho$ndim, function(j) c(-rho$inf.value, theta[[j]])[rho$y[, j]])
  theta.upper <- sapply(1:rho$ndim, function(j) c(theta[[j]], rho$inf.value)[rho$y[, j]])
  ## make matrix of exp.par.sd
  lev <- match(rho$error.structure$x, rho$error.structure$levels)
  #pred.lower <- (theta.lower - pred.fixed)/rho$sd.y
  #pred.upper <- (theta.lower - pred.fixed)/rho$sd.y
  #exp.par.sd.mat <- matrix(ncol = ncol(theta.upper), nrow = nrow(theta.upper))
  #for (l in 1:rho$ncor.levels){
  #  exp.par.sd.mat[lev == l, ] <-
  #  exp.par.sd[(l - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)]
  #}
  #for (i in 1:nrow(pred.lower)){
  #   pred.lower[i,] <- pred.lower[i,]/exp.par.sd[(lev[i] - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)]
  #   pred.upper[i,] <- pred.upper[i,]/exp.par.sd[(lev[i] - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)]
  #}
  pred.lower <- theta.lower - pred.fixed
  pred.upper <- theta.upper - pred.fixed
  for (l in 1:rho$ncor.levels){
  pred.lower[lev==l,] <-t(t(pred.lower[lev==l,]/rho$sd.y)/exp.par.sd[[l]])
  pred.upper[lev==l,] <-t(t((pred.upper[lev==l,])/rho$sd.y)/exp.par.sd[[l]])
  }
  #pred.lower <-t(t((theta.lower - pred.fixed)/rho$sd.y)/exp.par.sd)
  #pred.upper <-t(t((theta.upper - pred.fixed)/rho$sd.y)/exp.par.sd)

  # t(t(pred.upper)/exp.par.sd)
  # test speed up
  sigmas <- lapply(sigmas, cov2cor)#cov2cor(sigmas)
  list(U = pred.upper,
       L = pred.lower,
       sigmas = sigmas,
       std.dev = exp.par.sd)
}
#########################################################################
## transformation of the threshold parameters (to ensure monotonicity) ##
#########################################################################
transf.thresholds.fixall <- function(gamma, rho){
  rho$threshold.values
}


transf.thresholds.flexible <- function(gamma, rho){
  lapply(1:rho$ndim, function(j)
    if (anyNA(rho$threshold.values[[j]])){
    if (rho$ntheta[j] > 1) {
      cumsum(c(gamma[rho$ind.thresholds[[j]][1]],
               exp(gamma[rho$ind.thresholds[[j]][2:rho$ntheta[j]]])))
    } else if (rho$ntheta[j] == 1) gamma[rho$ind.thresholds[[j]]] else NULL
    } else rho$threshold.values[[j]]
  )
}

transf.thresholds.fix1.first <- function(gamma, rho){
  lapply(1:rho$ndim, function(j)
    if (length(rho$ind.thresholds[[j]]) >= 1) {
      cumsum(c(rho$threshold.values.fixed[[j]][1],
        exp(gamma[rho$ind.thresholds[[j]][seq_len(length(rho$ind.thresholds[[j]]))]])))
    } else if (length(rho$ind.thresholds[[j]]) == 0) rho$threshold.values.fixed[[j]] else NULL)
}

#gamma <- par[seq_len(rho$npar.thetas)]

transf.thresholds.fix2.first <- function(gamma, rho){
   lapply(1:rho$ndim, function(j)
     if (rho$npar.theta[j] >= 1) {
       c(rho$threshold.values.fixed[[j]][1],
        cumsum(c(rho$threshold.values.fixed[[j]][2] ,exp(gamma[rho$first.ind.theta[j] - 1 +
                                                     seq_len(rho$npar.theta[j])]))))
     } else {
       if (rho$npar.theta[j] == 0) rho$threshold.values.fixed[[j]]
     }
     )
}

# transf.thresholds.fix2.first <- function(gamma, rho){
#   lapply(1:rho$ndim, function(j)
#     if (rho$npar.theta[rho$threshold.constraints[j]] >= 1) {
#       c(rho$threshold.values.fixed[[j]][1],
#         cumsum(c(rho$threshold.values.fixed[[j]][2] ,exp(gamma[rho$first.ind.theta[rho$threshold.constraints[j]] - 1 +
#                                                                  seq_len(rho$npar.theta[rho$threshold.constraints[j]])]))))
#     } else {
#       if (rho$npar.theta[rho$threshold.constraints[j]] == 0) rho$threshold.values.fixed[[j]]
#     }
#   )
# }

transf.thresholds.fix2.firstlast <- function(gamma, rho){
  lapply(1:rho$ndim, function(j){
    gamma1 <- gamma[rho$first.ind.theta[j] - 1 + seq_len(rho$npar.theta[j])]
    recursive.theta <- function(i) {
      if (i == 0) 0
      else return ((exp(gamma1[i]) +
                      recursive.theta(i - 1))/(1 + exp(gamma1[i])))
    }
    if (rho$npar.theta[j] >= 1) {
      theta <- sapply(1:rho$npar.theta[j], function(i)
        recursive.theta(i))
      c(0, theta, 1) * (rho$threshold.values.fixed[[j]][2] - rho$threshold.values.fixed[[j]][1]) +
        rho$threshold.values.fixed[[j]][1]
    } else {
      rho$threshold.values.fixed[[j]]
    }
  })
}

# transf.thresholds.fix2.firstlast <- function(gamma, rho){
#   lapply(1:rho$ndim, function(j){
#     gamma1 <- gamma[rho$first.ind.theta[j] - 1 + seq_len(rho$npar.theta[rho$threshold.constraints[j]])]
#     recursive.theta <- function(i) {
#       if (i == 0) 0
#       else return ((exp(gamma1[i]) +
#                       recursive.theta(i - 1))/(1 + exp(gamma1[i])))
#     }
#     if (rho$npar.theta[rho$threshold.constraints[j]] >= 1) {
#       theta <- sapply(1:rho$npar.theta[rho$threshold.constraints[j]], function(i)
#         recursive.theta(i))
#       c(0, theta, 1) * (rho$threshold.values.fixed[[j]][2] - rho$threshold.values.fixed[[j]][1]) +
#         rho$threshold.values.fixed[[j]][1]
#     } else {
#       rho$threshold.values.fixed[[j]]
#     }
#   })
# }

##############################################################################
## transformation of the parameters of the correlation/covariance structure ##
##############################################################################

transf.sigmas.spheric <- function(par.nu, rho, exp.par.sd = NULL) {
  lapply(1:rho$ncor.levels, function(l) {
    nu <- par.nu[(l - 1) * rho$npar.cor + seq_len(rho$npar.cor)]
    angles <- pi * exp(nu)/(1 + exp(nu))
    cosmat <- diag(rho$ndim)
    cosmat[lower.tri(cosmat)] <- cos(angles)
    S1 <- matrix(0, nrow = rho$ndim, ncol = rho$ndim)
    S1[, 1L] <- 1
    S1[-1L, -1L][lower.tri(S1[-1L, -1L], diag = T)] <- sin(angles)
    tLmat <- sapply(1:rho$ndim,
                    function(j) cosmat[j, ] * cumprod(S1[j, ]))
    sigma <- crossprod(tLmat)
    if (is.null(exp.par.sd)) {
      return(sigma)
    } else {
     stdev <-  exp.par.sd[(l - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)]
     return(t(stdev * sigma) * stdev)
    }
  })
}

# transf.sigmas.corAR1 <- function(par.sigma, rho, exp.par.sd = NULL) {
#   lapply(1:rho$ncor.levels, function(l){
#   ###TODO include covariate dependent corAR1? (like corEqui)
#   r <- z2r(par.sigma[seq_len(rho$npar.cor) + l -1])
#   #r <- 0.9
#   sigma <- diag(rho$ndim)
#   sigma[lower.tri(sigma)]  <- r^sequence((rho$ndim-1):1)
#   ### CHECK: easier way?
#   sigma <- sigma + t(sigma) - diag(diag(sigma))
#   sigma
#   })
# }

transf.sigmas.corAR1 <- function(par.sigma, rho, exp.par.sd = NULL) {
  w <- par.sigma[seq_len(rho$npar.cor)]
  z <- rho$error.structure$x %*% w
  r <- z2r(z)
  lapply(seq_len(length(r)), function(i){
    sigma <- diag(rho$ndim)
    sigma[lower.tri(sigma)]  <- r[i]^sequence((rho$ndim-1):1)
    ### CHECK: easier way?
    sigma <- sigma + t(sigma) - diag(diag(sigma))
    sigma
  })
}

transf.sigmas.corEqui <- function(par.sigma, rho, exp.par.sd = NULL) {
  w <- par.sigma[seq_len(rho$npar.cor)]
  z <- rho$error.structure$x %*% w
  z2r(z)
}

z2r <- function (z) {
  ifelse(z > 354, 1, (exp(2 * z) - 1)/(1 + exp(2 * z)))
}


getInd.coef <- function(coef.constraints, coef.values) {
  ind <- matrix(NA,ncol = ncol(coef.constraints), nrow = nrow(coef.constraints))
  k <- 1
  i <- 1
  for(j in 1:ncol(coef.constraints)){
    if (is.na(coef.values[i,j])){
      ind[i,j] <- k
      k <- k + 1
    }
  }
  for(i in 2:nrow(coef.constraints)){
    for(j in 1:ncol(coef.constraints)){
      if (is.na(coef.values[i,j]) && !(coef.constraints[i,j] %in% coef.constraints[1:(i-1),j])){
        ind[i,j] <- k
        k <- k + 1
      }
      if (is.na(coef.values[i,j]) && coef.constraints[i,j] %in% coef.constraints[1:(i-1),j]){
        ind[i,j] <- ind[min(which(coef.constraints[i,j] == coef.constraints[1:(i-1),j])),j]
      }
    }
  }
  ind
}

getInd.thresholds.flexible <- function(threshold.constraints,rho){
  pointer <- lapply(1:rho$ndim, function(j) rep(0,rho$npar.theta[j]))
  pointer[[1]] <- seq_len(rho$npar.theta[1])
  k <- length(pointer[[1]])
  for(i in 2:rho$ndim){
    if (rho$npar.theta[i] == 0){
      pointer[[i]] <- integer(0)
    } else {
      if (threshold.constraints[i] %in% threshold.constraints[1:(i-1)]){
        min.ind <- min(which(threshold.constraints[i] == threshold.constraints[1:(i-1)]))
        if (rho$ntheta[i] != rho$ntheta[min.ind])
          stop("Constraints on threshold parameters are not valid (different number of categories)", call. = FALSE)
        pointer[[i]] <- pointer[[min.ind]]
      } else{
        pointer[[i]] <- (k+1):(k+rho$ntheta[i])
        k <- k + rho$ntheta[i]
      }
    }
    }
  pointer
}

getInd.thresholds.fix2 <- function(threshold.constraints,rho) {
  #pointer <- lapply((1:rho$ndim)[which(rho$npar.theta > 0)], function(j) 1:rho$npar.theta[j])
  pointer <- lapply((1:rho$ndim), function(j) seq_len(rho$npar.theta[j]))

  k <- length(pointer[[1]])
  #for(i in (1:rho$ndim)[which(rho$npar.theta > 0)][-1]){
  for(i in (1:rho$ndim)){
    if (rho$npar.theta[i] == 0){
      pointer[[i]] <- integer(0)
    } else {
    if (threshold.constraints[i] %in% threshold.constraints[1:(i - 1)]){
      min.ind <- min(which(threshold.constraints[i] == threshold.constraints[1:(i-1)]))
      if (rho$ntheta[i] != rho$ntheta[min.ind])
        stop("Constraints on threshold parameters are not valid (different number of categories)", call. = FALSE)
      pointer[[i]] <- pointer[[min.ind]]
    } else{
      pointer[[i]] <- (k+1):(k+rho$ntheta[i]-2)
      k <- k + rho$ntheta[i]-2
    }
    }
  }
  pointer
}

getInd.thresholds.fix1 <- function(threshold.constraints,rho) {
  pointer <- lapply(1:rho$ndim, function(j) rep(0,rho$npar.theta[j]))
  pointer[[1]] <- seq_len(rho$npar.theta[1])#1:(rho$ntheta[1]-1)
  k <- length(pointer[[1]])
  for(i in 2:rho$ndim){
    if (rho$npar.theta[i] == 0){
      pointer[[i]] <- integer(0)
    } else {
    if (threshold.constraints[i] %in% threshold.constraints[1:(i - 1)]){
      min.ind <- min(which(threshold.constraints[i] == threshold.constraints[1:(i - 1)]))
      if (rho$ntheta[i] != rho$ntheta[min.ind])
        stop("Constraints on threshold parameters are not valid (different number of categories)", call. = FALSE)
      pointer[[i]] <- pointer[[min.ind]]
    } else{
      pointer[[i]] <- (k+1):(k+rho$ntheta[i]-1)
      k <- k + rho$ntheta[i]-1
    }
    }
  }
  pointer
}

getInd.thresholds.fixall <- function(threshold.constraints,rho) {
  lapply(1:rho$ndim, function(j) integer())
}

get.labels.theta <- function(rho,j) {
  lev <- levels(rho$y[, j])
  sapply(1:(rho$ntheta[j]), function(i){
    paste(lev[i], lev[i + 1], sep = "|")
  })
}

#' @title Error Structures in Multord
#' @description Different \code{error.structures} are available in \pkg{MultOrd}:
#' \itemize{
#' \item general correlation structure (default) \code{corGeneral(~1)},
#' \item general covariance structure \code{covGeneral(~1)},
#' \item factor dependent correlation structure \code{covGeneral(~f)},
#' \item factor dependent covariance structure \code{covGeneral(~f)},
#' \item covariate dependent equicorrelation structure \code{corEqui(~X)},
#' \item AR(1) correlation structure \code{corAR1(~1)}, or
#' \item factor dependent AR(1) correlation structure \code{corAR1(~f)}.
#' }
#' See \code{\link{error.structures}} or 'Details'.
#' @details General (symmetric) correlation structure (cross-sectional)
#' @param formula \code{\link{formula}} object
#' @export
#' @name error.structures
corGeneral <- function(formula) UseMethod("corGeneral")
#' @rdname error.structures
#' @export
corGeneral.formula <- function(formula){
  return(list(type = "corGeneral", formula = formula))
}

#' @details Equicorelation structure (cross-sectional)
# #' @param formula \code{\link{formula}} object
#' @export
#' @rdname error.structures
corEqui <- function(formula) UseMethod("corEqui")
#' @rdname error.structures
#' @export
corEqui.formula <- function(formula){
  return(list(type = "corEqui", formula = formula))
}

#' @details General (symmetric) covariance structure (cross-sectional)
# #' @param formula \code{\link{formula}} object
#' @export
#' @rdname error.structures
covGeneral <- function(formula) UseMethod("covGeneral")
#' @rdname error.structures
#' @export
covGeneral.formula <- function(formula){
  return(list(type = "covGeneral", formula = formula))
}

#' @details AR(1) correlation structure (panel; only applicable in line with \code{\link{multord}})
# #' @param formula \code{\link{formula}} object
#' @export
#' @rdname error.structures
corAR1 <- function(formula) UseMethod("corAR1")
#' @rdname error.structures
#' @export
corAR1.formula <- function(formula){
  return(list(type = "corAR1", formula = formula))
}


backtransf.sigmas <- function(R){
  J <- nrow(R)
  l <- t(chol(R))
  angmat <- matrix(1,ncol=J,nrow=J)
  angmat[-1,1] <- acos(l[-1,1])
  for (j in 2:(J-1)){
    sinprod <- apply(sin(angmat[, seq_len(j-1), drop=F]), 1, prod) ## denominator in division
    angmat[-(1:j),j]<-acos((l/sinprod)[-(1:j),j])
  }
  angdivpi <- angmat[lower.tri(angmat)]/pi
  log(angdivpi/(1-angdivpi))
}

get.sigma.i <- function(object) {
  if (object$rho$error.structure$type == "corEqui") {
    lapply(1:object$rho$n, function(i) {
        tmp <- matrix(object$r[i], nrow = object$rho$ndim, ncol = object$rho$ndim)
        diag(tmp) <- 1
        tmp
      })
  } else {
  lev <- match(object$rho$error.structure$x, object$rho$error.structure$levels)
  lapply(1:object$rho$n, function(i) object$sigmas[lev[i]][[1]])
  }
}

#' @title Data preparation for multord
#'
#' @description
#' This function is an (internally) used to transforms the \code{data}, into a "multivariate setting",
#' where all repeated measurements are matched accordingly to their ID. A matrix of all ordinal responses with \code{J} columns
#' as well as a list of length \code{J} of matrices with all the covariates are created.
#' @param data a \code{data.frame}, where each row corresponds to a single measurement.
#' @param index is an (optional) argument that specifies the index for the subjects and the response index of the multiple measurement.
#' This is usually performed
#' by a character vector of length two specifying the column names of the subject index and
#' the multiple response index in data. The default value of index is NULL assuming that the
#' first column of data contains the subject index and the second column the multiple response index.
#' @param y.names column name of \code{data} where the ordinal observations are stored.
#' @param x.names column names of all the covariates in {data}.
#' @param y.levels (optional) list of length \code{J} that specifies the levels of each repeated measurement. If the categories
#' differ across repeated measurements (either the number of categories or the category labels) it is recommended to set them.
#' @param response.names (optional) vector of names of the repeated measurements in \code{data}
#'which specifies the ordering of the repeated measurements.
#' @export
#y.names <- "yy"
#y.names <- rho$response.name
#x.names <- rho$x.names

# data <- data_multord
# y.names <- rho$response.name
# x.names <- unique(c(rho$x.names,weights))
# y.levels = response.levels
# response.names = response.names = c("rater1", "rater2", "rater3")
# index = c("firmID", "rater")

multord.data <- function(data, index, y.names, x.names,
                         y.levels = NULL, response.names = NULL) {
  df <- list()
  if (any(duplicated(data[,index]))) stop("duplicated observation(s) for one index", call. = FALSE)
  index.levels <- levels(as.factor(data[, index[2]]))
  data.split <- split(data[,c(y.names, index[1])], data[, index[2]])
  data.split.y <- lapply(seq_len(length(data.split)), function(j) {
    colnames(data.split[[j]]) <- c(index.levels[j],index[1])
    data.split[[j]]})
  df$y <- Reduce(function(...) merge(..., by = index[1], all = TRUE), data.split.y )[, -1]
  response.names.NA <- c()
  if (!is.null(response.names)) {
    response.names.NA <- response.names[!response.names %in% colnames(df$y)]
    df$y <- cbind(df$y, matrix(NA, ncol = length(response.names.NA), nrow = nrow(df$y), dimnames = list(c(),response.names.NA)))
    df$y <- df$y[,as.character(response.names)]
  }
  colnames.y <- colnames(df$y)
  if (is.null(y.levels)) {
    df$y <- do.call(cbind.data.frame, lapply(1:ncol(df$y), function(j) ordered(df$y[, j])))
  } else {
    df$y <- do.call(cbind.data.frame, lapply(1:ncol(df$y), function(j) {
      if (!all(unique(df$y[!is.na(df$y[, j]), j]) %in% y.levels[[j]])) stop("levels of response do not match with y.levels", call. = FALSE)
      ordered(df$y[, j], levels = y.levels[[j]])}
      ))
  }
  colnames(df$y) <- colnames.y#paste0(y.names, 1:ncol(df$y))

  #X
  data.split.x <- split(data[, c(x.names, index[1])], data[, index[2]])
  names.x <- names(data.split.x)
  #REMOVE this
  #names.x <- names(data.split.x)[names(data.split.x) %in% response.names]

  #set colnames (otherwise warning due to identical colnames in reduce)
  data.split.x <- lapply(seq_len(length(data.split.x)), function(j) {
    colnames(data.split.x[[j]]) <- c(paste(x.names,j, sep = "."), index[1])
    data.split.x[[j]]})


  xdatadf <- Reduce(function(...) merge(...,by = index[1], all = TRUE), data.split.x)[, -1]
  xdatadf <- cbind(xdatadf, matrix(NA, ncol = length(response.names.NA) * length(x.names),
                                   nrow = nrow(xdatadf)))
  #xdatadf <- Reduce(function(x,y) merge(x,y,by = index[1], all = TRUE), split(data[,c(x.names, index[1])], data[, index[2]]) )[, -1]

  #REMOVE this
  # df$x <- lapply(1:ncol(df$y), function(i) {
  #   tmp <- xdatadf[,(i - 1) * length(x.names) + seq_len(length(x.names))]
  #   names(tmp) <- x.names
  #   tmp
  # })
  df$x <- lapply(1:(length(index.levels) + length(response.names.NA)), function(i) {
    tmp <- xdatadf[,(i - 1) * length(x.names) + seq_len(length(x.names))]
    names(tmp) <- x.names
    tmp
  })

  names(df$x) <- c(names.x, response.names.NA)
  if (!is.null(response.names)) df$x <- df$x[as.character(response.names)]
  df
}

check <- function(...){
  stopifnot(...)
}

#threshold.values <- list(c(NA,NA,NA),c(2,3,NA),c(1,2,NA,NA,NA))
# set.model <- function(rho){
#   if (rho$error.structure$type == "corAR1"){
#     m <- "PMOR"
#   } else{
#   if (rho$error.structure$type == "corGeneral"){
#     #if (length(all.vars(rho$error.structure$formula[[2]])) > 1) stop("only one factor is allowed in corGeneral", call. = FALSE)
#     rho$model <- "CMORcor"
#   }
#   if (rho$error.structure$type == "corEqui"){
#     rho$model <- "CMORcor"
#   }
#   if (rho$error.structure$type == "covGeneral"){
#     #if (length(all.vars(rho$error.structure$formula[[2]])) > 1) stop("only one factor is allowed in covGeneral", call. = FALSE)
#     rho$model <- "CMORcov"
#   }
#       if (is.null(rho$threshold.values)){
#       if ((rho$model == "CMORcor") && (attr(terms.formula(rho$formula), "intercept") == 1)){
#         stop("Remove intercept in formula (~ 0 + ...) or fix the first threshold in threshold.values", call. = FALSE)
#       }
#       if ((rho$model == "CMORcov") && (attr(terms.formula(rho$formula), "intercept") == 0)){
#         stop("Include intercept in formula (~ 1 + ...)", call. = FALSE)
#       }
#       m <- rho$model
#     } else{
#     if (rho$model == "CMORcor"){
#       if (all(sapply(1:rho$ndim, function(j) all(is.na(rho$threshold.values[[j]]))))){ #all thresholds NA
#         m <- "CMORcor"
#     } else if (all(sapply(1:rho$ndim, function(j) (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && all(which(!is.na(rho$threshold.values[[j]])) == 1))))){ #all first thresholds are not NA (and no additional)
#         #warning if no intercept
#               if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix a threshold in model with corGeneral", call. = FALSE)
#                   m <- "CMORcor2"
#     } else if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in all dimensions
#           if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix all thresholds in model with corGeneral", call. = FALSE)
#                   m <- "CMORcor"
#     } else if (any(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in at least one dimension
#           ##CHECK somewhere if all values are fixed in one dimension + ?warning(sense)
#               if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))  ||  all(
#                 (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && (which(!is.na(rho$threshold.values[[j]])) == 1)))))){ #all not NA or only first fixed for each j
#                   if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix a threshold in model with corGeneral", call. = FALSE)
#                     m <- "CMORcor2"
#             } else if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))  ||  all(is.na(rho$threshold.values[[j]]))))){ #all not NA or all NA for each j
#                   if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix a threshold in model with corGeneral", call. = FALSE)
#                     m <- "CMORcor"
#             } else stop("Either fix all thresholds in one or more outcome dimensions,
#                          or consistently in all other outcome dimensions, all first thresholds or none.
#                          In addition, intercept must be included in the model.", call. = FALSE)
#     } else{
#           stop("Model is not identifiable. One of the following could help:
#                1. Remove intercept in formula (~ 0 + ...) or
#                2. Fix all first thresholds in threshold.values or
#                3. Fix all thresholds in threshold.values or
#                4. Either fix all thresholds in one or more outcome dimensions,
#                   or consistently in all other outcome dimensions, all first thresholds or none", call. = FALSE)
#         }
#     }
#     if (rho$model == "CMORcov"){
#       if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix a threshold in model CMORcov", call. = FALSE)
#       if (all(sapply(1:rho$ndim, function(j){
#             (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,2)) #all first two thresholds are not NA
#             ) || ((length(rho$threshold.values[[j]]) == 1) &&  (which(!is.na(rho$threshold.values[[j]])) == 1))
#             }))){
#                 m <- "CMORcov"
#     } else if (all(sapply(1:rho$ndim, function(j){
#         (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j]))#all first and last two thresholds are not NA
#         ) || ((length(rho$threshold.values[[j]]) == 1) &&  (which(!is.na(rho$threshold.values[[j]])) == 1))
#          }))){
#         m <- "CMORcov2"
#     } else if (any(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in at least one dimensions
#         ##CHECK somewhere if all values are fixed in one dimensions + ?warning(sense)
#           if (all(sapply(1:rho$ndim, function(j) #
#                 #(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1))) &&  all((which(!is.na(rho$threshold.values[[j]])) == 1))) ||
#                 (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) &&  all((which(!is.na(rho$threshold.values[[j]])) == c(1,2)))) ||
#                 (all(!is.na(rho$threshold.values[[j]])))))){
#                     m <- "CMORcov"
#         } else if (all(sapply(1:rho$ndim, function(j)
#                 #(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1))) &&  all((which(!is.na(rho$threshold.values[[j]])) == 1))) ||
#                 (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) &&  all((which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j])))) ||
#                 (all(!is.na(rho$threshold.values[[j]])))))){
#                     m <- "CMORcov2"
#         } else stop("Either fix all thresholds in one or more outcome dimensions,
#                       or consistently in all other outcome dimensions, all first two thresholds or all first and last thresholds", call. = FALSE)
#     } else stop("Either fix first two or first and last threshold in threshold.values if you use a model with covGeneral", call. = FALSE)
#   }
#     }
#   }
#   m
# }

Try the MultOrd package in your browser

Any scripts or data that you put into this service are public.

MultOrd documentation built on May 2, 2019, 4:49 p.m.