R/family.binomial.R

Defines functions eimij.ebbinom rextbetabinom qextbetabinom pextbetabinom dextbetabinom extbetabinomial binom2.rho.ss binom2.Rho zipebcom seq2binomial betageometric betabinomialff betabinomialff.control Ebbin.ab rbetabinom pbetabinom dbetabinom rbetabinom.ab pbetabinom.ab dbetabinom.ab size.binomial my.dbinom pnorm2 pbinorm dnorm2 binom2.rho binom2.rho.control rbinom2.rho dbinom2.rho binom2.or rbinom2.or dbinom2.or betabinomial betabinomial.control

Documented in betabinomial betabinomialff betageometric binom2.or binom2.rho binom2.Rho binom2.rho.ss dbetabinom dbetabinom.ab dbinom2.or dbinom2.rho dextbetabinom dnorm2 dnorm2 Ebbin.ab eimij.ebbinom extbetabinomial my.dbinom pbetabinom pbetabinom.ab pbinorm pextbetabinom pnorm2 pnorm2 qextbetabinom rbetabinom rbetabinom.ab rbinom2.or rbinom2.rho rextbetabinom seq2binomial size.binomial zipebcom

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












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



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


  if (!is.matrix(y)) {
    yf <- as.factor(y)
    lev <- levels(yf)
    llev <- length(lev)
    if (llev != 4)
        stop("response must have 4 levels")
    nn <- length(yf)
    y <- matrix(0, nn, llev)
    y[cbind(1:nn, as.vector(unclass(yf)))] <- 1
    colnamesy <- paste(lev, ":", c("00", "01", "10", "11"), sep = "")
    dimnames(y) <- list(names(yf), colnamesy)
    input.type <- 1
  } else if (ncol(y) == 2) {
    if (!all(y == 0 | y == 1))
      stop("response must contains 0's and 1's only")
    col.index <- y[, 2] + 2*y[, 1] + 1    # 1:4
    nn <- nrow(y)
    y <- matrix(0, nn, 4)
    y[cbind(1:nn, col.index)] <- 1
    dimnames(y) <- list(dimnames(y)[[1]],
                        c("00", "01", "10", "11"))
    input.type <- 2
  } else if (ncol(y) == 4) {
    input.type <- 3
  } else
    stop("response unrecognized")



  nvec <- rowSums(y)

  w <- w * nvec
  y <- y / nvec   # Convert to proportions

  if (length(mustart) + length(etastart) == 0) {
    mu <- y + (1 / ncol(y) - y) / nvec
    dimnames(mu) <- dimnames(y)
    mustart <- mu
  }
})  # process.binomial2.data.VGAM






betabinomial.control <-
    function(save.weights = TRUE, ...) {
  list(save.weights = save.weights)
}






 betabinomial <-
  function(lmu = "logitlink",
           lrho = "logitlink",
           irho = NULL,
           imethod = 1, ishrinkage = 0.95,
           nsimEIM = NULL, zero = "rho") {
  lmu <- as.list(substitute(lmu))
  emu <- link2list(lmu)
  lmu <- attr(emu, "function.name")

  lrho <- as.list(substitute(lrho))
  erho <- link2list(lrho)
  lrho <- attr(erho, "function.name")



  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE,
                  positive = TRUE) ||
     imethod > 4)
    stop("argument 'imethod' must be 1, 2, 3 or 4")

  if (!is.Numeric(ishrinkage, length.arg = 1) ||
      ishrinkage < 0 ||
      ishrinkage > 1)
    stop("bad input for argument 'ishrinkage'")
  if (!is.null(nsimEIM)) {
    if (!is.Numeric(nsimEIM, length.arg = 1,
                    integer.valued = TRUE))
      stop("bad input for argument 'nsimEIM'")
    if (nsimEIM <= 10)
      warning("'nsimEIM' should be an integer",
              " greater than 10, say")
  }

  new("vglmff",
  blurb = c("Beta-binomial model\n",
            "Links:      ",
            namesof("mu",  lmu,  emu), ", ",
            namesof("rho", lrho, erho), "\n",
            "Mean:     mu", "\n",
            "Variance: mu*(1-mu)*(1+(w-1)*rho)/w"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints,
              x = x, .zero , M = M, M1 = 3,
              predictors.names = predictors.names)
  }), list( .zero = zero ))),


  infos = eval(substitute(function(...) {
    list(M1 = 3,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("mu", "rho"),
         imethod  = .imethod ,
         ishrinkage  = .ishrinkage ,
         nsimEIM  = .nsimEIM ,
         lmu  = .lmu ,
         lrho = .lrho ,
         zero = .zero )
  },
  list( .lmu = lmu, .lrho = lrho,
        .imethod = imethod,
        .ishrinkage = ishrinkage,
        .nsimEIM = nsimEIM, .zero = zero ))),


  initialize = eval(substitute(expression({
    if (!all(w == 1))
      extra$orig.w <- w

    if (is.null( .nsimEIM )) {
      save.weights <-
      control$save.weights <- FALSE
    }

    mustart.orig <- mustart
    eval(binomialff()@initialize)  # Note
    if (length(mustart.orig))
      mustart <- mustart.orig  # Retainif inputted


    ycounts <- if (is.numeric(extra$orig.w))
         y * w / extra$orig.w else
         y * w  # Convert proportions to counts
    if (max(abs(ycounts-round(ycounts))) > 1.0e-6)
      warning("the response (as counts) does not ",
              "appear to be integer-valued. ",
              "Am rounding to integer values.")
    ycounts <- round(ycounts)  # Now an integer
    predictors.names <-
    c(namesof("mu",  .lmu ,  .emu  , tag = FALSE),
      namesof("rho", .lrho , .erho , tag = FALSE))

    if (!length(etastart)) {
      betabinomial.Loglikfun <-
      function(rhoval, y, x, w, extraargs) {
        shape1 <-    extraargs$mustart  * (1-rhoval) / rhoval
        shape2 <- (1-extraargs$mustart) * (1-rhoval) / rhoval
        ycounts <- extraargs$ycounts  # Ought to be integer-valued
        nvec <- extraargs$nvec
        sum(dbetabinom.ab(ycounts,
               size = nvec, shape1 = shape1,
               shape2 = shape2, log = TRUE))
      }
      rho.grid <- seq(0.05, 0.95, len = 25)  # rvar =
      mustart.use <- if (length(mustart.orig)) {
        mustart.orig
      } else if ( .imethod == 1) {
        rep_len(weighted.mean(y, w), n)
      } else if ( .imethod == 2) {
        .ishrinkage * weighted.mean(y, w) +
        (1 - .ishrinkage ) * y
      } else if ( .imethod == 3) {
        ymat <- cbind(y)
        mat.temp <- matrix(colMeans(ymat),
                           nrow(ymat), ncol(ymat),
                           byrow = TRUE)
        0.5 * mustart + 0.5 * mat.temp
      } else {
        mustart
      }
        try.this <-
          grid.search(rho.grid,
            objfun = betabinomial.Loglikfun,
            y = y,  x = x, w = w,
            extraargs = list(
            ycounts = ycounts,
            nvec = if (is.numeric(extra$orig.w))
            round(w / extra$orig.w) else round(w),
            mustart = mustart.use))
      init.rho <- if (is.Numeric( .irho ))
                  rep_len( .irho , n) else
                  rep_len(try.this, n)
      etastart <-
    cbind(theta2eta(mustart.use, .lmu  ,  .emu ),
          theta2eta(init.rho,    .lrho , .erho ))
      mustart <- NULL  # As etastart been computed
    }
  }),
  list( .lmu = lmu, .lrho = lrho,
        .emu = emu, .erho = erho,
        .imethod = imethod,
        .ishrinkage = ishrinkage,
        .nsimEIM = nsimEIM, .irho = irho ))),
  linkinv = eval(substitute(
      function(eta, extra = NULL)
    eta2theta(eta[, 1], .lmu , earg = .emu ),
  list( .lmu = lmu, .emu = emu ))),
  last = eval(substitute(expression({
    misc$link <-    c(mu = .lmu , rho = .lrho)
    misc$earg <- list(mu = .emu , rho = .erho )

    misc$zero <- .zero
    misc$expected <- TRUE
    misc$nsimEIM <- .nsimEIM
    misc$rho <- 1 / (shape1 + shape2 + 1)
  }),
  list( .lmu = lmu, .lrho = lrho,
        .emu = emu, .erho = erho,
        .nsimEIM = nsimEIM, .zero = zero ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL, summation = TRUE) {
      ycounts <- if (is.numeric(extra$orig.w))
               y * w / extra$orig.w else
               y * w  # Convert propns to counts

    mymu <- eta2theta(eta[, 1], .lmu ,  .emu )
    rho  <- eta2theta(eta[, 2], .lrho , .erho )
    smallno <- 1.0e4 * .Machine$double.eps

    if (max(abs(ycounts-round(ycounts))) > smallno)
      warning("converting 'ycounts' to integer",
              " in @loglikelihood")
    ycounts <- round(ycounts)

    rho  <- pmax(rho,     smallno)
    rho  <- pmin(rho, 1 - smallno)
    shape1 <-      mymu  * (1 - rho) / rho
    shape2 <- (1 - mymu) * (1 - rho) / rho

    nvec <- if (is.numeric(extra$orig.w))
            round(w / extra$orig.w) else
            round(w)

    if (residuals) {
      stop("loglikelihood resid not implemented")
    } else {
      ll.elts <-
        (if (is.numeric(extra$orig.w))
           extra$orig.w else 1) *
        dbetabinom.ab(ycounts, size = nvec,
                      shape1 = shape1,
                      shape2 = shape2, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lmu = lmu, .lrho = lrho,
           .emu = emu, .erho = erho  ))),
  vfamily = c("betabinomial"),
  validparams = eval(substitute(
      function(eta, y, extra = NULL) {
    mymu <- eta2theta(eta[, 1], .lmu  , .emu  )
    rho  <- eta2theta(eta[, 2], .lrho , .erho )
    okay1 <-
      all(is.finite(mymu)) &&
      all(0 < mymu & mymu < 1) &&
      all(is.finite(rho )) &&
      all(0 < rho  & rho  < 1)
    okay1
  },
  list( .lmu = lmu, .lrho = lrho,
        .emu = emu, .erho = erho  ))),




  simslot = eval(substitute(
  function(object, nsim) {

    pwts <-
      if (length(pwts <- object@prior.weights) > 0)
          pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    w <- pwts
    eta <- predict(object)
    extra <- object@extra

    mymu <- eta2theta(eta[, 1],  .lmu ,  .emu )
    rho  <- eta2theta(eta[, 2], .lrho , .erho )
    nvec <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else
              round(w)

    rbetabinom(nsim * length(rho), size = nvec,
               prob = mymu, rho  = rho)
  },
  list( .lmu = lmu, .lrho = lrho,
        .emu = emu, .erho = erho  ))),





  deriv = eval(substitute(expression({
    nvec <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else
              round(w)
    ycounts <- if (is.numeric(extra$orig.w))
            y * w / extra$orig.w else
            y * w  # Convert propns to counts

    ycounts <- round(ycounts)
    mymu <- eta2theta(eta[, 1], .lmu  , .emu  )
    rho  <- eta2theta(eta[, 2], .lrho , .erho )
    smallno <- 100 * .Machine$double.eps
    rho  <- pmax(rho, smallno)
    rho  <- pmin(rho, 1-smallno)

    shape1 <-      mymu  * (1 - rho) / rho
    shape2 <- (1 - mymu) * (1 - rho) / rho
    dshape1.dmu <-  (1 - rho) / rho
    dshape2.dmu <- -(1 - rho) / rho
    dshape1.drho <-       -mymu  / rho^2
    dshape2.drho <-  -(1 - mymu) / rho^2

    dmu.deta  <- dtheta.deta(mymu, .lmu  , .emu  )
    drho.deta <- dtheta.deta(rho,  .lrho , .erho )

    dl.dmu <- dshape1.dmu *
        (digamma(shape1 + ycounts) -
         digamma(shape2 + nvec - ycounts) -
         digamma(shape1) + digamma(shape2))
    dl.drho <- (-1 / rho^2) *
        (mymu * digamma(shape1 + ycounts) +
        (1 - mymu) *
        digamma(shape2 + nvec - ycounts) -
        digamma(shape1 + shape2 + nvec) -
        mymu * digamma(shape1) -
        (1 - mymu) * digamma(shape2) +
        digamma(shape1 + shape2))
    ansd <- c(if (is.numeric(extra$orig.w))
        extra$orig.w else 1) *
        cbind(dl.dmu  * dmu.deta,
              dl.drho * drho.deta)
    ansd
  }),
  list( .lmu = lmu, .lrho = lrho,
        .emu = emu, .erho = erho  ))),
  weight = eval(substitute(expression({
    if (is.null( .nsimEIM )) {
      wz <- matrix(NA_real_, n, dimm(M))  # 3
      wz11 <-
        -(Ebbin.ab(nvec, shape1, shape2, TRUE) -
         trigamma(shape1 + shape2 + nvec) -
         trigamma(shape1) +
         trigamma(shape1 + shape2))
      wz22 <-
        -(Ebbin.ab(nvec, shape1, shape2, FALSE) -
        trigamma(shape1 + shape2 + nvec) -
        trigamma(shape2) +
        trigamma(shape1 + shape2))
      wz21 <- -(trigamma(shape1 + shape2) -
                trigamma(shape1 + shape2 + nvec))

      wz[, iam(1, 1, M)] <- dmu.deta^2 *
          (wz11 * dshape1.dmu^2 +
           wz22 * dshape2.dmu^2 +
           2 * wz21 * dshape1.dmu * dshape2.dmu)
      wz[, iam(2, 2, M)] <- drho.deta^2 *
          (wz11 * dshape1.drho^2 +
           wz22 * dshape2.drho^2 +
           2 * wz21 * dshape1.drho * dshape2.drho)
      wz[, iam(2, 1, M)] <- dmu.deta *
          drho.deta *
          (dshape1.dmu*(wz11 * dshape1.drho +
                        wz21 * dshape2.drho) +
           dshape2.dmu*(wz21 * dshape1.drho +
                        wz22 * dshape2.drho))

      wz * (if (is.numeric(extra$orig.w))
                extra$orig.w else 1)
    } else {
      run.varcov <- 0
      ind1 <- iam(NA, NA, M = M, both = TRUE,
                  diag = TRUE)
      dthetas.detas <- cbind(dmu.deta, drho.deta)

      for (ii in 1:( .nsimEIM )) {
        ysim <- rbetabinom.ab(n, size = nvec,
                              shape1 = shape1,
                              shape2 = shape2)
        dl.dmu <- dshape1.dmu *
            (digamma(shape1 + ysim) -
             digamma(shape2 + nvec - ysim) -
             digamma(shape1) + digamma(shape2))
        dl.drho <- (-1/rho^2) *
            (mymu * digamma(shape1 + ysim) +
            (1 - mymu) *
            digamma(shape2 + nvec - ysim) -
            digamma(shape1 + shape2 + nvec) -
            mymu * digamma(shape1) -
            (1 - mymu) * digamma(shape2) +
            digamma(shape1 + shape2))


        temp3 <- cbind(dl.dmu, dl.drho)  # n x M
        run.varcov <- run.varcov +
                      temp3[, ind1$row.index] *
                      temp3[, ind1$col.index]
      }
      run.varcov <- run.varcov / .nsimEIM


      wz <- if (intercept.only)
          matrix(colMeans(run.varcov),
                 n, ncol(run.varcov),
                 byrow = TRUE) else run.varcov

      wz <- wz * dthetas.detas[, ind1$row] *
                 dthetas.detas[, ind1$col]
      wz * (if (is.numeric(extra$orig.w))
              extra$orig.w else 1)
    }
  }),
  list( .lmu = lmu, .lrho = lrho,
        .emu = emu, .erho = erho,
        .nsimEIM = nsimEIM ))))
}  # betabinomial






 dbinom2.or <-
  function(mu1,
           mu2 = if (exchangeable) mu1 else
                 stop("'mu2' not specified"),
           oratio = 1,
           exchangeable = FALSE,
           tol = 0.001,
           colnames = c("00", "01", "10", "11"),
           ErrorCheck = TRUE) {
  if (ErrorCheck) {
    if (!is.Numeric(mu1, positive = TRUE) || max(mu1) >= 1)
      stop("bad input for argument 'mu1'")
    if (!is.Numeric(mu2, positive = TRUE) || max(mu2) >= 1)
      stop("bad input for argument 'mu2'")
    if (!is.Numeric(oratio, positive = TRUE))
      stop("bad input for argument 'oratio'")
    if (!is.Numeric(tol, positive = TRUE, length.arg = 1) ||
        tol > 0.1)
      stop("bad input for argument 'tol'")
    if (exchangeable && max(abs(mu1 - mu2)) > 0.00001)
      stop("argument 'exchangeable' is TRUE but 'mu1' and 'mu2' differ")
  }

  L <- max(length(mu1), length(mu2), length(oratio))
  if (length(oratio) != L) oratio <- rep_len(oratio, L)
  if (length(mu1   ) != L) mu1    <- rep_len(mu1,    L)
  if (length(mu2   ) != L) mu2    <- rep_len(mu2,    L)

  a.temp <- 1 + (mu1+mu2)*(oratio-1)
  b.temp <- -4 * oratio * (oratio-1) * mu1 * mu2
  temp <- sqrt(a.temp^2 + b.temp)
  p11 <- ifelse(abs(oratio-1) < tol,
                mu1*mu2,
               (a.temp-temp)/(2*(oratio-1)))
  p01 <- mu2 - p11
  p10 <- mu1 - p11
  p00 <- 1 - p11 - p01 - p10
  matrix(c(p00, p01, p10, p11), L, 4, dimnames = list(NULL, colnames))
}  # dbinom2.or




 rbinom2.or <-
  function(n, mu1,
           mu2 = if (exchangeable) mu1 else
                   stop("argument 'mu2' not specified"),
           oratio = 1,
           exchangeable = FALSE,
           tol = 0.001,
           twoCols = TRUE,
           colnames = if (twoCols) c("y1", "y2") else
                         c("00", "01", "10", "11"),
           ErrorCheck = TRUE) {
  use.n <- if ((length.n <- length(n)) > 1) length.n else
           if (!is.Numeric(n, integer.valued = TRUE,
                           length.arg = 1, positive = TRUE))
              stop("bad input for argument 'n'") else n

  if (ErrorCheck) {
    if (!is.Numeric(mu1, positive = TRUE) || max(mu1) >= 1)
      stop("bad input for argument 'mu1'")
    if (!is.Numeric(mu2, positive = TRUE) || max(mu2) >= 1)
      stop("bad input for argument 'mu2'")
    if (!is.Numeric(oratio, positive = TRUE))
      stop("bad input for argument 'oratio'")
    if (!is.Numeric(tol, positive = TRUE, length.arg = 1) ||
        tol > 0.1)
      stop("bad input for argument 'tol'")
    if (exchangeable && max(abs(mu1 - mu2)) > 0.00001)
      stop("argument 'exchangeable' is TRUE but 'mu1' and 'mu2' differ")
  }

  dmat <- dbinom2.or(mu1 = mu1, mu2 = mu2, oratio = oratio,
                     exchangeable = exchangeable,
                     tol = tol, ErrorCheck = ErrorCheck)

  answer <- matrix(0, use.n, 2,
                   dimnames = list(NULL,
                                   if (twoCols) colnames else NULL))
  yy <- runif(use.n)
  cs1 <- dmat[, "00"] + dmat[, "01"]
  cs2 <- cs1 + dmat[, "10"]
  index <- (dmat[, "00"] < yy) & (yy <= cs1)
  answer[index, 2] <- 1
  index <- (cs1 < yy) & (yy <= cs2)
  answer[index, 1] <- 1
  index <- (yy > cs2)
  answer[index,] <- 1
  if (twoCols) {
    answer
  } else {
    answer4 <- matrix(0, use.n, 4, dimnames = list(NULL, colnames))
    answer4[cbind(1:use.n, 1 + 2*answer[, 1] + answer[, 2])] <- 1
    answer4
  }
}  # rbinom2.or







 binom2.or <-
  function(lmu = "logitlink", lmu1 = lmu, lmu2 = lmu,
           loratio = "loglink",
           imu1 = NULL, imu2 = NULL, ioratio = NULL,
           zero = "oratio",  # zero = 3,
           exchangeable = FALSE,
           tol = 0.001,
           more.robust = FALSE) {

  lmu1 <- lmu1
  lmu2 <- lmu2


  lmu1 <- as.list(substitute(lmu1))
  emu1 <- link2list(lmu1)
  lmu1 <- attr(emu1, "function.name")

  lmu2 <- as.list(substitute(lmu2))
  emu2 <- link2list(lmu2)
  lmu2 <- attr(emu2, "function.name")


  loratio <- as.list(substitute(loratio))
  eoratio <- link2list(loratio)
  loratio <- attr(eoratio, "function.name")



  if (!is.logical(exchangeable))
    warning("argument 'exchangeable' should be a single logical")

  if (is.logical(exchangeable) && exchangeable &&
     ((lmu1 != lmu2) || !identical(emu1, emu2)))
    warning("exchangeable = TRUE but marginal links are not equal")



  if (!is.Numeric(tol, positive = TRUE, length.arg = 1) ||
      tol > 0.1)
    stop("bad input for argument 'tol'")


  new("vglmff",
  blurb = c("Bivariate binomial regression with an odds ratio\n",
            "Links:    ",
            namesof("mu1", lmu1, earg = emu1), ", ",
            namesof("mu2", lmu2, earg = emu2), "; ",
            namesof("oratio", loratio, earg = eoratio)),
  constraints = eval(substitute(expression({
    cm.intercept.default <- diag(3)
    constraints <- cm.VGAM(matrix(c(1, 1, 0, 0, 0, 1), 3, 2),
                     x = x,
                     bool = .exchangeable ,
                     constraints = constraints,
                     apply.int = TRUE,
                     cm.default           = cm.intercept.default,
                     cm.intercept.default = cm.intercept.default)
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 3,
                     predictors.names = predictors.names)
  }), list( .exchangeable = exchangeable, .zero = zero ))),
  deviance = Deviance.categorical.data.vgam,

  infos = eval(substitute(function(...) {
    list(M1 = 3,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("mu1", "mu2", "oratio"),
         exchangeable = .exchangeable ,
         lmu1 = .lmu1 ,
         lmu2 = .lmu2 ,
         loratio = .loratio ,
         zero = .zero )
  }, list( .lmu1 = lmu1,
           .lmu2 = lmu2,
           .loratio = loratio,
           .zero = zero,
           .exchangeable = exchangeable
         ))),


  initialize = eval(substitute(expression({
    mustart.orig <- mustart
    eval(process.binomial2.data.VGAM)
    if (length(mustart.orig))
      mustart <- mustart.orig  # Retain it if inputted


    predictors.names <-
       c(namesof("mu1",    .lmu1 ,    earg = .emu1 ,    short = TRUE),
         namesof("mu2",    .lmu2 ,    earg = .emu2 ,    short = TRUE),
         namesof("oratio", .loratio , earg = .eoratio , short = TRUE))


    if (!length(etastart)) {
        pmargin <- cbind(mustart[, 3] + mustart[, 4],
                         mustart[, 2] + mustart[, 4])
        ioratio <- if (length( .ioratio ))
                     rep_len( .ioratio , n) else
                     mustart[, 4] * mustart[, 1] / (mustart[, 2] *
                                                    mustart[, 3])
        if (length( .imu1 )) pmargin[, 1] <- .imu1
        if (length( .imu2 )) pmargin[, 2] <- .imu2
        etastart <-
          cbind(theta2eta(pmargin[, 1], .lmu1 , earg = .emu1 ),
                theta2eta(pmargin[, 2], .lmu2 , earg = .emu2 ),
                theta2eta(ioratio, .loratio , earg = .eoratio ))
    }
  }), list( .lmu1 = lmu1, .lmu2 = lmu2, .loratio = loratio,
            .emu1 = emu1, .emu2 = emu2, .eoratio = eoratio,
            .imu1 = imu1, .imu2 = imu2, .ioratio = ioratio ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    pmargin <- cbind(eta2theta(eta[, 1], .lmu1 , earg = .emu1 ),
                     eta2theta(eta[, 2], .lmu2 , earg = .emu2 ))
    oratio <- eta2theta(eta[, 3], .loratio , earg = .eoratio )
    a.temp <- 1 + (pmargin[, 1] + pmargin[, 2]) * (oratio - 1)
    b.temp <- -4 * oratio * (oratio-1) * pmargin[, 1] * pmargin[, 2]
    temp <- sqrt(a.temp^2 + b.temp)
    pj4 <- ifelse(abs(oratio-1) < .tol , pmargin[, 1] * pmargin[, 2],
                 (a.temp-temp)/(2*(oratio-1)))
    pj2 <- pmargin[, 2] - pj4
    pj3 <- pmargin[, 1] - pj4
    cbind("00" = 1-pj4-pj2-pj3,
          "01" = pj2,
          "10" = pj3,
          "11" = pj4)
  }, list( .lmu1 = lmu1, .lmu2 = lmu2, .loratio = loratio,
           .emu1 = emu1, .emu2 = emu2, .eoratio = eoratio,
           .tol = tol ))),
  last = eval(substitute(expression({
    misc$link <-    c(mu1 = .lmu1 , mu2 = .lmu2 , oratio = .loratio )

    misc$earg <- list(mu1 = .emu1 , mu2 = .emu2 , oratio = .eoratio )

    misc$tol <- .tol
    misc$expected <- TRUE
  }), list( .lmu1 = lmu1, .lmu2 = lmu2, .loratio = loratio,
            .emu1 = emu1, .emu2 = emu2, .eoratio = eoratio,
            .tol = tol ))),
  linkfun = eval(substitute(function(mu, extra = NULL) {
    pmargin <- cbind(mu[, 3] + mu[, 4], mu[, 2] + mu[, 4])
    oratio <- mu[, 4] * mu[, 1] / (mu[, 2] * mu[, 3])
    cbind(theta2eta(pmargin[, 1], .lmu1 , earg = .emu1),
          theta2eta(pmargin[, 2], .lmu2 , earg = .emu2),
          theta2eta(oratio,      .loratio, earg = .eoratio))
  },
  list( .lmu1 = lmu1, .lmu2 = lmu2, .loratio = loratio,
        .emu1 = emu1, .emu2 = emu2, .eoratio = eoratio ))),
  loglikelihood = eval(substitute(
      function(mu, y, w, residuals = FALSE, eta,
               extra = NULL, summation = TRUE) {
    if (residuals) {
      stop("loglikelihood resids not implemented")
    } else {
      if ( .more.robust) {
        vsmallno <-  1.0e4 * .Machine$double.xmin
        mu[mu < vsmallno] <- vsmallno
      }

      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
      }
    }
  }, list( .more.robust = more.robust ))),
  vfamily = c("binom2.or", "binom2"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    pmargin <- cbind(eta2theta(eta[, 1], .lmu1 , earg = .emu1 ),
                     eta2theta(eta[, 2], .lmu2 , earg = .emu2 ))
    oratio <- eta2theta(eta[, 3], .loratio , earg = .eoratio )
    okay1 <- all(is.finite(pmargin)) &&
             all(0 < pmargin & pmargin < 1) &&
             all(is.finite(oratio )) && all(0 < oratio)
    okay1
  }, list( .lmu1 = lmu1, .lmu2 = lmu2, .loratio = loratio,
           .emu1 = emu1, .emu2 = emu2, .eoratio = eoratio ))),








  deriv = eval(substitute(expression({
    smallno <- 1.0e4 * .Machine$double.eps
    mu.use <- mu
    mu.use[mu.use <     smallno] <-     smallno
    mu.use[mu.use > 1 - smallno] <- 1 - smallno
    pmargin <- cbind(mu.use[, 3] + mu.use[, 4],
                     mu.use[, 2] + mu.use[, 4])
    pmargin[, 1] <- pmax(    smallno, pmargin[, 1])
    pmargin[, 1] <- pmin(1 - smallno, pmargin[, 1])
    pmargin[, 2] <- pmax(    smallno, pmargin[, 2])
    pmargin[, 2] <- pmin(1 - smallno, pmargin[, 2])

    oratio <- mu.use[, 4]*mu.use[, 1] / (mu.use[, 2]*mu.use[, 3])
    use.oratio <- pmax(smallno, oratio)
    a.temp <- 1 + (pmargin[, 1]+pmargin[, 2])*(oratio-1)
    b.temp <- -4 * oratio * (oratio-1) * pmargin[, 1] * pmargin[, 2]
    temp9 <- sqrt(a.temp^2 + b.temp)

    coeff12 <- -0.5 + (2*oratio*pmargin - a.temp) / (2*temp9)
    dl.dmu1 <- coeff12[, 2] *
      (y[, 1] / mu.use[, 1] - y[, 3] / mu.use[, 3]) -
      (1+coeff12[, 2]) * (y[, 2] / mu.use[, 2] - y[, 4] / mu.use[, 4])

    dl.dmu2 <- coeff12[, 1] *
       (y[, 1] / mu.use[, 1] - y[, 2] / mu.use[, 2]) -
       (1+coeff12[, 1]) * (y[, 3] / mu.use[, 3] - y[, 4] / mu.use[, 4])

    coeff3 <- (y[, 1]/mu.use[, 1] - y[, 2]/mu.use[, 2] -
              y[, 3]/mu.use[, 3] + y[, 4]/mu.use[, 4])
    Vab <- pmax(smallno, 1 / (1 / mu.use[, 1] + 1 / mu.use[, 2] +
                              1 / mu.use[, 3] + 1 / mu.use[, 4]))
    dp11.doratio <- Vab / use.oratio
    dl.doratio <- coeff3 * dp11.doratio

    c(w) *
    cbind(dl.dmu1 * dtheta.deta(pmargin[, 1], .lmu1, earg = .emu1),
          dl.dmu2 * dtheta.deta(pmargin[, 2], .lmu2, earg = .emu2),
          dl.doratio * dtheta.deta(oratio, .loratio, earg = .eoratio))
  }), list( .lmu1 = lmu1, .lmu2 = lmu2, .loratio = loratio,
            .emu1 = emu1, .emu2 = emu2, .eoratio = eoratio ))),
  weight = eval(substitute(expression({
    Deltapi <- mu.use[, 3] * mu.use[, 2] - mu.use[, 4] * mu.use[, 1]
    myDelta  <- pmax(smallno, mu.use[, 1] * mu.use[, 2] *
                              mu.use[, 3] * mu.use[, 4])
    pqmargin <- pmargin * (1-pmargin)
    pqmargin[pqmargin < smallno] <- smallno

    wz <- matrix(0, n, 4)
    wz[, iam(1, 1, M)] <- (pqmargin[, 2] * Vab / myDelta) *
                      dtheta.deta(pmargin[, 1], .lmu1, earg = .emu1)^2
    wz[, iam(2, 2, M)] <- (pqmargin[, 1] * Vab / myDelta) *
                      dtheta.deta(pmargin[, 2], .lmu2, earg = .emu2)^2
    wz[, iam(3, 3, M)] <- (Vab / use.oratio^2) *
                 dtheta.deta(use.oratio, .loratio, earg = .eoratio)^2
    wz[, iam(1, 2, M)] <- (Vab * Deltapi / myDelta) *
                      dtheta.deta(pmargin[, 1], .lmu1, earg = .emu1) *
                      dtheta.deta(pmargin[, 2], .lmu2, earg = .emu2)
    c(w) * wz
  }), list( .lmu1 = lmu1, .lmu2 = lmu2, .loratio = loratio,
            .emu1 = emu1, .emu2 = emu2, .eoratio = eoratio ))))
}  # binom2.or






setClass("binom2",         contains = "vglmff")
setClass("binom2.or",      contains = "binom2")



setMethod("summaryvglmS4VGAM",  signature(VGAMff = "binom2.or"),
  function(object,
           VGAMff,
           ...) {


  cfit <- coefvlm(object, matrix.out = TRUE)

  if (rownames(cfit)[1] == "(Intercept)" &&
      all(cfit[-1, 3] == 0)) {
    object@post$oratio <- eta2theta(cfit[1, 3],
                                    link = object@misc$link[3],
                                    earg = object@misc$earg[[3]])
  }

  object@post
})



setMethod("showsummaryvglmS4VGAM",  signature(VGAMff = "binom2.or"),
  function(object,
           VGAMff,
           ...) {
 if (length(object@post$oratio) == 1 &&
      is.numeric(object@post$oratio)) {
    cat("\nOdds ratio: ", round(object@post$oratio, digits = 4), "\n")
  }
})














 dbinom2.rho <-
  function(mu1,
           mu2 = if (exchangeable) mu1 else stop("'mu2' not specified"),
           rho = 0,
           exchangeable = FALSE,
           colnames = c("00", "01", "10", "11"),
           ErrorCheck = TRUE) {
  if (ErrorCheck) {
    if (!is.Numeric(mu1, positive = TRUE) || max(mu1) >= 1)
      stop("bad input for argument 'mu1'")
    if (!is.Numeric(mu2, positive = TRUE) || max(mu2) >= 1)
      stop("bad input for argument 'mu2'")
    if (!is.Numeric(rho) || min(rho) <= -1 || max(rho) >= 1)
      stop("bad input for argument 'rho'")
    if (exchangeable && max(abs(mu1 - mu2)) > 0.00001)
      stop("argument 'exchangeable' is TRUE but 'mu1' and 'mu2' differ")
  }

  nn <- max(length(mu1), length(mu2), length(rho))
  rho <- rep_len(rho, nn)
  mu1 <- rep_len(mu1, nn)
  mu2 <- rep_len(mu2, nn)
  eta1 <- qnorm(mu1)
  eta2 <- qnorm(mu2)
  p11 <- pbinorm(eta1, eta2, cov12 = rho)
  p01 <- mu2 - p11
  p10 <- mu1 - p11
  p00 <- 1.0 - p01 - p10 - p11
  matrix(c(p00, p01, p10, p11), nn, 4,
         dimnames = list(NULL, colnames))
}  # dbinom2.rho



 rbinom2.rho <-
  function(n, mu1,
           mu2 = if (exchangeable) mu1 else
                   stop("argument 'mu2' not specified"),
           rho = 0,
           exchangeable = FALSE,
           twoCols = TRUE,
           colnames = if (twoCols) c("y1", "y2") else
                      c("00", "01", "10", "11"),
           ErrorCheck = TRUE) {
  use.n <- if ((length.n <- length(n)) > 1) length.n else
           if (!is.Numeric(n, integer.valued = TRUE,
                           length.arg = 1, positive = TRUE))
              stop("bad input for argument 'n'") else n

  if (ErrorCheck) {
    if (!is.Numeric(mu1, positive = TRUE) ||
        max(mu1) >= 1)
      stop("bad input for argument 'mu1'")
    if (!is.Numeric(mu2, positive = TRUE) ||
        max(mu2) >= 1)
      stop("bad input for argument 'mu2'")
    if (!is.Numeric(rho) || min(rho) <= -1 ||
        max(rho) >= 1)
      stop("bad input for argument 'rho'")


    if (exchangeable &&
        max(abs(mu1 - mu2)) > 0.00001)
      stop("argument 'exchangeable' is TRUE but 'mu1' and 'mu2' differ")
  }

  dmat <- dbinom2.rho(mu1 = mu1, mu2 = mu2, rho = rho,
                      exchangeable = exchangeable,
                      ErrorCheck = ErrorCheck)

  answer <- matrix(0, use.n, 2,
                   dimnames = list(NULL,
                                   if (twoCols) colnames else NULL))
  yy <- runif(use.n)
  cs1 <- dmat[, "00"] + dmat[, "01"]
  cs2 <- cs1 + dmat[, "10"]
  index <- (dmat[, "00"] < yy) & (yy <= cs1)
  answer[index, 2] <- 1
  index <- (cs1 < yy) & (yy <= cs2)
  answer[index, 1] <- 1
  index <- (yy > cs2)
  answer[index,] <- 1
  if (twoCols) {
    answer
  } else {
    answer4 <- matrix(0, use.n, 4, dimnames = list(NULL, colnames))
    answer4[cbind(1:use.n, 1 + 2*answer[, 1] + answer[, 2])] <- 1
    answer4
  }
}  # rbinom2.rho





binom2.rho.control <-
    function(save.weights = TRUE, ...) {
  list(save.weights = save.weights)
}





 binom2.rho <-
  function(lmu = "probitlink",  # added 20120817
           lrho = "rhobitlink",
           imu1 = NULL, imu2 = NULL, irho = NULL,
           imethod = 1,
           zero =  "rho",   # 3
           exchangeable = FALSE,
           grho = seq(-0.95, 0.95, by = 0.05),
           nsimEIM = NULL) {



  lrho <- as.list(substitute(lrho))
  erho <- link2list(lrho)
  lrho <- attr(erho, "function.name")

  lmu <- as.list(substitute(lmu))
  emu <- link2list(lmu)
  lmu <- attr(emu, "function.name")

  if (lmu != "probitlink")
    warning("arg 'lmu' should be 'probitlink'")

    lmu12 <- "probitlink"
    emu12 <- emu # list()



  if (is.Numeric(nsimEIM)) {
    if (!is.Numeric(nsimEIM, length.arg = 1,
                    integer.valued = TRUE))
      stop("bad input for argument 'nsimEIM'")
    if (nsimEIM <= 100)
      warning("'nsimEIM' should be an integer greater than 100")
  }

  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE, positive = TRUE) ||
      imethod > 2)
    stop("argument 'imethod' must be 1 or 2")




  new("vglmff",
  blurb = c("Bivariate probit model\n",
            "Links:    ",
            namesof("mu1", lmu12, earg = emu12), ", ",
            namesof("mu2", lmu12, earg = emu12), ", ",
            namesof("rho", lrho,  earg = erho)),
  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(c(1, 1, 0, 0, 0, 1), 3, 2),
                           x = x,
                           bool = .exchangeable ,
                           constraints = constraints,
                           apply.int = TRUE)
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 3,
                     predictors.names = predictors.names)
  }), list( .exchangeable = exchangeable, .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 3,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("mu1", "mu2", "rho"),
         lmu1 = .lmu12,
         lmu2 = .lmu12,
         lrho = .lrho ,
         zero = .zero )
  }, list( .lmu12 = lmu12, .lrho = lrho,
           .zero = zero ))),


  initialize = eval(substitute(expression({
    mustart.orig <- mustart
    eval(process.binomial2.data.VGAM)

    if (length(mustart.orig))
      mustart <- mustart.orig  # Retain it if inputted

    predictors.names <- c(
        namesof("mu1", .lmu12 , earg = .emu12 , short = TRUE),
        namesof("mu2", .lmu12 , earg = .emu12 , short = TRUE),
        namesof("rho", .lrho ,  earg = .erho,  short = TRUE))

    if (is.null( .nsimEIM )) {
      save.weights <- control$save.weights <- FALSE
    }


    ycounts <- if (is.numeric(extra$orig.w))
               y * c(w) / extra$orig.w else
               y * c(w)  # Convert proportions to counts
    if (max(abs(ycounts - round(ycounts))) > 1.0e-6)
       warning("the response (as counts) does not appear to ",
               "be integer-valued. Am rounding to integer values.")
    ycounts <- round(ycounts)  # Make sure it is an integer
    nvec <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else
              round(w)


    if (is.null(etastart)) {
      if (length(mustart.orig)) {
        mu1.init <- mustart.orig[, 3] + mustart.orig[, 4]
        mu2.init <- mustart.orig[, 2] + mustart.orig[, 4]
      } else if ( .imethod == 1) {
          glm1.fit <- glm(cbind(ycounts[, 3] + ycounts[, 4],
                                ycounts[, 1] + ycounts[, 2]) ~ x - 1,
                          fam = binomial("probit"))
          glm2.fit <- glm(cbind(ycounts[, 2] + ycounts[, 4],
                                ycounts[, 1] + ycounts[, 3]) ~ x - 1,
                          fam = binomial("probit"))
          mu1.init <- fitted(glm1.fit)
          mu2.init <- fitted(glm2.fit)
      } else if ( .imethod == 2) {
          mu1.init <- if (is.Numeric( .imu1 ))
                      rep_len( .imu1 , n) else
                      mu[, 3] + mu[, 4]
          mu2.init <- if (is.Numeric( .imu2 ))
                      rep_len( .imu2 , n) else
                      mu[, 2] + mu[, 4]
      } else {
        stop("bad value for argument 'imethod'")
      }



      binom2.rho.Loglikfun <-
          function(rhoval, y, x, w, extraargs) {
          init.mu1 <-    extraargs$initmu1
          init.mu2 <-    extraargs$initmu2
          ycounts  <-    extraargs$ycounts
          nvec     <-    extraargs$nvec
          eta1 <- qnorm(init.mu1)
          eta2 <- qnorm(init.mu2)
          p11 <- pbinorm(eta1, eta2, cov12 = rhoval)
          p01 <- pmin(init.mu2 - p11, init.mu2)
          p10 <- pmin(init.mu1 - p11, init.mu1)
          p00 <- 1.0 - p01 - p10 - p11
          mumat <- abs(cbind("00" = p00,
                             "01" = p01,
                             "10" = p10,
                             "11" = p11))
          mumat <- mumat / rowSums(mumat)
          mumat[mumat < 1.0e-100] <- 1.0e-100

          sum((if (is.numeric(extraargs$orig.w))
               extraargs$orig.w else 1) *
               dmultinomial(x = ycounts, size = nvec, prob = mumat,
                            log = TRUE, dochecking = FALSE))
        }
        rho.grid <- .grho # seq(-0.95, 0.95, len = 31)
        try.this <-
          grid.search(rho.grid, objfun = binom2.rho.Loglikfun,
                      y = y, x = x, w = w, extraargs = list(
                      orig.w = extra$orig.w,
                      ycounts = ycounts,
                      initmu1 = mu1.init,
                      initmu2 = mu2.init,
                      nvec = nvec ))


      rho.init <- if (is.Numeric( .irho ))
                   rep_len( .irho , n) else {
          try.this
      }

      etastart <- cbind(theta2eta(mu1.init, .lmu12 , earg = .emu12 ),
                        theta2eta(mu2.init, .lmu12 , earg = .emu12 ),
                        theta2eta(rho.init, .lrho ,  earg = .erho ))
      mustart <- NULL # Since etastart has been computed.
    }
  }), list( .lmu12 = lmu12, .lrho = lrho,
            .emu12 = emu12, .erho = erho,
                            .grho = grho,
                            .irho = irho,
            .imethod = imethod, .nsimEIM = nsimEIM,
            .imu1 = imu1, .imu2 = imu2 ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    pmargin <- cbind(eta2theta(eta[, 1], .lmu12 , earg = .emu12 ),
                     eta2theta(eta[, 2], .lmu12 , earg = .emu12 ))
    rho <- eta2theta(eta[, 3], .lrho , earg = .erho )
    p11 <- pbinorm(eta[, 1], eta[, 2], cov12 = rho)
    p01 <- pmin(pmargin[, 2] - p11, pmargin[, 2])
    p10 <- pmin(pmargin[, 1] - p11, pmargin[, 1])
    p00 <- 1.0 - p01 - p10 - p11
    ansmat <- abs(cbind("00" = p00,
                        "01" = p01,
                        "10" = p10,
                        "11" = p11))
    ansmat / rowSums(ansmat)
  }, list( .lmu12 = lmu12, .lrho = lrho,
           .emu12 = emu12, .erho = erho ))),
  last = eval(substitute(expression({
    misc$link <-    c(mu1 = .lmu12 , mu2 = .lmu12 , rho = .lrho )

    misc$earg <- list(mu1 = .emu12 , mu2 = .emu12 , rho = .erho )

    misc$nsimEIM <- .nsimEIM
    misc$expected <- TRUE
    misc$multipleResponses <- FALSE
  }),
  list( .lmu12 = lmu12, .lrho = lrho,
        .nsimEIM = nsimEIM,
        .emu12 = emu12, .erho = erho ))),

  loglikelihood = eval(substitute(
      function(mu, y, w, residuals = FALSE, eta,
               extra = NULL, summation = TRUE) {
    if (residuals) {
      stop("loglikelihood resids not implemented")
    } else {

      ycounts <- if (is.numeric(extra$orig.w))
                 y * c(w) / extra$orig.w else
                 y * c(w)  # propns 2 counts

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

      nvec <- if (is.numeric(extra$orig.w))
                round(w / extra$orig.w) else
                round(w)

      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
      }
    }
  }, list( .erho = erho ))),
  vfamily = c("binom2.rho", "binom2"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    pmargin <- cbind(eta2theta(eta[, 1], .lmu12 , earg = .emu12 ),
                     eta2theta(eta[, 2], .lmu12 , earg = .emu12 ))
    rhovec <- eta2theta(eta[, 3], .lrho , earg = .erho )
      okay1 <- all(is.finite(pmargin)) &&
               all( 0 < pmargin & pmargin < 1) &&
               all(is.finite(rhovec )) &&
               all(-1 < rhovec  & rhovec  < 1)
    okay1
  }, list( .lmu12 = lmu12, .lrho = lrho,
           .emu12 = emu12, .erho = erho ))),






  deriv = eval(substitute(expression({
    nvec <- if (is.numeric(extra$orig.w))
            round(w / extra$orig.w) else
            round(w)
    ycounts <- if (is.numeric(extra$orig.w))
               y * c(w) / extra$orig.w else
               y * c(w)  # Convert proportions to counts

    pmargin <- cbind(eta2theta(eta[, 1], .lmu12 , earg = .emu12 ),
                     eta2theta(eta[, 2], .lmu12 , earg = .emu12 ))
    rhovec <- eta2theta(eta[, 3], .lrho , earg = .erho )
    p11 <- pbinorm(eta[, 1], eta[, 2], cov12 = rhovec)
    p01 <- pmargin[, 2] - p11
    p10 <- pmargin[, 1] - p11
    p00 <- 1 - p01 - p10 - p11

    ABmat <- (eta[, 1:2] - rhovec * eta[, 2:1]) / (
      sqrt(pmax(1e5 * .Machine$double.eps, 1.0 - rhovec^2)))
    PhiA <- pnorm(ABmat[, 1])
    PhiB <- pnorm(ABmat[, 2])
    onemPhiA <- pnorm(ABmat[, 1], lower.tail = FALSE)
    onemPhiB <- pnorm(ABmat[, 2], lower.tail = FALSE)

    smallno <- 1000 * .Machine$double.eps
    p00[p00 < smallno] <- smallno
    p01[p01 < smallno] <- smallno
    p10[p10 < smallno] <- smallno
    p11[p11 < smallno] <- smallno

    dprob00 <- dbinorm(eta[, 1], eta[, 2], cov12 = rhovec)
    dl.dprob1 <-     PhiB * (ycounts[, 4]/p11 - ycounts[, 2]/p01) +
                 onemPhiB * (ycounts[, 3]/p10 - ycounts[, 1]/p00)
    dl.dprob2 <-     PhiA * (ycounts[, 4]/p11 - ycounts[, 3]/p10) +
                 onemPhiA * (ycounts[, 2]/p01 - ycounts[, 1]/p00)
    dl.drho <- (ycounts[, 4]/p11 - ycounts[, 3]/p10 -
                ycounts[, 2]/p01 + ycounts[, 1]/p00) * dprob00

    dprob1.deta <- dtheta.deta(pmargin[, 1], .lmu12 , earg = .emu12 )
    dprob2.deta <- dtheta.deta(pmargin[, 2], .lmu12 , earg = .emu12 )
    drho.deta <- dtheta.deta(rhovec, .lrho , earg = .erho )
    dthetas.detas <- cbind(dprob1.deta, dprob2.deta, drho.deta)

    (if (is.numeric(extra$orig.w)) extra$orig.w else 1) *
    cbind(dl.dprob1,
          dl.dprob2,
          dl.drho) * dthetas.detas
  }), list( .lmu12 = lmu12, .lrho = lrho,
            .emu12 = emu12, .erho = erho ))),
  weight = eval(substitute(expression({
    if (is.null( .nsimEIM )) {
      ned2l.dprob1prob1 <-      PhiB^2 * (1/p11 + 1/p01) +
                            onemPhiB^2 * (1/p10 + 1/p00)
      ned2l.dprob2prob2 <-      PhiA^2 * (1/p11 + 1/p10) +
                            onemPhiA^2 * (1/p01 + 1/p00)
      ned2l.dprob1prob2 <-      PhiA * (    PhiB/p11 - onemPhiB/p10) +
                            onemPhiA * (onemPhiB/p00 -     PhiB/p01)
      ned2l.dprob1rho <-     (PhiB * (1/p11 + 1/p01) -
                          onemPhiB * (1/p10 + 1/p00)) * dprob00
      ned2l.dprob2rho <-     (PhiA * (1/p11 + 1/p10) -
                          onemPhiA * (1/p01 + 1/p00)) * dprob00
      ned2l.drho2 <-  (1/p11 + 1/p01 + 1/p10 + 1/p00) * dprob00^2

      wz <- matrix(0, n, dimm(M))  # 6=dimm(M)
      wz[, iam(1, 1, M)] <- ned2l.dprob1prob1 * dprob1.deta^2
      wz[, iam(2, 2, M)] <- ned2l.dprob2prob2 * dprob2.deta^2
      wz[, iam(3, 3, M)] <- ned2l.drho2 * drho.deta^2
      wz[, iam(1, 2, M)] <- ned2l.dprob1prob2 * dprob1.deta * dprob2.deta
      wz[, iam(2, 3, M)] <- ned2l.dprob2rho * dprob2.deta * drho.deta
      wz[, iam(1, 3, M)] <- ned2l.dprob1rho * dprob1.deta * drho.deta
    } else {
      run.varcov <- 0
      ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)
      for (ii in 1:( .nsimEIM )) {
        ysim <- rbinom2.rho(n, mu1 = pmargin[, 1], mu2 = pmargin[, 2],
                            twoCols = FALSE, rho = rhovec)
        dl.dprob1 <-     PhiB * (ysim[, 4]/p11 - ysim[, 2]/p01) +
                     onemPhiB * (ysim[, 3]/p10 - ysim[, 1]/p00)
        dl.dprob2 <-     PhiA * (ysim[, 4]/p11 - ysim[, 3]/p10) +
                     onemPhiA * (ysim[, 2]/p01 - ysim[, 1]/p00)
        dl.drho <- (ysim[, 4]/p11 - ysim[, 3]/p10 -
                    ysim[, 2]/p01 + ysim[, 1]/p00) * dprob00

        rm(ysim)
        temp3 <- cbind(dl.dprob1, dl.dprob2, dl.drho)
          run.varcov <- ((ii-1) * run.varcov +
            temp3[, ind1$row.index] * temp3[, ind1$col.index]) / ii
      }
      wz <- if (intercept.only)
        matrix(colMeans(run.varcov),
               n, ncol(run.varcov), byrow = TRUE) else run.varcov

      wz <- wz * dthetas.detas[, ind1$row] * dthetas.detas[, ind1$col]
    }
    c(w) * wz
  }), list( .nsimEIM = nsimEIM ))))
}  # binom2.rho





 dnorm2 <- function(x, y, rho = 0, log = FALSE) {
  warning("decommissioning dnorm2() soon; use ",
          "dbinorm(..., cov12 = rho) instead")

  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  logdnorm2 <-
    (-0.5*(x * (x - 2*y*rho) + y^2) / (1.0 - rho^2)) - log(2 * pi) -
      0.5 * log1p(-rho^2)

  if (log.arg) {
    logdnorm2
  } else {
    exp(logdnorm2)
  }
}  # dnorm2




 pbinorm <-
 function(q1, q2,
          mean1 = 0, mean2 = 0,
          var1 = 1, var2 = 1,
          cov12 = 0) {



  sd1 <- sqrt(var1)
  sd2 <- sqrt(var2)
  rho <- cov12 / (sd1 * sd2)

  if (anyNA(q1)    || anyNA(q2)    ||
      anyNA(sd1)   || anyNA(sd2)   ||
      anyNA(mean1) || anyNA(mean2) || anyNA(rho))
    stop("no NAs allowed in arguments or variables 'q1', 'q2',",
         " 'mean1', 'mean2', 'sd1', 'sd2', 'cov12'")
  if (min(rho) < -1 || max(rho) > +1)
    stop("correlation 'rho' is out of range")


  if (length(mean1) > 1 && length(mean2) == 1 &&
      length(var1) == 1 && length(var2)  == 1 && length(cov12) == 1)
    warning("the call to pnorm2() seems based on the old version ",
            "of the arguments")

  LLL <- max(length(q1), length(q2),
             length(mean1), length(mean2),
             length(sd1), length(sd2),
             length(rho))
  if (length(q1)    != LLL) q1    <- rep_len(q1,     LLL)
  if (length(q2)    != LLL) q2    <- rep_len(q2,     LLL)
  if (length(mean1) != LLL) mean1 <- rep_len(mean1,  LLL)
  if (length(mean2) != LLL) mean2 <- rep_len(mean2,  LLL)
  if (length(sd1)   != LLL) sd1   <- rep_len(sd1,    LLL)
  if (length(sd2)   != LLL) sd2   <- rep_len(sd2,    LLL)
  if (length(rho)   != LLL) rho   <- rep_len(rho,    LLL)

  Zedd1 <- Z1 <- (q1 - mean1) / sd1
  Zedd2 <- Z2 <- (q2 - mean2) / sd2

  is.inf1.neg <- is.infinite(Z1) & Z1 < 0  # -Inf
  is.inf1.pos <- is.infinite(Z1) & Z1 > 0  # +Inf
  is.inf2.neg <- is.infinite(Z2) & Z2 < 0  # -Inf
  is.inf2.pos <- is.infinite(Z2) & Z2 > 0  # +Inf
  Zedd1[is.inf1.neg] <- 0
  Zedd1[is.inf1.pos] <- 0
  Zedd2[is.inf2.neg] <- 0
  Zedd2[is.inf2.pos] <- 0

  ans <- Zedd1
  singler <- ifelse(length(rho) == 1, 1, 0)
  answer <- .C("pnorm2ccc",
       ah = as.double(-Zedd1), ak = as.double(-Zedd2),
       r = as.double(rho), size = as.integer(LLL),
       singler = as.integer(singler),
       ans = as.double(ans))$ans
  if (any(answer < 0.0))
    warning("some negative values returned")

  answer[is.inf1.neg] <- 0
  answer[is.inf1.pos] <- pnorm(Z2[is.inf1.pos])  # pnorm(Z2[is.inf1.neg])
  answer[is.inf2.neg] <- 0
  answer[is.inf2.pos] <- pnorm(Z1[is.inf2.pos])  # pnorm(Z1[is.inf2.neg])

  answer
}  # pbinorm




 pnorm2 <- function(x1, x2,
                    mean1 = 0, mean2 = 0,
                    var1 = 1, var2 = 1,
                    cov12 = 0) {


  warning("decommissioning pnorm2() soon; use ",
          "pbinorm() instead")


  sd1 <- sqrt(var1)
  sd2 <- sqrt(var2)
  rho <- cov12 / (sd1 * sd2)

  if (anyNA(x1)    || anyNA(x2)    ||
      anyNA(sd1)   || anyNA(sd2)   ||
      anyNA(mean1) || anyNA(mean2) || anyNA(rho))
    stop("no NAs allowed in arguments or variables 'x1', 'x2',",
         " 'mean1', 'mean2', 'sd1', 'sd2', 'cov12'")
  if (min(rho) < -1 || max(rho) > +1)
    stop("correlation 'rho' is out of range")


  if (length(mean1) > 1 && length(mean2) == 1 &&
      length(var1) == 1 && length(var2)  == 1 && length(cov12) == 1)
    warning("the call to pnorm2() seems based on the old version ",
            "of the arguments")

  LLL <- max(length(x1), length(x2),
             length(mean1), length(mean2),
             length(sd1), length(sd2),
             length(rho))
  if (length(x1)    != LLL) x1    <- rep_len(x1,     LLL)
  if (length(x2)    != LLL) x2    <- rep_len(x2,     LLL)
  if (length(mean1) != LLL) mean1 <- rep_len(mean1,  LLL)
  if (length(mean2) != LLL) mean2 <- rep_len(mean2,  LLL)
  if (length(sd1)   != LLL) sd1   <- rep_len(sd1,    LLL)
  if (length(sd2)   != LLL) sd2   <- rep_len(sd2,    LLL)
  if (length(rho)   != LLL) rho   <- rep_len(rho,    LLL)

  Z1 <- (x1 - mean1) / sd1
  Z2 <- (x2 - mean2) / sd2

  ans <- Z1
  singler <- ifelse(length(rho) == 1, 1, 0)
  answer <- .C("pnorm2ccc",
       ah = as.double(-Z1), ak = as.double(-Z2), r = as.double(rho),
       size = as.integer(LLL), singler = as.integer(singler),
       ans = as.double(ans))$ans
  if (any(answer < 0.0))
    warning("some negative values returned")
  answer
}  # pnorm2






my.dbinom <- function(x,
                      size = stop("no 'size' argument"),
                      prob = stop("no 'prob' argument")) {

  exp(lfactorial(size) - lfactorial(size - x) - lfactorial(x) +
      x * logitlink(prob) + size * log1p(-prob))
}



 size.binomial <- function(prob = 0.5, link = "loglink") {
  if (any(prob <= 0 | prob >= 1))
    stop("some values of prob out of range")


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



  new("vglmff",
    blurb = c("Binomial with n unknown, prob known (prob = ",
            prob, ")\n",
            "Links:    ",
            namesof("size", link, tag = TRUE),
            " (treated as real-valued)\n",
            "Variance:  Var(Y) = size * prob * (1-prob);",
            " Var(size) is intractable"),
  initialize = eval(substitute(expression({
    predictors.names <- "size"
    extra$temp2 <- rep_len( .prob , n)

    if (is.null(etastart)) {
      nvec <- (y + 0.1) / extra$temp2
      etastart <- theta2eta(nvec, .link )
    }
  }), list( .prob = prob, .link = link ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    nvec <- eta2theta(eta, .link )
    nvec * extra$temp2
  }, list( .link = link ))),
  last = eval(substitute(expression({
    misc$link <- c(size = .link )

    misc$prob <- extra$temp2

  }), list( .link = link ))),
  linkfun = eval(substitute(
      function(mu, extra = NULL) {
    nvec <- mu / extra$temp2
    theta2eta(nvec, .link )
  }, list( .link = link ))),
  loglikelihood = eval(substitute(
      function(mu, y, w, residuals = FALSE, eta,
               extra = NULL, summation = TRUE) {
    nvec <- mu / extra$temp2
    if (residuals) {
      stop("loglikelihood resids not implemented")
    } else {

      ll.elts <-
        c(w) * (lfactorial(nvec) - lfactorial(y) -
                lfactorial(nvec - y) + y * logitlink( .prob ) +
                nvec * log1p(- ( .prob )))
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .prob = prob ))),
  vfamily = c("size.binomial"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    nvec <- eta2theta(eta, .link )
    okay1 <- all(is.finite(nvec)) && all( 0 < nvec)
    okay1
  }, list( .link = link ))),
  deriv = eval(substitute(expression({
    nvec <- mu/extra$temp2
    dldnvec <- digamma(nvec + 1) -
               digamma(nvec - y + 1) +
               log1p(-extra$temp2)
    dnvecdeta <- dtheta.deta(nvec, .link )
    c(w) * cbind(dldnvec * dnvecdeta)
  }), list( .link = link ))),
  weight = eval(substitute(expression({
      d2ldnvec2 <- trigamma(nvec + 1) -
          trigamma(nvec - y + 1)
    d2ldnvec2[y == 0] <- -sqrt( .Machine$double.eps )
    wz <- -c(w) * dnvecdeta^2 * d2ldnvec2
    wz
  }), list( .link = link ))))
}  # size.binomial






 dbetabinom.ab <-
  function(x, size, shape1, shape2, log = FALSE,
           Inf.shape = exp(20),  # 1e6, originally
           limit.prob = 0.5  # Strictly should be NaN
          ) {


  Bigg  <- Inf.shape
  Bigg2 <- Inf.shape  # big.shape  # exp(34)  # Found empirically
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)


  LLL <- max(length(x), length(size), length(shape1),
             length(shape2), length(limit.prob))
  if (length(x)      != LLL) x      <- rep_len(x,      LLL)
  if (length(size)   != LLL) size   <- rep_len(size,   LLL)
  if (length(shape1) != LLL) shape1 <- rep_len(shape1, LLL)
  if (length(shape2) != LLL) shape2 <- rep_len(shape2, LLL)
  if (length(limit.prob) != LLL)
    limit.prob <- rep_len(limit.prob, LLL)
  is.infinite.shape1 <- is.infinite(shape1)  # Includes -Inf !!
  is.infinite.shape2 <- is.infinite(shape2)

  ans <- x
  ans[TRUE] <- log(0)
  ans[is.na(x)]  <- NA
  ans[is.nan(x)] <- NaN


  ok0 <- !is.na(shape1) & !is.na(shape2) & !is.na(x) & !is.na(size)
  okk <- (round(x) == x) & (x >= 0) & (x <= size) &
         !is.infinite.shape1 & !is.infinite.shape2 & ok0
  if (any(okk)) {
    ans[okk] <- lchoose(size[okk], x[okk]) +
                lbeta(shape1[okk]             + x[okk],
                      shape2[okk] + size[okk] - x[okk]) -
                lbeta(shape1[okk], shape2[okk])


    endpt1 <- (x == size) &
              ((shape1 < 1/Bigg) | (shape2 < 1/Bigg)) & ok0
    if (any(endpt1)) {
      ans[endpt1] <-
        lgamma(  size[endpt1] + shape1[endpt1]) +
        lgamma(shape1[endpt1] + shape2[endpt1]) -
       (lgamma(  size[endpt1] + shape1[endpt1] + shape2[endpt1]) +
        lgamma(shape1[endpt1]))
    }  # endpt1




    endpt2 <-
        (x == 0) &
        ((shape1 < 1/Bigg) | (shape2 < 1/Bigg)) &
       ok0
    if (any(endpt2)) {
      ans[endpt2] <-
        lgamma(  size[endpt2] + shape2[endpt2]) +
        lgamma(shape1[endpt2] + shape2[endpt2]) -
       (lgamma(  size[endpt2] + shape1[endpt2] +
                                shape2[endpt2]) +
        lgamma(shape2[endpt2]))
    }  # endpt2





    endpt3 <- ((Bigg < shape1) | (Bigg < shape2)) & ok0
    if (any(endpt3)) {



      ans[endpt3] <- lchoose(size[endpt3], x[endpt3]) +
      lbeta(shape1[endpt3] + x[endpt3],
            shape2[endpt3] + size[endpt3] - x[endpt3]) -
      lbeta(shape1[endpt3], shape2[endpt3])
    }  # endpt3
  }  # if (any(okk))


  if (!log.arg) {
    ans <- exp(ans)
  }


  ok1 <- !is.infinite.shape1 &  is.infinite.shape2  # rho==0 & prob==0
  ok2 <-  is.infinite.shape1 & !is.infinite.shape2  # rho==0 & prob==1
  ok3 <-  Bigg2 < shape1     &  Bigg2 < shape2
  ok4 <-  is.infinite.shape1 &  is.infinite.shape2  # prob undefined



  if (any(ok3, na.rm = TRUE)) {
    ok33 <- !is.na(ok3) & ok3
    prob1 <- shape1[ok33] / (shape1[ok33] + shape2[ok33])
    ans[ok33] <- dbinom(x = x[ok33], size = size[ok33],
                        prob = prob1, log = log.arg)

    if (any(ok4)) {
      ans[ok4] <- dbinom(x = x[ok4], size = size[ok4],
                         prob = limit.prob[ok4],
                         log = log.arg)
    }
  }  # ok3

  if (any(ok1))
    ans[ok1] <- dbinom(x = x[ok1], size = size[ok1],
                       prob = 0,  # finite / (finite + Inf) == 0
                       log = log.arg)
  if (any(ok2))
    ans[ok2] <- dbinom(x = x[ok2], size = size[ok2],
                       prob = 1,  # Inf / (finite + Inf) == 1
                       log = log.arg)


  ans[is.na(shape1) | shape1 < 0] <- NaN
  ans[is.na(shape2) | shape2 < 0] <- NaN
  a.NA <- is.na(x) | is.na(size) |
          is.na(shape1) | is.na(shape2)
  ans[a.NA] <- NA  # 20180217

  ans
}  # dbetabinom.ab



 pbetabinom.ab <-
  function(q, size, shape1, shape2,
           limit.prob = 0.5,  # Should be NaN
           log.p = FALSE) {



  LLL <- max(length(q), length(size), length(shape1),
             length(shape2), length(limit.prob))

  if (length(q)       != LLL) q      <- rep_len(q,      LLL)
  if (length(shape1)  != LLL) shape1 <- rep_len(shape1, LLL)
  if (length(shape2)  != LLL) shape2 <- rep_len(shape2, LLL)
  if (length(size)    != LLL) size   <- rep_len(size,   LLL)
  if (length(limit.prob) != LLL)
    limit.prob <- rep_len(limit.prob, LLL)


  ind3 <- !is.na(size) & !is.na(q) & q > size
  q[ind3] <- size[ind3]  # Useful if q == Inf as it makes q finite

  ans <- q   # Retains names(q)
  ans[] <- NA_real_  # Handles NAs in size, shape1, etc. hopefully

  if (all(size       == size[1])   &&
      all(shape1     == shape1[1]) &&  # Cant handle NAs in comparison
      all(shape2     == shape2[1]) &&
      all(limit.prob == limit.prob[1]) &&
      !any(is.na(c(size, shape1, shape2, limit.prob)))) {
    qstar <- floor(q)
    maxqstar <- max(qstar, na.rm = TRUE)
    tmp2 <-
      dbetabinom.ab(if (maxqstar >= 0)
                      0:maxqstar else -1,
                    size = size[1],
                    shape1 = shape1[1],
                    shape2 = shape2[1],
                    limit.prob = limit.prob[1])
    unq <- unique(qstar[!is.na(qstar)])
    for (ii in unq) {
      index <- !is.na(qstar) & (qstar == ii)
      ans[index] <- if (ii >= 0) sum(tmp2[1:(1+ii)]) else tmp2
    }
  } else {
    for (ii in 1:LLL) {
      qstar <- floor(q[ii])
      qvec <- if (!is.na(qstar) && qstar >= 0)
                  0:qstar else NA
      ans[ii] <-
        sum(dbetabinom.ab(x = qvec,
                size = size[ii],
                shape1 = shape1[ii],
                shape2 = shape2[ii],
                limit.prob = limit.prob[ii]))
    }
  }



  ind4 <- !is.na(size) & !is.na(q) &
          !is.na(shape1) & !is.na(shape2) & q < 0
  ans[ind4] <- 0




  if (log.p) log(ans) else ans
}  # pbetabinom.ab



rbetabinom.ab <-
  function(n, size, shape1, shape2,
    limit.prob = 0.5,  # Strictly should be NaN
    .dontuse.prob = NULL   # 20180814 temporary!!
    ) {


  use.n <- if ((length.n <- length(n)) > 1)
     length.n else
  if (!is.Numeric(n, integer.valued = TRUE,
        length.arg = 1, positive = TRUE))
    stop("bad input for argument 'n'") else n
  if (length(size)   != use.n)
      size   <- rep_len(size,   use.n)
  if (length(shape1) != use.n)
      shape1 <- rep_len(shape1, use.n)
  if (length(shape2) != use.n)
      shape2 <- rep_len(shape2, use.n)
  if (length(limit.prob) != use.n)
    limit.prob <- rep_len(limit.prob, use.n)

  ans <- rep_len(NA_real_, use.n)
  ind3 <- !is.na(shape1) & !is.na(shape2) &
((is.infinite(shape1) & is.infinite(shape2))) # |

  if (sum.ind3 <- sum(ind3))
    ans[ind3]<- rbinom(sum.ind3, size = size[ind3],
                       prob = limit.prob[ind3])

  if (ssum.ind3 <- sum(!ind3))
    ans[!ind3] <-
      rbinom(ssum.ind3, size = size[!ind3],
             prob = rbeta(n = ssum.ind3,
                          shape1 = shape1[!ind3],
                          shape2 = shape2[!ind3]))

  ans[is.na(shape1) | shape1 < 0] <- NaN
  ans[is.na(shape2) | shape2 < 0] <- NaN

  ans
}  # rbetabinom.ab






dbetabinom <-
  function(x, size, prob, rho = 0,
           log = FALSE) {
  dbetabinom.ab(x = x, size = size,
     shape1 = prob * (1 - rho) / rho,
     shape2 = (1 - prob) * (1 - rho) / rho,
     limit.prob = prob,  # 20180216, as rho = 0.
     log = log)
}



pbetabinom <-
  function(q, size, prob, rho = 0,
           log.p = FALSE) {
  pbetabinom.ab(q, size = size,
     shape1 = prob * (1 - rho) / rho,
     shape2 = (1 - prob) * (1 - rho) / rho,
     limit.prob = prob,  # 20180217, as rho = 0.
     log.p = log.p)
}



rbetabinom <-
  function(n, size, prob, rho = 0
          ) {
  rbetabinom.ab(n, size = size,
     shape1 = prob *(1 - rho) / rho,
     shape2 = (1 - prob)*(1 - rho) / rho,
     limit.prob = prob  # 20180217, as rho = 0.
               )
}






 Ebbin.ab <-
    function(nvec, shape1, shape2, first) {



  NN <- length(nvec)
  ans <- rep_len(0.0, NN)
  if (first) {
    for (ii in 1:NN) {
      temp639 <- lbeta(shape1[ii], shape2[ii])
      yy <- 0:nvec[ii]
      ans[ii] <- ans[ii] +
          sum(trigamma(shape1[ii] + yy) *
          exp(lchoose(nvec[ii], yy) +
              lbeta(shape1[ii] + yy,
                    shape2[ii] + nvec[ii] - yy) -
              temp639))
    }  # ii
  } else {
    for (ii in 1:NN) {
      temp639 <- lbeta(shape1[ii], shape2[ii])
      yy <- 0:nvec[ii]
      ans[ii] <- ans[ii] +
        sum(trigamma(nvec[ii] + shape2[ii] - yy) *
            exp(lchoose(nvec[ii], yy) +
                lbeta(shape1[ii] + yy,
                      shape2[ii] + nvec[ii] - yy) -
                temp639))
    }
  }
  ans
}  # Ebbin.ab 



betabinomialff.control <-
    function(save.weights = TRUE, ...) {
    list(save.weights = save.weights)
}




 betabinomialff <-
   function(lshape1 = "loglink",
            lshape2 = "loglink",
            ishape1 = 1, ishape2 = NULL,
            imethod = 1,
            ishrinkage = 0.95,
            nsimEIM = NULL,
            zero = NULL) {




  lshape1 <- as.list(substitute(lshape1))
  eshape1 <- link2list(lshape1)
  lshape1 <- attr(eshape1, "function.name")


  lshape2 <- as.list(substitute(lshape2))
  eshape2 <- link2list(lshape2)
  lshape2 <- attr(eshape2, "function.name")


  if (!is.Numeric(ishape1, positive = TRUE))
    stop("bad input for argument 'ishape1'")
  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE,
                  positive = TRUE) ||
     imethod > 3)
    stop("argument 'imethod' must be 1, 2 or 3")

   if (length(ishape2) &&
       !is.Numeric(ishape2, positive = TRUE))
    stop("bad input for argument 'ishape2'")

  if (!is.null(nsimEIM)) {
    if (!is.Numeric(nsimEIM, length.arg = 1,
                    integer.valued = TRUE))
      stop("bad input for argument 'nsimEIM'")
    if (nsimEIM <= 10)
      warning("'nsimEIM' should be an integer ",
              "greater than 10, say")
  }


  new("vglmff",
  blurb = c("Beta-binomial model\n",
            "Links:    ",
   namesof("shape1", lshape1, eshape1), ", ",
   namesof("shape2", lshape2, eshape2), "\n",
   "Mean:     mu = shape1 / (shape1+shape2)", "\n",
   "Variance: mu * (1-mu) * (1+(w-1)*rho) / w, ",
   "where rho = 1 / (shape1+shape2+1)"),
  constraints = eval(substitute(expression({
    constraints <-
      cm.zero.VGAM(constraints, x = x, .zero ,
             M = M, M1 = 2,
             predictors.names = predictors.names)
  }), list( .zero = zero ))),


  infos = eval(substitute(function(...) {
    list(M1 = 2,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("shape1", "shape2"),
         lshape1 = .lshape1 ,
         lshape2 = .lshape2 ,
         zero = .zero )
  },
  list( .zero = zero,
        .lshape1 = lshape1,
        .lshape2 = lshape2 ))),


  initialize = eval(substitute(expression({
    if (!all(w == 1))
      extra$orig.w <- w

    if (is.null( .nsimEIM )) {
      save.weights <-
      control$save.weights <- FALSE
    }

    mustart.orig <- mustart
    eval(binomialff()@initialize)
    if (length(mustart.orig))
      mustart <- mustart.orig  # Retainif inputted
    predictors.names <-
      c(namesof("shape1", .lshape1 , .eshape1 ,
                tag = FALSE),
        namesof("shape2", .lshape2 , .eshape2 ,
                tag = FALSE))

    if (!length(etastart)) {

      mustart.use <- if (length(mustart.orig))
                     mustart.orig else
                     mustart

      shape1 <- rep_len( .ishape1 , n)
      shape2 <- if (length( .ishape2 )) {
                  rep_len( .ishape2 , n)
        } else if (length(mustart.orig)) {
          shape1 * (1 / mustart.use - 1)
        } else if ( .imethod == 1) {
          shape1 * (1 / weighted.mean(y, w) - 1)
        } else if ( .imethod == 2) {
          tmp7 <- (1 - .ishrinkage ) * y +
              .ishrinkage * weighted.mean(y, w)
           shape1 * (1 / tmp7 - 1)
        } else {
          shape1 *
          (1 / weighted.mean(mustart.use, w) - 1)
        }
      ycounts <- if (is.numeric(extra$orig.w))
                 y * w / extra$orig.w else
                 y * w # Convert propns to counts
      if (max(abs(ycounts-round(ycounts))) > 1e-6)
         warning("the response (as counts) ",
         "does not appear to be integer-valued.",
         " Am rounding to integer values.")
      ycounts <- round(ycounts)  # Now an integer
    etastart <-
    cbind(theta2eta(shape1, .lshape1 , .eshape1 ),
          theta2eta(shape2, .lshape2 , .eshape2 ))
      mustart <- NULL  # Since etastart computed.
    }
  }),
  list( .lshape1 = lshape1, .lshape2 = lshape2,
        .eshape1 = eshape1, .eshape2 = eshape2,
        .ishape1 = ishape1, .ishape2 = ishape2,
        .nsimEIM = nsimEIM,
        .imethod = imethod,
        .ishrinkage = ishrinkage ))),
  linkinv = eval(substitute(
    function(eta, extra = NULL) {
      shape1 <- eta2theta(eta[, 1], .lshape1 ,
                                    .eshape1 )
      shape2 <- eta2theta(eta[, 2], .lshape2 ,
                                    .eshape2 )
      shape1 / (shape1 + shape2)
  },
  list( .lshape1 = lshape1, .lshape2 = lshape2,
        .eshape1 = eshape1, .eshape2 = eshape2 ))),
  last = eval(substitute(expression({
    misc$link <-    c("shape1" = .lshape1 ,
                      "shape2" = .lshape2 )
    misc$earg <- list("shape1" = .eshape1 ,
                      "shape2" = .eshape2 )

    shape1 <- eta2theta(eta[, 1], .lshape1 ,
                                  .eshape1 )
    shape2 <- eta2theta(eta[, 2], .lshape2 ,
                                  .eshape2 )
    misc$rho <- 1 / (shape1 + shape2 + 1)
    misc$expected <- TRUE
    misc$nsimEIM <- c( .nsimEIM )
    misc$zero <- c( .zero )
  }),
  list( .lshape1 = lshape1, .lshape2 = lshape2,
        .eshape1 = eshape1, .eshape2 = eshape2,
        .nsimEIM = nsimEIM, .zero = zero ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL, summation = TRUE) {
      ycounts <- if (is.numeric(extra$orig.w))
                 y * w / extra$orig.w else
                 y * w # Convert propns to counts

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

    shape1 <- eta2theta(eta[, 1], .lshape1 ,
                                  .eshape1 )
    shape2 <- eta2theta(eta[, 2], .lshape2 ,
                                  .eshape2 )
    nvec <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else
              round(w)
    if (residuals) {
      stop("loglikelihood resids not implemented")
    } else {
      ll.elts <-
        (if (is.numeric(extra$orig.w))
           extra$orig.w else 1) *
        dbetabinom.ab(ycounts, size = nvec,
                      shape1 = shape1,
                      shape2 = shape2, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  },
  list( .lshape1 = lshape1, .lshape2 = lshape2,
        .eshape1 = eshape1, .eshape2 = eshape2 ))),
  vfamily = c("betabinomialff"),
  validparams = eval(substitute(
      function(eta, y, extra = NULL) {
    shape1 <- eta2theta(eta[, 1], .lshape1 ,
                                  .eshape1 )
    shape2 <- eta2theta(eta[, 2], .lshape2 ,
                                  .eshape2 )
    okay1 <-
    all(is.finite(shape1)) && all(0 < shape1) &&
    all(is.finite(shape2)) && all(0 < shape2)
    okay1
  },
  list( .lshape1 = lshape1, .lshape2 = lshape2,
        .eshape1 = eshape1, .eshape2 = eshape2 ))),





  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <-
              object@prior.weights) > 0)
              pwts else
              weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    w <- pwts
    eta <- predict(object)
    extra <- object@extra
    shape1 <- eta2theta(eta[, 1], .lshape1 ,
                                  .eshape1 )
    shape2 <- eta2theta(eta[, 2], .lshape2 ,
                                  .eshape2 )
    nvec <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else
              round(w)
    rbetabinom.ab(nsim * length(shape1),
                  size = nvec,
                  shape1 = shape1,
                  shape2 = shape2)
  },
  list( .lshape1 = lshape1, .lshape2 = lshape2,
        .eshape1 = eshape1, .eshape2 = eshape2 ))),



  deriv = eval(substitute(expression({
    nvec <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else
              round(w)
    ycounts <- if (is.numeric(extra$orig.w))
              y * w / extra$orig.w else
              y * w # Convert propns to counts
    shape1 <- eta2theta(eta[, 1], .lshape1 ,
                                  .eshape1 )
    shape2 <- eta2theta(eta[, 2], .lshape2 ,
                                  .eshape2 )

    dshape1.deta <- dtheta.deta(shape1, .lshape1 ,
                                        .eshape1 )
    dshape2.deta <- dtheta.deta(shape2, .lshape2 ,
                                        .eshape2 )

    dl.dshape1 <-
        digamma(shape1 + ycounts) -
        digamma(shape1 + shape2 + nvec) -
        digamma(shape1) +
        digamma(shape1 + shape2)
    dl.dshape2 <-
        digamma(nvec + shape2 - ycounts) -
        digamma(shape1 + shape2 + nvec) -
        digamma(shape2) +
        digamma(shape1 + shape2)

    (if (is.numeric(extra$orig.w))
       extra$orig.w else 1) *
    cbind(dl.dshape1 * dshape1.deta,
          dl.dshape2 * dshape2.deta)
  }),
  list( .lshape1 = lshape1, .lshape2 = lshape2,
        .eshape1 = eshape1, .eshape2 = eshape2 ))),
  weight = eval(substitute(expression({
    if (is.null( .nsimEIM)) {
      wz <- matrix(NA_real_, n, dimm(M))  # 3
      wz[, iam(1, 1, M)] <-
        -(Ebbin.ab(nvec, shape1, shape2, TRUE) -
        trigamma(shape1 + shape2 + nvec) -
        trigamma(shape1) +
        trigamma(shape1 + shape2)) *
        dshape1.deta^2
      wz[, iam(2, 2, M)] <-
        -(Ebbin.ab(nvec, shape1, shape2, FALSE) -
        trigamma(shape1 + shape2 + nvec) -
        trigamma(shape2) +
        trigamma(shape1 + shape2)) *
        dshape2.deta^2
      wz[, iam(2, 1, M)] <-
          -(trigamma(shape1 + shape2) -
            trigamma(shape1 + shape2 + nvec)) *
            dshape1.deta * dshape2.deta
      wz * (if (is.numeric(extra$orig.w))
                extra$orig.w else 1)
    } else {
      run.varcov <- 0
      ind1 <- iam(NA, NA, M = M, both = TRUE,
                  diag = TRUE)
      dthetas.detas <-
        cbind(dshape1.deta, dshape2.deta)

      for (ii in 1:( .nsimEIM )) {
        ysim <- rbetabinom.ab(n, size = nvec,
                              shape1 = shape1,
                              shape2 = shape2)
        checkargs = ( .checkargs )
        dl.dshape1 <- digamma(shape1 + ysim) -
           digamma(shape1 + shape2 + nvec) -
           digamma(shape1) +
           digamma(shape1 + shape2)
        dl.dshape2 <-
            digamma(nvec+shape2 - ysim) -
            digamma(shape1 + shape2 + nvec) -
            digamma(shape2) +
            digamma(shape1 + shape2)
        rm(ysim)
        temp3 <- cbind(dl.dshape1,
                       dl.dshape2)  # n x M
        run.varcov <- ((ii-1) * run.varcov +
                     temp3[, ind1$row.index]*
                     temp3[, ind1$col.index]) / ii
      }
      wz <- if (intercept.only)
          matrix(colMeans(run.varcov),
                 n, ncol(run.varcov),
                 byrow = TRUE) else run.varcov

      wz <- wz * dthetas.detas[, ind1$row] *
                 dthetas.detas[, ind1$col]
      wz * (if (is.numeric(extra$orig.w))
              extra$orig.w else 1)
    }
  }),
  list( .lshape1 = lshape1, .lshape2 = lshape2,
        .eshape1 = eshape1, .eshape2 = eshape2,
        .nsimEIM = nsimEIM ))))
}  # betabinomialff






 betageometric <-
  function(lprob = "logitlink", lshape = "loglink",
           iprob = NULL, ishape = 0.1,
           moreSummation = c(2, 100),
           tolerance = 1.0e-10,
           zero = NULL) {
  lprob <- as.list(substitute(lprob))
  eprob <- link2list(lprob)
  lprob <- attr(eprob, "function.name")

  lshape <- as.list(substitute(lshape))
  eshape <- link2list(lshape)
  lshape <- attr(eshape, "function.name")


  if (!is.Numeric(ishape, positive = TRUE))
    stop("bad input for argument 'ishape'")
  if (!is.Numeric(moreSummation, positive = TRUE,
                  length.arg = 2, integer.valued = TRUE))
    stop("bad input for argument 'moreSummation'")
  if (!is.Numeric(tolerance, positive = TRUE, length.arg = 1) ||
      1.0 - tolerance >= 1.0)
      stop("bad input for argument 'tolerance'")



  new("vglmff",
  blurb = c("Beta-geometric distribution\n",
            "Links:    ",
   namesof("prob",  lprob,  earg = eprob), ", ",
   namesof("shape", lshape, earg = eshape)),

  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints,
             x = x, .zero , M = M, M1 = 2,
             predictors.names = predictors.names)
  }), list( .zero = zero ))),


  infos = eval(substitute(function(...) {
    list(M1 = 2,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("prob", "shape"),
         lprob  = .lprob ,
         lshape = .lshape ,
         zero = .zero )
  }, list( .lprob = lprob, .lshape = lshape,
           .zero = zero ))),


  initialize = eval(substitute(expression({
    eval(geometric()@initialize)

    predictors.names <-
         c(namesof("prob",  .lprob  , .eprob  , tag = FALSE),
           namesof("shape", .lshape , .eshape , tag = FALSE))

    if (length( .iprob ))
      prob.init <- rep_len( .iprob , n)

    if (!length(etastart) ||
      NCOL(etastart) != 2) {
      shape.init <- rep_len( .ishape , n)
      etastart <-
        cbind(theta2eta(prob.init,  .lprob ,  earg = .eprob ),
              theta2eta(shape.init, .lshape , earg = .eshape ))
      }
  }), list( .iprob = iprob, .ishape = ishape, .lprob = lprob,
            .eprob = eprob, .eshape = eshape,
            .lshape = lshape ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    prob  <- eta2theta(eta[, 1], .lprob ,  earg = .eprob )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )
    mymu <- (1-prob) / (prob - shape)
    ifelse(mymu >= 0, mymu, NA)
  }, list( .lprob = lprob, .lshape = lshape,
           .eprob = eprob, .eshape = eshape ))),
  last = eval(substitute(expression({
    misc$link <-    c("prob" = .lprob , "shape" = .lshape )

    misc$earg <- list("prob" = .eprob , "shape" = .eshape )

    if (intercept.only) {
      misc$shape1 <- shape1[1]  # These quantities computed in @deriv
      misc$shape2 <- shape2[1]
    }
    misc$expected <- TRUE
    misc$tolerance <- .tolerance
    misc$zero <- .zero
    misc$moreSummation = .moreSummation
  }),
  list( .lprob = lprob, .lshape = lshape,
        .eprob = eprob, .eshape = eshape,
        .tolerance = tolerance,
        .moreSummation = moreSummation,
        .zero = zero ))),
  loglikelihood = eval(substitute(
      function(mu, y, w, residuals = FALSE, eta,
               extra = NULL, summation = TRUE) {
   prob  <- eta2theta(eta[, 1], .lprob  , .eprob  )
   shape <- eta2theta(eta[, 2], .lshape , .eshape )
    ans <- log(prob)
    maxy <- max(y)
    if (residuals) {
      stop("loglikelihood resids not implemented")
    } else {
    for (ii in 1:maxy) {
      index <- (ii <= y)
      ans[index] <- ans[index] +
    log1p(-prob[index] + (ii-1) * shape[index]) -
    log1p((ii-1) * shape[index])
    }
    ans <- ans - log1p((y+1-1) * shape)




      ll.elts <- w * ans
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lprob = lprob, .lshape = lshape,
           .eprob = eprob, .eshape = eshape ))),
  vfamily = c("betageometric"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    prob  <- eta2theta(eta[, 1], .lprob ,  earg = .eprob )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )
    okay1 <- all(is.finite(prob )) && all(0 < prob  & prob < 1) &&
             all(is.finite(shape)) && all(0 < shape)
    okay1
  },
  list( .lprob = lprob, .lshape = lshape,
        .eprob = eprob, .eshape = eshape ))),




  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    prob  <- eta2theta(eta[, 1], .lprob  , earg = .eprob  )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )
    rbetageom(nsim * length(shape),
              shape1 = shape, shape2 = shape)
  },
  list( .lprob = lprob, .lshape = lshape,
        .eprob = eprob, .eshape = eshape ))),




  deriv = eval(substitute(expression({
    prob  <- eta2theta(eta[, 1], .lprob ,  earg = .eprob )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )
    shape1 <-      prob  / shape
    shape2 <- (1 - prob) / shape
    dprob.deta  <- dtheta.deta(prob , .lprob  , earg = .eprob  )
    dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
    dl.dprob <- 1 / prob
    dl.dshape <- 0 * y
    maxy <- max(y)
    for (ii in 1:maxy) {
      index <- (ii <= y)
      dl.dprob[index] <- dl.dprob[index] -
                         1/(1-prob[index]+(ii-1) * shape[index])
      dl.dshape[index] <- dl.dshape[index] +
                         (ii-1)/(1-prob[index]+(ii-1) * shape[index]) -
                         (ii-1)/(1+(ii-1) * shape[index])
    }
    dl.dshape <- dl.dshape - (y+1 -1)/(1+(y+1 -1) * shape)
    c(w) * cbind(dl.dprob * dprob.deta,
                 dl.dshape * dshape.deta)
  }),
  list( .lprob = lprob, .lshape = lshape,
        .eprob = eprob, .eshape = eshape ))),
  weight = eval(substitute(expression({
    wz <- matrix(0, n, dimm(M))  #3=dimm(2)
    wz[, iam(1, 1, M)] <- 1 / prob^2
    moresum <- .moreSummation
    maxsummation <- round(maxy * moresum[1] + moresum[2])
    for (ii in 3:maxsummation) {
      temp7 <- 1 - pbetageom(q = ii-1-1, shape1 = shape1,
                                         shape2 = shape2)
      denom1 <- (1-prob+(ii-2)*shape)^2
      denom2 <- (1+(ii-2)*shape)^2
      wz[, iam(1, 1, M)] <- wz[, iam(1, 1, M)] + temp7 / denom1
      wz[, iam(1, 2, M)] <- wz[, iam(1, 2, M)] - (ii-2) * temp7 / denom1
      wz[, iam(2, 2, M)] <- wz[, iam(2, 2, M)] +
                        (ii-2)^2 * temp7 / denom1 -
                        (ii-1)^2 * temp7 / denom2
      if (max(temp7) < .tolerance ) break
    }
    ii <- 2
    temp7 <- 1 - pbetageom(q=ii-1-1, shape1 = shape1, shape2 = shape2)
    denom1 <- (1-prob+(ii-2)*shape)^2
    denom2 <- (1+(ii-2)*shape)^2

    wz[, iam(1, 1, M)] <- wz[, iam(1, 1, M)] + temp7 / denom1
    wz[, iam(2, 2, M)] <- wz[, iam(2, 2, M)] - (ii-1)^2 * temp7 / denom2
    wz[, iam(1, 1, M)] <- wz[, iam(1, 1, M)] * dprob.deta^2
    wz[, iam(2, 2, M)] <- wz[, iam(2, 2, M)] * dshape.deta^2
    wz[, iam(2, 1, M)] <- wz[, iam(2, 1, M)] * dprob.deta * dshape.deta
    c(w) * wz
  }),
  list( .lprob = lprob, .lshape = lshape,
        .eprob = eprob, .eshape = eshape,
        .moreSummation = moreSummation,
        .tolerance = tolerance ))))
}  # betageometric





 seq2binomial <-
   function(lprob1 = "logitlink",
            lprob2 = "logitlink",
            iprob1 = NULL,    iprob2 = NULL,
            parallel = FALSE,
            zero = NULL) {
  apply.parint <- TRUE

  lprob1 <- as.list(substitute(lprob1))
  eprob1 <- link2list(lprob1)
  lprob1 <- attr(eprob1, "function.name")

  lprob2 <- as.list(substitute(lprob2))
  eprob2 <- link2list(lprob2)
  lprob2 <- attr(eprob2, "function.name")


  if (length(iprob1) &&
     (!is.Numeric(iprob1, positive = TRUE) ||
     max(iprob1) >= 1))
    stop("bad input for argument 'iprob1'")
  if (length(iprob2) &&
     (!is.Numeric(iprob2, positive = TRUE) ||
     max(iprob2) >= 1))
    stop("bad input for argument 'iprob2'")


  new("vglmff",
  blurb = c("Sequential binomial distribution ",
            "(Crowder and Sweeting, 1989)\n",
            "Links:    ",
  namesof("prob1", lprob1, earg = eprob1), ", ",
  namesof("prob2", lprob2, earg = eprob2)),

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

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("prob1", "prob2"),
         lprob1 = .lprob1 ,
         lprob2 = .lprob2 ,
         zero = .zero )
  }, list( .zero = zero ))),


  initialize = eval(substitute(expression({
    if (!is.vector(w))
      stop("the 'weights' argument must be a vector")

    if (any(abs(w - round(w)) > 1e-6))
      stop("the 'weights' argument does not seem to be integer-valued")


    if (ncol(y <- cbind(y)) != 2)
      stop("the response must be a 2-column matrix")

    if (any(y < 0 | y > 1))
      stop("the response must have values between 0 and 1")


    w <- round(w)
    rvector <- w * y[, 1]
    if (any(abs(rvector - round(rvector)) > 1.0e-8))
      warning("number of successes in column one ",
              "should be integer-valued")
    svector <- rvector * y[, 2]
    if (any(abs(svector - round(svector)) > 1.0e-8))
      warning("number of successes in column",
              " two should be integer-valued")

    predictors.names <-
        c(namesof("prob1", .lprob1 , .eprob1 , tag = FALSE),
          namesof("prob2", .lprob2 , .eprob2 , tag = FALSE))

    prob1.init <- if (is.Numeric( .iprob1))
                   rep_len( .iprob1 , n) else
                   rep_len(weighted.mean(y[, 1], w = w), n)
    prob2.init <- if (is.Numeric( .iprob2 ))
                   rep_len( .iprob2 , n) else
                   rep_len(weighted.mean(y[, 2], w = w*y[, 1]), n)

    if (!length(etastart)) {
      etastart <-
        cbind(theta2eta(prob1.init, .lprob1 , earg = .eprob1 ),
              theta2eta(prob2.init, .lprob2 , earg = .eprob2 ))
    }
  }),
  list( .iprob1 = iprob1, .iprob2 = iprob2,
        .lprob1 = lprob1, .lprob2 = lprob2,
        .eprob1 = eprob1, .eprob2 = eprob2 ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    prob1 <- eta2theta(eta[, 1], .lprob1 , .eprob1 )
    prob2 <- eta2theta(eta[, 2], .lprob2 , .eprob2 )
    cbind(prob1, prob2)
  }, list( .lprob1 = lprob1, .lprob2 = lprob2,
           .eprob1 = eprob1, .eprob2 = eprob2 ))),
  last = eval(substitute(expression({
    misc$link <-    c("prob1" = .lprob1 , "prob2" = .lprob2 )

    misc$earg <- list("prob1" = .eprob1 , "prob2" = .eprob2 )

    misc$expected <- TRUE
    misc$zero <- .zero
    misc$parallel <- .parallel
    misc$apply.parint <- .apply.parint
  }), list( .lprob1 = lprob1, .lprob2 = lprob2,
            .eprob1 = eprob1, .eprob2 = eprob2,
            .parallel = parallel,
            .apply.parint = apply.parint,
            .zero = zero ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta, extra = NULL,
             summation = TRUE) {
    prob1 <- eta2theta(eta[, 1], .lprob1 , earg = .eprob1 )
    prob2 <- eta2theta(eta[, 2], .lprob2 , earg = .eprob2 )

    smallno <- 100 * .Machine$double.eps
    prob1 <- pmax(prob1,   smallno)
    prob1 <- pmin(prob1, 1-smallno)
    prob2 <- pmax(prob2,   smallno)
    prob2 <- pmin(prob2, 1-smallno)
    mvector <- w
    rvector <- w * y[, 1]
    svector <- rvector * y[, 2]
    if (residuals) {
      stop("loglikelihood resids not implemented")
    } else {
      ans1 <-
          dbinom(rvector, size = mvector,
                 prob = prob1, log = TRUE) +
          dbinom(svector, size = rvector,
                 prob = prob2, log = TRUE)

      ll.elts <- ans1
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  },
  list( .lprob1 = lprob1, .lprob2 = lprob2,
        .eprob1 = eprob1, .eprob2 = eprob2 ))),
  vfamily = c("seq2binomial"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    prob1 <- eta2theta(eta[, 1], .lprob1 , earg = .eprob1 )
    prob2 <- eta2theta(eta[, 2], .lprob2 , earg = .eprob2 )
    okay1 <- all(is.finite(prob1)) && all(0 < prob1 & prob1 < 1) &&
             all(is.finite(prob2)) && all(0 < prob2 & prob2 < 1)
    okay1
  }, list( .lprob1 = lprob1, .lprob2 = lprob2,
           .eprob1 = eprob1, .eprob2 = eprob2 ))),
  deriv = eval(substitute(expression({
    prob1 <- eta2theta(eta[, 1], .lprob1 , .eprob1 )
    prob2 <- eta2theta(eta[, 2], .lprob2 , .eprob2 )
    smallno <- 100 * .Machine$double.eps
    prob1 <- pmax(prob1, smallno)
    prob1 <- pmin(prob1, 1-smallno)
    prob2 <- pmax(prob2, smallno)
    prob2 <- pmin(prob2, 1-smallno)

    dprob1.deta <- dtheta.deta(prob1, .lprob1, .eprob1)
    dprob2.deta <- dtheta.deta(prob2, .lprob2, .eprob2)

    mvector <- w
    rvector <- w * y[, 1]
    svector <- rvector * y[, 2]

    dl.dprob1 <- rvector / prob1 - (mvector-rvector) / (1-prob1)
    dl.dprob2 <- svector / prob2 - (rvector-svector) / (1-prob2)

    cbind(dl.dprob1 * dprob1.deta, dl.dprob2 * dprob2.deta)
  }), list( .lprob1 = lprob1, .lprob2 = lprob2,
            .eprob1 = eprob1, .eprob2 = eprob2 ))),
  weight = eval(substitute(expression({
    wz <- matrix(0, n, M)
    wz[, iam(1, 1, M)] <- (dprob1.deta^2) / (prob1 * (1-prob1))
    wz[, iam(2, 2, M)] <- (dprob2.deta^2) * prob1 / (prob2 * (1-prob2))
    c(w) * wz
  }), list( .lprob1 = lprob1, .lprob2 = lprob2,
            .eprob1 = eprob1, .eprob2 = eprob2 ))))
}  # seq2binomial



 zipebcom   <-
    function(lmu12 = "clogloglink",
             lphi12 = "logitlink",
             loratio = "loglink",
             imu12 = NULL, iphi12 = NULL,
             ioratio = NULL,
             zero = c("phi12", "oratio"),
             tol = 0.001, addRidge = 0.001) {


  lmu12 <- as.list(substitute(lmu12))
  emu12 <- link2list(lmu12)
  lmu12 <- attr(emu12, "function.name")

  lphi12 <- as.list(substitute(lphi12))
  ephi12 <- link2list(lphi12)
  lphi12 <- attr(ephi12, "function.name")

  loratio <- as.list(substitute(loratio))
  eoratio <- link2list(loratio)
  loratio <- attr(eoratio, "function.name")



  if (!is.Numeric(tol, positive = TRUE, length.arg = 1) ||
      tol > 0.1)
      stop("bad input for argument 'tol'")
  if (!is.Numeric(addRidge, length.arg = 1, positive = TRUE) ||
      addRidge > 0.5)
    stop("bad input for argument 'addRidge'")

  if (lmu12 != "clogloglink")
    warning("argument 'lmu12' should be 'clogloglink'")


  new("vglmff",
  blurb = c("Exchangeable bivariate ", lmu12,
            " odds-ratio model based on\n",
            "a zero-inflated Poisson distribution\n\n",
            "Links:    ",
  namesof("mu12",   lmu12,   earg = emu12), ", ",
  namesof("phi12",  lphi12,  earg = ephi12), ", ",
  namesof("oratio", loratio, earg = eoratio)),

  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints,
            x = x, .zero , M = M, M1 = 3,
            predictors.names = predictors.names)
    }), list( .zero = zero ))),


  infos = eval(substitute(function(...) {
    list(M1 = 3,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("mu12", "phi12", "oratio"),
         lmu12   = .lmu12 ,
         lphi12  = .lphi12 ,
         loratio = .loratio ,
         zero = .zero )
  }, list( .zero = zero,
           .lmu12 = lmu12, .lphi12 = lphi12, .loratio = loratio
         ))),



  initialize = eval(substitute(expression({
    eval(process.binomial2.data.VGAM)

    predictors.names <- c(
             namesof("mu12",   .lmu12   , earg = .emu12   , tag = FALSE),
             namesof("phi12",  .lphi12  , earg = .ephi12  , tag = FALSE),
             namesof("oratio", .loratio , earg = .eoratio , tag = FALSE))

    propY1.eq.0 <- weighted.mean(y[,'00'], w) + weighted.mean(y[,'01'], w)
    propY2.eq.0 <- weighted.mean(y[,'00'], w) + weighted.mean(y[,'10'], w)
    if (length( .iphi12) && any( .iphi12 > propY1.eq.0))
      warning("iphi12 must be less than the sample proportion of Y1==0")
    if (length( .iphi12) && any( .iphi12 > propY2.eq.0))
      warning("iphi12 must be less than the sample proportion of Y2==0")

    if (!length(etastart)) {
        pstar.init <- ((mu[, 3]+mu[, 4]) + (mu[, 2]+mu[, 4])) / 2
        phi.init <- if (length(.iphi12)) rep_len( .iphi12 , n) else
            min(propY1.eq.0 * 0.95, propY2.eq.0 * 0.95, pstar.init/1.5)
        oratio.init <- if (length( .ioratio)) rep_len( .ioratio , n) else
                  mu[, 4]*mu[, 1]/(mu[, 2]*mu[, 3])
        mu12.init <- if (length(.imu12)) rep_len( .imu12 , n) else
            pstar.init / (1-phi.init)

        etastart <- cbind(
            theta2eta(mu12.init,   .lmu12 ,   earg = .emu12 ),
            theta2eta(phi.init,    .lphi12,  earg = .ephi12),
            theta2eta(oratio.init, .loratio, earg = .eoratio))
        mustart <- NULL
      }
  }), list( .lmu12 = lmu12, .lphi12 = lphi12, .loratio = loratio,
            .emu12 = emu12, .ephi12 = ephi12, .eoratio = eoratio,
            .imu12 = imu12, .iphi12 = iphi12, .ioratio = ioratio ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    A1vec  <- eta2theta(eta[, 1], .lmu12  , earg = .emu12  )
    phivec <- eta2theta(eta[, 2], .lphi12 , earg = .ephi12 )
    pmargin <- matrix((1 - phivec) * A1vec, nrow(eta), 2)
    oratio <- eta2theta(eta[, 3], .loratio , earg = .eoratio )
    a.temp <- 1 + (pmargin[, 1]+pmargin[, 2])*(oratio-1)
    b.temp <- -4 * oratio * (oratio-1) * pmargin[, 1] * pmargin[, 2]
    temp <- sqrt(a.temp^2 + b.temp)
    pj4 <- ifelse(abs(oratio-1) < .tol, pmargin[, 1]*pmargin[, 2],
                 (a.temp-temp)/(2*(oratio-1)))
    pj2 <- pmargin[, 2] - pj4
    pj3 <- pmargin[, 1] - pj4
    cbind("00" = 1-pj4-pj2-pj3, "01" = pj2, "10" = pj3, "11" = pj4)
  }, list( .tol = tol,
           .lmu12 = lmu12, .lphi12 = lphi12, .loratio = loratio,
           .emu12 = emu12, .ephi12 = ephi12, .eoratio = eoratio ))),
  last = eval(substitute(expression({
    misc$link <-    c("mu12"   = .lmu12 ,
                      "phi12"  = .lphi12 ,
                      "oratio" = .loratio )

    misc$earg <- list("mu12"   = .emu12 ,
                      "phi12"  = .ephi12 ,
                      "oratio" = .eoratio )

    misc$tol <- .tol
    misc$expected <- TRUE
    misc$addRidge <- .addRidge
  }),
  list( .tol = tol, .addRidge = addRidge,
        .lmu12 = lmu12, .lphi12 = lphi12,
        .loratio = loratio,
        .emu12 = emu12, .ephi12 = ephi12,
        .eoratio = eoratio ))),
  loglikelihood =
      function(mu, y, w, residuals = FALSE, eta,
               extra = NULL, summation = TRUE) {
    if (residuals) {
      stop("loglikelihood resids not implemented")
    } 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("zipebcom"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    A1vec  <- eta2theta(eta[, 1], .lmu12 ,  earg = .emu12 )
    smallno <- .Machine$double.eps^(2/4)
    A1vec[A1vec > 1.0 - smallno] <- 1.0 - smallno

    phivec <- eta2theta(eta[, 2], .lphi12 , earg = .ephi12 )
    oratio <- eta2theta(eta[, 3], .loratio , earg = .eoratio )
    okay1 <- all(is.finite(A1vec )) && all(0 < A1vec  & A1vec  < 1) &&
             all(is.finite(phivec)) && all(0 < phivec & phivec < 1) &&
             all(is.finite(oratio)) && all(0 < oratio)
    okay1
  },
  list( .lmu12 = lmu12, .lphi12 = lphi12,
        .loratio = loratio,
        .emu12 = emu12, .ephi12 = ephi12,
        .eoratio = eoratio ))),
  deriv = eval(substitute(expression({
    A1vec  <- eta2theta(eta[, 1], .lmu12 ,  earg = .emu12 )
    smallno <- .Machine$double.eps^(2/4)
    A1vec[A1vec > 1.0 -smallno] <- 1.0 - smallno

    phivec <- eta2theta(eta[, 2], .lphi12, earg = .ephi12)
    pmargin <- matrix((1 - phivec) * A1vec, n, 2)
    oratio <- eta2theta(eta[, 3], .loratio, earg = .eoratio)

    Vab <- 1 / (1/mu[, 1] + 1/mu[, 2] + 1/mu[, 3] + 1/mu[, 4])
    Vabc <- 1/mu[, 1] + 1/mu[, 2]
    denom3 <- 2 * oratio * mu[, 2] + mu[, 1] + mu[, 4]
    temp1 <- oratio * mu[, 2] + mu[, 4]
    dp11star.dp1unstar <- 2*(1-phivec)*Vab * Vabc
    dp11star.dphi1 <- -2 * A1vec * Vab * Vabc
    dp11star.doratio <- Vab / oratio
    yandmu <- (y[, 1]/mu[, 1] - y[, 2]/mu[, 2] - y[, 3]/mu[, 3] +
               y[, 4]/mu[, 4])
    dp11.doratio <- Vab / oratio
    check.dl.doratio <- yandmu * dp11.doratio

    cyandmu <- (y[, 2]+y[, 3])/mu[, 2] - 2 * y[, 1]/mu[, 1]
    dl.dmu1 <- dp11star.dp1unstar * yandmu + (1-phivec) * cyandmu
    dl.dphi1 <- dp11star.dphi1 * yandmu - A1vec * cyandmu
    dl.doratio <- check.dl.doratio
    dthetas.detas =
      cbind(dtheta.deta(A1vec,  .lmu12 ,   earg = .emu12 ),
            dtheta.deta(phivec, .lphi12,  earg = .ephi12),
            dtheta.deta(oratio, .loratio, earg = .eoratio))
    c(w) * cbind(dl.dmu1,
                 dl.dphi1,
                 dl.doratio) * dthetas.detas
  }), list( .lmu12 = lmu12, .lphi12 = lphi12, .loratio = loratio,
            .emu12 = emu12, .ephi12 = ephi12, .eoratio = eoratio ))),
  weight = eval(substitute(expression({
    wz <- matrix(0, n, 4)
    alternwz11 <- 2 * (1-phivec)^2 *
                  (2/mu[, 1] + 1/mu[, 2] - 2*Vab*Vabc^2) *
                  (dthetas.detas[, 1])^2
    wz[, iam(1, 1, M)] <- alternwz11

    alternwz22 <- 2* A1vec^2 *(2/mu[, 1] + 1/mu[, 2] - 2*Vab*Vabc^2) *
                  (dthetas.detas[, 2])^2
    wz[, iam(2, 2, M)] <- alternwz22

    alternwz12 <- -2*A1vec*(1-phivec)*
                  (2/mu[, 1] + 1/mu[, 2] - 2*Vab*Vabc^2) *
                   dthetas.detas[, 1] * dthetas.detas[, 2]
    wz[, iam(1, 2, M)] <- alternwz12

    alternwz33 <- (Vab / oratio^2) * dthetas.detas[, 3]^2
    wz[, iam(3, 3, M)] <- alternwz33

    wz[, 1:2] <- wz[, 1:2] * (1 + .addRidge)
    c(w) * wz
  }), list( .addRidge = addRidge ))))
}  # zipebcom






 binom2.Rho <-
  function(rho = 0, imu1 = NULL, imu2 = NULL,
           exchangeable = FALSE, nsimEIM = NULL) {
  lmu12 <- "probitlink"
  emu12 <- list()

  if (is.Numeric(nsimEIM)) {
    if (!is.Numeric(nsimEIM, length.arg = 1,
                    integer.valued = TRUE))
      stop("bad input for argument 'nsimEIM'")
    if (nsimEIM <= 100)
      warning("'nsimEIM' should be an integer greater than 100")
  }
  if (min(rho) <= -1 || 1 <= max(rho))
    stop("argument 'rho' should lie in (-1, 1)")


  new("vglmff",
  blurb = c("Bivariate probit model with rho = ", format(rho), "\n",
            "Links:    ",
            namesof("mu1", lmu12, earg = emu12), ", ",
            namesof("mu2", lmu12, earg = emu12)),
  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(c(1, 1), 2, 1), x = x,
                           bool = .exchangeable ,
                           constraints = constraints,
                           apply.int = TRUE)
  }), list( .exchangeable = exchangeable ))),
  deviance = Deviance.categorical.data.vgam,
  infos = eval(substitute(function(...) {
    list(M1 = 3,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("mu1", "mu2"),
         lmu12   = .lmu12 )
  }, list( .lmu12 = lmu12 ))),


  initialize = eval(substitute(expression({
    eval(process.binomial2.data.VGAM)
    predictors.names <- c(
                  namesof("mu1", .lmu12 , earg = .emu12 , short = TRUE),
                  namesof("mu2", .lmu12 , earg = .emu12 , short = TRUE))

    if (is.null( .nsimEIM )) {
      save.weights <- control$save.weights <- FALSE
    }
    if (is.null(etastart)) {
      mu1.init <- if (is.Numeric( .imu1 ))
                    rep_len( .imu1 , n) else
                    mu[, 3] + mu[, 4]
      mu2.init <- if (is.Numeric( .imu2 ))
                    rep_len( .imu2 , n) else
                    mu[, 2] + mu[, 4]
      etastart <-
        cbind(theta2eta(mu1.init, .lmu12 , earg = .emu12 ),
              theta2eta(mu2.init, .lmu12 , earg = .emu12 ))
      mustart <- NULL
    }
  }), list( .lmu12 = lmu12, .emu12 = emu12, .nsimEIM = nsimEIM,
            .imu1 = imu1, .imu2 = imu2 ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    pmargin <- cbind(eta2theta(eta[, 1], .lmu12 , earg = .emu12 ),
                     eta2theta(eta[, 2], .lmu12 , earg = .emu12 ))
    rhovec <- rep_len( .rho , nrow(eta))
    p11 <- pbinorm(eta[, 1], eta[, 2], cov12 = rhovec)
    p01 <- pmin(pmargin[, 2] - p11, pmargin[, 2])
    p10 <- pmin(pmargin[, 1] - p11, pmargin[, 1])
    p00 <- 1 - p01 - p10 - p11
    ansmat <- abs(cbind("00"=p00, "01"=p01, "10"=p10, "11"=p11))
    ansmat / rowSums(ansmat)
  }, list( .lmu12 = lmu12,
           .emu12 = emu12, .rho = rho ))),
  last = eval(substitute(expression({
      misc$link <-    c(mu1 = .lmu12 , mu2 = .lmu12 )

      misc$earg <- list(mu1 = .emu12 , mu2 = .emu12 )

      misc$nsimEIM <- .nsimEIM
      misc$expected <- TRUE
      misc$rho <- .rho
  }),
  list( .lmu12 = lmu12,
        .emu12 = emu12,
        .rho = rho, .nsimEIM = nsimEIM ))),
  loglikelihood = eval(substitute(
      function(mu, y, w, residuals = FALSE, eta,
               extra = NULL, summation = TRUE) {
    if (residuals) {
      stop("loglikelihood resids not implemented")
    } else {
      ycounts <- if (is.numeric(extra$orig.w))
                  y * w / extra$orig.w else
                  y * w # 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
      }
    }
  }, list( .rho = rho ))),
  vfamily = c("binom2.Rho", "binom2"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    pmargin <- cbind(eta2theta(eta[, 1], .lmu12 , earg = .emu12 ),
                     eta2theta(eta[, 2], .lmu12 , earg = .emu12 ))
    okay1 <- all(is.finite(pmargin)) && all(0 < pmargin & pmargin < 1)
    okay1
  },
  list( .lmu12 = lmu12, .emu12 = emu12,
        .rho = rho ))),
  deriv = eval(substitute(expression({
    pmargin <- cbind(eta2theta(eta[, 1], .lmu12 , earg = .emu12 ),
                     eta2theta(eta[, 2], .lmu12 , earg = .emu12 ))
    rhovec <- rep_len( .rho , nrow(eta))
    p11 <- pbinorm(eta[, 1], eta[, 2], cov12 = rhovec)
    p01 <- pmargin[, 2]-p11
    p10 <- pmargin[, 1]-p11
    p00 <- 1-p01-p10-p11

    ABmat <- (eta[, 1:2] -
              rhovec * eta[, 2:1]) /  sqrt(pmax(1e5 * .Machine$double.eps,
                                                1.0 - rhovec^2))
    PhiA <- pnorm(ABmat[, 1])
    PhiB <- pnorm(ABmat[, 2])
    onemPhiA <- pnorm(ABmat[, 1], lower.tail = FALSE)
    onemPhiB <- pnorm(ABmat[, 2], lower.tail = FALSE)

    smallno <- 1000 * .Machine$double.eps
    p00[p00 < smallno] <- smallno
    p01[p01 < smallno] <- smallno
    p10[p10 < smallno] <- smallno
    p11[p11 < smallno] <- smallno

    dprob00 <- dibinorm(eta[, 1], eta[, 2], cov12 = rhovec)
    dl.dprob1 <- PhiB*(y[, 4]/p11-y[, 2]/p01) +
                onemPhiB*(y[, 3]/p10-y[, 1]/p00)
    dl.dprob2 <- PhiA*(y[, 4]/p11-y[, 3]/p10) +
                onemPhiA*(y[, 2]/p01-y[, 1]/p00)
    dprob1.deta <- dtheta.deta(pmargin[, 1], .lmu12 , earg = .emu12 )
    dprob2.deta <- dtheta.deta(pmargin[, 2], .lmu12 , earg = .emu12 )
    dthetas.detas <- cbind(dprob1.deta, dprob2.deta)

    c(w) * cbind(dl.dprob1, dl.dprob2) * dthetas.detas
  }), list( .lmu12 = lmu12, .emu12 = emu12, .rho = rho ))),
  weight = eval(substitute(expression({
    if (is.null( .nsimEIM)) {
      ned2l.dprob1prob1 <- PhiB^2 *(1/p11+1/p01) +
                           onemPhiB^2 *(1/p10+1/p00)
      ned2l.dprob2prob2 <- PhiA^2 *(1/p11+1/p10) +
                           onemPhiA^2 *(1/p01+1/p00)
      ned2l.dprob1prob2 <- PhiA * (PhiB/p11 - onemPhiB/p10) +
                       onemPhiA * (onemPhiB/p00 - PhiB/p01)
      wz <- matrix(0, n, dimm(M))  # 6=dimm(M)
      wz[, iam(1, 1, M)] <- ned2l.dprob1prob1 * dprob1.deta^2
      wz[, iam(2, 2, M)] <- ned2l.dprob2prob2 * dprob2.deta^2
      wz[, iam(1, 2, M)] <- ned2l.dprob1prob2 * dprob1.deta * dprob2.deta
    } else {
      run.varcov <- 0
      ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)
      for (ii in 1:( .nsimEIM )) {
        ysim <- rbinom2.rho(n = n, mu1 = pmargin[, 1],
                                   mu2 = pmargin[, 2],
                         twoCols = FALSE, rho = rhovec)
        dl.dprob1 <- PhiB * (ysim[, 4]/p11-ysim[, 2]/p01) +
                             onemPhiB * (ysim[, 3]/p10-ysim[, 1]/p00)
        dl.dprob2 <- PhiA * (ysim[, 4]/p11-ysim[, 3]/p10) +
                             onemPhiA * (ysim[, 2]/p01-ysim[, 1]/p00)

        rm(ysim)
        temp3 <- cbind(dl.dprob1, dl.dprob2)
        run.varcov <- ((ii-1) * run.varcov +
                      temp3[, ind1$row.index] *
                      temp3[, ind1$col.index]) / ii
      }
      wz <- if (intercept.only)
        matrix(colMeans(run.varcov),
                 n, ncol(run.varcov), byrow = TRUE) else run.varcov

      wz <- wz * dthetas.detas[, ind1$row] *
                 dthetas.detas[, ind1$col]
    }
    c(w) * wz
  }), list( .nsimEIM = nsimEIM ))))
}  # binom2.Rho





 binom2.rho.ss <-
  function(lrho = "rhobitlink",
           lmu = "probitlink",  # added 20120817
           imu1 = NULL, imu2 = NULL, irho = NULL,
           imethod = 1,
           zero = 3,
           exchangeable = FALSE,
           grho = seq(-0.95, 0.95, by = 0.05)) {



  lrho <- as.list(substitute(lrho))
  e.rho <- link2list(lrho)
  l.rho <- attr(e.rho, "function.name")

  lmu <- as.list(substitute(lmu))
  emu <- link2list(lmu)
  lmu <- attr(emu, "function.name")

  if (lmu != "probitlink")
    warning("argument 'lmu' should be 'probitlink'. Changing it.")

    lmu12 <- "probitlink"  # But emu may contain some arguments.
    emu12 <- emu  # list()


  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE, positive = TRUE) ||
      imethod > 2)
    stop("argument 'imethod' must be 1 or 2")


  new("vglmff",
  blurb = c("Bivariate probit model with sample selection\n",
            "Links:    ",
            namesof("mu1", lmu12, earg = emu12), ", ",
            namesof("mu2", lmu12, earg = emu12), ", ",
            namesof("rho", l.rho, earg = e.rho)),
  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(c(1, 1, 0, 0, 0, 1), 3, 2),
                           x = x,
                           bool = .exchangeable ,
                           constraints = constraints,
                           apply.int = TRUE)
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 3,
                     predictors.names = predictors.names)
  }), list( .exchangeable = exchangeable, .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 3,
         multipleResponses = FALSE,
         parameters.names = c("mu1", "mu2", "rho"),
         zero = .zero )
  }, list( .zero = zero ))),


  initialize = eval(substitute(expression({

    if (!is.matrix(y))
      stop("response must be a 2- or 3-column matrix")
    ncoly <- ncol(y)

    temp5 <-
    w.y.check(w = w, y = y,
              ncol.w.min = 1,
              ncol.w.max = 1,
              ncol.y.min = 2,
              ncol.y.max = 3,
              Is.integer.y = TRUE,
              Is.nonnegative.y = TRUE,
              out.wy = TRUE,
              colsyperw = ncoly,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y
    if (!all(c(y) == 0 | c(y) == 1))
      stop("response matrix must have values 0 and 1 only")


    if (ncoly == 2) {
      extra$ymat2col <- y
      y <- cbind("0"  = 1 - y[, 1],
                 "10" = y[, 1] * (1 - y[, 2]),
                 "11" = y[, 1] *      y[, 2])
    } else {
      if (!all(rowSums(y) == 1))
        stop("response matrix must have two 0s and one 1 in each row")
      y1vec <- 1 - y[, 1]  # Not a 0 means a 1.
      y2vec <- ifelse(y1vec == 1, y[, 3], 0)
      extra$ymat2col <- cbind(y1vec, y2vec)
    }


    predictors.names <- c(
    namesof("mu1", .lmu12 , .emu12 , short = TRUE),
    namesof("mu2", .lmu12 , .emu12 , short = TRUE),
    namesof("rho", .l.rho , .e.rho , short = TRUE))


    ycounts <- y
    nvec <- 1



    if (!length(etastart)) {
      if (length(mustart)) {
        mu1.init <- mustart[, 1]
        mu2.init <- mustart[, 2]
      } else if ( .imethod == 1) {
        mu1.init <- weighted.mean(extra$ymat2col[, 1], c(w))
        index1 <- (extra$ymat2col[, 1] == 1)
        mu2.init <- weighted.mean(extra$ymat2col[index1, 2],
                                  w[index1, 1])
        mu1.init <- rep_len(mu1.init, n)
        mu2.init <- rep_len(mu2.init, n)

      } else if ( .imethod == 2) {
 warning("not working yet2")
          glm1.fit <- glm(ycounts ~ x - 1,
                          weights = c(w),
                          fam = binomial("probit"))
          glm2.fit <- glm(ycounts[, 2:1] ~ x - 1,
                          weights = c(w),
                          fam = binomial("probit"))
          mu1.init <- fitted(glm1.fit)
          mu2.init <- fitted(glm2.fit)
      } else {
        stop("bad value for argument 'imethod'")
      }

      if (length( .imu1 ))
        mu1.init <- rep_len( .imu1 , n)
      if (length( .imu2 ))
        mu2.init <- rep_len( .imu2 , n)



      binom2.rho.ss.Loglikfun <-
          function(rhoval, y, x, w, extraargs) {
          init.mu1 <-    extraargs$initmu1
          init.mu2 <-    extraargs$initmu2
          ymat2col <-    extraargs$ymat2col
          nvec     <-    extraargs$nvec
          eta1 <- qnorm(init.mu1)
          eta2 <- qnorm(init.mu2)

          smallno <- 1000 * .Machine$double.eps
          p11 <- pmax(smallno, pbinorm(eta1, eta2, cov12 = rhoval))
          p10 <- pmax(smallno, pnorm( eta1) - p11)
          p0  <- pmax(smallno, pnorm(-eta1))

          mumat <- abs(cbind("0"  = p0,
                             "10" = p10,
                             "11" = p11))  # rows sum to unity

          smallpos <- 1.0e-100
          mumat[mumat < smallpos] <- smallpos
          ycounts <- y  # n x 3
          use.mu <- mumat  # cbind(p0, p10, p11)

          retval <-
          sum(c(w) *
              dmultinomial(x = ycounts,
                           size = nvec, prob = use.mu,  # mumat,
                           log = TRUE, dochecking = FALSE))
          retval
        }
        rho.grid <- .grho  # seq(-0.95, 0.95, len = 31)
      try.this <-
        grid.search(rho.grid, objfun = binom2.rho.ss.Loglikfun,
                    y = y, x = x, w = w, extraargs = list(
                    ymat2col = extra$ymat2col,
                    initmu1  = mu1.init,
                    initmu2  = mu2.init,
                    nvec = nvec ))

      rho.init <- if (is.Numeric( .irho )) rep_len( .irho , n) else {
          try.this
      }



      etastart <- cbind(theta2eta(mu1.init, .lmu12 , earg = .emu12 ),
                        theta2eta(mu2.init, .lmu12 , earg = .emu12 ),
                        theta2eta(rho.init, .l.rho , earg = .e.rho ))
    }
    mustart <- NULL
  }), list( .lmu12 = lmu12, .l.rho = l.rho,
            .emu12 = emu12, .e.rho = e.rho,
                            .grho = grho,
                            .irho = irho,
            .imethod = imethod,
            .imu1 = imu1, .imu2 = imu2 ))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
    rhovec <- eta2theta(eta[, 3], .l.rho , earg = .e.rho )

    smallno <- 1000 * .Machine$double.eps
    p11 <- pmax(smallno, pbinorm(eta[, 1], eta[, 2], cov12 = rhovec))
    p10 <- pmax(smallno, pnorm( eta[, 1]) - p11)
    p0  <- pmax(smallno, pnorm(-eta[, 1]))
    sumprob <- p11 + p10 + p0
    p11 <- p11 / sumprob
    p10 <- p10 / sumprob
    p0  <- p0  / sumprob

    ansmat <- abs(cbind("0"  = p0,  # p0 == P(Y_1 = 0)
                        "10" = p10,
                        "11" = p11))
    ansmat
  }, list( .lmu12 = lmu12, .l.rho = l.rho,
           .emu12 = emu12, .e.rho = e.rho ))),
  last = eval(substitute(expression({
    misc$link <-    c(mu1 = .lmu12 , mu2 = .lmu12 , rho = .l.rho )

    misc$earg <- list(mu1 = .emu12 , mu2 = .emu12 , rho = .e.rho )

    misc$expected <- TRUE
    misc$multipleResponses <- FALSE
  }),
  list( .lmu12 = lmu12, .l.rho = l.rho,
        .emu12 = emu12, .e.rho = e.rho ))),

  loglikelihood = eval(substitute(
      function(mu, y, w, residuals = FALSE, eta,
               extra = NULL, summation = TRUE) {
    if (residuals) {
      stop("loglikelihood resids not implemented")
    } else {

      ycounts <- y  # n x 3
      nvec <- 1

      smallno <- 1000 * .Machine$double.eps
      rhovec <- eta2theta(eta[, 3], .l.rho , earg = .e.rho )
      p11 <- pmax(smallno, pbinorm(eta[, 1], eta[, 2], cov12 = rhovec))
      p10 <- pmax(smallno, pnorm( eta[, 1]) - p11)
      p0  <- pmax(smallno, pnorm(-eta[, 1]))
      sumprob <- p11 + p10 + p0
      p11 <- p11 / sumprob
      p10 <- p10 / sumprob
      p0  <- p0  / sumprob


      ll.elts <-
        c(w) * dmultinomial(x = ycounts, size = nvec,
                            prob = mu,  # use.mu,
                            log = TRUE, dochecking = FALSE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .l.rho = l.rho, .e.rho = e.rho ))),
  vfamily = c("binom2.rho.ss", "binom2"),
 validparams = eval(substitute(
     function(eta, y, extra = NULL) {
    pmargin <- cbind(eta2theta(eta[, 1], .lmu12 , earg = .emu12 ),
                     eta2theta(eta[, 2], .lmu12 , earg = .emu12 ))
    rhovec <-        eta2theta(eta[, 3], .l.rho , earg = .e.rho )
    okay1 <- all(is.finite(pmargin)) &&
             all( 0 < pmargin & pmargin < 1) &&
             all(is.finite(rhovec )) &&
             all(-1 < rhovec  & rhovec  < 1)
    okay1
  },
  list( .lmu12 = lmu12, .l.rho = l.rho,
        .emu12 = emu12, .e.rho = e.rho ))),

  deriv = eval(substitute(expression({
    nvec <- 1
    ycounts <- extra$ymat2col

    pmargin <-
      cbind(eta2theta(eta[, 1], .lmu12 , .emu12 ),
            eta2theta(eta[, 2], .lmu12 , .emu12 ))
    rhovec <- eta2theta(eta[, 3], .l.rho , .e.rho )


    smallno <- 1000 * .Machine$double.eps
    p11 <- pmax(smallno, pbinorm(eta[, 1], eta[, 2], cov12 = rhovec))
    p10 <- pmax(smallno, pnorm( eta[, 1]) - p11)
    p0  <- pmax(smallno, pnorm(-eta[, 1]))
    sumprob <- p11 + p10 + p0
    p11 <- p11 / sumprob
    p10 <- p10 / sumprob
    p0  <- p0  / sumprob


    BAmat <- (eta[, 1:2] - rhovec * eta[, 2:1]) / (
      sqrt(pmax(1e5 * .Machine$double.eps,
                1.0 - rhovec^2)))


    PhiA     <- pnorm(BAmat[, 2])
    PhiB     <- pnorm(BAmat[, 1])
    onemPhiA <- pnorm(BAmat[, 2], lower.tail = FALSE)
    onemPhiB <- pnorm(BAmat[, 1], lower.tail = FALSE)


  mycode <- FALSE  # zz
  mycode <- TRUE   # zz

 if (mycode) {
   dprob00 <- dibinorm(eta[, 1], eta[, 2], cov12 = rhovec)
   dl.dprob1 <-     PhiA *  ycounts[, 1] *      ycounts[, 2]  / p11 +
                onemPhiA *  ycounts[, 1] * (1 - ycounts[, 2]) / p10 -
                       (1 - ycounts[, 1]) / p0
   dl.dprob2 <-     PhiB * (ycounts[, 1] *      ycounts[, 2]  / p11 -
                            ycounts[, 1] * (1 - ycounts[, 2]) / p10)
   dl.drho   <-  dprob00 * (ycounts[, 1] *      ycounts[, 2]  / p11 -
                            ycounts[, 1] * (1 - ycounts[, 2]) / p10)

    dprob1.deta <- dtheta.deta(pmargin[, 1], .lmu12 , earg = .emu12 )
    dprob2.deta <- dtheta.deta(pmargin[, 2], .lmu12 , earg = .emu12 )
    drho...deta <- dtheta.deta(rhovec,       .l.rho , earg = .e.rho )

    ans.deriv <- c(w) * cbind(dl.dprob1 * dprob1.deta,
                              dl.dprob2 * dprob2.deta,
                              dl.drho   * drho...deta)
 }  # else {
    eta1 <- eta[, 1]  # dat1 %*% params[1:X1.d2]
    eta2 <- eta[, 2]  # dat2 %*% params[(X1.d2 + 1):(X1.d2 + X2.d2)]
    corr.st <- eta[, 3]  # params[(X1.d2 + X2.d2 + 1)]
    corr <- rhovec # tanh(corr.st)

    dat <- ycounts

    y1.y2  <-      dat[, 1] *      dat[, 2]
    y1.cy2 <-      dat[, 1] * (1 - dat[, 2])
    cy1    <- (1 - dat[, 1])

    d.r <- 1/sqrt(pmax(10000 * .Machine$double.eps, 1 - corr^2))
    A <- pnorm((eta2 - corr * eta1) * d.r)
    A.c <- 1 - A
    B <- pnorm((eta1 - corr * eta2) * d.r)
    p11 <- pmax(pbinorm(eta1, eta2, cov12 = corr),
                1000 * .Machine$double.eps)
    p10 <- pmax(pnorm( eta1) - p11, 1000 * .Machine$double.eps)
    p0  <- pmax(pnorm(-eta1), 1000 * .Machine$double.eps)
    d.n1 <- dnorm(eta1)
    d.n2 <- dnorm(eta2)
    d.n1n2 <- dibinorm(eta1, eta2, cov12 = corr)
    drh.drh.st <- 4 * exp(2 * corr.st)/(exp(2 * corr.st) + 1)^2

    dl.dbe1 <- d.n1 * (y1.y2/p11 * A + y1.cy2/p10 * A.c - cy1/p0)
    dl.dbe2 <- d.n2 * B * (y1.y2/p11 - y1.cy2/p10)
    dl.drho <- d.n1n2 * (y1.y2/p11 - y1.cy2/p10) * drh.drh.st

    ans.deriv2 <- c(w) * cbind(dl.dbe1, dl.dbe2, dl.drho)
 # }










    ans.deriv
  }), list( .lmu12 = lmu12, .l.rho = l.rho,
            .emu12 = emu12, .e.rho = e.rho ))),

  weight = eval(substitute(expression({



 if (mycode) {
    ned2l.dprob1prob1 <-      PhiA^2 / p11 +
                          onemPhiA^2 / p10 +
                                   1 / p0
    ned2l.dprob2prob2 <-   (1/p11 + 1/p10) * PhiB^2
    ned2l.drho2       <-   (1/p11 + 1/p10) * dprob00^2

    ned2l.dprob1prob2 <-    PhiA * PhiB / p11 - onemPhiA * PhiB / p10
    ned2l.dprob1rho   <-   (PhiA/p11 -  onemPhiA/p10) * dprob00
    ned2l.dprob2rho   <-   (1/p11 + 1/p10) * PhiB * dprob00

    wz <- matrix(0, n, dimm(M))  # 6=dimm(M)
    wz[, iam(1, 1, M)] <- ned2l.dprob1prob1 * dprob1.deta^2
    wz[, iam(2, 2, M)] <- ned2l.dprob2prob2 * dprob2.deta^2
    wz[, iam(3, 3, M)] <- ned2l.drho2       * drho...deta^2
    wz[, iam(1, 2, M)] <- ned2l.dprob1prob2 * dprob1.deta * dprob2.deta
    wz[, iam(1, 3, M)] <- ned2l.dprob1rho   * dprob1.deta * drho...deta
    wz[, iam(2, 3, M)] <- ned2l.dprob2rho   * dprob2.deta * drho...deta
  }  # else {

    ned2l.be1.be1 <- (A^2/p11 + A.c^2/p10 + 1/p0)      * d.n1^2
    ned2l.be2.be2 <- (  1/p11 +     1/p10) * B^2       * d.n2^2
    ned2l.rho.rho <- (  1/p11 +     1/p10) * d.n1n2^2  * drh.drh.st^2

    ned2l.be1.be2 <- (A *  B/p11  - A.c *  B/p10)  * d.n1   * d.n2
ned2l.be1.rho <-(A *(1/p11) - A.c * (1/p10)) * d.n1n2 * d.n1 * drh.drh.st
ned2l.be2.rho <- B *(1/p11  +        1/p10)  * d.n1n2 * d.n2 * drh.drh.st




    WZ <- matrix(0, n, dimm(M))  # 6=dimm(M)
    WZ[, iam(1, 1, M)] <- ned2l.be1.be1
    WZ[, iam(2, 2, M)] <- ned2l.be2.be2
    WZ[, iam(3, 3, M)] <- ned2l.rho.rho
    WZ[, iam(1, 2, M)] <- ned2l.be1.be2
    WZ[, iam(1, 3, M)] <- ned2l.be1.rho
    WZ[, iam(2, 3, M)] <- ned2l.be2.rho

    c(w) * wz
  }), list( .zero = zero ))))
}  # binom2.rho.ss







 extbetabinomial <-
  function(lmu = "logitlink",
           lrho = "cloglink",  # (-eps, 1)
           zero = "rho",
           irho = 0,  # NULL, 0.01,
           grho = c(0, 0.05, 0.1, 0.2),
           vfl = FALSE, Form2 = NULL,
           imethod = 1, ishrinkage = 0.95
           ) {
  lmu <- as.list(substitute(lmu))
  emu <- link2list(lmu)
  lmu <- attr(emu, "function.name")

  lrho <- as.list(substitute(lrho))
  erho <- link2list(lrho)
  lrho <- attr(erho, "function.name")

  if (!is.logical(vfl) || length(vfl) != 1)
    stop("argument 'vfl' must be TRUE or FALSE")

  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE,
                  positive = TRUE) ||
     imethod > 4)
    stop("argument 'imethod' must be 1, 2, 3 or 4")

  if (!is.Numeric(ishrinkage, length.arg = 1) ||
      ishrinkage < 0 || ishrinkage > 1)
    stop("bad input for argument 'ishrinkage'")

  new("vglmff",
  blurb = c("Extended beta-binomial model (Prentice 1986)\n",
            "Links:      ",
            namesof("mu",  lmu,  emu), ", ",
            namesof("rho", lrho, erho), "\n",
            "Mean:     mu", "\n",
            "Variance: mu*(1-mu)*(1+(w-1)*rho)/w"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints,
              x = x, .zero , M = M, M1 = 2,
              predictors.names = predictors.names)

    if ( .vfl && M != 2)
      stop("vfl = TRUE only allowed when M == 2")
    LC <- length(constraints)
    if ( .vfl && LC <= 2)
      stop("vfl = T only allowed if ncol(x) > 2")
    if ( .vfl && !is.zero( .zero ))
      stop("Need zero = NULL when vfl = TRUE")
    if ( .vfl ) {
      constraints <- cm.VGAM(rbind(0, 1), x = x,
                       bool = .Form2 ,
                       constraints = constraints)
      mterms <- 0
      for (jay in 1:LC) {  # Include the intercept
        if (!all(c(constraints[[jay]]) == 0:1)) {
          mterms <- mterms + 1
          constraints[[jay]] <- rbind(1, 0)
        }
      }  # jay
      if (mterms == 0)
        warning("no terms for 'mu'... something",
                " looks awry")
      if (mterms == LC)
        warning("no terms for 'rho'... something",
                " looks awry")
    }  # vfl
  }), list( .zero = zero,
            .vfl = vfl ,
            .Form2 = Form2))),


  infos = eval(substitute(function(...) {
    list(M1 = 2,
         dpqrfun = "extbetabinom",
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("mu", "rho"),
         imethod  = .imethod ,
         ishrinkage  = .ishrinkage ,
         vfl = .vfl , Form2 = .Form2 ,
         lmu  = .lmu ,
         lrho = .lrho ,
         vfl  = .vfl ,
         zero = .zero )
  },
  list( .lmu = lmu, .lrho = lrho,
        .imethod = imethod, .zero = zero,
        .vfl = vfl, .Form2 = Form2,
        .ishrinkage = ishrinkage ))),


  initialize = eval(substitute(expression({
    if (!all(w == 1))
      extra$orig.w <- w
    
    mustart.orig <- mustart
    eval(binomialff()@initialize)  # Note
    if (length(mustart.orig)) # Retain if 
      mustart <- mustart.orig  # inputted

    ycounts <- if (is.numeric(extra$orig.w))
         y * w / extra$orig.w else
         y * w  # Convert proportions to counts
    if (max(abs(ycounts - round(ycounts))) > 1e-6)
      warning("the response (as counts) does not ",
              "appear to be integer-valued. ",
              "Am rounding to integer values.")
    ycounts <- round(ycounts)  # Now an integer
    nvec <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else
              round(w)
    extra$size <- nvec  # For @validparams
    predictors.names <-
    c(namesof("mu",  .lmu ,  .emu  , tag = FALSE),
      namesof("rho", .lrho , .erho , tag = FALSE))

    if (!length(etastart)) {
      ebbinom.Loglikfun <-
        function(rhoval, y, x, w, extraargs) {
        muvec <- rep(extraargs$mustart, length(w))
        rhovec <- rep(rhoval, length(w))
        ycounts <- extraargs$ycounts  # integer
        nvec <- extraargs$nvec  
        if (extraargs$trace) {
          cat(".")
          flush.console()
        }
        ans1 <- sum(dextbetabinom(ycounts,
                nvec, muvec, rhovec, log = TRUE))
        ans1
      }  # ebbinom.Loglikfun
      rho.grid <- c( .grho )
      mustart.use <- if (length(mustart.orig)) {
        mustart.orig
      } else if ( .imethod == 1) {
        rep_len(weighted.mean(y, w), n)
      } else if ( .imethod == 2) {
        .ishrinkage * weighted.mean(y, w) +
        (1 - .ishrinkage ) * y
      } else if ( .imethod == 3) {
        ymat <- cbind(y)
        mat.temp <- matrix(colMeans(ymat),
                           nrow(ymat), ncol(ymat),
                           byrow = TRUE)
        0.5 * mustart + 0.5 * mat.temp
      } else {  # imethod == 4
        mustart
      }

      if (!is.Numeric( .irho )) {
        if (trace) {
          cat("Starting grid search")
          flush.console()
        }
        try.this <-
        grid.search(rho.grid,
            objfun = ebbinom.Loglikfun,
            y = y,  x = x, w = w,
            extraargs = list(
              trace = trace,
              ycounts = ycounts,
              nvec = if (is.numeric(extra$orig.w))
                round(w / extra$orig.w) else
                round(w),
              mustart = mustart.use))
        if (trace) {
          cat("\n")
          flush.console()
        }
      }  # irho is not Numerical
      init.rho <- if (is.Numeric( .irho ))
                  rep_len( .irho , n) else
                  rep_len(try.this, n)
      etastart <-
    cbind(theta2eta(mustart.use, .lmu  ,  .emu ),
          theta2eta(init.rho,    .lrho , .erho ))
      mustart <- NULL  # As etastart been computed
    }
  }),
  list( .lmu = lmu, .lrho = lrho, .irho = irho,
        .emu = emu, .erho = erho, .grho = grho,
        .imethod = imethod,
        .ishrinkage = ishrinkage ))),
  linkinv = eval(substitute(
      function(eta, extra = NULL)
    eta2theta(eta[, 1], .lmu , earg = .emu ),
  list( .lmu = lmu, .emu = emu ))),
  last = eval(substitute(expression({
    misc$link <-    c(mu = .lmu , rho = .lrho)
    misc$earg <- list(mu = .emu , rho = .erho )
  }),
  list( .lmu = lmu, .lrho = lrho,
        .emu = emu, .erho = erho, .zero = zero ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL, summation = TRUE) {
      ycounts <- if (is.numeric(extra$orig.w))
               y * w / extra$orig.w else
               y * w  # Convert propns to counts

    mymu <- eta2theta(eta[, 1], .lmu ,  .emu )
    rho  <- eta2theta(eta[, 2], .lrho , .erho )

    if (max(abs(ycounts - round(ycounts))) > 1e-6)
      warning("converting 'ycounts' to integer",
              " in @loglikelihood")
    ycounts <- round(ycounts)

    nvec <- if (is.numeric(extra$orig.w))
            round(w / extra$orig.w) else round(w)

    if (residuals) {
      stop("loglikelihood resid not implemented")
    } else {
      ll.elts <- (if (is.numeric(extra$orig.w))
                 extra$orig.w else 1) *
        dextbetabinom(ycounts, nvec, mymu,
                      rho, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lmu = lmu, .lrho = lrho,
           .emu = emu, .erho = erho  ))),
  vfamily = c("extbetabinomial"),
  validparams = eval(substitute(
      function(eta, y, extra = NULL) {
    mymu <- eta2theta(eta[, 1], .lmu  , .emu  )
    rho  <- eta2theta(eta[, 2], .lrho , .erho )
    size <- extra$size
    lowlim <- pmax(-mymu / (size - mymu - 1),
       -(1 - mymu) / (size - (1 - mymu) - 1))
    okay1 <-
      all(is.finite(mymu)) &&
      all(is.finite(rho)) &&
      all(0 < mymu & mymu < 1) &&
      all(lowlim < rho & rho < 1)
    okay1
  },
  list( .lmu = lmu, .lrho = lrho,
        .emu = emu, .erho = erho  ))),




  simslot = eval(substitute(
  function(object, nsim) {

    pwts <-
      if (length(pwts <- object@prior.weights) > 0)
          pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    w <- pwts
    eta <- predict(object)
    extra <- object@extra

    mymu <- eta2theta(eta[, 1],  .lmu ,  .emu )
    rho  <- eta2theta(eta[, 2], .lrho , .erho )
    nvec <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else
              round(w)

    rextbetabinom(nsim * length(rho), nvec,
                  prob = mymu, rho  = rho)
  },
  list( .lmu = lmu, .lrho = lrho,
        .emu = emu, .erho = erho  ))),



  deriv = eval(substitute(expression({
    size <- if (is.numeric(extra$orig.w))
              round(w / extra$orig.w) else
              round(w)
    ycounts <- if (is.numeric(extra$orig.w))
            y * w / extra$orig.w else
            y * w  # Convert propns to counts
    ycounts <- round(ycounts)
    prob <-
    mymu <- eta2theta(eta[, 1], .lmu  , .emu  )
    rho  <- eta2theta(eta[, 2], .lrho , .erho )

    dmu.deta  <- dtheta.deta(mymu, .lmu  , .emu  )
    drho.deta <- dtheta.deta(rho,  .lrho , .erho )
    tmp4 <- ned2l.ebbinom(ycounts,
               size, prob, rho)
    dl.dmu  <- tmp4$deriv1[, 1]
    dl.drho <- tmp4$deriv1[, 2]

    orig.w <- c(if (is.numeric(extra$orig.w))
                extra$orig.w else 1)
    ansd <- orig.w *
    cbind(dl.dmu  * dmu.deta,
          dl.drho * drho.deta)
    ansd
  }),
  list( .lmu = lmu, .lrho = lrho,
        .emu = emu, .erho = erho  ))),
  weight = eval(substitute(expression({
    wz <- matrix(NA_real_, n, dimm(M))  # 3
    wz11 <- tmp4$Eim[, 1]
    wz22 <- tmp4$Eim[, 2]
    wz12 <- tmp4$Eim[, 3]
    wz[, iam(1, 1, M)] <- wz11 * dmu.deta^2
    wz[, iam(2, 2, M)] <- wz22 * drho.deta^2
    wz[, iam(1, 2, M)] <- wz12 * dmu.deta * 
                                 drho.deta
    wz * orig.w
  }),
  list( .lmu = lmu, .lrho = lrho,
        .emu = emu, .erho = erho ))))
}  # extbetabinomial





dextbetabinom <-
    function(x, size, prob, rho = 0, log = FALSE,
             forbycol = TRUE) {
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)
  if (!is.logical(forbycol) || length(forbycol) != 1)
    stop("bad input for argument 'forbycol'")
  if (!forbycol)
    stop("argument 'forbycol' must be TRUE")

  L <- max(length(x), length(size), length(prob),
           length(rho))
  if (length(x)    != L) x    <- rep_len(x,    L)
  if (length(size) != L) size <- rep_len(size, L)
  if (length(prob) != L) prob <- rep_len(prob, L)
  if (length(rho)  != L) rho  <- rep_len(rho,  L)


  if (all(is.finite(rho)) && all(rho == 0))
    return(dbinom(x, size, prob, log = log.arg))


  bad0 <- !is.finite(size) | size != round(size) |
          !is.finite(prob) |
          !is.finite(rho) |
          size < 3 | 1 <= rho |
          prob < 0 | 1 < prob
  bad1 <- bad0 | x != round(x)
  bad2 <- bad0 | !is.finite(x)
  bad  <- bad0 | bad1 | bad2
  logpdf <- lchoose(size, x)
  logpdf[bad] <- NA  # x + size + prob + rho
  L0s <- numeric(L)
  maxs <- max(size, na.rm = TRUE)
  Gama <- rho / (1 - rho)
  if (any(!bad)) {
    if (forbycol)  # Loop over the cols
    for (i in 0:maxs) {
      logpdf <- logpdf + 
        ifelse(!bad & i < x,
               log(prob + Gama * i), L0s) +
        ifelse(!bad & i < size - x,
               log1p(-prob + Gama * i), L0s) -
        ifelse(!bad & i < size,
               log1p(Gama * i), L0s)
    }  # for i
  }  # any(!bad)
  
  
  logpdf[!bad0 & x < 0   ] <- log(0)
  logpdf[!bad0 & size < x] <- log(0)
  vecTF <- bad & !bad2   # & x != round(x)
  if (any(vecTF)) {
    warning("non-integer x. Returning 0 values")
    logpdf[vecTF] <- log(0)
  }
  logpdf[ bad0] <- NaN
  if (log.arg) logpdf else exp(logpdf)
}  # dextbetabinom

    


pextbetabinom <-
  function(q, size, prob, rho = 0,
           lower.tail = TRUE,
           forbycol = TRUE) {
  q <- floor(q)
  L <- max(length(q), length(size), length(prob),
           length(rho))
  if (length(q)    != L) q    <- rep_len(q,    L)
  if (length(size) != L) size <- rep_len(size, L)
  if (length(prob) != L) prob <- rep_len(prob, L)
  if (length(rho)  != L) rho  <- rep_len(rho,  L)

  bad0 <- !is.finite(size) | size != round(size) |
          !is.finite(prob) |
          !is.finite(rho) |
          size < 3 | 1 <= rho |
          prob < 0 | 1 < prob
  bad1 <- bad0 | q != round(q)
  bad2 <- bad0 | !is.finite(q)
  bad  <- bad0 | bad1 | bad2

      
  if (all(is.finite(rho)) && all(rho == 0))
    return(pbinom(q, size, prob,  # log.p = log.p,
                  lower.tail = lower.tail))

  ans <- q + size + prob + rho
  ok3 <- !bad & 0 <= q
  if (any(ok3)) {
    ans[ok3] <- mapply(function(q, size, prob, rho)
      sum(dextbetabinom(0:q, size, prob, rho)),
      q      =      q[ok3],
      size   =   size[ok3],
      prob   =   prob[ok3],
      rho    =    rho[ok3])
  }
  ans[!bad0 & size < q      ] <- 1
  ans[!bad0 & q < 0         ] <- 0
  ans[ bad0] <- NaN

  if (!lower.tail)
    ans <- 1 - ans
  ans
}  # pextbetabinom



qextbetabinom <-
    function(p, size, prob, rho = 0,
             forbycol = TRUE) {

  L <- max(length(p), length(size), length(prob),
           length(rho))
  if (length(p)    != L) p    <- rep_len(p,    L)
  if (length(size) != L) size <- rep_len(size, L)
  if (length(prob) != L) prob <- rep_len(prob, L)
  if (length(rho)  != L) rho  <- rep_len(rho,  L)

  bad0 <- !is.finite(size) | size != round(size) |
          !is.finite(prob) |
          !is.finite(rho)  |
          is.na(p) | is.na(size) |
          is.na(prob) | is.na(rho) |
          size < 3 | 1 <= rho |
          prob < 0 | 1 < prob
  bad1 <- bad0 | !is.finite(p) | p <= 0 | 1 <= p
  bad  <- bad0 | bad1

      
  if (all(is.finite(rho)) && all(rho == 0))
    return(qbinom(p, size, prob))

  ans <- p + size + prob + rho

  lo <- numeric(L) - 0.25
  hi <- size + 0.25
  approx.ans <- lo  # True at lhs
  dont.iterate <- bad

  foo <- function(q, size, prob, rho, p)
    pextbetabinom(q, size, prob, rho) - p

  lhs <- dont.iterate |
         p <= dextbetabinom(0, size, prob, rho)


  if (any(!lhs)) {
  approx.ans[!lhs] <-
    bisection.basic(foo, lo[!lhs], hi[!lhs],
                    tol = 1/16,
                    size = size[!lhs],
                    prob = prob[!lhs],
                    rho  =  rho[!lhs],
                    p = p[!lhs])
  faa <- floor(approx.ans[!lhs])
  tmp <-
    ifelse(pextbetabinom(faa, rho = rho[!lhs],
             prob = prob[!lhs],
             size = size[!lhs]) < p[!lhs] &
           p[!lhs] <= pextbetabinom(faa+1,
                 prob = prob[!lhs],
                 size = size[!lhs],
                 rho  =  rho[!lhs]),
           faa+1, faa)
  ans[!lhs] <- tmp
  }  # any(!lhs)
        
  vecTF <- !bad0 & !is.na(p) &
           p <= dextbetabinom(0, size, prob, rho)
  ans[vecTF] <- 0

  ans[!bad0 & !is.na(p) & p == 0] <- 0
  vecTF <- !bad0 & !is.na(p) & p == 1
  ans[vecTF] <- size[vecTF]
  ans[!bad0 & !is.na(p) & p <  0] <- NaN
  ans[!bad0 & !is.na(p) & p >  1] <- NaN
  ans[ bad0] <- NaN
  ans
}  # qextbetabinom






rextbetabinom <-
  function(n, size, prob, rho = 0) {
     
  if (all(is.finite(rho)) && all(rho == 0))
    return(rbinom(n, size, prob))

  qextbetabinom(runif(n), size, prob, rho)
}  # rextbetabinom







eimij.ebbinom <-
  function(x,  # x is a scalar
           size, prob, rho) {  # These r vectors
  L0s <- numeric(L <- length(rho))  # not x!
  xvec <- rep_len(x, L)  # For vectorized ifelse()
  ned2l.dprob2 <- ned2l.dGama2 <-
  ned2l.dp.dG <- L0s
  maxs <- max(size, na.rm = TRUE)
  Gama <- rho / (1 - rho)
  for (i in 0:maxs) {
    ned2l.dprob2 <- ned2l.dprob2 +
      ifelse(i < xvec,
             1 / (prob + Gama * i)^2,  L0s) +
      ifelse(i < size - xvec,
             1 / (1 - prob + Gama * i)^2, L0s)
    ned2l.dGama2 <- ned2l.dGama2 +
      ifelse(i < xvec,
             (i / (prob + Gama * i))^2,  L0s) +
      ifelse(i < size - xvec,
             (i / (1 - prob + Gama * i))^2, L0s) -
      ifelse(i < size,
             (i / (1 + Gama * i))^2, L0s)
    ned2l.dp.dG <- ned2l.dp.dG +
      ifelse(i < xvec,
             i / (prob + Gama * i)^2,  L0s) -
      ifelse(i < size - xvec,
             i / (1 - prob + Gama * i)^2, L0s)
  }  # for i
  cbind(ned2l.dprob2,  # In matrix-band format
        ned2l.dGama2,
        ned2l.dp.dG)
}  # eimij.ebbinom






ned2l.ebbinom <-  # x and others are
    function(x, size, prob, rho,  # n-vectors
             forbycol = TRUE) {
  L0s <- numeric(L <- length(rho))
  ned2l.dprob2 <- ned2l.dGama2 <- 
  ned2l.dp.dG <- L0s
  dl.dprob <- dl.dGama <- L0s
  maxs <- max(size, na.rm = TRUE)
  Gama <- rho / (1 - rho)
  dGama.drho <- (1 + Gama)^2
    for (i in 0:maxs) {
      dl.dprob <- dl.dprob + 
        ifelse(i < x,
               1 / (prob + Gama * i), L0s) -
        ifelse(i < size - x,
               1 / (1 - prob + Gama * i), L0s)
      dl.dGama <- dl.dGama + 
        ifelse(i < x,
               i / (prob + Gama * i), L0s) +
        ifelse(i < size - x,
               i / (1 - prob + Gama * i), L0s) -
        ifelse(i < size,
               i / (1 + Gama * i), L0s)
      tmpmat <- eimij.ebbinom(i, size, prob, rho)
      pmfvec <- dextbetabinom(i, size, prob, rho)
      ned2l.dprob2 <- ned2l.dprob2 +
        pmfvec * tmpmat[, 1]  # Expectation
      ned2l.dGama2 <- ned2l.dGama2 +
        pmfvec * tmpmat[, 2]
      ned2l.dp.dG <- ned2l.dp.dG +
        pmfvec * tmpmat[, 3]
    }  # for i

  list(deriv1 = cbind(dl.dprob,
                      dl.dGama * dGama.drho),
       Eim = cbind(ned2l.dprob2,
                   ned2l.dGama2 * dGama.drho^2,
                   ned2l.dp.dG  * dGama.drho))
}  # ned2l.ebbinom

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.