R/family.categorical.R

Defines functions ordsup.vglm R2latvar is.zero.vglm is.zero.matrix is.parallel.vglm is.parallel.matrix prplot.control prplot margefff cratio.derivs subsetarray3 margeff findFirstMethod ordpoissonProbs ordpoisson InverseBrat Brat .brat.indices .brat.alpha bratt brat acat.deriv acat propodds cumulative vglm.VGAMcategorical.control vglm.multinomial.control vglm.multinomial.deviance.control cratio sratio dmultinomial Deviance.categorical.data.vgam multinomial CM.equid CM.symm0 CM.symm1 CM.qlogis CM.qnorm CM.ones CM.free rdiag simslotVGAMcat

Documented in acat acat.deriv brat Brat bratt CM.equid CM.equid CM.free CM.ones CM.qlogis CM.qnorm CM.qnorm CM.symm0 CM.symm0 CM.symm1 CM.symm1 CM.symm1 CM.symm1 cratio cratio.derivs cumulative Deviance.categorical.data.vgam dmultinomial findFirstMethod InverseBrat is.parallel.matrix is.parallel.vglm is.zero.matrix is.zero.vglm margeff multinomial ordpoisson ordsup.vglm propodds prplot prplot.control R2latvar simslotVGAMcat sratio subsetarray3 vglm.multinomial.control vglm.multinomial.deviance.control vglm.VGAMcategorical.control

# These functions are
# Copyright (C) 1998-2024 T.W. Yee, University of Auckland.
# All rights reserved.




















 simslotVGAMcat <- function(object, nsim) {
    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (max(abs(pwts - round(pwts))) > 0.001 ||
        any(pwts < 0.5))
      stop("prior weights seem not to be ",
           "positive or integer-valued")
    pwts <- c(round(pwts))  # Make sure, a vector
    nn <- nobs(object)
    M <- npred(object)
    fv <- fitted(object)
    label.it <- length(colnames(fv)) &&
                length(rownames(fv))
    ansarr <- array(NA_real_, c(M + 1, nsim, nn))
    if (label.it)
      dimnames(ansarr) <-
        list(colnames(fv),
             paste("sim", seq_len(nsim), sep = "_"),
             rownames(fv))
    if (object@misc$intercept.only &&
        all(pwts == pwts[1]) &&
        all(object@offset == 0)) {
      tmp2 <- c(rmultinom(nsim * nn, pwts[1], fv[1, ]))
      ansarr[] <- tmp2
    } else {
      for (i in seq(nn))
        ansarr[,, i] <- rmultinom(nsim, pwts[i], fv[i, ])
    }
    ansarr <- aperm(ansarr, c(3, 1, 2))
    attr(ansarr, "Verbatim") <- TRUE  # Removed later
    ansarr
}  # simslotVGAMcat










rdiag <- function(...) {
  a <- as.vector(unlist(list(...)))
  if (length(a) == 1) {
      m <- matrix(0, a, a)
      m[cbind(a:1, 1:a)] <- 1
  } else {
      m <- matrix(0, length(a), length(a))
      m[cbind(length(a):1, 1:length(a))] <- a
  }
  m
}  # rdiag



 CM.free <- function(M, Trev = FALSE, Tref = 1) {
  diag(M)
}  # CM.free



 CM.ones <- function(M, Trev = FALSE, Tref = 1) {
  matrix(1, M, 1)
}  # CM.ones



 CM.qnorm <- function(M, Trev = FALSE, Tref = 1) {
  matrix(qnorm(seq(M) / (M + 1)) *
         ifelse(Trev, -1, 1), M, 1) 
}  # CM.qnorm



 CM.qlogis <- function(M, Trev = FALSE, Tref = 1) {
  matrix(qlogis(seq(M) / (M + 1)) *
         ifelse(Trev, -1, 1), M, 1) 
}  # CM.qlogis







 CM.symm1 <- function(M, Trev = FALSE, Tref = 1) {
  if (M <  1)
    stop("argument 'M' should never be < 1")
  if (M == 1)
    return(cbind(1))
  if (M == 2)
    return(cbind(1, c(-1, 1)))
  Modd <- 1 + floor(M / 2)
  H1 <- matrix(0, M, Modd)
  if (M %% 2 == 0)  # Add another row
    H1 <- rbind(0, H1)
  H1[, 1] <- 1
  for (jay in 1:(Modd - 1)) {
    H1[Modd - jay, 1 + jay] <- -1
    H1[Modd + jay, 1 + jay] <-  1
  }
  if (M %% 2 == 0)
    H1 <- H1[-Modd, , drop = FALSE]
  H1
}  # CM.symm1



 CM.symm0 <- function(M, Trev = FALSE, Tref = 1) {
  H1 <- CM.symm1(M)
  if (M == 1)
    stop("cannot have 'symm0' when M == 1")
  H1[, -1, drop = FALSE]
}  # CM.symm0










 CM.equid <-
    function(M,
             Trev = FALSE,
             Tref = 1) {
  if (!is.logical(Trev) && length(Trev) != 1)
    stop("bad input for argument 'Trev'")
  if (is.character(Tref) && Tref == "M")
    Tref <- M
  if (!is.numeric(Tref))
    Tref <- as.numeric(Tref)
  if (M <  1)
    stop("argument 'M' should never be less than 1")
  if (M == 1)
    return(cbind(1))
  H1 <- matrix(1, M, 2)
  H1[, 2] <- if (Trev) rev(seq(M)) else seq(M)
  H1[, 2] <- H1[, 2] - H1[Tref, 2]
  H1
}  # CM.equid














 multinomial <-
  function(zero = NULL, parallel = FALSE,
      nointercept = NULL, refLevel = "(Last)",
      ynames = FALSE,  # 20240216
      imethod = 1,
      imu = NULL, byrow.arg = FALSE,  # 20211105
      Thresh = NULL,  # "free",
      Trev = FALSE,  # reverse (propodds())
      Tref = if (Trev) "M" else 1,
      whitespace = FALSE) {


      

  if (length( Thresh ) &&
      !exists(paste0("CM.", Thresh),
              mode = "function"))
    stop("No function ", paste0("CM.", Thresh))

  if (length(refLevel) != 1)
    stop("the length of 'refLevel' must be one")

  if ( is.numeric(refLevel) &&
      !is.Numeric(refLevel, integer.valued = TRUE,
                  positive = TRUE))
    stop("arg 'refLevel' not a positive integer")

  if (is.character(refLevel)) {
    if (refLevel == "(Last)")
      refLevel <- -1
  }
  if (is.factor(refLevel)) {
    if (is.ordered(refLevel))
      warning("argument 'refLevel' is from an ",
              "ordered factor")

    refLevel <- as.character(refLevel) == levels(refLevel)
    refLevel <- (seq_along(refLevel))[refLevel]
    if (!is.Numeric(refLevel, length.arg = 1,
        integer.valued = TRUE, positive = TRUE))
      stop("could not coerce 'refLevel' into a",
           " single positive integer")
  }

  if (!is.logical(ynames) || length(ynames) != 1)
    stop("bad input for 'ynames'")



  stopifnot(is.logical(whitespace) &&
            length(whitespace) == 1)
  fillerChar <- ifelse(whitespace, " ", "")


  new("vglmff",
      blurb = c("Multinomial logit model\n\n",
            "Links:    ",
         if (is.numeric(refLevel)) {
         if (refLevel < 0) {
           ifelse(whitespace,
             "log(mu[,j] / mu[,M+1]), j = 1:M,\n",
             "log(mu[,j]/mu[,M+1]), j=1:M,\n")
         } else {
            if (refLevel == 1) {
              paste0("log(mu[,", "j]",
                     fillerChar, "/", fillerChar,
                     "mu[,", refLevel, "]), j",
                     fillerChar, "=", fillerChar,
                     "2:(M+1),\n")
            } else {
              paste0("log(mu[,", "j]", fillerChar,
                     "/",
                     "mu[,", refLevel, "]), j",
                     fillerChar, "=",
                     fillerChar,
                     "c(1:", refLevel-1,
                     ",", fillerChar, refLevel+1,
                     ":(M+1)),\n")
            }
         }
         } else {  # refLevel is character
             paste0("log(mu[,", "j]", fillerChar,
                    "/",
                    "mu[,'", refLevel, "']), j",
                    fillerChar, " != '",
                    fillerChar,
                    refLevel, "',\n")
         },
         "Variance: ",
           ifelse(whitespace,
        "mu[,j] * (1 - mu[,j]); -mu[,j] * mu[,k]",
        "mu[,j]*(1-mu[,j]); -mu[,j]*mu[,k]")),

  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(1, M, 1), x = x,
                           bool = .parallel ,
                       apply.int = TRUE,
                       constraints = constraints)
    constraints <- cm.zero.VGAM(constraints, x = x,
              .zero , M = M, M1 = M,
              predictors.names = predictors.names)
 constraints <- cm.nointercept.VGAM(constraints, x,
                             .nointercept , M)

    if (length( .Thresh ))
      constraints[["(Intercept)"]] <-
        do.call(paste0("CM.", .Thresh ),
                list(M = M, Trev = .Trev ,
                     Tref = .Tref ))
  }),
  list( .parallel = parallel, .zero = zero,
        .Thresh = Thresh, .Tref = Tref,
        .Trev = Trev,
        .nointercept = nointercept ))),

  deviance = Deviance.categorical.data.vgam,

  infos = eval(substitute(function(...) {
    list(parallel = .parallel ,
         refLevel = .refLevel ,  # original
         M1 = -1,
         link = "multilogitlink",
         link1parameter = FALSE,  # multiparameter
         mixture.links = FALSE,  # Not a mixture.
         expected = TRUE,
         imethod = .imethod ,
         multipleResponses = FALSE,
         parameters.names = as.character(NA),
         Thresh = .Thresh ,
         Tref = .Tref ,
         Trev = .Trev ,
         ynames = .ynames ,
         zero = .zero )
  },
  list( .zero = zero, .ynames = ynames,
        .refLevel = refLevel, .imethod = imethod,
        .Thresh = Thresh,
        .Tref = Tref ,
        .Trev = Trev ,
        .parallel = parallel ))),
  initialize = eval(substitute(expression({

    if (is.factor(y) && is.ordered(y))
      warning("response should be nominal, not ordinal")


    delete.zero.colns <- TRUE
    eval(process.categorical.data.VGAM)

    if ( .imethod == 2) {  # MLE for intercept-only.
      mustart <- matrix(colMeans(y), n, NCOL(y), byrow = TRUE)
    }
    if ( .imethod > 2) stop("argument 'imethod' unmatched")
    

    M <- ncol(y) - 1
    use.refLevel <- if (is.numeric( .refLevel )) {
      if ( .refLevel < 0) M + 1 else ( .refLevel )
    } else {  # Is character. Match it with the response levels.
      tmp6 <- match( .refLevel , colnames(y))
      if (is.na(tmp6))
        stop("could not match 'refLevel' with ",
             "any columns of the response matrix")
      tmp6
    }


    if (use.refLevel > (M + 1))
      stop("argument 'refLevel' has a value that is too high")
    extra$use.refLevel <- use.refLevel  # Used in
    pmone <- rep_len(1, M + 1)


    allbut.refLevel <- seq(M + 1)[-use.refLevel]
    if ( .ynames )
      if (is.null(cnsy <- colnames(y)))
        stop("response matrix has no colnames")
    choice1 <- if ( .ynames )
      cnsy[allbut.refLevel] else allbut.refLevel
    choice2 <- if ( .ynames )
      cnsy[use.refLevel] else use.refLevel
    choice3 <- if ( .ynames ) "P[" else "mu[,"
    predictors.names <-
        paste0("log(", choice3, choice1, "]",
               .fillerChar, "/", .fillerChar,
               choice3, choice2, "])")






    extra$colnames.y <- colnames(y)


    imu <- as.vector( .imu )
    if (length(imu)) {
      mustart <- matrix(imu, n, NCOL(y), byrow = .byrow.arg )
    }
  }),
  list( .refLevel = refLevel,
        .ynames = ynames,
        .fillerChar = fillerChar,
        .imu = imu, .byrow.arg = byrow.arg,
        .imethod = imethod,
        .whitespace = whitespace ))),
  linkinv = eval(substitute( function(eta, extra = NULL) {
    if (anyNA(eta))
      warning("there are NAs in eta in slot inverse")
    ans <-
      multilogitlink(eta,  # .refLevel ,
                     refLevel = extra$use.refLevel,
                     inverse = TRUE)
    if (anyNA(ans))
      warning("there are NAs here in slot linkinv")
    if (min(ans) == 0 || max(ans) == 1)
      warning("fitted probabilities numerically ",
              "0 or 1 occurred")

    label.cols.y(ans, colnames.y = extra$colnames.y, NOS = 1)
  }), list( .refLevel = refLevel )),
  last = eval(substitute(expression({
    misc$link <- "multilogitlink"

    misc$earg <- list(multilogitlink = list(
      M = M,
      refLevel = use.refLevel
    ))

    dy <- dimnames(y)
    if (!is.null(dy[[2]]))
      dimnames(fit$fitted.values) <- dy

    misc$nointercept <- ( .nointercept )
    misc$parallel <- ( .parallel )
    misc$refLevel <- use.refLevel
    misc$refLevel.orig <- ( .refLevel )
    misc$zero <- ( .zero )
  }),
  list( .refLevel = refLevel, .nointercept = nointercept,
        .parallel = parallel,
        .zero = zero ))),

  linkfun = eval(substitute( function(mu, extra = NULL) {

    multilogitlink(mu,
                   refLevel = extra$use.refLevel)
  }), list( .refLevel = refLevel )),

  loglikelihood =
    function(mu, y, w, residuals = FALSE, eta, extra = NULL,
             summation = TRUE) {
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ycounts <- if (is.numeric(extra$orig.w))
                   y * w / extra$orig.w else
                   y * w  # Convert proportions to counts
      nvec <- if (is.numeric(extra$orig.w))
                round(w / extra$orig.w) else round(w)

      smallno <- 1.0e4 * .Machine$double.eps
      if (max(abs(ycounts - round(ycounts))) > smallno)
        warning("converting 'ycounts' to integer in ",
                "@loglikelihood")
      ycounts <- round(ycounts)
      ll.elts <-
        (if (is.numeric(extra$orig.w)) extra$orig.w else 1) *
         dmultinomial(x = ycounts, size = nvec, prob = mu,
                      log = TRUE, dochecking = FALSE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  },
  vfamily = c("multinomial", "VGAMcategorical"),




  simslot = eval(substitute(
  function(object, nsim) {
    simslotVGAMcat(object, nsim)
  }, list( .refLevel = refLevel ))),



  hadof = eval(substitute(
  function(eta, extra = list(),
           linpred.index = 1, w = 1,
           dim.wz = c(NROW(eta), NCOL(eta) * (NCOL(eta)+1)/2), 
           deriv = 1, ...) {


  multinomial.eim.deriv1 <- function(mu, jay, w) {
    M <- ncol(mu)  # Not ncol(mu)-1 coz 1 coln has been deleted
    MM12 <- M * (M + 1) / 2  # Full matrix
    wz1 <- matrix(0, NROW(mu), MM12)
    ind5 <- iam(NA, NA, M = M, both = TRUE)
    for (i in 1:MM12) {
      i1 <- ind5$row.index[i]
      j1 <- ind5$col.index[i]
      if (i1 == jay && j1 == jay) {
        wz1[, iam(i1, j1, M = M)] <-  (1 - 2 * mu[, i1])  # *
      }
      if (i1 == jay && j1 != jay) {
        wz1[, iam(i1, j1, M = M)] <- -mu[, j1]  # -
      }
      if (i1 != jay && j1 == jay) {
        wz1[, iam(i1, j1, M = M)] <- -mu[, i1]  # -(1 - ...) *
      }
    }  # for (i)
    c(w) * wz1
  }  # multinomial.eim.deriv1



  multinomial.eim.deriv2 <- function(mu, jay, w) {
    M <- ncol(mu)  # Not ncol(mu)-1 coz 1 coln has been deleted
    zz.yettodo <- 0  # NA_real_
    MM12 <- M * (M + 1) / 2  # Full matrix
    wz1 <- matrix(NA_real_, NROW(mu), MM12)
    ind5 <- iam(NA, NA, M = M, both = TRUE)
    for (i in 1:MM12) {
      i1 <- ind5$row.index[i]
      j1 <- ind5$col.index[i]
      if (i1 == jay && j1 == jay) {
        wz1[, iam(i1, j1, M = M)] <-  zz.yettodo
      }
      if (i1 != jay && j1 == i1) {
        wz1[, iam(i1, j1, M = M)] <- -mu[, i1] * mu[, jay] *
                           (1 - 2 * mu[, i1] - 2 * mu[, jay] +
                                6 * mu[, i1] * mu[, jay])
                                            
      }
      if (i1 == jay && j1 != jay) {
        wz1[, iam(i1, j1, M = M)] <- zz.yettodo
      }
      if (i1 != jay && j1 == jay) {
        wz1[, iam(i1, j1, M = M)] <- zz.yettodo
      }
      if (any(is.na(wz1[, iam(i1, j1, M = M)]))) {
        wz1[, iam(i1, j1, M = M)] <-   2 * mu[, i1]  * 
                                       mu[, j1] * mu[, jay] *
                                       (1 - 3 * mu[, jay])
      }
    }  # for (i)
 cat("\n\n\n\n")
    c(w) * wz1
  }  # multinomial.eim.deriv2


    M <- NCOL(eta)
    use.refLevel <- extra$use.refLevel  # Restore its value
    if (!is.numeric(use.refLevel)) {
      warning("variable 'use.refLevel' cannot be found. ",
              "Trying the original value.")
      use.refLevel <- .refLevel  # Only if numeric...
      if (use.refLevel == "(Last)")
        use.refLevel <- M+1
    }
    mu.use <- multilogitlink(eta, refLevel = use.refLevel,
                             inverse = TRUE)
    mu.use <- pmax(mu.use, .Machine$double.eps * 1.0e-0)


    index <- iam(NA, NA, M, both = TRUE, diag = TRUE)
    myinc <- (index$row.index >= use.refLevel)
    index$row.index[myinc] <- index$row.index[myinc] + 1
    myinc <- (index$col.index >= use.refLevel)
    index$col.index[myinc] <- index$col.index[myinc] + 1


    switch(as.character(deriv),
       "0" =  {
         wz <- -mu.use[, index$row] * mu.use[, index$col]
         wz[, 1:M] <- wz[, 1:M] + mu.use[, -use.refLevel ]
         c(w) * wz
       },
       "1" = {
         multinomial.eim.deriv1(mu.use[, -use.refLevel,
                                       drop = FALSE],
                                jay = linpred.index, w = w)
       },
       "2" = {
         multinomial.eim.deriv2(mu.use[, -use.refLevel,
                                       drop = FALSE],
                                jay = linpred.index, w = w)
       },
       stop("argument 'deriv' must be 0 or 1 or 2"))
  }, list( .refLevel = refLevel  # End of @hadof
         ))),




 validparams =
   eval(substitute(function(eta, y, extra = NULL) {
    probs <-
      multilogitlink(eta, refLevel = extra$use.refLevel,
                     inverse = TRUE)  # ( .refLevel )
     okay1 <- all(is.finite(probs)) &&
              all(0 < probs & probs < 1)
    okay1
  }, list( .refLevel = refLevel ))),
  deriv = eval(substitute(expression({
    use.refLevel <- extra$use.refLevel  # Restore its value
    ansd <- ( y[, -use.refLevel, drop = FALSE] -
             mu[, -use.refLevel, drop = FALSE])  # AMLM


    c(w) * ansd
  }), list( .refLevel = refLevel ))),
  weight = eval(substitute(expression({
    mytiny <- (mu <       sqrt(.Machine$double.eps)) |
              (mu > 1.0 - sqrt(.Machine$double.eps))

    if (M == 1) {
      wz <- mu[, 3 - use.refLevel] * (1 - mu[, 3 - use.refLevel])
    } else {  # M > 1
      index <- iam(NA, NA, M, both = TRUE, diag = TRUE)
      myinc <- (index$row.index >= use.refLevel)
      index$row.index[myinc] <- index$row.index[myinc] + 1
      myinc <- (index$col.index >= use.refLevel)
      index$col.index[myinc] <- index$col.index[myinc] + 1
      wz <- -mu[, index$row, drop = FALSE] *
             mu[, index$col, drop = FALSE]
      wz[, 1:M] <- wz[, 1:M] + mu[, -use.refLevel ]  # AMLM
    }





    atiny <- (mytiny %*% rep(1, ncol(mu))) > 0
    if (any(atiny)) {
      if (M == 1)
        wz[atiny] <- wz[atiny] * (1 + .Machine$double.eps^0.5) +
                                      .Machine$double.eps else
        wz[atiny, 1:M] <- .Machine$double.eps +
        wz[atiny, 1:M] * (1 + .Machine$double.eps^0.5)
    }  # atiny
    c(w) * wz
  }), list( .refLevel = refLevel ))))
}  # multinomial()







 process.categorical.data.VGAM <- expression({






    extra$y.integer <- TRUE


  if (!all(w == 1))
    extra$orig.w <- w

  if (!is.matrix(y)) {
    yf <- as.factor(y)
    lev <- levels(yf)
    llev <- length(lev)
    nn <- length(yf)
    y <- matrix(0, nn, llev)
    y[cbind(1:nn,as.vector(unclass(yf)))] <- 1
    dimnames(y) <- list(names(yf), lev)

    if (llev <= 1)
      stop("the response matrix does not have ",
           "2 or more columns")
  } else {
    nn <- nrow(y)
  }


  nvec <- rowSums(y)

  if (min(y) < 0 || any(round(y) != y))
    stop("the response must be non-negative ",
         "counts (integers)")

  if (!exists("delete.zero.colns") ||
      (exists("delete.zero.colns") && delete.zero.colns)) {
    sumy2 <- colSums(y)
    if (any(index <- sumy2 == 0)) {
      y <- y[, !index, drop = FALSE]
      sumy2 <- sumy2[!index]
      if (all(index) || ncol(y) <= 1)
        stop("'y' matrix has 0 or 1 columns")
      warning("Deleted ", sum(!index), " columns of the", 
              " response matrix due to zero counts")
    }
  }


  if (any(miss <- (nvec == 0))) {
    smiss <- sum(miss)
    warning("Deleted ", smiss,
            " rows of the response matrix due to zero counts")
    x <- x[!miss,, drop = FALSE]
    y <- y[!miss,, drop = FALSE]
    w <- cbind(w)
    w <- w[!miss,, drop = FALSE]

    nvec <- nvec[!miss]
    nn <- nn - smiss
  }

  w <- w * nvec

  nvec[nvec == 0] <- 1
  y <- prop.table(y, 1)   # Convert to proportions


  if (length(mustart) + length(etastart) == 0) {
      mustart <- y + (1 / ncol(y) - y) / nvec
  }
})





Deviance.categorical.data.vgam <-
  function(mu, y, w, residuals = FALSE, eta, extra = NULL,
           summation = TRUE) {




  if (ncol(y) == 1 || ncol(mu) == 1)
    stop("arguments 'y' and 'mu' must have at ",
         "least 2 columns")



  double.eps <- sqrt( .Machine$double.xmin )
  devy <- y
  nonz <- (y != 0)
  devy[nonz] <- y[nonz] * log(y[nonz])

  devmu <- 0 * y  # filler; y*log(mu) gives a warning
  if (any(smallmu <- (mu * (1 - mu) < double.eps))) {
    warning("fitted values close to 0 or 1")
    smu <- mu[smallmu]
    smy <-  y[smallmu]
    smu <- ifelse(smu < double.eps, double.eps, smu)
    devmu[smallmu] <- smy * log(smu)
  }
  devmu[!smallmu] <- y[!smallmu] * log(mu[!smallmu])

  devi <- 2 * (devy - devmu)

  if (residuals) {
    M <- if (is.matrix(eta)) ncol(eta) else 1
    if (M > 1)
      return(NULL)
    devi <- devi %*% rep_len(1, ncol(devi))  #Dev=\sum_i devi[i]
    return(c(sign(y[, 1] - mu[, 1]) * sqrt(abs(devi) * w)))
  } else {
    dev.elts <- c(w) * devi
    if (summation) {
      sum(dev.elts)
    } else {
      dev.elts
    }
  }
}







 dmultinomial <-
    function(x, size = NULL, prob, log = FALSE,
             dochecking = TRUE, smallno = 1.0e-7) {
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)



  x <- as.matrix(x)
  prob <- as.matrix(prob)
  if (((K <- ncol(x)) <= 1) ||
       ncol(prob) != K)
    stop("arguments 'x' and 'prob' must be matrices with ",
         "two or more columns")
  if (dochecking) {
    if (min(prob) < 0)
      stop("argument 'prob' contains some negative values")
    if (any(abs((rsprob <- rowSums(prob)) - 1) > smallno))
      stop("some rows of 'prob' do not add to unity")
    if (any(abs(x - round(x)) > smallno))
      stop("argument 'x' should be integer-valued")
    if (length(size)) {
      if (any(abs(size - rowSums(x)) > smallno))
        stop("rowSums(x) does not agree with argument 'size'")
    } else {
      size <- round(rowSums(x))
    }
  } else {
    if (!length(size))
      size <- round(rowSums(prob))
  }
  logdensity <- lfactorial(size) +
                rowSums(x * log(prob) - lfactorial(x))
  if (log.arg) logdensity else exp(logdensity)
}  # dmultinomial()











 sratio <-
  function(link = "logitlink",
           parallel = FALSE, reverse = FALSE,
           zero = NULL,
           ynames = FALSE,
    Thresh = NULL,   # "free",
    Trev = reverse,
    Tref = if (Trev) "M" else 1,
    whitespace = FALSE) {
  link <- as.list(substitute(link))
  earg  <- link2list(link)
  link <- attr(earg, "function.name")

  if (length( Thresh ) &&
      !exists(paste0("CM.", Thresh),
              mode = "function"))
    stop("No function ", paste0("CM.", Thresh))

  if (!is.logical(reverse) || length(reverse) != 1)
    stop("arg 'reverse' not a single logical")

  stopifnot(is.logical(whitespace) &&
            length(whitespace) == 1)
  fillerChar <- ifelse(whitespace, " ", "")

  if (!is.logical(ynames) || length(ynames) != 1)
    stop("bad input for 'ynames'")

  new("vglmff",
  blurb = c("Stopping ratio model\n\n",
            "Links:    ",
            namesof(if (reverse)
            ifelse(whitespace, "P[Y = j+1|Y <= j+1]",
                               "P[Y=j+1|Y<=j+1]") else
            ifelse(whitespace, "P[Y = j|Y >= j]",
                               "P[Y=j|Y>=j]"),
                   link, earg = earg), "\n",
            "Variance: ",
            ifelse(whitespace,
                   "mu[,j] * (1 - mu[,j]); -mu[,j] * mu[,k]",
                   "mu[,j]*(1-mu[,j]); -mu[,j]*mu[,k]")),
  infos = eval(substitute(function(...) {
    list(M1 = NA,  # zz -1?
         Q1 = NA,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = as.character(NA),
         parallel = .parallel ,
         reverse = .reverse ,
         Thresh = .Thresh ,
         Tref = .Tref ,
         Trev = .Trev ,
         whitespace = .whitespace ,
         ynames = .ynames ,
         zero = .zero ,
         link = .link )
  },
  list( .link = link, .ynames = ynames,
        .zero = zero,
        .parallel = parallel,
        .reverse = reverse,
        .Thresh = Thresh,
        .Tref = Tref,
        .Trev = Trev,
        .whitespace = whitespace ))),

  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(1, M, 1), x = x,
                     bool = .parallel ,
                     constraints = constraints)
    constraints <- cm.zero.VGAM(constraints, x = x,
             .zero , M = M, M1 = M,
             predictors.names = predictors.names)




    if (names(constraints)[1] != "(Intercept)")
      stop("the model has no intercept term")



    if (length( .Thresh ))
      constraints[["(Intercept)"]] <-
        do.call(paste0("CM.", .Thresh ),
                list(M = M, Trev = .Trev ,
                     Tref = .Tref ))
  }),
  list( .parallel = parallel,
        .Thresh = Thresh,
        .Tref = Tref ,
        .Trev = Trev ,
        .zero = zero ))),
  deviance = Deviance.categorical.data.vgam,

  initialize = eval(substitute(expression({

    if (is.factor(y) && !is.ordered(y))
      warning("response should be ordinal---see ordered()")



    delete.zero.colns <- TRUE
    eval(process.categorical.data.VGAM)
    extra$wy.prod <- TRUE
    M <- ncol(y) - 1


    if ( .ynames )  # overwrite mynames
      if (is.null(cnsy <- colnames(y)))
        stop("response matrix has no colnames")
    choice1 <- if ( .ynames ) cnsy[1:M] else
               1:M
    choice2 <- if ( .ynames ) cnsy[2:(M+1)] else
               2:(M+1)
    choice3 <- ifelse( .ynames , "P[", "P[Y")
    choice3 <- if ( .ynames ) paste0(choice3,
      if ( .reverse ) choice2 else choice1) else
      paste0(choice3,
             .fillerChar, "=", .fillerChar ,
             if ( .reverse ) choice2 else choice1)
    mynames <- if ( .reverse )
       paste0(choice3,
              ifelse( .ynames , "|", "|Y"),
              .fillerChar , "<=", .fillerChar ,
              choice2, "]") else
       paste0(choice3,
              ifelse( .ynames , "|", "|Y"),
              .fillerChar , ">=", .fillerChar ,
              choice1, "]")


    predictors.names <-
      namesof(mynames, .link , short = TRUE, earg = .earg )
    y.names <- param.names("mu", M+1)

    extra$mymat <- if ( .reverse )
       tapplymat1(y, "cumsum") else
       tapplymat1(y[, ncol(y):1, drop = FALSE],
                  "cumsum")[, ncol(y):1, drop = FALSE]
    extra$colnames.y  <- colnames(y)
  }),
  list( .earg = earg, .link = link,
       .reverse = reverse, .ynames = ynames,
        .fillerChar = fillerChar,
        .whitespace = whitespace ))),

  linkinv = eval(substitute( function(eta, extra = NULL) {
    if (!is.matrix(eta))
      eta <- as.matrix(eta)
    fv.mat <- if ( .reverse ) {
      M <- NCOL(eta)
      djr <- eta2theta(eta, .link , earg = .earg )
      temp <- tapplymat1(1 - djr[, M:1, drop = FALSE],
                         "cumprod")[, M:1, drop = FALSE]
      cbind(1, djr) * cbind(temp, 1)
    } else {
      dj <- eta2theta(eta, .link , earg = .earg )
      temp <- tapplymat1(1 - dj, "cumprod")
      cbind(dj, 1) * cbind(1, temp)
    }
    label.cols.y(fv.mat, colnames.y = extra$colnames.y,
                 NOS = 1)
  },
  list( .earg = earg, .link = link,
        .reverse = reverse) )),
  last = eval(substitute(expression({
    misc$link <- rep_len( .link , M)
    names(misc$link) <- mynames

    misc$earg <- vector("list", M)
    names(misc$earg) <- names(misc$link)
    for (ii in 1:M)
      misc$earg[[ii]] <- .earg

    misc$parameters <- mynames
    misc$reverse <- .reverse
    misc$fillerChar <- .fillerChar
    misc$whitespace <- .whitespace

    extra <- list()  # kill what was used
  }),
  list( .earg = earg, .link = link, .reverse = reverse,
        .fillerChar = fillerChar,
        .whitespace = whitespace ))),
  linkfun = eval(substitute( function(mu, extra = NULL) {
    cump <- tapplymat1(mu, "cumsum")
    if ( .reverse ) {
      djr <- mu[, -1, drop = FALSE] / cump[, -1, drop = FALSE]
      theta2eta(djr, .link , earg = .earg )
    } else {
      M <- ncol(mu) - 1
      dj <- if (M == 1)
              mu[, 1, drop = FALSE] else
              mu[, 1:M, drop = FALSE] / (
              1 - cbind(0, cump[, 1:(M-1), drop = FALSE]))
      theta2eta(dj, .link , earg = .earg )
    }
  },
  list( .earg = earg, .link = link, .reverse = reverse) )),
  loglikelihood =
    function(mu, y, w, residuals = FALSE, eta, extra = NULL,
             summation = TRUE)
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ycounts <- if (is.numeric(extra$orig.w))
                 y * w / extra$orig.w else
                 y * w  # Convert proportions to counts
      nvec <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else round(w)

      smallno <- 1.0e4 * .Machine$double.eps
      if (max(abs(ycounts - round(ycounts))) > smallno)
        warning("converting 'ycounts' to integer ",
                "in @loglikelihood")
      ycounts <- round(ycounts)

      ll.elts <-
        (if (is.numeric(extra$orig.w)) extra$orig.w else 1) *
         dmultinomial(x = ycounts, size = nvec, prob = mu,
                      log = TRUE, dochecking = FALSE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    },
  vfamily = c("sratio", "VGAMordinal", "VGAMcategorical"),

  simslot = eval(substitute(
  function(object, nsim) {
    simslotVGAMcat(object, nsim)
  }, list(.link = link, .earg = earg ))),


  validparams =
      eval(substitute(function(eta, y, extra = NULL) {
    djr <- eta2theta(eta, .link , earg = .earg )  # dj or djr
    okay1 <- all(is.finite(djr)) && all(0 < djr & djr < 1)
    okay1
  }, list( .earg = earg, .link = link, .reverse = reverse) )),
  deriv = eval(substitute(expression({
    if (!length(extra$mymat)) {
      extra$mymat <-
        if ( .reverse )
          tapplymat1(y, "cumsum") else
          tapplymat1(y[, ncol(y):1, drop = FALSE],
             "cumsum")[, ncol(y):1, drop = FALSE]
    }

    ans <- if ( .reverse ) {
      djr <- eta2theta(eta, .link , earg = .earg )
      ddjr.deta <- dtheta.deta(djr, .link , earg = .earg )
      Mp1 <- ncol(extra$mymat)
      if ( .link == "logitlink") {
        ans2 <- c(w) *  # New
                      (y[, -1, drop = FALSE] * (1 - djr) -
           extra$mymat[, -Mp1, drop = FALSE] * djr)
        ans2
      } else {  # Not the default link
        c(w) * ddjr.deta *
        (y[, -1, drop = FALSE] / djr -
         extra$mymat[, -Mp1, drop = FALSE] / (1 - djr))
      }
    } else {  # Not reverse
      dj <- eta2theta(eta, .link , earg = .earg )
      ddj.deta <- dtheta.deta(dj, .link , earg = .earg )
      if ( .link == "logitlink") {
        ans2 <- c(w) *  # New
              (y[, -ncol(y), drop = FALSE] * (1 - dj) -
           extra$mymat[, -1, drop = FALSE] * dj)
        ans2
      } else {  # Not "logitlink"
        c(w) * ddj.deta *  # Orig.
        (y[, -ncol(y), drop = FALSE] / dj -
         extra$mymat[, -1, drop = FALSE] / (1 - dj))
      }
    }
    ans
  }),
  list( .earg = earg, .link = link, .reverse = reverse) )),
  weight = eval(substitute(expression({
    if ( .reverse ) {
      cump <- tapplymat1(mu, "cumsum")
      wz <- if ( .link == "logitlink") {
        ans2 <- c(w) *  # New
                   (mu[, -1, drop = FALSE] * (1 - djr)^2 +
                 cump[, 1:M, drop = FALSE] * djr^2)
        ans2
      } else {  # Not "logitlink"
        c(w) * ddjr.deta^2 *
           (mu[, -1, drop = FALSE] / djr^2 +
         cump[, 1:M, drop = FALSE] / (1 - djr)^2)
      }
    } else {  # Not reversed
      ccump <- tapplymat1(mu[, ncol(mu):1, drop = FALSE],
                   "cumsum")[, ncol(mu):1, drop = FALSE]
      wz <- if ( .link == "logitlink") {
        ans2 <- c(w) *  # New
                 (mu[, 1:M, drop = FALSE] * (1 - dj)^2 +
                ccump[, -1, drop = FALSE] * dj^2)
        ans2
      } else {  # Not "logitlink"
        c(w) * ddj.deta^2 *
          (mu[, 1:M, drop = FALSE] / dj^2 +
         ccump[, -1, drop = FALSE] / (1 - dj)^2)
      }
    }

    wz
  }),
  list( .earg = earg, .link = link, .reverse = reverse ))))
}  # sratio()






 cratio <-
  function(link = "logitlink",
           parallel = FALSE, reverse = FALSE,
           zero = NULL,
           ynames = FALSE,
    Thresh = NULL,   # "free",
    Trev = reverse,
    Tref = if (Trev) "M" else 1,
    whitespace = FALSE) {
  link <- as.list(substitute(link))
  earg  <- link2list(link)
  link <- attr(earg, "function.name")

  if (length( Thresh ) &&
      !exists(paste0("CM.", Thresh),
              mode = "function"))
    stop("No function ", paste0("CM.", Thresh))

  if (!is.logical(reverse) || length(reverse) != 1)
    stop("argument 'reverse' not a single logical")

  stopifnot(is.logical(whitespace) &&
            length(whitespace) == 1)
  fillerChar <- ifelse(whitespace, " ", "")

  if (!is.logical(ynames) || length(ynames) != 1)
    stop("bad input for 'ynames'")


  new("vglmff",
  blurb = c("Continuation ratio model\n\n",
            "Links:    ",
            namesof(if (reverse)
            ifelse(whitespace, "P[Y < j+1|Y <= j+1]",
                               "P[Y<j+1|Y<=j+1]") else
            ifelse(whitespace, "P[Y > j|Y >= j]",
                               "P[Y>j|Y>=j]"),
                   link, earg = earg),
            "\n",
            "Variance: ",
            ifelse(whitespace,
                   "mu[,j] * (1 - mu[,j]); -mu[,j] * mu[,k]",
                   "mu[,j]*(1-mu[,j]); -mu[,j]*mu[,k]")),

  infos = eval(substitute(function(...) {
    list(M1 = NA,  # zz -1?
         Q1 = NA,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = as.character(NA),
         parallel = .parallel ,
         reverse = .reverse ,
         Thresh = .Thresh ,
         Tref = .Tref ,
         Trev = .Trev ,
         whitespace = .whitespace ,
         ynames = .ynames ,
         zero = .zero ,
         link = .link )
  },
  list( .link = link, .ynames = ynames,
        .zero = zero,
        .parallel = parallel,
        .reverse = reverse,
        .Thresh = Thresh,
        .Tref = Tref,
        .Trev = Trev,
        .whitespace = whitespace ))),


  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(1, M, 1), x = x,
                       bool = .parallel ,
                       constraints = constraints)
    constraints <- cm.zero.VGAM(constraints, x = x,
            .zero , M = M, M1 = M,
            predictors.names = predictors.names)


    if (names(constraints)[1] != "(Intercept)")
      stop("the model needs an intercept term")



    if (length( .Thresh ))
      constraints[["(Intercept)"]] <-
        do.call(paste0("CM.", .Thresh ),
                list(M = M, Trev = .Trev ,
                     Tref = .Tref ))
  }),
  list( .parallel = parallel,
        .Thresh = Thresh,
        .Tref = Tref,
        .Trev = Trev,
        .zero = zero ))),

  deviance = Deviance.categorical.data.vgam,

  initialize = eval(substitute(expression({

    if (is.factor(y) && !is.ordered(y))
      warning("response should be ordinal---",
              "see ordered()")


    delete.zero.colns <- TRUE
    eval(process.categorical.data.VGAM)
    M <- ncol(y) - 1




    if ( .ynames )  # overwrite mynames
      if (is.null(cnsy <- colnames(y)))
        stop("response matrix has no colnames")
    choice1 <- if ( .ynames ) cnsy[1:M] else
               1:M
    choice2 <- if ( .ynames ) cnsy[2:(M+1)] else
               2:(M+1)
    choice3 <- ifelse( .ynames , "P[", "P[Y")
    choice3 <- if ( .ynames ) paste0(choice3,
     .fillerChar ,
     if ( .reverse ) "<" else ">",
     .fillerChar ) else
      paste0(choice3,
             .fillerChar,
             if ( .reverse ) "<" else ">",
             .fillerChar )
    choice3 <-
      paste0(choice3,
             if ( .reverse ) choice2 else choice1,
             "|")
    choice3 <- if ( .ynames ) paste0(choice3,
     .fillerChar ,
      if ( .reverse ) "<=" else ">=",
     .fillerChar ,
      if ( .reverse ) choice2 else choice1) else
      paste0(choice3, "Y",
             .fillerChar,
             if ( .reverse ) "<=" else ">=",
             .fillerChar,
             if ( .reverse ) choice2 else choice1)
    mynames <- paste0(choice3, "]")


  if (F)  # orig.
    mynames <- if ( .reverse )
      paste0("P[Y", .fillerChar , "<",
             .fillerChar , choice2,
             "|Y",   .fillerChar , "<=",
             .fillerChar , choice2, "]") else
      paste0("P[Y", .fillerChar , ">",
             .fillerChar , choice1,
             "|Y",   .fillerChar , ">=",
             .fillerChar , choice1, "]")
    predictors.names <-
      namesof(mynames, .link , .earg , short = TRUE)
    y.names <-  param.names("mu", M+1)

    extra$mymat <- if ( .reverse )
       tapplymat1(y, "cumsum") else
       tapplymat1(y[, ncol(y):1, drop = FALSE],
                  "cumsum")[, ncol(y):1, drop = FALSE]

    extra$colnames.y  <- colnames(y)
  }),
  list( .earg = earg, .link = link,
        .reverse = reverse, .ynames = ynames,
        .fillerChar = fillerChar,
        .whitespace = whitespace ))),

  linkinv = eval(substitute( function(eta, extra = NULL) {
    if (!is.matrix(eta))
      eta <- as.matrix(eta)
    fv.mat <- if ( .reverse ) {
      M <- ncol(eta)
      djrs <- eta2theta(eta, .link , earg = .earg )
      temp <- tapplymat1(djrs[, M:1, drop = FALSE],
                         "cumprod")[, M:1, drop = FALSE]
      cbind(1, 1 - djrs) * cbind(temp, 1)
    } else {
      djs <- eta2theta(eta, .link , earg = .earg )
      temp <- tapplymat1(djs, "cumprod")
      cbind(1 - djs, 1) * cbind(1, temp)
    }
    label.cols.y(fv.mat, NOS = 1,
                 colnames.y = extra$colnames.y)
  },
  list( .earg = earg, .link = link, .reverse = reverse) )),
  last = eval(substitute(expression({

    misc$link <- rep_len( .link , M)
    names(misc$link) <- mynames

    misc$earg <- vector("list", M)
    names(misc$earg) <- names(misc$link)
    for (ii in 1:M)
      misc$earg[[ii]] <- .earg

    misc$parameters <- mynames
    misc$reverse <- ( .reverse )
    misc$fillerChar <- ( .fillerChar )
    misc$whitespace <- ( .whitespace )


    extra <- list()  # kill what was used
  }),
  list( .earg = earg, .reverse = reverse,
        .fillerChar = fillerChar, .link = link,
        .whitespace = whitespace ))),
  linkfun = eval(substitute( function(mu, extra = NULL) {
    cump <- tapplymat1(mu, "cumsum")
    if ( .reverse ) {
      djrs <- 1 - mu[, -1, drop = FALSE] / (
                cump[, -1, drop = FALSE])
      theta2eta(djrs, .link , earg = .earg )
    } else {
      M <- ncol(mu) - 1
      djs <- if (M == 1)
             1 - mu[, 1,   drop = FALSE] else
             1 - mu[, 1:M, drop = FALSE] / (
             1 - cbind(0, cump[, 1:(M-1), drop = FALSE]))
      theta2eta(djs, .link , earg = .earg )
    }
  },
  list( .earg = earg, .link = link, .reverse = reverse) )),
  loglikelihood =
    function(mu, y, w, residuals = FALSE, eta, extra = NULL,
             summation = TRUE) {
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ycounts <- if (is.numeric(extra$orig.w))
                 y * w / extra$orig.w else
                 y * w  # Convert proportions to counts
      nvec <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else round(w)

      smallno <- 1.0e4 * .Machine$double.eps
      if (max(abs(ycounts - round(ycounts))) > smallno)
        warning("converting 'ycounts' to integer ",
                "in @loglikelihood")
      ycounts <- round(ycounts)

      ll.elts <-
        (if (is.numeric(extra$orig.w)) extra$orig.w else 1) *
         dmultinomial(x = ycounts, size = nvec, prob = mu,
                      log = TRUE, dochecking = FALSE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  },
  vfamily = c("cratio", "VGAMordinal", "VGAMcategorical"),

  simslot = eval(substitute(
  function(object, nsim) {
    simslotVGAMcat(object, nsim)
  }, list(.link = link, .earg = earg ))),


  validparams = eval(substitute(
      function(eta, y, extra = NULL) {
    djrs <- eta2theta(eta, .link , .earg )  # djs or djrs
    okay1 <- all(is.finite(djrs)) && all(0 < djrs & djrs < 1)
    okay1
  }, list( .earg = earg, .link = link, .reverse = reverse) )),

  deriv = eval(substitute(expression({
    if (!length(extra$mymat)) {
      extra$mymat <-
        if ( .reverse )
          tapplymat1(y, "cumsum") else
          tapplymat1(y[, ncol(y):1, drop = FALSE],
                     "cumsum")[, ncol(y):1, drop = FALSE]
    }

    if ( .reverse ) {
      djrs <- eta2theta(eta, .link , earg = .earg )
      ddjrs.deta <- dtheta.deta(djrs, .link , earg = .earg )
      Mp1 <- ncol(extra$mymat)
      if ( .link == "logitlink") {
        ans2 <- -c(w) *  # New
                     (y[, -1, drop = FALSE] * djrs -
          extra$mymat[, -Mp1, drop = FALSE] * (1 - djrs))
        ans2
      } else {  # Not the default link
        -c(w) * ddjrs.deta *  # Orig.
          (y[, -1, drop = FALSE] / (1 - djrs) -
          extra$mymat[, -Mp1, drop = FALSE] / djrs)
      }
    } else {  # Not reverse
      djs <- eta2theta(eta, .link , earg = .earg )
      ddjs.deta <- dtheta.deta(djs, .link , earg = .earg )
      if ( .link == "logitlink") {
        ans2 <- -c(w) *   # New
             (y[, -ncol(y), drop = FALSE] * djs -
          extra$mymat[, -1, drop = FALSE] * (1 - djs))
        ans2
      } else {  # Not "logitlink"
        -c(w) * (y[, -ncol(y), drop = FALSE] / (1 - djs) -
        extra$mymat[, -1, drop = FALSE] / djs) *
        ddjs.deta
      }
    }
  }),
  list( .earg = earg, .link = link, .reverse = reverse) )),
  weight = eval(substitute(expression({
    if ( .reverse ) {
      cump <- tapplymat1(mu, "cumsum")
      wz <- if ( .link == "logitlink") {
        ans2 <- c(w) *  # New
               (mu[, -1, drop = FALSE] * djrs^2 +
             cump[, 1:M, drop = FALSE] * (1 - djrs)^2)
        ans2
      } else {  # Not "logitlink"
        c(w) * ddjrs.deta^2 *
          (mu[, -1, drop = FALSE] / (1 - djrs)^2 +
        cump[, 1:M, drop = FALSE] / djrs^2)
      }
    } else {  # Not reversed
      ccump <- tapplymat1(mu[, ncol(mu):1, drop = FALSE],
                   "cumsum")[, ncol(mu):1, drop = FALSE]
      wz <- if ( .link == "logitlink") {
        ans2 <- c(w) *  # New
           (mu[, 1:M, drop = FALSE] * djs^2 +
          ccump[, -1, drop = FALSE] * (1 - djs)^2)
        ans2
      } else {  # Not "logitlink"
        c(w) * ddjs.deta^2 *
         (mu[, 1:M, drop = FALSE] / (1 - djs)^2 +
        ccump[, -1, drop = FALSE] / djs^2)
      }
    }
    wz
  }),
  list( .earg = earg, .link = link, .reverse = reverse ))))
}  # cratio()






 vglm.multinomial.deviance.control <-
  function(maxit = 31, panic = FALSE, ...) {
  if (maxit < 1) {
      warning("bad value of maxit; using 31 instead")
      maxit <- 31
  }
  list(maxit = maxit, panic = as.logical(panic)[1])
}


 vglm.multinomial.control <-
  function(maxit = 21, panic = FALSE,
           criterion = c("aic1", "aic2",
                         names( .min.criterion.VGAM )),
           ...) {
  if (mode(criterion) != "character" &&
      mode(criterion) != "name")
    criterion <- as.character(substitute(criterion))
  criterion <- match.arg(criterion,
      c("aic1", "aic2", names( .min.criterion.VGAM )))[1]

  if (maxit < 1) {
    warning("bad value of maxit; using 21 instead")
    maxit <- 21
  }
  list(maxit = maxit,
       panic = as.logical(panic)[1],
       criterion = criterion,
       min.criterion = c("aic1" = FALSE, "aic2" = TRUE,
                         .min.criterion.VGAM))
}


 vglm.VGAMcategorical.control <-
  function(maxit = 30,
           trace = FALSE,
           panic = TRUE, ...) {
  if (maxit < 1) {
    warning("bad value of maxit; using 200 instead")
    maxit <- 200
  }
  list(maxit = maxit,
       trace = as.logical(trace)[1],
       panic = as.logical(panic)[1])
}  # vglm.VGAMcategorical.control()








 cumulative <-
  function(link = "logitlink",
    parallel = FALSE,  # Doesnt apply 2 intercept
    reverse = FALSE,
    multiple.responses = FALSE,
    ynames = FALSE,  # 20240216
    Thresh = NULL,  # "free",
    Trev = reverse,
    Tref = if (Trev) "M" else 1,
    whitespace = FALSE) {

  if (length( Thresh ) &&
      !exists(paste0("CM.", Thresh),
              mode = "function"))
    stop("No function ", paste0("CM.", Thresh))
  if (multiple.responses && !(Thresh == "free" ||
      is.null(Thresh)))
    stop("'Thresh' must be NULL or 'free'")

  if (!is.logical(ynames) || length(ynames) != 1)
    stop("bad input for 'ynames'")

  apply.parint <- FALSE


  link <- as.list(substitute(link))
  earg  <- link2list(link)
  link <- attr(earg, "function.name")



  stopifnot(is.logical(whitespace) &&
            length(whitespace) == 1)
  fillerChar <- ifelse(whitespace, " ", "")


  if (!is.logical(multiple.responses) ||
      length(multiple.responses) != 1)
    stop("argument 'multiple.responses' must be ",
         "a single logical")
  if (!is.logical(reverse) || length(reverse) != 1)
    stop("'reverse' must be a single logical")


  new("vglmff",
  blurb = if ( multiple.responses )
          c(paste("Multivariate cumulative", link, "model\n\n"),
          "Links:   ",
          namesof(if (reverse)
          ifelse(whitespace, "P[Y1 >= j+1]", "P[Y1>=j+1]") else
          ifelse(whitespace, "P[Y1 <= j]",   "P[Y1<=j]"),
          link, earg = earg),
          ", ...") else
          c(paste("Cumulative", link, "model\n\n"),
          "Links:   ",
          namesof(if (reverse)
          ifelse(whitespace, "P[Y >= j+1]", "P[Y>=j+1]") else
          ifelse(whitespace, "P[Y <= j]",   "P[Y<=j]"),
          link, earg = earg)),
  infos = eval(substitute(function(...) {
    list(M1 = NA,  # zz -1?
         Q1 = NA,
         expected = TRUE,
         hadof = TRUE,
         multipleResponses = .multiple.responses ,
         parameters.names = as.character(NA),
         parallel = .parallel ,
         reverse = .reverse ,
         Thresh = .Thresh ,
         Tref = .Tref ,
         Trev = .Trev ,
         whitespace = .whitespace ,
         ynames = .ynames ,
         link =  .link )
  },
  list( .link = link, .ynames = ynames,
        .parallel = parallel,
        .multiple.responses = multiple.responses,
        .reverse = reverse,
        .Thresh = Thresh,
        .Tref = Tref,
        .Trev = Trev,
        .whitespace = whitespace ))),

  constraints = eval(substitute(expression({
    if ( .multiple.responses ) {
      if ( !length(constraints) ) {
        Llevels <- extra$Llevels
        NOS <- extra$NOS
        Hk.matrix <- kronecker(diag(NOS),
                     matrix(1, Llevels - 1, 1))
      constraints <- cm.VGAM(Hk.matrix, x = x,
                       bool = .parallel ,
                       apply.int = .apply.parint ,
                       constraints = constraints)
      }
    } else {
    constraints <- cm.VGAM(matrix(1, M, 1), x = x,
                       bool = .parallel ,
                       apply.int = .apply.parint ,
                       constraints = constraints)
    }


    if (names(constraints)[1] != "(Intercept)")
      stop("the model needs an intercept term")


    if (length( .Thresh ))
      constraints[["(Intercept)"]] <-
        do.call(paste0("CM.", .Thresh ),
                list(M = M, Trev = .Trev ,
                     Tref = .Tref ))
  }),
  list( .parallel = parallel,
        .multiple.responses = multiple.responses,
        .Thresh = Thresh,
        .Tref = Tref,
        .Trev = Trev,
        .apply.parint = apply.parint ))),
  deviance = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta, extra = NULL) {

    answer <-
    if ( .multiple.responses ) {
      totdev <- 0
      NOS <- extra$NOS
      Llevels <- extra$Llevels
      for (iii in 1:NOS) {
        cindex <- (iii-1)*(Llevels-1) + 1:(Llevels-1)
        aindex <- (iii-1)*(Llevels  ) + 1:(Llevels)
        totdev <- totdev +
                  Deviance.categorical.data.vgam(
                    mu = mu[, aindex, drop = FALSE],
                    y = y[, aindex, drop = FALSE], w = w,
                    residuals = residuals,
                    eta = eta[, cindex, drop = FALSE],
                    extra = extra,
                    summation = TRUE)
      }
      totdev
    } else {
      Deviance.categorical.data.vgam(mu = mu,
          y = y, w = w, residuals = residuals,
          eta = eta, extra = extra,
          summation = TRUE)
    }
    answer
  },
  list( .earg = earg, .link = link,
        .multiple.responses = multiple.responses
       ) )),
  initialize = eval(substitute(expression({

    if (colnames(x)[1] != "(Intercept)")
      warning("there seems to be no intercept term!")


    if (is.factor(y) && !is.ordered(y))
      warning("response should be ordinal---see ordered()")


    extra$multiple.responses <- .multiple.responses
    if ( .multiple.responses ) {
      if ( .ynames )
          stop("multiple.responses = T: cannot ",
               "have 'ynames = T' too")
      checkCut(y)  # Check input; stops if an error
      if (any(w != 1) || NCOL(w) != 1)
        stop("the 'weights' argument must be a ",
             "vector of all 1s")
      Llevels <- max(y)
      delete.zero.colns <- FALSE
      orig.y <- cbind(y)  # Convert y into a matrix if necsry
      NOS <- NCOL(orig.y)
      use.y <- use.mustart <- NULL
      for (iii in 1:NOS) {
        y <- as.factor(orig.y[,iii])
        eval(process.categorical.data.VGAM)
        use.y <- cbind(use.y, y)
        use.mustart <- cbind(use.mustart, mustart)
      }
      mustart <- use.mustart
      y <- use.y  # n x (Llevels*NOS)
      M <- NOS * (Llevels - 1)
      mynames <- y.names <- NULL
      for (iii in 1:NOS) {
        Y.names <- paste0("Y", iii)
        mu.names <- paste0("mu", iii, ".")
        mynames <- c(mynames, if ( .reverse )
          paste0("P[", Y.names, ">=",
                 2:Llevels, "]") else
          paste0("P[", Y.names, "<=",
                 1:(Llevels - 1), "]"))
        y.names <- c(y.names,
                     param.names(mu.names,
                                 Llevels))
      }

      predictors.names <-
        namesof(mynames, .link , short = TRUE, earg = .earg )

      extra$NOS <- NOS
      extra$Llevels <- Llevels
  } else {


      delete.zero.colns <- TRUE

      eval(process.categorical.data.VGAM)
      M <- ncol(y) - 1



      if ( .ynames )  # overwrite mynames
        if (is.null(cnsy <- colnames(y)))
          stop("response matrix has no colnames")
      choice1 <- if ( .ynames ) cnsy[1:M] else
                 1:M
      choice2 <- if ( .ynames ) cnsy[2:(M+1)] else
               2:(M+1)
      choice3 <- ifelse( .ynames , "P[", "P[Y")
      choice4 <- ifelse( .ynames , "", "=")

      mynames <- if ( .reverse )
        paste0(choice3,   # "P[Y",
               .fillerChar ,
               ">=",
               .fillerChar, choice2, "]") else
        paste0(choice3,  # "P[Y",
               .fillerChar ,
               "<=",
               .fillerChar, choice1,  "]")
      

      predictors.names <-
        namesof(mynames, .link , short = TRUE,
                earg = .earg )

      if (NCOL(w) == 1) {
          if (length(mustart) && all(c(y) %in% c(0, 1)))
            for (iii in 1:ncol(y))
                mustart[,iii] <- weighted.mean(y[,iii], w)
      }

    extra$colnames.y  <- colnames(y)
  }
  }),
  list( .reverse = reverse, .ynames = ynames,
        .multiple.responses = multiple.responses,
        .link = link, .earg = earg,
        .fillerChar = fillerChar,
        .whitespace = whitespace ))),


    linkinv = eval(substitute( function(eta, extra = NULL) {
        answer <-
        if ( .multiple.responses ) {
          NOS <- extra$NOS
          Llevels <- extra$Llevels
          fv.mat <- matrix(0, nrow(eta), NOS * Llevels)
          for (iii in 1:NOS) {
            cindex <- (iii-1)*(Llevels-1) + 1:(Llevels-1)
            aindex <- (iii-1)*(Llevels) + 1:(Llevels)
            if ( .reverse ) {
              ccump <- cbind(1,
                      eta2theta(eta[, cindex, drop = FALSE],
                                .link , earg = .earg ))
              fv.mat[, aindex] <-
                  cbind(-tapplymat1(ccump, "diff"),
                        ccump[, ncol(ccump)])
            } else {
        cump <- cbind(eta2theta(eta[, cindex, drop = FALSE],
                                      .link ,
                                      earg = .earg ),
                            1)
              fv.mat[, aindex] <-
                  cbind(cump[, 1], tapplymat1(cump, "diff"))
            }
          }
          label.cols.y(fv.mat, NOS = NOS,
                colnames.y = if (is.null(extra$colnames.y))
                                    NULL else
                rep_len(extra$colnames.y, ncol(fv.mat)))
        } else {
          fv.mat <- if ( .reverse ) {
                    ccump <- cbind(1, eta2theta(eta, .link ,
                                         earg = .earg ))
                    cbind(-tapplymat1(ccump, "diff"),
                          ccump[, ncol(ccump)])
          } else {
              cump <- cbind(eta2theta(eta, .link ,
                                      earg = .earg ), 1)
            cbind(cump[, 1], tapplymat1(cump, "diff"))
          }
          label.cols.y(fv.mat, NOS = 1,
                  colnames.y = extra$colnames.y)
        }
        answer
  },
  list( .reverse = reverse,
        .multiple.responses = multiple.responses,
        .link = link, .earg = earg))),

  last = eval(substitute(expression({
    if ( .multiple.responses ) {
      misc$link <- .link
      misc$earg <- list( .earg )

    } else {
      misc$link <- rep_len( .link , M)
      names(misc$link) <- mynames

      misc$earg <- vector("list", M)
      names(misc$earg) <- names(misc$link)
      for (ii in 1:M)
        misc$earg[[ii]] <- .earg

    }

    misc$fillerChar <- .fillerChar
    misc$whitespace <- .whitespace

    misc$parameters <- mynames
    misc$reverse <- .reverse
    misc$parallel <- .parallel
    misc$multiple.responses <- .multiple.responses
  }),
  list( .reverse = reverse, .parallel = parallel,
        .link = link, .earg = earg,
        .fillerChar = fillerChar,
        .multiple.responses = multiple.responses,
        .whitespace = whitespace ))),

  linkfun = eval(substitute( function(mu, extra = NULL) {
    answer <-
    if ( .multiple.responses ) {
      NOS <- extra$NOS
      Llevels <- extra$Llevels
      eta.matrix <- matrix(0, nrow(mu), NOS*(Llevels-1))
      for (iii in 1:NOS) {
        cindex <- (iii-1)*(Llevels-1) + 1:(Llevels-1)
        aindex <- (iii-1)*(Llevels) + 1:(Llevels)
        cump <- tapplymat1(as.matrix(mu[, aindex]), "cumsum")
        eta.matrix[, cindex] <-
            theta2eta(if ( .reverse )
                          1-cump[, 1:(Llevels-1)] else
                  cump[, 1:(Llevels-1)], .link , earg = .earg )
      }
      eta.matrix
    } else {
      cump <- tapplymat1(as.matrix(mu), "cumsum")
      M <- NCOL(mu) - 1
      theta2eta(if ( .reverse ) 1-cump[, 1:M] else cump[, 1:M],
                .link , earg = .earg )
    }
  if (nrow(mu) == 1)
    answer <- rbind(answer)
    answer
  },
 list( .link = link, .earg = earg,
       .multiple.responses = multiple.responses,
       .reverse = reverse ))),
  loglikelihood =
    function(mu, y, w, residuals = FALSE, eta, extra = NULL,
             summation = TRUE) {
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ycounts <- if (is.numeric(extra$orig.w))
                 y * w / extra$orig.w else
                 y * w # Convert proportions to counts
      nvec <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else
              round(w)

      smallno <- 1.0e4 * .Machine$double.eps
      if (max(abs(ycounts - round(ycounts))) > smallno)
        warning("converting 'ycounts' to integer ",
                "in @loglikelihood")
      ycounts <- round(ycounts)

      ll.elts <-
        (if (is.numeric(extra$orig.w)) extra$orig.w else 1) *
         dmultinomial(x = ycounts, size = nvec, prob = mu,
                      log = TRUE, dochecking = FALSE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  },
  vfamily = c("cumulative", "VGAMordinal", "VGAMcategorical"),


  simslot = eval(substitute(
  function(object, nsim) {
    simslotVGAMcat(object, nsim)
  }, list(.link = link, .earg = earg ))),


  hadof = eval(substitute(
  function(eta, extra = list(),
           linpred.index = 1, w = 1,
           dim.wz = c(NROW(eta), NCOL(eta) * (NCOL(eta)+1)/2), 
           deriv = 1, ...) {
  if ( .multiple.responses )
    return(NA_real_, dim.wz[1], dim.wz[2])



  cumulative.eim.deriv1 <-
    function(mu, jay, w, reverse = FALSE) {
    M <- ncol(mu) - 1
    wz1 <- matrix(0, NROW(mu), M + M-1)  # Tridiagonal
    wz1[, iam(jay, jay, M = M)] <-  1 / (mu[, jay+1])^2 -
                                    1 / (mu[, jay  ])^2
    if (1 <= jay-1)
      wz1[, iam(jay-1, jay-1, M = M)] <-  -1 / (mu[, jay  ])^2
    if (jay+1 <= M)
      wz1[, iam(jay+1, jay+1, M = M)] <-   1 / (mu[, jay+1])^2
    if (1 < M && jay+1 <= M)
      wz1[, iam(jay, jay+1, M = M)] <-    -1 / (mu[, jay+1])^2
    if (1 < M && 1 <= jay-1)
      wz1[, iam(jay-1, jay, M = M)] <-     1 / (mu[, jay  ])^2
    (if (reverse) -c(w) else c(w)) * wz1
  }  # cumulative.eim.deriv1



  cumulative.eim.deriv2 <- function(mu, jay, w) {
    M <- ncol(mu) - 1
    wz2 <- matrix(0, NROW(mu), M + M-1)  # Tridiagonal
    wz2[, iam(jay, jay, M = M)] <-   1 / (mu[, jay+1])^3 +
                                     1 / (mu[, jay  ])^3
    if (1 <= jay-1)
      wz2[, iam(jay-1, jay-1, M = M)] <-  1 / (mu[, jay  ])^3
    if (jay+1 <= M)
      wz2[, iam(jay+1, jay+1, M = M)] <-  1 / (mu[, jay+1])^3
    if (1 < M && jay+1 <= M)
      wz2[, iam(jay, jay+1, M = M)] <-   -1 / (mu[, jay+1])^3
    if (1 < M && 1 <= jay-1)
      wz2[, iam(jay-1, jay, M = M)] <-   -1 / (mu[, jay  ])^3
    2 * c(w) * wz2
  }  # cumulative.eim.deriv2



    probs <- if ( .reverse ) {
        ccump <- cbind(1, eta2theta(eta, .link , earg = .earg ))
        cbind(-tapplymat1(ccump, "diff"), ccump[, ncol(ccump)])
      } else {
        cump <- cbind(eta2theta(eta, .link , earg = .earg ), 1)
        cbind(cump[, 1], tapplymat1(cump, "diff"))
      }
    mu.use <- pmax(probs, .Machine$double.eps * 1.0e-0)

    switch(as.character(deriv),
       "0" =  {
         M <- ncol(eta)
         wz <- c(w) * (1 / mu.use[, 1:M] + 1 / mu.use[, -1])
         if (M > 1)
           wz <- cbind(wz, -c(w) / mu.use[, 2:M])
         wz
       },
       "1" = {
           cumulative.eim.deriv1(mu = mu.use,
                                 jay = linpred.index, w = w,
                                 reverse = .reverse )
       },
       "2" = {
           cumulative.eim.deriv2(mu = mu.use,
                                 jay = linpred.index, w = w)
       },
       stop("argument 'deriv' must be 0 or 1 or 2"))
  },
  list( .link = link, .earg = earg,
        .multiple.responses = multiple.responses,
        .reverse = reverse ))),




  validparams = eval(substitute(function(eta, y, extra = NULL) {
    if ( .multiple.responses ) {
      return(TRUE)
    }
    probs <-
      if ( .reverse ) {
        ccump <- cbind(1, eta2theta(eta, .link , .earg ))
        cbind(-tapplymat1(ccump, "diff"),
              ccump[, ncol(ccump)])
      } else {
        cump <- cbind(eta2theta(eta, .link , .earg ), 1)
        cbind(cump[, 1], tapplymat1(cump, "diff"))
      }
    okay1 <- all(is.finite(probs)) &&
             all(0 < probs & probs < 1)
    if (!okay1)
      warning("It seems that the nonparallelism assumption ",
    "has resulted in intersecting linear/additive ",
    "predictors. Try propodds() or fitting a partial ",
    "nonproportional odds model or choosing ",
    "some other link function, etc.")
    okay1
  },
  list( .link = link, .earg = earg,
        .multiple.responses = multiple.responses,
        .reverse = reverse ))),

  deriv = eval(substitute(expression({
    prob.lo <- .Machine$double.eps  * 1.0e-0
    prob.hi <- 1 - prob.lo
    mu.use <- pmax(mu, prob.lo)
    mu.use <- pmin(mu.use, prob.hi)  # 20230701
    range.mu <- range(mu)
    if (range.mu[1] < prob.lo)
      warning("some probabilities are very close to 0")
    if (range.mu[2] > prob.hi)
      warning("some probabilities are very close to 1")

    deriv.answer <-
    if ( .multiple.responses ) {
      NOS <- extra$NOS
      Llevels <- extra$Llevels
      dcump.deta <- resmat <- matrix(0, n, NOS * (Llevels-1))
      for (iii in 1:NOS) {
        cindex <- (iii-1)*(Llevels-1) + 1:(Llevels-1)
        aindex <- (iii-1)*(Llevels)   + 1:(Llevels-1)
        cump <- eta2theta(eta[,cindex, drop = FALSE],
                         .link , earg = .earg )
        dcump.deta[,cindex] <- dtheta.deta(cump, .link ,
                                           earg = .earg )
        resmat[,cindex] <-
            (y[,   aindex, drop = FALSE] / mu.use[,   aindex,
                                             drop = FALSE] -
             y[, 1+aindex, drop = FALSE] / mu.use[, 1+aindex,
                                             drop = FALSE])
      }
      (if ( .reverse ) -c(w)  else c(w)) * dcump.deta * resmat
    } else {
      cump <- eta2theta(eta    , .link , earg = .earg )
      dcump.deta <- dtheta.deta(cump, .link , earg = .earg )
      c(if ( .reverse ) -c(w) else c(w)) *
      dcump.deta *
      (     y[, -(M+1), drop = FALSE] /
       mu.use[, -(M+1), drop = FALSE] -
            y[, -1, drop = FALSE] /
       mu.use[, -1, drop = FALSE])
    }
    deriv.answer
  }),
  list( .link = link, .earg = earg,
        .multiple.responses = multiple.responses,
        .reverse = reverse ))),
  weight = eval(substitute(expression({
    if ( .multiple.responses ) {
      NOS <- extra$NOS
      Llevels <- extra$Llevels
      wz <- matrix(0, n,
           NOS*(Llevels-1))  # Diag elts only for a start
      for (iii in 1:NOS) {
        cindex <- (iii-1)*(Llevels-1) + 1:(Llevels-1)
        aindex <- (iii-1)*(Llevels)   + 1:(Llevels-1)
        wz[, cindex] <-
            c(w) * dcump.deta[, cindex, drop = FALSE]^2 *
            (1 / mu.use[,   aindex, drop = FALSE] +
             1 / mu.use[, 1+aindex, drop = FALSE])
      }
      if (Llevels - 1 > 1) {
        iii <- 1
        oindex <- (iii-1) * (Llevels-1) + 1:(Llevels-2)
        wz <- cbind(wz, -c(w) *
              dcump.deta[, oindex] * dcump.deta[, 1+oindex])


        if (NOS > 1) {
          cptrwz <- ncol(wz)  # Like a pointer
          wz <- cbind(wz,
                      matrix(0, nrow(wz),
                            (NOS-1) * (Llevels-1)))
          for (iii in 2:NOS) {
            oindex <- (iii-1) * (Llevels-1) + 1:(Llevels-2)
            wz[,cptrwz + 1 + (1:(Llevels-2))] <-
                  -c(w) * dcump.deta[, oindex] *
                       dcump.deta[, 1 + oindex]
     cptrwz <- cptrwz + Llevels - 1  # Move it along a bit
          }
        }



        }
    } else {  # Not .multiple.responses here
      wz <- c(w) * dcump.deta^2 *
            (1 / mu.use[, 1:M, drop = FALSE] +
             1 / mu.use[,  -1, drop = FALSE])
      if (M > 1)
        wz <- cbind(wz,
                    -c(w) *
                    dcump.deta[,  -M, drop = FALSE] *
                    dcump.deta[, 2:M, drop = FALSE] / (
                        mu.use[, 2:M, drop = FALSE]))
    }  # End of not .multiple.responses here
    wz
  }),
  list( .multiple.responses = multiple.responses,
        .earg = earg, .link = link ))))
}  # cumulative()





 propodds <-
     function(reverse = TRUE,
              whitespace = FALSE,
              ynames = FALSE,
       Thresh = NULL,  # "free",
       Trev = reverse,
       Tref = if (Trev) "M" else 1) {
  if (!is.logical(reverse) ||
      length(reverse) != 1)
    stop("'reverse' not a single logical")
  if (length( Thresh ) &&
      !exists(paste0("CM.", Thresh),
              mode = "function"))
    stop("No function ", paste0("CM.", Thresh))

  cumulative(parallel = TRUE, reverse = reverse,
             ynames = ynames,
             Thresh = Thresh,
             Trev = Trev, Tref = Tref,
             whitespace = whitespace)
}





 acat <-
  function(link = "loglink", parallel = FALSE,
           reverse = FALSE, zero = NULL,
           ynames = FALSE,  # 20240216
       Thresh = NULL,  # "free",
       Trev = reverse,
       Tref = if (Trev) "M" else 1,
       whitespace = FALSE) {


  link <- as.list(substitute(link))
  earg  <- link2list(link)
  link <- attr(earg, "function.name")

  if (length( Thresh ) &&
      !exists(paste0("CM.", Thresh),
              mode = "function"))
    stop("No function ", paste0("CM.", Thresh))

  if (!is.logical(reverse) || length(reverse) != 1)
    stop("'reverse' not a single logical")

  stopifnot(is.logical(whitespace) &&
            length(whitespace) == 1)
  fillerChar <- ifelse(whitespace, " ", "")

  if (!is.logical(ynames) || length(ynames) != 1)
    stop("bad input for 'ynames'")

  new("vglmff",
  blurb = c("Adjacent-categories model\n\n",
            "Links:    ",
    namesof(if (reverse)
                ifelse(whitespace, "P[Y = j] / P[Y = j + 1]",
                       "P[Y=j]/P[Y=j+1]") else
                ifelse(whitespace, "P[Y = j + 1] / P[Y = j]",
                                   "P[Y=j+1]/P[Y=j]"),
            link, earg = earg), "\n",
            "Variance: ",
            ifelse(whitespace,
            "mu[,j] * (1 - mu[,j]); -mu[,j] * mu[,k]",
            "mu[,j]*(1-mu[,j]); -mu[,j]*mu[,k]")),
  infos = eval(substitute(function(...) {
    list(M1 = NA,  # zz -1?
         Q1 = NA,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = as.character(NA),
         parallel = .parallel ,
         reverse = .reverse ,
         Thresh = .Thresh ,
         Tref = .Tref ,
         Trev = .Trev ,
         whitespace = .whitespace ,
         ynames = .ynames ,
         zero = .zero ,
         link = .link )
  },
  list( .link = link,
        .zero = zero, .ynames = ynames,
        .parallel = parallel,
        .reverse = reverse,
        .Thresh = Thresh,
        .Tref = Tref,
        .Trev = Trev,
        .whitespace = whitespace ))),

  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(1, M, 1), x = x,
                       bool = .parallel ,
                       constraints = constraints)
    constraints <- cm.zero.VGAM(constraints, x = x,
              .zero , M = M, M1 = M,
              predictors.names = predictors.names)


    if (names(constraints)[1] != "(Intercept)")
      stop("the model needs an intercept term")


    if (length( .Thresh ))
      constraints[["(Intercept)"]] <-
        do.call(paste0("CM.", .Thresh ),
                list(M = M, Trev = .Trev ,
                     Tref = .Tref ))
  }),
  list( .parallel = parallel,
        .Thresh = Thresh,
        .Tref = Tref,
        .Trev = Trev,
        .zero = zero ))),

  deviance = Deviance.categorical.data.vgam,

  initialize = eval(substitute(expression({

    if (is.factor(y) && !is.ordered(y))
      warning("response should be ordinal---see ordered()")



    delete.zero.colns <- TRUE
    eval(process.categorical.data.VGAM)
    M <- ncol(y) - 1


    if ( .ynames )   # overwrite mynames
      if (is.null(cnsy <- colnames(y)))
        stop("response matrix has no colnames")
    choice1 <- if ( .ynames ) cnsy[1:M] else
               1:M
    choice2 <- if ( .ynames ) cnsy[2:(M+1)] else
               2:(M+1)
    choice3 <- ifelse( .ynames , "P[", "P[Y")
    choice4 <- ifelse( .ynames , "", "=")
    mynames <- if ( .reverse )
      paste0(choice3,  # .fillerChar ,
             choice4, .fillerChar ,
             choice1, .fillerChar , "]",
             .fillerChar , "/", .fillerChar ,
             choice3, .fillerChar ,
             choice4, .fillerChar ,
             choice2, .fillerChar , "]") else
      paste0(choice3,  # .fillerChar ,
             choice4, .fillerChar ,
             choice2, .fillerChar , "]",
             .fillerChar , "/", .fillerChar ,
             choice3,  # .fillerChar ,
             choice4, .fillerChar ,
             choice1, .fillerChar ,  "]")

    predictors.names <-
      namesof(mynames, .link , short = TRUE, earg = .earg )
    y.names <- param.names("mu", M+1)

    extra$colnames.y  <- colnames(y)
  }),
  list( .earg = earg, .link = link,
        .reverse = reverse, .ynames = ynames,
        .fillerChar = fillerChar,
        .whitespace = whitespace ))),

  linkinv = eval(substitute( function(eta, extra = NULL) {
    if (!is.matrix(eta))
      eta <- as.matrix(eta)
    M <- ncol(eta)
    fv.mat <- if ( .reverse ) {
      zetar <- eta2theta(eta, .link , earg = .earg )
      temp <- tapplymat1(zetar[, M:1],
                         "cumprod")[, M:1, drop = FALSE]
      cbind(temp, 1) / drop(1 + temp %*% rep(1, ncol(temp)))
    } else {
      zeta <- eta2theta(eta, .link , earg = .earg )
      temp <- tapplymat1(zeta, "cumprod")
      cbind(1, temp) / drop(1 + temp %*% rep(1, ncol(temp)))
    }
    label.cols.y(fv.mat, NOS = 1,
                 colnames.y = extra$colnames.y)
  },
  list( .earg = earg, .link = link,
        .reverse = reverse) )),

  last = eval(substitute(expression({
    misc$link <- rep_len( .link , M)
    names(misc$link) <- mynames

    misc$earg <- vector("list", M)
    names(misc$earg) <- names(misc$link)
    for (ii in 1:M)
      misc$earg[[ii]] <- .earg

    misc$parameters <- mynames
    misc$reverse <- .reverse
    misc$fillerChar <- .fillerChar
    misc$whitespace <- .whitespace
  }),
  list( .earg = earg, .link = link,
        .reverse = reverse,
        .fillerChar = fillerChar,
        .whitespace = whitespace ))),
  linkfun = eval(substitute( function(mu, extra = NULL) {
    M <- ncol(mu) - 1
    theta2eta(if ( .reverse )
        mu[, 1:M, drop = FALSE] / mu[,  -1, drop = FALSE] else
        mu[,  -1, drop = FALSE] / mu[, 1:M, drop = FALSE],
        .link , earg = .earg )
  },
  list( .earg = earg, .link = link,
        .reverse = reverse) )),
  loglikelihood =
    function(mu, y, w, residuals = FALSE, eta, extra = NULL,
             summation = TRUE) {
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ycounts <- if (is.numeric(extra$orig.w))
                 y * w / extra$orig.w else
                 y * w  # Convert proportions to counts
      nvec <- if (is.numeric(extra$orig.w))
                  round(w / extra$orig.w) else round(w)

      smallno <- 1.0e4 * .Machine$double.eps
      if (max(abs(ycounts - round(ycounts))) > smallno)
        warning("converting 'ycounts' to ",
                "integer in @loglikelihood")
      ycounts <- round(ycounts)

      ll.elts <-
        (if (is.numeric(extra$orig.w)) extra$orig.w else 1) *
         dmultinomial(x = ycounts, size = nvec, prob = mu,
                      log = TRUE, dochecking = FALSE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  },
  vfamily = c("acat", "VGAMordinal", "VGAMcategorical"),


  simslot = eval(substitute(
  function(object, nsim) {
    simslotVGAMcat(object, nsim)
  }, list(.link = link, .earg = earg ))),






  validparams = eval(substitute(function(eta, y, extra = NULL) {
    if (!is.matrix(eta))
      eta <- as.matrix(eta)
    M <- ncol(eta)
    probs <- if ( .reverse ) {
      zetar <- eta2theta(eta, .link , earg = .earg )
      temp <- tapplymat1(zetar[, M:1],
                         "cumprod")[, M:1, drop = FALSE]
      cbind(temp, 1) / drop(1 + temp %*% rep(1, ncol(temp)))
    } else {
      zeta <- eta2theta(eta, .link , earg = .earg )
      temp <- tapplymat1(zeta, "cumprod")
      cbind(1, temp) / drop(1 + temp %*% rep(1, ncol(temp)))
    }
    okay1 <- all(is.finite(probs)) && all(0 < probs & probs < 1)
    okay1
  },
  list( .earg = earg, .link = link,
        .reverse = reverse) )),
  deriv = eval(substitute(expression({
    zeta <- eta2theta(eta, .link , .earg )  # May be zetar
    dzeta.deta <- dtheta.deta(zeta, .link , earg = .earg )


  
    
    d1 <- acat.deriv(zeta, M = M, n = n, reverse = .reverse )
    score <- attr(d1, "gradient") / d1
 

    answer <- if ( .reverse ) {
      cumy <- tapplymat1(y, "cumsum")
      c(w) * dzeta.deta * (cumy[, 1:M] / zeta - score)
    } else {
      ccumy <- tapplymat1(y[, ncol(y):1, drop = FALSE],
                          "cumsum")[, ncol(y):1, drop = FALSE]
      c(w) * dzeta.deta * (ccumy[, -1, drop = FALSE] / zeta -
             score)
    }


    answer
  }),
  list( .earg = earg, .link = link,
        .reverse = reverse) )),
  weight = eval(substitute(expression({
    wz <- matrix(NA_real_, n, dimm(M))

    hess <- attr(d1, "hessian") / d1

    if (M > 1)
      for (jay in 1:(M - 1))
        for (kay in (jay + 1):M)
          wz[, iam(jay, kay, M)] <-
            (hess[, jay, kay] -
            score[, jay] * score[, kay]) *
            dzeta.deta[, jay] * dzeta.deta[, kay]
    if ( .reverse ) {
      cump <- tapplymat1(mu, "cumsum")
      wz[, 1:M] <- (cump[, 1:M, drop = FALSE] / zeta^2 -
                    score^2) * dzeta.deta^2
    } else {
      ccump <- tapplymat1(mu[, ncol(mu):1, drop = FALSE],
                          "cumsum")[, ncol(mu):1, drop = FALSE]
      wz[, 1:M] <- (ccump[, -1, drop = FALSE] / zeta^2 -
                    score^2) * dzeta.deta^2
    }
    c(w) * wz
  }),
  list( .earg = earg, .link = link,
        .reverse = reverse ))))
}  # acat()




acat.deriv <- function(zeta, reverse, M, n) {

  alltxt <- NULL
  for (ii in 1:M) {
    index <- if (reverse) ii:M else 1:ii
    vars <- paste0("zeta", index)
    txt <- paste(vars, collapse = "*")
    alltxt <- c(alltxt, txt)
  }
  alltxt <- paste(alltxt, collapse = " + ")
  alltxt <- paste(" ~ 1 +", alltxt)
  txt.f <- as.formula(alltxt)

  allvars <- param.names("zeta", M)
  d1 <- deriv3(txt.f, allvars, hessian = TRUE)

  zeta <- as.matrix(zeta)
  for (ii in 1:M)
    assign(paste0("zeta", ii), zeta[, ii])

  ans <- eval(d1)
  ans
}  # acat.deriv()






 brat <- function(refgp = "last",
                  refvalue = 1,
                  ialpha = 1) {
  if (!is.Numeric(ialpha, positive = TRUE))
    stop("'ialpha' must contain positive values only")
  if (!is.Numeric(refvalue, length.arg = 1, positive = TRUE))
    stop("'refvalue' must be a single positive value")

  if (!is.character(refgp) &&
      !is.Numeric(refgp, length.arg = 1,
                  integer.valued = TRUE, positive = TRUE))
    stop("'refgp' must be a single positive integer")


  new("vglmff",
  blurb = c(paste("Bradley-Terry model (without ties)\n\n"),
            "Links:   ",
            namesof("alpha's", "loglink")),
  infos = eval(substitute(function(...) {
    list(M1 = NA,  # zz -1?
         Q1 = NA,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = as.character(NA),
         refvalue = .refvalue ,
         refgp = .refgp ,
         ialpha = .ialpha )
  }, list( .ialpha = ialpha,
           .refgp = refgp,
           .refvalue = refvalue ))),

  initialize = eval(substitute(expression({
    are.ties <- attr(y, "are.ties")  # If Brat() was used
    if (is.logical(are.ties) && are.ties)
        stop("use bratt(), not brat(), when there are ties")

    try.index <- 1:400
    M <- (seq_along(try.index))[(try.index+1)*(try.index) ==
                                 ncol(y)]
    if (!is.finite(M))
      stop("cannot determine 'M'")
    ialpha <- matrix(rep_len( .ialpha , M), n, M, byrow = TRUE)
    etastart <- matrix(theta2eta(ialpha, "loglink",
                                 earg = list(theta = NULL)),
                       n, M, byrow = TRUE)
    refgp <- .refgp
    if (!intercept.only)
      warning("this function only works with ",
              "intercept-only models")
    extra$ybrat.indices <- .brat.indices(NCo = M+1,
                                         are.ties = FALSE)
    uindex <- if ( .refgp == "last")
                  1:M else (1:(M+1))[-( .refgp ) ]

    predictors.names <-
      namesof(paste0("alpha", uindex), "loglink",
              short = TRUE)

  }),
  list( .refgp = refgp, .ialpha = ialpha ))),

  linkinv = eval(substitute( function(eta, extra = NULL) {
    probs <- NULL
    eta <- as.matrix(eta)  # in case M = 1; prior to 20171227
    for (ii in 1:nrow(eta)) {
      alpha <- .brat.alpha(eta2theta(eta[ii, ], "loglink",
                                     earg = list(theta = NULL)),
                           .refvalue , .refgp )
      alpha1 <- alpha[extra$ybrat.indices[, "rindex"]]
      alpha2 <- alpha[extra$ybrat.indices[, "cindex"]]
      probs <- rbind(probs, alpha1 / (alpha1 + alpha2))
    }
    if (NROW(probs) == NROW(eta) &&
        NCOL(probs) == NCOL(eta))
      dimnames(probs) <- dimnames(eta)
    probs
  },
  list( .refgp = refgp, .refvalue = refvalue) )),

  last = eval(substitute(expression({
    misc$link <- rep_len("loglink", M)
    names(misc$link) <- paste0("alpha", uindex)

    misc$earg <- vector("list", M)
    names(misc$earg) <- names(misc$link)
    for (ii in 1:M)
      misc$earg[[ii]] <- list(theta = NULL)

    misc$refgp <- .refgp
    misc$refvalue <- .refvalue
  }),
  list( .refgp = refgp, .refvalue = refvalue ))),

  loglikelihood =
    function(mu, y, w, residuals = FALSE, eta, extra = NULL,
             summation = TRUE) {
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ycounts <- if (is.numeric(extra$orig.w))
                 y * w / extra$orig.w else
                 y * w  # Convert proportions to counts
      nvec <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else round(w)

      smallno <- 1.0e4 * .Machine$double.eps
      if (max(abs(ycounts - round(ycounts))) > smallno)
        warning("converting 'ycounts' to integer ",
                "in @loglikelihood")
      ycounts <- round(ycounts)

      ll.elts <-
        (if (is.numeric(extra$orig.w)) extra$orig.w else 1) *
         dmultinomial(x = ycounts, size = nvec, prob = mu,
                      log = TRUE, dochecking = FALSE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  },
  vfamily = c("brat", "VGAMcategorical"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    probs <- NULL
    eta <- as.matrix(eta)  # in case M = 1
    for (ii in 1:nrow(eta)) {
      alpha <- .brat.alpha(eta2theta(eta[ii, ], "loglink",
                                     earg = list(theta = NULL)),
                           .refvalue , .refgp )
      alpha1 <- alpha[extra$ybrat.indices[, "rindex"]]
      alpha2 <- alpha[extra$ybrat.indices[, "cindex"]]
      probs <- rbind(probs, alpha1 / (alpha1 + alpha2))
    }


    okay1 <- all(is.finite(probs)) && all(0 < probs & probs < 1)
    okay1
  },
  list( .refvalue = refvalue, .refgp = refgp) )),
  deriv = eval(substitute(expression({
    ans <- NULL
    uindex <- if ( .refgp == "last")
                  1:M else (1:(M+1))[-( .refgp ) ]
    eta <- as.matrix(eta)  # in case M = 1
    for (ii in 1:nrow(eta)) {
      alpha <- .brat.alpha(eta2theta(eta[ii, ], "loglink",
                             earg = list(theta = NULL)),
                           .refvalue, .refgp )
      ymat <- InverseBrat(y[ii, ], NCo = M+1, diag = 0)
      answer <- rep_len(0, M)
      for (aa in 1:(M+1)) {
        answer <- answer + (1 - (aa == uindex)) *
          (ymat[uindex, aa] * alpha[aa] - ymat[aa, uindex] *
          alpha[uindex]) / (alpha[aa] + alpha[uindex])
      }
      ans <- rbind(ans, w[ii] * answer)
    }
    dimnames(ans) <- dimnames(eta)
    ans
  }),
  list( .refvalue = refvalue, .refgp = refgp) )),
  weight = eval(substitute(expression({
    wz <- matrix(0, n, dimm(M))
    for (ii in 1:nrow(eta)) {
      alpha <- .brat.alpha(eta2theta(eta[ii, ], "loglink",
                                     earg = list(theta = NULL)),
                          .refvalue, .refgp)
      ymat <- InverseBrat(y[ii, ], NCo = M+1, diag = 0)
      for (aa in 1:(M+1)) {
        wz[ii, 1:M] <- wz[ii, 1:M] + (1 - (aa == uindex)) *
                       (ymat[aa, uindex] + ymat[uindex, aa]) *
            alpha[aa] * alpha[uindex] / (
            alpha[aa] + alpha[uindex])^2
      }
      if (M > 1) {
        ind5 <- iam(1, 1, M, both = TRUE, diag = FALSE)
        wz[ii, (M+1):ncol(wz)] <-
          -(ymat[cbind(uindex[ind5$row], uindex[ind5$col])] +
            ymat[cbind(uindex[ind5$col], uindex[ind5$row])]) *
         alpha[uindex[ind5$col]] * alpha[uindex[ind5$row]] / (
         alpha[uindex[ind5$row]] + alpha[uindex[ind5$col]])^2
      }
    }
    wz <- c(w) * wz
    wz
  }),
  list( .refvalue = refvalue, .refgp = refgp ))))
}  # brat()







 bratt <- function(refgp = "last",
                   refvalue = 1,
                   ialpha = 1,
                   i0 = 0.01) {
  if (!is.Numeric(i0, length.arg = 1, positive = TRUE))
    stop("'i0' must be a single positive value")
  if (!is.Numeric(ialpha, positive = TRUE))
    stop("'ialpha' must contain positive values only")
  if (!is.Numeric(refvalue, length.arg = 1, positive = TRUE))
    stop("'refvalue' must be a single positive value")

  if (!is.character(refgp) &&
     !is.Numeric(refgp, length.arg = 1,
                 integer.valued = TRUE, positive = TRUE))
    stop("'refgp' must be a single positive integer")


  new("vglmff",
  blurb = c(paste("Bradley-Terry model (with ties)\n\n"),
            "Links:   ",
            namesof("alpha's", "loglink"), ", log(alpha0)"),
  infos = eval(substitute(function(...) {
    list(M1 = NA,  # zz -1?
         Q1 = NA,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = as.character(NA),
         refvalue = .refvalue ,
         refgp = .refgp ,
         i0 = .i0 ,
         ialpha = .ialpha )
  }, list( .ialpha = ialpha,
           .i0 = i0,
           .refgp = refgp,
           .refvalue = refvalue ))),

  initialize = eval(substitute(expression({
    try.index <- 1:400
    M <- (seq_along(try.index))[(try.index*(try.index-1)) ==
                                ncol(y)]
    if (!is.Numeric(M, length.arg = 1, integer.valued = TRUE))
      stop("cannot determine 'M'")
    NCo <- M  # Number of contestants

    are.ties <- attr(y, "are.ties")  # If Brat() was used
    if (is.logical(are.ties)) {
      if (!are.ties)
        stop("use brat(), not bratt(), when there are no ties")
      ties <- attr(y, "ties")
    } else {
      are.ties <- FALSE
      ties <- 0 * y
    }

    ialpha <- rep_len( .ialpha, NCo-1)
    ialpha0 <- .i0
    etastart <-
      cbind(matrix(theta2eta(ialpha,
                             "loglink",
                             list(theta = NULL)),
                   n, NCo-1, byrow = TRUE),
            theta2eta(rep_len(ialpha0, n), "loglink",
                      list(theta = NULL)))
    refgp <- .refgp
    if (!intercept.only)
      warning("this function only works with ",
              "intercept-only models")
    extra$ties <- ties # Flat (1-row) matrix
    extra$ybrat.indices <-
        .brat.indices(NCo = NCo, are.ties = FALSE)
    extra$tbrat.indices <- .brat.indices(NCo = NCo,
                                         are.ties = TRUE)
    extra$dnties <- dimnames(ties)
    uindex <- if (refgp == "last")
                1:(NCo-1) else (1:(NCo))[-refgp ]

    predictors.names <- c(
      namesof(paste0("alpha", uindex), "loglink",
              short = TRUE),
      namesof("alpha0", "loglink", short = TRUE))
  }),
  list( .refgp = refgp,
        .i0 = i0,
        .ialpha = ialpha ))),

  linkinv = eval(substitute( function(eta, extra = NULL) {
    probs <- qprobs <- NULL
    M <- ncol(eta)
    for (ii in 1:nrow(eta)) {
      alpha <- .brat.alpha(eta2theta(eta[ii, -M],
                                     "loglink"),
                           .refvalue , .refgp )
      alpha0 <- loglink(eta[ii, M], inverse = TRUE)
      alpha1 <- alpha[extra$ybrat.indices[, "rindex"]]
      alpha2 <- alpha[extra$ybrat.indices[, "cindex"]]
         probs <- rbind( probs,
                         alpha1 / (alpha1 + alpha2 + alpha0))
        qprobs <- rbind(qprobs,
                         alpha0 / (alpha1 + alpha2 + alpha0))
    }
    if (length(extra$dnties))
      dimnames(qprobs) <- extra$dnties
    attr(probs, "probtie") <- qprobs
    probs
  },
  list( .refgp = refgp, .refvalue = refvalue) )),
  last = eval(substitute(expression({
    misc$link <- rep_len("loglink", M)
    names(misc$link) <- c(paste0("alpha", uindex),
                          "alpha0")


    misc$earg <- vector("list", M)
    names(misc$earg) <- names(misc$link)
    for (ii in 1:M)
      misc$earg[[ii]] <- list(theta = NULL)


    misc$refgp <- .refgp
    misc$refvalue <- .refvalue
    misc$alpha  <- alpha
    misc$alpha0 <- alpha0
  }),
  list( .refgp = refgp, .refvalue = refvalue ))),
  loglikelihood =
    function(mu, y, w, residuals = FALSE, eta, extra = NULL,
             summation = TRUE) {
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * (y * log(mu) +
                 0.5 * extra$ties * log(attr(mu, "probtie")))
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  },
  vfamily = c("bratt", "VGAMcategorical"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    probs <- qprobs <- NULL
    M <- ncol(eta)
    for (ii in 1:nrow(eta)) {
      alpha <- .brat.alpha(eta2theta(eta[ii, -M],
                                     "loglink"),
                           .refvalue , .refgp )
      alpha0 <- loglink(eta[ii, M], inverse = TRUE)
      alpha1 <- alpha[extra$ybrat.indices[, "rindex"]]
      alpha2 <- alpha[extra$ybrat.indices[, "cindex"]]
       probs <- rbind( probs,
                      alpha1 / (alpha1 + alpha2 + alpha0))
      qprobs <- rbind(qprobs,
                      alpha0 / (alpha1 + alpha2 + alpha0))
    }
    okay1 <-
    all(is.finite( probs)) && all(0 <  probs &  probs < 1) &&
    all(is.finite(qprobs)) && all(0 < qprobs & qprobs < 1)
    okay1
  },
  list( .refvalue = refvalue, .refgp = refgp) )),
  deriv = eval(substitute(expression({
    ans <- NULL
    ties <- extra$ties
    NCo <- M
    uindex <- if ( .refgp == "last")
                1:(M-1) else (1:(M))[-( .refgp )]
    eta <- as.matrix(eta)
    for (ii in 1:nrow(eta)) {
      alpha <- .brat.alpha(eta2theta(eta[ii, -M], "loglink",
                                     list(theta = NULL)),
                           .refvalue, .refgp )
      alpha0 <- loglink(eta[ii, M], inverse = TRUE)
      ymat <- InverseBrat(   y[ii, ], NCo = M, diag = 0)
      tmat <- InverseBrat(ties[ii, ], NCo = M, diag = 0)
      answer <- rep_len(0, NCo-1)  # deriv wrt eta[-M]
      for (aa in 1:NCo) {
        Daj <- alpha[aa] + alpha[uindex] + alpha0
        pja <- alpha[uindex] / Daj
        answer <- answer + alpha[uindex] *
            (-ymat[aa, uindex] + ymat[uindex, aa] *
             (1 - pja) / pja - tmat[uindex, aa]) / Daj
      }
      deriv0 <- 0 # deriv wrt eta[M]
      for (aa in 1:(NCo-1))
        for (bb in (aa+1):NCo) {
          Dab <- alpha[aa] + alpha[bb] + alpha0
          qab <- alpha0 / Dab
          deriv0 <- deriv0 + alpha0 *
                    (-ymat[aa, bb] - ymat[bb,aa] +
                    tmat[aa, bb] * (1 - qab) / qab) / Dab
          }
        ans <- rbind(ans, w[ii] * c(answer, deriv0))
    }
    dimnames(ans) <- dimnames(eta)
    ans
  }),
  list( .refvalue = refvalue, .refgp = refgp) )),
  weight = eval(substitute(expression({
    wz <- matrix(0, n, dimm(M))   # includes diagonal
    for (ii in 1:nrow(eta)) {
      alpha <- .brat.alpha(eta2theta(eta[ii, -M], "loglink",
                           earg = list(theta = NULL)),
                          .refvalue, .refgp)
      alpha0 <- loglink(eta[ii, M], inverse = TRUE)
      ymat <- InverseBrat(   y[ii, ], NCo = M, diag = 0)
      tmat <- InverseBrat(ties[ii, ], NCo = M, diag = 0)

      for (aa in 1:(NCo)) {
        Daj <- alpha[aa] + alpha[uindex] + alpha0
        pja <- alpha[uindex] / Daj
        nja <- ymat[aa, uindex] + ymat[uindex, aa] +
               tmat[uindex, aa]
        wz[ii, 1:(NCo-1)] <- wz[ii, 1:(NCo - 1)] +
                             alpha[uindex]^2 * nja *
                             (1 - pja) / (pja * Daj^2)
        if (aa < NCo)
          for (bb in (aa+1):(NCo)) {
            nab <- ymat[aa,bb] + ymat[bb,aa] + tmat[bb,aa]
            Dab <- alpha[aa] + alpha[bb] + alpha0
            qab <- alpha0 / Dab
            wz[ii, NCo] <- wz[ii,NCo] + alpha0^2 * nab *
                     (1-qab) / (qab * Dab^2)
          }
      }

      if (NCo > 2) {
        ind5 <- iam(1, 1, M = NCo, both = TRUE, diag = FALSE)
        alphajunk <- c(alpha, junk = NA)
        mat4 <- cbind(uindex[ind5$row],uindex[ind5$col])
        wz[ii,(M+1):ncol(wz)] <-
            -(ymat[mat4] + ymat[mat4[, 2:1]] +
           tmat[mat4]) * alphajunk[uindex[ind5$col]] *
           alphajunk[uindex[ind5$row]] / (alpha0 +
           alphajunk[uindex[ind5$row]] +
           alphajunk[uindex[ind5$col]])^2
      }
      for (sss in seq_along(uindex)) {
        jay <- uindex[sss]
        naj <- ymat[, jay] + ymat[jay, ] + tmat[, jay]
        Daj <- alpha[jay] + alpha + alpha0
        wz[ii, iam(sss, NCo, M = NCo, diag = TRUE)] <-
            -alpha[jay] * alpha0 * sum(naj / Daj^2)
      }
    }
    wz <- c(w) * wz
    wz
  }),
  list( .refvalue = refvalue, .refgp = refgp ))))
}  # bratt()



.brat.alpha <- function(vec, value, posn) {
  if (is.character(posn))
    if (posn != "last")
      stop("can only handle \"last\"") else
      return(c(vec, value))
  c(if (posn == 1) NULL else vec[1:(posn-1)], value,
    if (posn == length(vec) + 1) NULL else
      vec[posn:length(vec)])
}


.brat.indices <- function(NCo, are.ties = FALSE) {
  if (!is.Numeric(NCo, length.arg = 1,
                  integer.valued = TRUE) ||
      NCo < 2)
    stop("bad input for 'NCo'")
  m <- diag(NCo)
  if (are.ties) {
    cbind(rindex = row(m)[col(m) <  row(m)],
          cindex = col(m)[col(m) <  row(m)])
  } else
    cbind(rindex = row(m)[col(m) != row(m)],
          cindex = col(m)[col(m) != row(m)])
}



 Brat <- function(mat,
                  ties = 0 * mat,
                  string = c(">", "=="),
                  whitespace = FALSE) {


  stopifnot(is.logical(whitespace) &&
            length(whitespace) == 1)
  fillerChar <- ifelse(whitespace, " ", "")
  string <- paste0(fillerChar, string, fillerChar)


  allargs <- list(mat)  # ,...
  callit <- if (length(names(allargs))) names(allargs) else
            as.character(seq_along(allargs))
  ans <- ans.ties <- NULL
  for (ii in seq_along(allargs)) {
    m <- allargs[[ii]]
    if (!is.matrix(m) || dim(m)[1] != dim(m)[2])
      stop("m must be a square matrix")

    diag(ties) <- 0
    if (!all(ties == t(ties)))
      stop("ties must be a symmetric matrix")
    are.ties <- any(ties > 0)
    diag(ties) <- NA

    diag(m) <- 0  # Could have been NAs
    if (anyNA(m))
      stop("missing values not allowed ",
           "(except on the diagonal)")
    diag(m) <- NA

    dm <- as.data.frame.table(m)
    dt <- as.data.frame.table(ties)
    dm <- dm[!is.na(dm$Freq), ]
    dt <- dt[!is.na(dt$Freq), ]
    usethis1 <- paste0(dm[, 1], string[1], dm[, 2])
    usethis2 <- paste0(dm[, 1], string[2], dm[, 2])
    ans <- rbind(ans, matrix(dm$Freq, nrow = 1))
    ans.ties <- rbind(ans.ties, matrix(dt$Freq, nrow = 1))
  }
  dimnames(ans) <- list(callit, usethis1)
  dimnames(ans.ties) <- list(callit, usethis2)
  attr(ans, "ties") <- ans.ties
  attr(ans, "are.ties") <- are.ties
  ans
}  # Brat()




InverseBrat <-
  function(yvec,
  NCo = (1:900)[(1:900)*((1:900)-1) == ncol(rbind(yvec))],
  multiplicity = if (is.matrix(yvec)) nrow(yvec) else 1,
           diag = NA,
           string = c(">", "=="),
           whitespace = FALSE) {



  stopifnot(is.logical(whitespace) &&
            length(whitespace) == 1)
  fillerChar <- ifelse(whitespace, " ", "")
  string <- paste0(fillerChar, string, fillerChar)


  ans <- array(diag, c(NCo, NCo, multiplicity))
  yvec.orig <- yvec
  yvec <- c(yvec)
  ptr <- 1
  for (mul in 1:multiplicity)
    for (i1 in 1:(NCo))
      for (i2 in 1:(NCo))
        if (i1 != i2) {
          ans[i2, i1, mul] <- yvec[ptr]
          ptr <- ptr + 1
        }
  ans <- if (multiplicity > 1) ans else matrix(ans, NCo, NCo)

  if (is.array(yvec.orig) ||
      is.matrix(yvec.orig)) {
    names.yvec <- dimnames(yvec.orig)[[2]]
    ii <- strsplit(names.yvec, string[1])
    cal <- NULL
    for (kk in c(NCo, 1:(NCo-1)))
      cal <- c(cal, (ii[[kk]])[1])
    if (multiplicity>1) {
      dimnames(ans) <- list(cal, cal, dimnames(yvec.orig)[[1]])
    } else {
      dimnames(ans) <- list(cal, cal)
    }
  }
  ans
}  # InverseBrat()






 ordpoisson <-
  function(cutpoints,
           countdata = FALSE, NOS = NULL, Levels = NULL,
           init.mu = NULL, parallel = FALSE,
           zero = NULL,
           link = "loglink") {

  link <- as.list(substitute(link))
  earg  <- link2list(link)
  link <- attr(earg, "function.name")




  fcutpoints <- cutpoints[is.finite(cutpoints)]
  if (!is.Numeric(fcutpoints, integer.valued = TRUE) ||
      any(fcutpoints < 0))
    stop("'cutpoints' must have non-negative integer or Inf ",
         "values only")
  if (is.finite(cutpoints[length(cutpoints)]))
    cutpoints <- c(cutpoints, Inf)

  if (!is.logical(countdata) || length(countdata) != 1)
    stop("argument 'countdata' must be a single logical")
  if (countdata) {
    if (!is.Numeric(NOS, integer.valued = TRUE, positive = TRUE))
      stop("'NOS' must have integer values only")
    if (!is.Numeric(Levels, integer.valued = TRUE,
                    positive = TRUE) ||
        any(Levels < 2))
      stop("'Levels' must have integer values (>= 2) only")
    Levels <- rep_len(Levels, NOS)
  }


  new("vglmff",
  blurb = c(paste("Ordinal Poisson model\n\n"),
            "Link:     ", namesof("mu", link, earg)),


  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(1, M, 1), x = x,
                        bool = .parallel ,
                        apply.int = TRUE,
                        constraints = constraints)
    constraints <- cm.zero.VGAM(constraints,
             x = x, .zero , M = M, M1 = 1,
             predictors.names = predictors.names)
  }),
  list( .parallel = parallel, .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 1,
         expected = TRUE,
         multipleResponses = TRUE,
         parameters.names = c("mu"),
         lmu = .link ,
         zero = .zero )
  }, list( .zero = zero, .link = link ))),

  initialize = eval(substitute(expression({
    orig.y <- cbind(y)  # Convert y into a matrix if necessary
    if ( .countdata ) {
      extra$NOS <- M <- NOS <- .NOS
      extra$Levels <- Levels <- .Levels
      y.names <- dimnames(y)[[2]]  # Hopefully
    } else {
      if (any(w != 1) || NCOL(w) != 1)
        stop("the 'weights' argument must be ",
             "a vector of all ones")
        extra$NOS <- M <- NOS <- if (is.Numeric( .NOS ))
                                     .NOS else
          ncol(orig.y)
      Levels <- rep_len(if (is.Numeric( .Levels ))
                        .Levels else 0, NOS)
      if (!is.Numeric( .Levels ))
        for (iii in 1:NOS) {
          Levels[iii] <- length(unique(sort(orig.y[,iii])))
        }
      extra$Levels <- Levels
    }


    initmu <- if (is.Numeric( .init.mu ))
      rep_len( .init.mu , NOS) else NULL
    cutpoints <- rep_len( .cutpoints, sum(Levels))
    delete.zero.colns <- FALSE
    use.y <- if ( .countdata ) y else matrix(0, n, sum(Levels))
    use.etastart <- matrix(0, n, M)
    cptr <- 1
    for (iii in 1:NOS) {
        y <- factor(orig.y[,iii], levels=(1:Levels[iii]))
        if ( !( .countdata )) {
  eval(process.categorical.data.VGAM)  # Creates mustart & y
          use.y[,cptr:(cptr+Levels[iii]-1)] <- y
        }
        use.etastart[,iii] <- if (is.Numeric(initmu))
            initmu[iii] else
            median(cutpoints[cptr:(cptr+Levels[iii]-1-1)])
        cptr <- cptr + Levels[iii]
    }
    mustart <- NULL # Overwrite it
    etastart <- theta2eta(use.etastart, .link , earg = .earg )
    y <- use.y  # n x sum(Levels)
    M <- NOS
    for (iii in 1:NOS) {
        mu.names <- paste0("mu", iii, ".")
    }

    ncoly <- extra$ncoly <- sum(Levels)
    cp.vector <- rep_len( .cutpoints , ncoly)
    extra$countdata <- .countdata
    extra$cutpoints <- cp.vector
    extra$n <- n

    mynames <- param.names("mu", M, skip1 = TRUE)
    predictors.names <-
      namesof(mynames, .link , earg = .earg , tag = FALSE)
  }),
  list( .link = link, .countdata = countdata,
        .earg = earg,
        .cutpoints = cutpoints, .NOS = NOS,
        .Levels = Levels,
        .init.mu = init.mu))),
  linkinv = eval(substitute( function(eta, extra = NULL) {
    mu <- eta2theta(eta, .link , earg = .earg )  # Pois means
    mu <- cbind(mu)
    mu
  },
  list( .link = link,
       .earg = earg,
       .countdata = countdata ))),
  last = eval(substitute(expression({
    if ( .countdata ) {
      misc$link <- .link
      misc$earg <- list( .earg )
    } else {
      misc$link <- rep_len( .link , M)
      names(misc$link) <- mynames

      misc$earg <- vector("list", M)
      names(misc$earg) <- names(misc$link)
      for (ii in 1:M)
        misc$earg[[ii]] <- .earg
    }
    misc$parameters <- mynames
    misc$countdata <- .countdata
    misc$true.mu = FALSE    # $fitted is not a true mu
  }),
  list( .link = link,
       .countdata = countdata, .earg = earg ))),
    loglikelihood =
      function(mu, y, w, residuals = FALSE, eta, extra = NULL,
               summation = TRUE) {
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      probs <- ordpoissonProbs(extra, mu)
      index0 <- y == 0
      probs[index0] <- 1
      pindex0 <- probs == 0
      probs[pindex0] <- 1
      if (summation) {
        sum(pindex0) * (-1.0e+10) + sum(w * y * log(probs))
      } else {
        stop("20140311; 'summation=F' not done yet")
      }
      }
  },
  vfamily = c("ordpoisson", "VGAMcategorical"),
  deriv = eval(substitute(expression({
    probs <- ordpoissonProbs(extra, mu)
    probs.use <- pmax(probs, .Machine$double.eps * 1.0e-0)

    cp.vector <- extra$cutpoints
    NOS <- extra$NOS
    Levels <- extra$Levels
    resmat <- matrix(0, n, M)
    dl.dprob <- y / probs.use
    dmu.deta <- dtheta.deta(mu, .link , earg = .earg )
    dprob.dmu <- ordpoissonProbs(extra, mu, deriv = 1)
    cptr <- 1
    for (iii in 1:NOS) {
      for (kkk in 1:Levels[iii]) {
       resmat[,iii] <- resmat[,iii] +
                      dl.dprob[,cptr] * dprob.dmu[,cptr]
       cptr <- cptr + 1
      }
    }
    resmat <- c(w) * resmat * dmu.deta
    resmat
  }),
  list( .link = link,
       .earg = earg, .countdata = countdata ))),
  weight = eval(substitute(expression({
    d2l.dmu2 <- matrix(0, n, M)  # Diagonal matrix
    cptr <- 1
    for (iii in 1:NOS) {
      for (kkk in 1:Levels[iii]) {
        d2l.dmu2[,iii] <- d2l.dmu2[,iii] +
            dprob.dmu[,cptr]^2 / probs.use[,cptr]
        cptr <- cptr + 1
      }
    }
    wz <- c(w) * d2l.dmu2 * dmu.deta^2
    wz
  }),
  list( .earg = earg, .link = link,
       .countdata = countdata ))))
}



ordpoissonProbs <- function(extra, mu, deriv = 0) {
  cp.vector <- extra$cutpoints
  NOS <- extra$NOS
  if (deriv == 1) {
    dprob.dmu <- matrix(0, extra$n, extra$ncoly)
  } else {
    probs <- matrix(0, extra$n, extra$ncoly)
  }
  mu <- cbind(mu)
  cptr <- 1
  for (iii in 1:NOS) {
    if (deriv == 1) {
      dprob.dmu[,cptr] <- -dpois(x = cp.vector[cptr],
                                 lambda = mu[,iii])
    } else {
      probs[,cptr] <- ppois(q = cp.vector[cptr],
                            lambda = mu[,iii])
    }
    cptr <- cptr + 1
    while (is.finite(cp.vector[cptr])) {
      if (deriv == 1) {
        dprob.dmu[,cptr] <-
            dpois(x = cp.vector[cptr-1], lambda = mu[,iii]) -
            dpois(x = cp.vector[cptr  ], lambda = mu[,iii])
      } else {
        probs[,cptr] <-
            ppois(q = cp.vector[cptr  ], lambda = mu[,iii]) -
            ppois(q = cp.vector[cptr-1], lambda = mu[,iii])
      }
      cptr <- cptr + 1
    }
    if (deriv == 1) {
        dprob.dmu[,cptr] <-
            dpois(x = cp.vector[cptr-1], lambda = mu[,iii]) -
            dpois(x = cp.vector[cptr  ], lambda = mu[,iii])
    } else {
        probs[,cptr] <-
            ppois(q = cp.vector[cptr  ], lambda = mu[,iii]) -
            ppois(q = cp.vector[cptr-1], lambda = mu[,iii])
    }
    cptr <- cptr + 1
  }
  if (deriv == 1) dprob.dmu else probs
}









findFirstMethod <- function(methodsfn, charvec) {
  answer <- NULL
  for (ii in seq_along(charvec)) {
    if (existsMethod(methodsfn,
                     signature(VGAMff = charvec[ii]))) {
      answer <- charvec[ii]
      break
    }
  }
  answer
}



margeff <- function(object, subset = NULL, ...) {


  try.this <- findFirstMethod("margeffS4VGAM",
                              object@family@vfamily)
  if (length(try.this)) {
    margeffS4VGAM(object = object,
            subset = subset,
            VGAMff = new(try.this),
        ...)
  } else {
    stop("Could not find a methods function for ",
         "'margeffS4VGAM' emanating ",
         "from '", object@family@vfamily[1], "'")
  }
}






subsetarray3 <- function(array3, subset = NULL) {
  if (is.null(subset)) {
    return(array3)
  } else
    if (is.numeric(subset) && (length(subset) == 1)) {
      return(array3[, , subset])
  } else {
    return(array3[, , subset])
  }
  warning("argument 'subset' unmatched. Doing nothing")
  array3
}





setClass("VGAMcategorical",     contains = "vglmff")

setClass("VGAMordinal",         contains = "VGAMcategorical")
setClass("multinomial",         contains = "VGAMcategorical")

setClass("acat",                contains = "VGAMordinal")
setClass("cumulative",          contains = "VGAMordinal")
setClass("cratio",              contains = "VGAMordinal")
setClass("sratio",              contains = "VGAMordinal")








setMethod("margeffS4VGAM",
          signature(VGAMff = "VGAMcategorical"),
  function(object,
           subset = NULL,
           VGAMff,
           ...) {
  object@post$M <- M   <- object@misc$M
  object@post$n <- nnn <- object@misc$n
  invisible(object)
  })  # "VGAMcategorical"






setMethod("margeffS4VGAM",  signature(VGAMff ="multinomial"),
  function(object,
           subset = NULL,
           VGAMff,
           ...) {

  object <- callNextMethod(VGAMff = VGAMff,
                           object = object,
                           subset = subset,
                           ...)

  M   <- object@misc$M
  nnn <- object@misc$n
    cfit <- coefvlm(object, matrix.out = TRUE)
    rlev <- object@misc$refLevel
    if (!length(rlev))
      relev <- M+1  # Default
    Bmat <- matrix(0, nrow(cfit), 1 + ncol(cfit))
    Bmat[, -rlev] <- cfit
    ppp   <- nrow(Bmat)
    pvec1 <- fitted(object)[1, ]
    rownames(Bmat) <- rownames(cfit)
    colnames(Bmat) <-
      if (length(names(pvec1))) names(pvec1) else
                 param.names("mu", M+1)


    BB <- array(Bmat, c(ppp, M+1, nnn))
    pvec  <- c(t(fitted(object)))
    pvec  <- rep(pvec, each = ppp)
    temp1 <- array(BB * pvec, c(ppp, M+1, nnn))
    temp2 <- aperm(temp1, c(2, 1, 3))  # (M+1) x ppp x nnn
    temp2 <- colSums(temp2)  # ppp x nnn
    temp2 <- array(rep(temp2, each = M+1), c(M+1, ppp, nnn))
    temp2 <- aperm(temp2, c(2, 1, 3))  # ppp x (M+1) x nnn
    temp3 <- pvec
    ans.mlm <- array((BB - temp2) * temp3, c(ppp, M+1, nnn),
               dimnames = list(dimnames(Bmat)[[1]],
                               dimnames(Bmat)[[2]],
                               dimnames(fitted(object))[[1]]))
    return(subsetarray3(ans.mlm, subset = subset))
  })  # "multinomial"






setMethod("margeffS4VGAM",  signature(VGAMff = "VGAMordinal"),
  function(object,
           subset = NULL,
           VGAMff,
           ...) {
  M   <- object@misc$M
  nnn <- object@misc$n

  object@post$reverse <- object@misc$reverse
  object@post$linkfunctions <-
      linkfunctions <- object@misc$link
  object@post$all.eargs <- all.eargs <- object@misc$earg
  object@post$Bmat <- Bmat <- coefvlm(object,
                                      matrix.out = TRUE)
  object@post$ppp <- nrow(Bmat)
  etamat <- predict(object)

  hdot <- Thetamat <- etamat
  for (jlocal in 1:M) {
    Thetamat[, jlocal] <- eta2theta(etamat[, jlocal],
                                    linkfunctions[jlocal],
                                    all.eargs[[jlocal]])
    hdot[, jlocal] <- dtheta.deta(Thetamat[, jlocal],
                                  linkfunctions[jlocal],
                                  all.eargs[[jlocal]])
  }  # jlocal


  object@post$hdot <- hdot
  object@post$Thetamat <- Thetamat
  object
  })  # "VGAMordinal"





setClass("poissonff",   contains = "VGAMcategorical")



setMethod("margeffS4VGAM",  signature(VGAMff = "poissonff"),
  function(object,
           subset = NULL,
           VGAMff,
           ...) {



  etamat <- predict(object)  # nnn x M
  nnn <- nrow(etamat)
  M   <- ncol(etamat)
  if (!all(linkfunvlm(object) == "loglink"))
    stop("Only objects using 'loglink' are permitted here")

  mustar <- fitted(object)  # nnn x NOS
  M1 <- npred(object, type = "one.response")
  if (M1 != 1) stop("M1 should be one")
  NOS <- M / M1
  ymat <- depvar(object)
  cnymat <- colnames(ymat)
  rnymat <- rownames(ymat)
  coefmat <- coef(object, matrix = TRUE)  # ppp x NOS
  ppp <- nrow(coefmat)

  betamat <- matrix(c(coefmat), nnn, ppp * NOS, byrow = TRUE)
  ans <- betamat * kronecker(mustar, matrix(1, 1, ppp))

  dim(ans) <- c(nnn, ppp, NOS)
  ans <- aperm(ans, c(2, 3, 1))  # ppp x NOS x nnn
  dimnames(ans) <- list(rownames(coefmat), cnymat, rnymat)

  ans
})  # poissonff








setClass("negbinomial",   contains = "VGAMcategorical")




setMethod("margeffS4VGAM",  signature(VGAMff = "negbinomial"),
  function(object,
           subset = NULL,
           VGAMff,
           ...) {


  etamat <- predict(object)  # nnn x M
  nnn <- nrow(etamat)
  M   <- ncol(etamat)
  if (!all(linkfunvlm(object)[c(TRUE, FALSE)] == "loglink"))
    stop("Only objects whose mean is modelled using ",
         "'loglink' are permitted here")

  mustar <- fitted(object)  # nnn x NOS
  M1 <- npred(object, type = "one.response")
  if (M1 != 2) stop("M1 should be two")
  NOS <- M / M1
  ymat <- depvar(object)
  cnymat <- colnames(ymat)
  rnymat <- rownames(ymat)
  coefmat <- coef(object, matrix = TRUE)[, c(TRUE, FALSE),
                          drop = FALSE]  # ppp x NOS
  ppp <- nrow(coefmat)

  betamat <- matrix(c(coefmat), nnn, ppp * NOS, byrow = TRUE)
  ans <- betamat * kronecker(mustar, matrix(1, 1, ppp))

  dim(ans) <- c(nnn, ppp, NOS)
  ans <- aperm(ans, c(2, 3, 1))  # ppp x NOS x nnn
  dimnames(ans) <- list(rownames(coefmat), cnymat, rnymat)

  ans
})  # negbinomial






setClass("posnegbinomial",   contains = "VGAMcategorical")




setMethod("margeffS4VGAM",
          signature(VGAMff = "posnegbinomial"),
  function(object,
           subset = NULL,
           VGAMff,
           ...) {


  etamat <- predict(object)  # nnn x M
  nnn <- nrow(etamat)
  M   <- ncol(etamat)
  if (!all(linkfunvlm(object)[c(TRUE, FALSE)] == "loglink"))
    stop("Only objects whose mean is modelled using ",
         "'loglink' are permitted here")

  mustar <- fitted(object)  # nnn x NOS
  M1 <- npred(object, type = "one.response")
  if (M1 != 2) stop("M1 should be two")
  NOS <- M / M1
  ymat <- depvar(object)
  cnymat <- colnames(ymat)
  rnymat <- rownames(ymat)
  coefmat <- coef(object, matrix = TRUE)[, c(TRUE, FALSE),
                          drop = FALSE]  # ppp x NOS
  ppp <- nrow(coefmat)

  betamat <- matrix(c(coefmat), nnn, ppp * NOS, byrow = TRUE)
  ans <- betamat * kronecker(mustar, matrix(1, 1, ppp))

  dim(ans) <- c(nnn, ppp, NOS)
  ans <- aperm(ans, c(2, 3, 1))  # ppp x NOS x nnn
  dimnames(ans) <- list(rownames(coefmat), cnymat, rnymat)

  ans
})  # posnegbinomial









setClass("tobit",          contains = "VGAMcategorical")



setMethod("margeffS4VGAM",  signature(VGAMff = "tobit"),
  function(object,
           subset = NULL,
           VGAMff,
           ...) {


  object <- callNextMethod(VGAMff = VGAMff,
                           object = object,
                           subset = subset,
                           ...)

  etamat <- predict(object)  # nnn x M
  nnn <- nrow(etamat)
  M   <- ncol(etamat)
  Lowmat <- object@misc$Lower
  Uppmat <- object@misc$Upper
  censoL <- object@extra$censoredL
  censoU <- object@extra$censoredU

  mustar <- fitted(object, type.fitted = "uncensored")
  M1 <- npred(object, type = "one.response")
  NOS <- M / M1
  sigmahat <- matrix(NA_real_, nnn, NOS)
  for (jay in seq(NOS))
    sigmahat[, jay] <-
      eta2theta(etamat[, M1 * jay],
                link = object@misc$link[[M1 * jay]],
                earg = object@misc$earg[[M1 * jay]])

  ymat <- depvar(object)
  cnymat <- colnames(ymat)
  rnymat <- rownames(ymat)
  ymat[censoL] <- Lowmat[censoL]
  ymat[censoU] <- Uppmat[censoU]
  zedd <- (ymat - mustar) / sigmahat  # nnn x NOS
  coefmat <- coef(object, matrix = TRUE)[, c(TRUE, FALSE),
                                         drop = FALSE]
  ppp <- nrow(coefmat)

  betamat <- matrix(c(coefmat), nnn, ppp * NOS, byrow = TRUE)
  ans <- betamat  *
         kronecker(dnorm(zedd) / sigmahat, matrix(1, 1, ppp))

  uncens <- !censoL && !censoU
  ans[uncens] <- ans[uncens] * zedd[uncens] / sigmahat[uncens]

  ans[censoL] <- -ans[censoL]


  dim(ans) <- c(nnn, ppp, NOS)
  ans <- aperm(ans, c(2, 3, 1))  # ppp x NOS x nnn
  dimnames(ans) <- list(rownames(coefmat), cnymat, rnymat)

  ans
})  # "tobit"











setMethod("margeffS4VGAM",  signature(VGAMff = "cumulative"),
  function(object,
           subset = NULL,
           VGAMff,
           ...) {


  object <- callNextMethod(VGAMff = VGAMff,
                           object = object,
                           subset = subset,
                           ...)
  reverse <- object@post$reverse
  linkfunctions <- object@post$linkfunctions
  all.eargs <- object@post$all.eargs
  Bmat <- cfit <- object@post$Bmat
  ppp <- object@post$ppp
  etamat <- predict(object)  # nnn x M
  fitmat <- fitted(object)   # nnn x (M + 1)
  nnn <- nrow(etamat)
  M   <- ncol(etamat)
  hdot <- object@post$hdot
  Thetamat <- object@post$Thetamat




    hdot.big <- kronecker(hdot, matrix(1, ppp, 1))  # Enlarged
    resmat <- cbind(hdot.big, 1)
    resmat[, 1] <- ifelse(reverse, -1, 1) *
                   hdot.big[, 1] * cfit[, 1]

    if (M > 1) {
      for (jlocal in 2:M) {
        resmat[, jlocal] <- ifelse(reverse, -1, 1) *
          (hdot.big[, jlocal    ] * cfit[, jlocal    ] -
           hdot.big[, jlocal - 1] * cfit[, jlocal - 1])
      }  # jlocal

    }  # if

    resmat[, M+1] <- ifelse(reverse, 1, -1) *
                     hdot.big[, M] * cfit[, M]

    ans.cum <- array(resmat, c(ppp, nnn, M+1),
               dimnames = list(dimnames(Bmat)[[1]],
                               dimnames(fitted(object))[[1]],
                               dimnames(fitted(object))[[2]]))
    ans.cum <- aperm(ans.cum, c(1, 3, 2))  # ppp x (M+1) x nnn

    subsetarray3(ans.cum, subset = subset)
  })  # "cumulative"







setMethod("margeffS4VGAM",  signature(VGAMff = "acat"),
  function(object,
           subset = NULL,
           VGAMff,
           ...) {


  object <- callNextMethod(VGAMff = VGAMff,
                           object = object,
                           subset = subset,
                           ...)
  reverse <- object@post$reverse
  linkfunctions <- object@post$linkfunctions
  all.eargs <- object@post$all.eargs
  Bmat <- cfit <- object@post$Bmat
  ppp <- object@post$ppp
  etamat <- predict(object)  # nnn x M
  fitmat <- fitted(object)   # nnn x (M + 1)
  nnn <- nrow(etamat)
  M   <- ncol(etamat)
  hdot <- object@post$hdot
  Thetamat <- object@post$Thetamat




    expcs.etamat <- if (reverse)
      exp(tapplymat1(etamat[, M:1, drop = FALSE],
                     "cumsum")[, M:1, drop = FALSE]) else
      exp(tapplymat1(etamat, "cumsum"))
    csexpcs.etavec <- rowSums(expcs.etamat)



    if (!all(object@misc$link == "loglink"))
      stop("currently only the 'loglink' link is supported")


  acat.derivs <- function(jay, tee,
                          M, expcs.etamat, Thetamat,
                          prob1, probMplus1,
                          reverse = FALSE) {

    if (jay > M+1) stop("argument 'jay' out of range")
    if (M   < tee) stop("argument 'tee' out of range")

    if (reverse) {  # ,,,,,,,,,,,,,,,,,,,,,,,,,

      dpMplus1.detat <- -(probMplus1^2) *
                        rowSums(expcs.etamat[, 1:tee,
                                drop = FALSE])
      if (jay == M+1) {
        return(dpMplus1.detat)
      }
      if (jay <= tee) {
          return((probMplus1 + dpMplus1.detat) *
                 expcs.etamat[, jay])
      }
      if (tee < jay) {
        return(dpMplus1.detat * expcs.etamat[, jay])
      }
    } else {  # reverse = FALSE ,,,,,,,,,,,,,,,,,,,,,,,,,,,

      dp1.detat <-
          -(prob1^2) * rowSums(expcs.etamat[, tee:M,
                                            drop = FALSE])
      if (jay == 1) {
        return(dp1.detat)
      }
      if (jay <= tee) {
        return(dp1.detat * expcs.etamat[, jay-1])
      }
      if (tee < jay) {
        return((prob1 + dp1.detat) * expcs.etamat[, jay-1])
      }
    } # reverse ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
  }  # acat.derivs



  A        <- array(0, c(i = nnn, vars = ppp, probs = M + 1,
                         etas = M))
  ansarray <- array(0, c(vars = ppp, i = nnn, probs = M + 1))
  if (reverse) {
    probMplus1 <- 1 / (1 +  csexpcs.etavec) # Last level of Y
  } else {
    prob1 <- 1 / (1 +  csexpcs.etavec)  # First level of Y
  }

    for (jlocal in 1:(M+1)) {
      for (tlocal in 1:M) {
        A[, , jlocal, tlocal] <-
          acat.derivs(jay = jlocal, tee = tlocal,
                      M = M, expcs.etamat = expcs.etamat,
                      Thetamat = Thetamat,
                      prob1 = prob1, probMplus1 = probMplus1,
                      reverse = reverse)
      }
    }


    A <- aperm(A, c(2, 1, 3, 4))  # c(ppp, nnn, M+1, M)
    for (jlocal in 1:(M + 1)) {
      for (tlocal in 1:M) {
        ansarray[,, jlocal]  <- ansarray[,, jlocal] +
                     A[,, jlocal, tlocal] * Bmat[, tlocal]
      }
    }
    ans.acat <- aperm(ansarray, c(1, 3, 2))  # c(ppp,M+1,nnn)
    dimnames(ans.acat) <- list(rownames(Bmat),
                               colnames(fitmat),
                               rownames(etamat))
    subsetarray3(ans.acat, subset = subset)
  })  # "acat"





 cratio.derivs <-
    function(jay, tee,
             hdot, M, cpThetamat, Thetamat,
             reverse = FALSE) {

    if (jay >= M+1) stop("argument 'jay' out of range")
    if (M   <  tee) stop("argument 'tee' out of range")

    if (reverse) {  # ,,,,,,,,,,,,,,,,,,,,,,,,,,
      if (jay == 1) {
        return(hdot[, tee] * cpThetamat[, 1] / Thetamat[, tee])
      }

      if (jay-1 == tee) {
        return(-hdot[, jay-1] * cpThetamat[, jay])
      }
      if (jay <= tee) {
        return((1 - Thetamat[, jay-1]) *
               hdot[, tee] *
               cpThetamat[, jay] / Thetamat[, tee])
      }
      return(rep_len(0, nrow(Thetamat)))  # Since jay-1 > tee
    } else {  # reverse = FALSE ,,,,,,,,,,,,,,,,,,,,,,,,,,,,

      if (jay == 1 && tee == 1) {
        return(-hdot[, 1])
      }

      if (jay == tee) {
        return(-hdot[, jay] * cpThetamat[, jay-1])
      }
      if (tee < jay) {
        return((1 - Thetamat[, jay]) *
               hdot[, tee] *
               cpThetamat[, jay-1] / Thetamat[, tee])
      }
      return(rep_len(0, nrow(Thetamat)))  # Since jay < tee
    } # reverse ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
  }  # cratio.derivs






setMethod("margeffS4VGAM",  signature(VGAMff = "cratio"),
  function(object,
           subset = NULL,
           VGAMff,
           ...) {


  object <- callNextMethod(VGAMff = VGAMff,
                           object = object,
                           subset = subset,
                           ...)
  reverse <- object@post$reverse
  linkfunctions <- object@post$linkfunctions
  all.eargs <- object@post$all.eargs
  Bmat <- cfit <- object@post$Bmat
  ppp <- object@post$ppp
  etamat <- predict(object)  # nnn x M
  fitmat <- fitted(object)   # nnn x (M + 1)
  nnn <- nrow(etamat)
  M   <- ncol(etamat)
  hdot <- object@post$hdot
  Thetamat <- object@post$Thetamat







  vfamily <- object@family@vfamily
  c.nots <- any(vfamily == "cratio")

  if (any(vfamily == "cratio")) {
    cpThetamat <- if (reverse)
      tapplymat1(    Thetamat[, M:1, drop = FALSE],
                 "cumprod")[, M:1, drop = FALSE] else
      tapplymat1(    Thetamat, "cumprod")
  }



  A        <- array(0, c(i = nnn, vars = ppp, probs = M + 1,
                         etas = M))
  ansarray <- array(0, c(vars = ppp, i = nnn, probs = M + 1))


  choosemat <- if (c.nots) Thetamat else 1 - Thetamat
  if (min(choosemat) <= 0)
    warning("division by 0 may occur")






    if (reverse) {
      for (tlocal in 1:M) {
        for (jlocal in 1:tlocal) {
          A[, , jlocal, tlocal] <-
            cratio.derivs(jay = jlocal, tee = tlocal,
                          hdot = ifelse(c.nots, 1, -1) * hdot,
                          M = M, cpThetamat = cpThetamat,
                          Thetamat = choosemat,
                          reverse = reverse)
        }
      }
      if (M > 1)
        for (jlocal in 2:M) {
          A[, , jlocal, jlocal-1] <-
            cratio.derivs(jay = jlocal, tee = jlocal-1,
                          hdot = ifelse(c.nots, 1, -1) * hdot,
                          M = M, cpThetamat = cpThetamat,
                          Thetamat = choosemat,
                          reverse = reverse)
        }
    } else {
     for (jlocal in 1:M) {
        for (tlocal in 1:jlocal) {
          A[, , jlocal, tlocal] <-
            cratio.derivs(jay = jlocal, tee = tlocal,
                          hdot = ifelse(c.nots, 1, -1) * hdot,
                          M = M, cpThetamat = cpThetamat,
                          Thetamat = choosemat,
                          reverse = reverse)
        }
      }
    }

    if (reverse) {
      A[, , M+1, M] <- ifelse(c.nots, -1, 1) * hdot[, M]
    } else {
      for (jlocal in 1:M) {
        for (tlocal in 1:jlocal) {
          A[, , M+1, tlocal] <- if (c.nots) {
            A[, , M+1, tlocal] - A[, , jlocal, tlocal]
          } else {
              -hdot[, tlocal] *
                  cpThetamat[, M] / choosemat[, tlocal]
          }
        }
      }
    }

    A <- aperm(A, c(2, 1, 3, 4))  # c(ppp, nnn, M+1, M)
    for (jlocal in 1:(M + 1)) {
      for (tlocal in 1:M) {
        ansarray[,, jlocal]  <- ansarray[,, jlocal] +
                          A[,, jlocal, tlocal] * Bmat[, tlocal]
      }
    }
    ans.csratio <- aperm(ansarray, c(1, 3, 2))  # c(ppp,M+1,nnn)
    dimnames(ans.csratio) <- list(rownames(Bmat),
                                  colnames(fitmat),
                                  rownames(etamat))
    subsetarray3(ans.csratio,
                 subset = subset)  # "cratio" & "sratio"
  })  # "cratio"






setMethod("margeffS4VGAM",  signature(VGAMff = "sratio"),
  function(object,
           subset = NULL,
           VGAMff,
           ...) {








  object <- callNextMethod(VGAMff = VGAMff,
                           object = object,
                           subset = subset,
                           ...)
  reverse <- object@post$reverse
  linkfunctions <- object@post$linkfunctions
  all.eargs <- object@post$all.eargs
  Bmat <- cfit <- object@post$Bmat
  ppp <- object@post$ppp
  etamat <- predict(object)  # nnn x M
  fitmat <- fitted(object)   # nnn x (M + 1)
  nnn <- nrow(etamat)
  M   <- ncol(etamat)
  hdot <- object@post$hdot
  Thetamat <- object@post$Thetamat




  vfamily <- object@family@vfamily
  c.nots <- any(vfamily == "cratio")
  if (any(vfamily == "sratio")) {
    cpThetamat <- if (reverse)
      tapplymat1(1 - Thetamat[, M:1, drop = FALSE],
                 "cumprod")[, M:1, drop = FALSE] else
      tapplymat1(1 - Thetamat, "cumprod")
  }



  A        <- array(0, c(i = nnn, vars = ppp, probs = M + 1,
                         etas = M))
  ansarray <- array(0, c(vars = ppp, i = nnn, probs = M + 1))


  choosemat <- if (c.nots) Thetamat else 1 - Thetamat
  if (min(choosemat) <= 0)
    warning("division by 0 may occur")











    if (reverse) {
      for (tlocal in 1:M) {
        for (jlocal in 1:tlocal) {
          A[, , jlocal, tlocal] <-
            cratio.derivs(jay = jlocal, tee = tlocal,
                          hdot = ifelse(c.nots, 1, -1) * hdot,
                          M = M, cpThetamat = cpThetamat,
                          Thetamat = choosemat,
                          reverse = reverse)
        }
      }
      if (M > 1)
        for (jlocal in 2:M) {
          A[, , jlocal, jlocal-1] <-
            cratio.derivs(jay = jlocal, tee = jlocal-1,
                          hdot = ifelse(c.nots, 1, -1) * hdot,
                          M = M, cpThetamat = cpThetamat,
                          Thetamat = choosemat,
                          reverse = reverse)
        }
    } else {
     for (jlocal in 1:M) {
        for (tlocal in 1:jlocal) {
          A[, , jlocal, tlocal] <-
            cratio.derivs(jay = jlocal, tee = tlocal,
                          hdot = ifelse(c.nots, 1, -1) * hdot,
                          M = M, cpThetamat = cpThetamat,
                          Thetamat = choosemat,
                          reverse = reverse)
        }
      }
    }

    if (reverse) {
      A[, , M+1, M] <- ifelse(c.nots, -1, 1) * hdot[, M]
    } else {
      for (jlocal in 1:M) {
        for (tlocal in 1:jlocal) {
          A[, , M+1, tlocal] <- if (c.nots) {
            A[, , M+1, tlocal] - A[, , jlocal, tlocal]
          } else {
              -hdot[, tlocal] *
                  cpThetamat[, M] / choosemat[, tlocal]
          }
        }
      }
    }

    A <- aperm(A, c(2, 1, 3, 4))  # c(ppp, nnn, M+1, M)
    for (jlocal in 1:(M + 1)) {
      for (tlocal in 1:M) {
        ansarray[,, jlocal]  <- ansarray[,, jlocal] +
                        A[,, jlocal, tlocal] * Bmat[, tlocal]
      }
    }
    ans.csratio <- aperm(ansarray, c(1, 3, 2))  # c(ppp,M+1,nnn)
    dimnames(ans.csratio) <- list(rownames(Bmat),
                                  colnames(fitmat),
                                  rownames(etamat))
    subsetarray3(ans.csratio,
                 subset = subset)  # "cratio" & "sratio"
  })  # "sratio"








 margefff <- function(object, subset = NULL) {


  ii <- subset
  if (!is(object, "vglm"))
    stop("'object' is not a vglm() object")
  if (!any(temp.logical <-
    is.element(c("multinomial", "cumulative", "acat",
                 "cratio", "sratio"),
               object@family@vfamily)))
    stop("'object' is not a 'multinomial' or 'acat' ",
         "or 'cumulative'  or 'cratio' or 'sratio' VGLM!")
  vfamily <- object@family@vfamily
  if (is(object, "vgam"))
    stop("'object' is a vgam() object")
  if (length(object@control$xij))
    stop("'object' contains 'xij' terms")
  if (length(object@misc$form2))
    stop("'object' contains 'form2' terms")

  oassign <- object@misc$orig.assign
  if (any(unlist(lapply(oassign, length)) > 1))
    warning("some terms in 'object' create more than one ",
            "column of the LM design matrix")

  nnn <- object@misc$n
  M <- object@misc$M  # ncol(B)  # length(pvec) - 1


    if (any(vfamily == "multinomial")) {
    rlev <- object@misc$refLevel
    cfit <- coefvlm(object, matrix.out = TRUE)
    B <- if (!length(rlev)) {
      cbind(cfit, 0)
    } else {
      if (rlev == M+1) {  # Default
        cbind(cfit, 0)
      } else if (rlev == 1) {
        cbind(0, cfit)
      } else {
        cbind(cfit[, 1:(rlev-1)], 0, cfit[, rlev:M])
      }
    }
    ppp   <- nrow(B)
    pvec1 <- fitted(object)[1, ]
    colnames(B) <- if (length(names(pvec1))) names(pvec1) else
                   param.names("mu", M+1)

    if (is.null(ii)) {
      BB <- array(B, c(ppp, M+1, nnn))
      pvec  <- c(t(fitted(object)))
      pvec  <- rep(pvec, each = ppp)
      temp1 <- array(BB * pvec, c(ppp, M+1, nnn))
      temp2 <- aperm(temp1, c(2, 1, 3))  # (M+1) x ppp x nnn
      temp2 <- colSums(temp2)  # ppp x nnn
      temp2 <- array(rep(temp2, each = M+1), c(M+1, ppp, nnn))
      temp2 <- aperm(temp2, c(2, 1, 3))  # ppp x (M+1) x nnn
      temp3 <- pvec
      ans <- array((BB - temp2) * temp3, c(ppp, M+1, nnn),
             dimnames = list(dimnames(B)[[1]],
                             dimnames(B)[[2]],
                             dimnames(fitted(object))[[1]]))
      return(ans)
    } else if (is.numeric(ii) && length(ii) == 1) {
        pvec  <- fitted(object)[ii, ]
        temp1 <- B * matrix(pvec, ppp, M+1, byrow = TRUE)
        temp2 <- matrix(rowSums(temp1), ppp, M+1)
        temp3 <- matrix(pvec, nrow(B), M+1, byrow = TRUE)
        return((B - temp2) * temp3)
    } else {
        if (is.logical(ii))
          ii <- (1:nnn)[ii]

        ans <- array(0, c(ppp, M+1, length(ii)),
               dimnames = list(dimnames(B)[[1]],
                               dimnames(B)[[2]],
                    dimnames(fitted(object)[ii, ])[[1]]))
        for (ilocal in seq_along(ii)) {
          pvec  <- fitted(object)[ii[ilocal], ]
          temp1 <- B * matrix(pvec, ppp, M+1, byrow = TRUE)
          temp2 <- matrix(rowSums(temp1), ppp, M+1)
          temp3 <- matrix(pvec, nrow(B), M+1, byrow = TRUE)
          ans[ , , ilocal] <- (B - temp2) * temp3
        }
        return(ans)
      }
  }  # "multinomial"





  reverse <- object@misc$reverse
  linkfunctions <- object@misc$link
  all.eargs <- object@misc$earg
  B <- cfit <- coefvlm(object, matrix.out = TRUE)
  ppp <- nrow(B)
  etamat <- predict(object)  # nnn x M
  fitmat <- fitted(object)   # nnn x (M + 1)
  nnn <- nrow(etamat)


  hdot <- Thetamat <- etamat
  for (jlocal in 1:M) {
    Thetamat[, jlocal] <- eta2theta(etamat[, jlocal],
                                    linkfunctions[jlocal],
                                    all.eargs[[jlocal]])
    hdot[, jlocal] <- dtheta.deta(Thetamat[, jlocal],
                                  linkfunctions[jlocal],
                                  all.eargs[[jlocal]])
  }  # jlocal



  if (any(vfamily == "acat")) {
    expcs.etamat <- if (reverse)
      exp(tapplymat1(etamat[, M:1, drop = FALSE],
                     "cumsum")[, M:1, drop = FALSE]) else
      exp(tapplymat1(etamat, "cumsum"))
    csexpcs.etavec <- rowSums(expcs.etamat)
  }
  if (any(vfamily == "cratio")) {
    cpThetamat <- if (reverse)
      tapplymat1(    Thetamat[, M:1, drop = FALSE],
                 "cumprod")[, M:1, drop = FALSE] else
      tapplymat1(    Thetamat, "cumprod")
  }
  if (any(vfamily == "sratio")) {
    cpThetamat <- if (reverse)
      tapplymat1(1 - Thetamat[, M:1, drop = FALSE],
                 "cumprod")[, M:1, drop = FALSE] else
      tapplymat1(1 - Thetamat, "cumprod")
  }





     if (is.logical(is.multivariateY <-
       object@misc$multiple.responses) && is.multivariateY)
    stop("cannot handle ",
         "cumulative(multiple.responses = TRUE)")




  if (any(vfamily == "cumulative")) {
    hdot.big <- kronecker(hdot, matrix(1, ppp, 1))  # Enlarged
    resmat <- cbind(hdot.big, 1)
    resmat[, 1] <- ifelse(reverse, -1, 1) *
                   hdot.big[, 1] * cfit[, 1]

    if (M > 1) {
      for (jlocal in 2:M)
        resmat[, jlocal] <- ifelse(reverse, -1, 1) *
          (hdot.big[, jlocal    ] * cfit[, jlocal    ] -
           hdot.big[, jlocal - 1] * cfit[, jlocal - 1])

    }  # jlocal

    resmat[, M+1] <- ifelse(reverse, 1, -1) *
        hdot.big[, M] * cfit[, M]

    temp1 <- array(resmat, c(ppp, nnn, M+1),
               dimnames = list(dimnames(B)[[1]],
                               dimnames(fitted(object))[[1]],
                               dimnames(fitted(object))[[2]]))
    temp1 <- aperm(temp1, c(1, 3, 2))  # ppp x (M+1) x nnn

    if (is.null(ii)) {
      return(temp1)
    } else
    if (is.numeric(ii) && (length(ii) == 1)) {
      return(temp1[, , ii])
    } else {
      return(temp1[, , ii])
    }
  }  # "cumulative"









  if (any(vfamily == "acat")) {
    if (!all(object@misc$link == "loglink"))
      stop("currently only the 'loglink' link is supported")


  acat.derivs <- function(jay, tee,
                          M, expcs.etamat, Thetamat,
                          prob1, probMplus1,
                          reverse = FALSE) {

    if (jay > M+1) stop("argument 'jay' out of range")
    if (M   < tee) stop("argument 'tee' out of range")

    if (reverse) {  # ,,,,,,,,,,,,,,,,,,,,,,,,,

      dpMplus1.detat <- -(probMplus1^2) *
              rowSums(expcs.etamat[, 1:tee, drop = FALSE])
      if (jay == M+1) {
        return(dpMplus1.detat)
      }
      if (jay <= tee) {
          return((probMplus1 + dpMplus1.detat) *
                 expcs.etamat[, jay])
      }
      if (tee < jay) {
        return(dpMplus1.detat * expcs.etamat[, jay])
      }
    } else {  # reverse = FALSE ,,,,,,,,,,,,,,,,,,

      dp1.detat <-
          -(prob1^2) * rowSums(expcs.etamat[, tee:M,
                                            drop = FALSE])
      if (jay == 1) {
        return(dp1.detat)
      }
      if (jay <= tee) {
        return(dp1.detat * expcs.etamat[, jay-1])
      }
      if (tee < jay) {
        return((prob1 + dp1.detat) * expcs.etamat[, jay-1])
      }
    } # reverse ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
  }  # acat.derivs


  A        <- array(0, c(i = nnn, vars = ppp, probs = M + 1,
                         etas = M))
  ansarray <- array(0, c(vars = ppp, i = nnn, probs = M + 1))
  if (reverse) {
    probMplus1 <- 1 / (1 +  csexpcs.etavec) # Last level of Y
  } else {
    prob1 <- 1 / (1 +  csexpcs.etavec)  # First level of Y
  }

    for (jlocal in 1:(M+1)) {
      for (tlocal in 1:M) {
        A[, , jlocal, tlocal] <-
          acat.derivs(jay = jlocal, tee = tlocal,
                      M = M, expcs.etamat = expcs.etamat,
                      Thetamat = Thetamat,
                      prob1 = prob1, probMplus1 = probMplus1,
                      reverse = reverse)
      }
    }


    A <- aperm(A, c(2, 1, 3, 4))  # c(ppp, nnn, M+1, M)
    for (jlocal in 1:(M + 1)) {
      for (tlocal in 1:M) {
        ansarray[,, jlocal]  <- ansarray[,, jlocal] +
                  A[,, jlocal, tlocal] * B[, tlocal]
      }
    }
    ans.acat <- aperm(ansarray, c(1, 3, 2))  # c(ppp, M+1, nnn)
    dimnames(ans.acat) <- list(rownames(B),
                               colnames(fitmat),
                               rownames(etamat))
    return(ans.acat)
  }  # "acat"






  c.nots <- any(vfamily == "cratio")

     cratio.derivs <-
         function(jay, tee,
                  hdot, M, cpThetamat, Thetamat,
                  reverse = FALSE) {

    if (jay >= M+1) stop("argument 'jay' out of range")
    if (M   <  tee) stop("argument 'tee' out of range")

    if (reverse) {  # ,,,,,,,,,,,,,,,,,,,,,,,,,
      if (jay == 1) {
          return(hdot[, tee] *
                 cpThetamat[, 1] / Thetamat[, tee])
      }

      if (jay-1 == tee) {
        return(-hdot[, jay-1] * cpThetamat[, jay])
      }
      if (jay <= tee) {
        return((1 - Thetamat[, jay-1]) *
               hdot[, tee] *
               cpThetamat[, jay] / Thetamat[, tee])
      }
      return(rep_len(0, nrow(Thetamat)))  # Since jay-1 > tee
    } else {  # reverse = FALSE ,,,,,,,,,,,,,,,,,

      if (jay == 1 && tee == 1) {
        return(-hdot[, 1])
      }

      if (jay == tee) {
        return(-hdot[, jay] * cpThetamat[, jay-1])
      }
      if (tee < jay) {
        return((1 - Thetamat[, jay]) *
        hdot[, tee] * cpThetamat[, jay-1] / Thetamat[, tee])
      }
      return(rep_len(0, nrow(Thetamat)))  # Since jay < tee
    } # reverse ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
  }  # cratio.derivs


  A        <- array(0, c(i = nnn, vars = ppp, probs = M + 1,
                         etas = M))
  ansarray <- array(0, c(vars = ppp, i = nnn, probs = M + 1))


  choosemat <- if (c.nots) Thetamat else 1 - Thetamat
  if (min(choosemat) <= 0)
    warning("division by 0 may occur")




  if (any(vfamily == "cratio" | vfamily == "sratio")) {


    if (reverse) {
      for (tlocal in 1:M) {
        for (jlocal in 1:tlocal) {
          A[, , jlocal, tlocal] <-
            cratio.derivs(jay = jlocal, tee = tlocal,
                          hdot = ifelse(c.nots, 1, -1) * hdot,
                          M = M, cpThetamat = cpThetamat,
                          Thetamat = choosemat,
                          reverse = reverse)
        }
      }
      if (M > 1)
        for (jlocal in 2:M) {
          A[, , jlocal, jlocal-1] <-
            cratio.derivs(jay = jlocal, tee = jlocal-1,
                          hdot = ifelse(c.nots, 1, -1) * hdot,
                          M = M, cpThetamat = cpThetamat,
                          Thetamat = choosemat,
                          reverse = reverse)
        }
    } else {
     for (jlocal in 1:M) {
        for (tlocal in 1:jlocal) {
          A[, , jlocal, tlocal] <-
            cratio.derivs(jay = jlocal, tee = tlocal,
                          hdot = ifelse(c.nots, 1, -1) * hdot,
                          M = M, cpThetamat = cpThetamat,
                          Thetamat = choosemat,
                          reverse = reverse)
        }
      }
    }

    if (reverse) {
      A[, , M+1, M] <- ifelse(c.nots, -1, 1) * hdot[, M]
    } else {
      for (jlocal in 1:M) {
        for (tlocal in 1:jlocal) {
          A[, , M+1, tlocal] <- if (c.nots) {
            A[, , M+1, tlocal] - A[, , jlocal, tlocal]
          } else {
    -hdot[, tlocal] * cpThetamat[, M] / choosemat[, tlocal]
          }
        }
      }
    }

    A <- aperm(A, c(2, 1, 3, 4))  # c(ppp, nnn, M+1, M)
    for (jlocal in 1:(M + 1)) {
      for (tlocal in 1:M) {
        ansarray[,, jlocal]  <- ansarray[,, jlocal] +
                       A[,, jlocal, tlocal] * B[, tlocal]
      }
    }
    ans.csratio <- aperm(ansarray, c(1, 3, 2))  # c(ppp,M+1,nnn)
    dimnames(ans.csratio) <- list(rownames(B),
                                  colnames(fitmat),
                                  rownames(etamat))
    return(ans.csratio)
  }  # "cratio" and "sratio"


}  # margefff








prplot <- function(object,
                  control = prplot.control(...), ...) {


  if (!any(slotNames(object) == "family") ||
      !any(object@family@vfamily == "VGAMcategorical"))
    stop("'object' does not seem to be a VGAM categorical ",
         "model object")

  if (!any(object@family@vfamily == "cumulative"))
    stop("'object' does not seem to be a VGAM categorical ",
         "model object")

    control <- prplot.control(...)


  object <- plot.vgam(object, plot.arg = FALSE,
                      raw = FALSE)  # , ...

  if (length(names(object@preplot)) != 1)
      stop("object needs to have only one term")


  MM <- object@misc$M
  use.y <- cbind((object@preplot[[1]])$y)
  Constant <- attr(object@preplot, "Constant")
  if (is.numeric(Constant) && length(Constant) == ncol(use.y))
    use.y <- use.y + matrix(Constant, nrow(use.y), ncol(use.y),
                            byrow = TRUE)
  for (ii in 1:MM) {
    use.y[, ii] <- eta2theta(use.y[, ii],
                             link = object@misc$link[[ii]],
                             earg = object@misc$earg[[ii]])
  }
  if (ncol(use.y) != MM) use.y = use.y[, 1:MM, drop = FALSE]

  use.x <- (object@preplot[[1]])$x
  myxlab <- if (length(control$xlab))
           control$xlab else (object@preplot[[1]])$xlab
  mymain <- if (MM <= 3)
    paste(object@misc$parameters,
          collapse = ", ") else
    paste(object@misc$parameters[c(1, MM)],
          collapse = ",...,")
  if (length(control$main)) mymain = control$main
  if (length(control$ylab)) myylab = control$ylab

  matplot(use.x, use.y, type = "l",
          xlab = myxlab, ylab = myylab,
          lty = control$lty,
          col = control$col, las = control$las,
          xlim = if (is.Numeric(control$xlim))
                 control$xlim else range(use.x),
          ylim = if (is.Numeric(control$ylim))
                 control$ylim else range(use.y),
          main=mymain)
  if (control$rug.arg)
    rug(use.x, col=control$rcol, lwd=control$rlwd)

  invisible(object)
}  # prplot




 prplot.control <-
  function(xlab = NULL, ylab = "Probability",
           main = NULL,
           xlim = NULL, ylim = NULL,
           lty = par()$lty,
           col = par()$col,
           rcol = par()$col,
           lwd = par()$lwd,
           rlwd = par()$lwd,
           las = par()$las,
           rug.arg  = FALSE,
           ...) {

    list(xlab = xlab, ylab = ylab,
         xlim = xlim, ylim = ylim,
         lty = lty, col = col, rcol = rcol,
         lwd = lwd, rlwd = rlwd, rug.arg = rug.arg,
         las = las, main = main)
}










is.parallel.matrix <- function(object, ...)
  is.matrix(object) && all(!is.na(object)) &&
  all(c(object) == 1) && ncol(object) == 1


is.parallel.vglm <-
  function(object, type = c("term", "lm"), ...) {

  type <- match.arg(type, c("term", "lm"))[1]
  Hlist <- constraints(object, type = type)

  unlist(lapply(Hlist, is.parallel.matrix))
}


if (!isGeneric("is.parallel"))
  setGeneric("is.parallel", function(object, ...)
             standardGeneric("is.parallel"),
             package = "VGAM")


setMethod("is.parallel",  "matrix", function(object, ...)
    is.parallel.matrix(object, ...))


setMethod("is.parallel",  "vglm", function(object, ...)
    is.parallel.vglm(object, ...))






is.zero.matrix <- function(object, ...) {

  rnames <- rownames(object)
  intercept.index <- if (length(rnames)) {
    if (any(rnames == "(Intercept)")) {
      (seq_along(rnames))[rnames == "(Intercept)"]
    } else {
      stop("the matrix doesnt have an intercept")
      NULL
    }
  } else {
      stop("the matrix doesnt have an intercept")
      NULL
  }

  if (nrow(object) <= 1)
    stop("the matrix needs to have more than one row",
         ", i.e., more than ",
         "an intercept on the RHS of the formula")

  cfit <- object[-intercept.index, , drop = FALSE]

  foo <- function(conmat.col)
    all(!is.na(conmat.col)) &&
    all(c(conmat.col) == 0)

  unlist(apply(cfit, 2, foo))
}  # is.zero.matrix



is.zero.vglm <- function(object, ...) {
  is.zero.matrix(coef(object, matrix = TRUE))
}



if (!isGeneric("is.zero"))
  setGeneric("is.zero", function(object, ...)
             standardGeneric("is.zero"),
             package = "VGAM")


setMethod("is.zero", "matrix", function(object, ...)
          is.zero.matrix(object, ...))



setMethod("is.zero", "NULL",
          function(object, ...)
          is.null(object))
setMethod("is.zero", "character",
          function(object, ...)
              is.character(object) &&
              length(object) == 1 &&
              object == "")
setMethod("is.zero", "logical",
          function(object, ...)
    length(object) == 1 &&
    is.na(object))  # Based on class(NA)





setMethod("is.zero",  "vglm", function(object, ...)
          is.zero.vglm(object, ...))






setMethod("showvglmS4VGAM",
          signature(VGAMff = "acat"),
  function(object,
           VGAMff,
           ...) {
  cat("\nThis is an adjacent categories model with",
      1 + object@misc$M, "levels\n")
  invisible(object)
  })


setMethod("showvgamS4VGAM",
          signature(VGAMff = "acat"),
  function(object,
           VGAMff,
           ...) {
  cat("\nThis is an adjacent categories model with",
      1 + object@misc$M, "levels\n")
  invisible(object)
  })



setMethod("showvglmS4VGAM",
          signature(VGAMff = "multinomial"),
  function(object,
           VGAMff,
           ...) {
  cat("\nThis is a multinomial logit model with",
      1 + object@misc$M, "levels\n")
  invisible(object)
  })


setMethod("showvgamS4VGAM",
          signature(VGAMff = "multinomial"),
  function(object,
           VGAMff,
           ...) {
  cat("\nThis is a multinomial logit model with",
      1 + object@misc$M, "levels\n")
  invisible(object)
  })






R2latvar <- function(object) {


  if (!is(object, "vglm"))
      stop("argument 'object' is not a vglm() fit")
 

  vfam <- object@family@vfamily
  fam.permitted <- c("cumulative", "propodds",
                     "binomialff")
  if (!(vfam[1] %in% fam.permitted))
      stop("the family function must be one of ",
           fam.permitted)


  linkfn <- linkfun(object)[1]
  link.permitted <- c("logitlink", "probitlink",
                      "clogloglink")  # "cauchitlink"
  if (!(linkfn %in% link.permitted))
      stop("allowable link functions supported are ",
           link.permitted)

  infos <- object@family@infos()
  if (length(infos$parallel) == 1 &&
      is.logical(infos$parallel))
    if (!infos$parallel)
      stop("the linear predictors are not parallel")
  if (!all(unlist(constraints(object)[-1]) == 1))
      stop("the linear predictors are not parallel")



  eta1 <- predict(object)[, 1]
  offset <- switch(linkfn,
                   logitlink   = (pi^2)/3,
                   probitlink  = 1,
                   clogloglink = (pi^2) / 6,
                   stop("link unrecognized"))

  veta1 <- var(eta1)
  veta1 / (veta1 + offset) 
}  # R2latvar








ordsup.vglm <-
  function(object, all.vars = FALSE, confint = FALSE, ...) {


  if (!is(object, "vglm"))
    stop("argument 'object' is not a vglm() fit")
 

  vfam <- object@family@vfamily
  fam.permitted <- c("uninormal", "propodds", "cumulative")
  if (!(vfam[1] %in% fam.permitted))
    stop("the family function must be one of ",
         paste(fam.permitted, collapse = ", "))


  linkfns <- linkfun(object)
  link.permitted <- c("identitylink", "logitlink", "probitlink")
  if (!any(linkfns %in% link.permitted))
    stop("allowable link functions supported are ",
         paste("'", link.permitted, "'",
               sep = "", collapse = ", "))


  R <- npred(object) / npred(object, type = "one.response")      
  if (R > 1)
    stop("currently cannot handle multiple responses")


  infos <- object@family@infos()
  cobj <- coef(object)
  cobj.mat <- coef(object, matrix = TRUE)
  x.LM  <- model.matrix(object, type =  "lm")
  is.binary <- apply(if (has.intercept(object))
                     x.LM[, -1, drop = FALSE] else x.LM,
                     2, function(coln)
                         length(unique(coln)) == 2)

  if (all.vars) {
    is.binary[] <- TRUE
  } else {
    if (sum(is.binary) == 0) {
        warning("no binary explanatory variables; ",
                "returning a NULL")
      return(NULL)
    }
  }


  is.binary <- names(is.binary)[is.binary]


      

  switch(vfam[1],
  uninormal = {
    if (!all(linkfns[c(TRUE, FALSE)] == "identitylink"))
      stop("must use 'identitylink' for the mean parameter")
    if (!all(cobj.mat[-1, c(FALSE, TRUE)] == 0))
      stop("the sdev or var parameter must be intercept-only")
    sdev <- numeric(ncol(cobj.mat) / 2)
    for (jay in 1:(ncol(cobj.mat) / 2))
      sdev[jay] <- eta2theta(cobj.mat["(Intercept)", 2*jay],
                             link = object@misc$link[2*jay],
                             earg = object@misc$earg[2*jay])
    if (infos$parameters.names[2] == "var")
      sdev <- sqrt(sdev)  # It was actually the variance before
  },
  cumulative = {

  if (length(infos$parallel) == 1 &&
      is.logical(infos$parallel))
    if (!infos$parallel)
      stop("the linear predictors are not parallel")
  if (!all(unlist(constraints(object)[-1]) == 1))
      stop("the linear predictors are not parallel")


  reverse <- infos$reverse
  if (FALSE)
  if (!reverse) {
    stop("the 'reverse' argument must be TRUE")
  }
  },
  {
  stop("family unrecognized")
  })



  gamma.fun <- function(cobj, is.binary, vfam,
                        reverse, linkfns = NULL) {
  switch(vfam[1],
  uninormal = {
    gamma <- pnorm(cobj[ is.binary ] / (sqrt(2) * sdev))
  },
  cumulative = {
    gamma <-
    switch(linkfns[1],
  clogloglink = clogloglink(ifelse(reverse, 1, -1) *  # Unsure
                            cobj[ is.binary ], inverse = TRUE),
    logitlink   =   logitlink(ifelse(reverse, 1, -1) *
                              cobj[ is.binary ] / sqrt(2),
                              inverse = TRUE),
    probitlink  =  probitlink(ifelse(reverse, 1, -1) *
                              cobj[ is.binary ] / sqrt(2),
                              inverse = TRUE))
  },
  zzzz = {
  })
  gamma
  }


      

  gamma <- gamma.fun(cobj, is.binary, vfam, reverse, linkfns)
  Delta <- 2 * gamma - 1

      

  

  if (confint) {
    ans2 <- confint(object, parm = is.binary, ...)
    if (!is.matrix(ans2))
      ans2 <- matrix(ans2, 1, 2,
                     dimnames = list(is.binary, names(ans2)))
    ans2.low <- gamma.fun(ans2[, 1], TRUE,
                          vfam, reverse, linkfns)
    ans2.upp <- gamma.fun(ans2[, 2], TRUE,
                          vfam, reverse, linkfns)
    Delta.low <- 2 * ans2.low - 1
    Delta.upp <- 2 * ans2.upp - 1
    names(ans2.low) <- names(ans2.upp) <-
    names(Delta.low) <- names(Delta.upp) <- is.binary
  }


  c(
  list(gamma = gamma,
       Delta = Delta),
       if (confint) {
         list(lower.gamma = ans2.low,
              upper.gamma = ans2.upp,
              Lower.Delta = Delta.low,
              Upper.Delta = Delta.upp)
       } else NULL)
}  # ordsup.vglm



 if (!isGeneric("ordsup"))
     setGeneric("ordsup",
  function(object, ...) standardGeneric("ordsup"))

setMethod("ordsup",  "vglm",
  function(object, ...) ordsup.vglm(object, ...))

Try the VGAM package in your browser

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

VGAM documentation built on Sept. 18, 2024, 9:09 a.m.