R/summary.aodml.R

In aods3: Analysis of Overdispersed Data using S3 Methods

summary.aodml <- function(object, ...) {
  ## Vector b
  # function to insert NAs in se when there are NAs
  insNA <- function(b, se) {
    if (any(is.na(b))) {
      nb <- length(b)
      SE <- rep(NA, nb)
      j <- 1
      for (i in seq(nb)) {
        if (!is.na(b[i])) {
          SE[i] <- se[j]
          j <- j + 1
        }
      }
      se <- SE
    }
    se
  }
  
  # check whether any coef was set to a fixed value
  nb <- length(coef(object))
  param <- object$param
  vpar <- object$varparam
  
  # check whether any estimated var was negative and replaces these elements with NA's
  diagv <- diag(vpar)
  if (!all(is.na(diagv))) {
    if (any(diagv[!is.na(diagv)] < 0)) {
      diagv[diagv < 0] <- NA
      warning("At least one variance was < 0 in the var-cov matrix. Any such element was replaced with NA.\n")
    }
  }
  
  # position of coef
  pos1 <- seq(nb)
  
  # position of coef set to a fixed value, if any
  fp <- match("fixpar", table = names(object$call))
  pos2 <- NA
  if (!is.na(fp))
    pos2 <- eval(object$call$fixpar[[2]])
  
  # remove coef set to a fixed value, if any
  pos3 <- setdiff(pos1, pos2)
  BCoef <- data.frame()
  
  # compute new var-cov mat, coef vector and position of term(s) to be tested
  if (length(pos3) > 0) {
    b3 <- param[pos3]
    v3 <- if (object$singular.hessian == 0 && !all(is.na(diagv))) diagv[pos3] else diagv
    
    # coef, se, z and t test
    se3 <- sqrt(v3)
    se3 <- insNA(b3, se3)
    BCoef <- data.frame(b = b3, se = se3, z = b3 / se3,
                        P = 2 * (1 - pnorm(abs(b3) / se3)))
    nam <- names(b3)
    rownames(BCoef) <- nam
    colnames(BCoef) <- c("Estimate", "Std. Error", "z value", "Pr(> |z|)")
  }
  
  # coef set to a fixed value, if any
  pos4 <- setdiff(pos1, pos3)
  FixedBCoef <- data.frame()
  if (length(pos4) > 0)
    FixedBCoef <- data.frame(Value = param[pos4])
  
  ## Vector phi
  
  # position of scoef
  pos1 <- (nb + 1):length(param)
  
  # position of coef set to a fixed value, if any
  fp <- match("fixpar", table = names(object$call))
  pos2 <- NA
  if (!is.na(fp)) pos2 <- eval(object$call$fixpar[[2]])
  
  # remove coef set to a fixed value, if any
  pos3 <- setdiff(pos1, pos2)
  Phi <- data.frame()
  
  # compute new var-cov mat, coef vector and position of term(s) to be tested
  if (length(pos3) > 0) {
    b3 <- param[pos3]
    if (object$singular.hessian == 0 && !all(is.na(diagv))) {
      va3 <- diagv[pos3]
      va3[va3 < 0] <- NA
      se3 <- sqrt(va3)
      se3 <- insNA(b3, se3)
    } else {
      se3 <- rep(NA, length(b3))
    }
    
    Phi <- data.frame(b  = b3, se = se3)
    nam <- names(b3)
    rownames(Phi) <- nam
    colnames(Phi) <- c("Estimate", "Std. Error") #, "z value", "Pr(> z)"
  }
  
  # print coef set to a fixed value, if any
  pos4 <- setdiff(pos1, pos3)
  FixedPhi <- data.frame()
  if (length(pos4) > 0)
    FixedPhi <- data.frame(Value = param[pos4])
  
  structure(
    list(object = object, BCoef = BCoef, FixedBCoef = FixedBCoef,
         Phi = Phi, FixedPhi = FixedPhi),
    class = "summary.aodml")
}