R/quasi_func.R

In assessor: Assessment Tools for Regression Models with Discrete and Semicontinuous Outcomes

#' @keywords internal
.resid_quasi_finalize <- function(x,
                                  y,
                                  main,
                                  xlab = "s",
                                  ylab = expression(hat(U) * "(s)"),
                                  xlim = range(x, na.rm = TRUE),
                                  line_args = list(),
                                  ...) {

  defaults <- list(
    type     = "l",
    main     = main,
    ylab     = ylab,
    xlab     = xlab,
    cex.lab  = 2,
    cex.axis = 2,
    cex.main = 2,
    lwd      = 2,
    xlim     = xlim
  )

  qq_args <- modifyList(defaults, list(...))

  do.call(plot, c(list(x = x, y = y), qq_args))

  default_abline <- list(
    a   = 0,
    b   = 1,
    col = "red",
    lty = 5,
    lwd = 2
  )
  ab_args <- modifyList(default_abline, line_args)
  do.call(abline, ab_args)

  invisible(NULL)
}


#' @keywords internal
listvec <- function(x) {
  x[1]:x[2]
}

## Binary
margin01 <- function(u, y, q0, h) {
  wei <- 1 * ((q0 - u)^2 < 5 * h^2) * (1 - ((q0 - u)^2) / h^2 / 5)
  l <- sum(wei * 1 * (y == 0)) / sum(wei)
  l
}
margin01 <- Vectorize(margin01, "u")

## bandwidth selector
bandwidth01 <- function(y, q0) {
  bw <- npregbw(ydat = 1 * (y == 0), xdat = q0, ckertype = "epanechnikov")
  return(bw$bw)
}

#' @keywords internal
#' @importFrom utils modifyList
resid.bin_quasi <- function(model, line_args, ...){
  y   <- model$y
  q10 <- 1 - model$fitted.values

  h       <- bandwidth01(y = y, q0 = q10)
  x.input <- seq(0, 1, length.out = 101)
  y_curve <- margin01(x.input, y = y, q0 = q10, h = h)

  xlim_vec <- c(
    min(pbinom(0, size = 1, prob = q10)),
    max(pbinom(0, size = 1, prob = q10))
  )

  .resid_quasi_finalize(
    x         = x.input,
    y         = y_curve,
    main      = "Quasi, Binary",
    xlim      = xlim_vec,
    line_args = line_args,
    ...
  )
}


## Poisson
marginesti <- function(u, y, lambdaf, h) {
  n <- length(y)
  # F^{-1}(u)
  inv1 <- qpois(u, lambdaf)
  inv1m <- ifelse(inv1 > 0, inv1 - 1, 0)
  n1 <- cbind(inv1, inv1m)
  # H^+ and H^_
  p1 <- ppois(n1, lambdaf)
  # find out which one is closer to u
  ind1 <- apply(abs(p1 - u), 1, which.min)
  # weight function
  wei <- 1 * ((p1[cbind(1:n, ind1)] - u)^2 < 5 * h^2) *
    (1 - ((p1[cbind(1:n, ind1)] - u)^2) / h^2 / 5)
  l <- sum(wei * 1 * (y <= n1[cbind(1:n, ind1)])) / sum(wei)
  l
}
marginesti <- Vectorize(marginesti, "u")

####### bandwidth selector
bandwidthp <- function(y, lambdaf){
  n <- length(y)
  newout <- unlist(sapply(split(cbind(qpois(0.1, lambdaf), qpois(0.9, lambdaf)), 1:n), listvec))
  newlambda <- rep(lambdaf, times = qpois(0.9, lambdaf) - qpois(0.1, lambdaf) + 1)
  newx <- ppois(newout, newlambda)
  newy <- 1 * (rep(y, times = qpois(0.9, lambdaf) - qpois(0.1, lambdaf) + 1) <= newout)
  bws <- npregbw(ydat = newy[which(newx <= 0.9)], xdat = newx[which(newx <= 0.9)], ckertype = "epanechnikov")
  return(bws$bw)
}

#' @keywords internal
resid.pois_quasi <- function(model, line_args, ...){
  y        <- model$y
  lambda1f <- model$fitted.values

  h       <- bandwidthp(y = y, lambdaf = lambda1f)
  x.input <- seq(0, 1, length.out = 101)
  y_curve <- marginesti(u = x.input, y = y, lambdaf = lambda1f, h = h)

  xlim_vec <- c(min(ppois(0, lambda = lambda1f)), 1)

  .resid_quasi_finalize(
    x         = x.input,
    y         = y_curve,
    main      = "Quasi, Poisson",
    xlim      = xlim_vec,
    line_args = line_args,
    ...
  )
}

############
##   NB   ##
############
marginnb <- function(u, y, lambdaf, sizef,h) {
  n <- length(y)
  inv1 <- qnbinom(u, mu = lambdaf, size = sizef)
  inv1m <- ifelse(inv1 > 0, inv1 - 1, 0)
  n1 <- cbind(inv1, inv1m)
  p1 <- pnbinom(n1, mu = lambdaf, size = sizef)
  ind1 <- apply(abs(p1 - u), 1, which.min)
  wei <- 1 * ((p1[cbind(1:n, ind1)] - u)^2 < 5 * h^2) *
    (1 - ((p1[cbind(1:n, ind1)] - u)^2) / h^2 / 5)
  l <- sum(wei * 1 * (y <= n1[cbind(1:n, ind1)])) / sum(wei)
  l
}

marginnb <- Vectorize(marginnb, "u")

## bandwidth selector
bandwidthnb <- function(y, lambdaf, sizef) {
  n <- length(y)
  newout <- unlist(sapply(split(cbind(qnbinom(0.1, mu = lambdaf, size = sizef), qnbinom(0.9, mu = lambdaf, size = sizef)), 1:n), listvec))
  newlambda <- rep(lambdaf, times = qnbinom(0.9, mu = lambdaf, size = sizef) - qnbinom(0.1, mu = lambdaf, size = sizef) + 1)
  newx <- pnbinom(newout, mu = newlambda, size = sizef)
  newy <- 1 * (rep(y, times = qnbinom(0.9, mu = lambdaf, size = sizef) - qnbinom(0.1, mu = lambdaf, size = sizef) + 1) <= newout)

  newy1 <- newy[which(newx <= 0.9)]
  newx1 <- newx[which(newx <= 0.9)]
  bws <- npregbw(ydat = newy1, xdat = newx1, ckertype = "epanechnikov")
  return(bws$bw)
}

#' @keywords internal
resid.nb_quasi <- function(model, line_args, ...){
  y        <- model$y
  lambda1f <- model$fitted.values
  size1f   <- summary(model)$theta

  h       <- bandwidthnb(y = y, lambdaf = lambda1f, sizef = size1f)
  x.input <- seq(0, 1, length.out = 101)
  y_curve <- marginnb(u = x.input, y = y,
                      lambdaf = lambda1f, sizef = size1f, h = h)

  xlim_vec <- c(min(pnbinom(0, size = size1f, mu = lambda1f)), 1)

  .resid_quasi_finalize(
    x         = x.input,
    y         = y_curve,
    main      = "Quasi, NB",
    xlim      = xlim_vec,
    line_args = line_args,
    ...
  )
}

###########
## Ord  ###
###########
### \hat{U} y is the outcome, q0=P(y<=1), q1=P(y<=2)
marginm <- function(x, y, p1, h) { # p1 = t(apply(model$fitted.values,1 ,cumsum))[,-n]
  n <- length(y)
  ind1 <- apply(abs(p1 - x), 1, which.min)
  wei <- 1 * ((p1[cbind(1:n, ind1)] - x)^2 < 5 * h^2) *
    (1 - ((p1[cbind(1:n, ind1)] - x)^2) / h^2 / 5)
  l <- sum(wei * 1 * (y <= (ind1 - 1))) / sum(wei)
  l
}

marginm <- Vectorize(marginm, "x")

### Bandwidth selection
bandwidthord <- function(y, p1) {
  k <- max(y)
  n <- length(y)
  xdat. <- rep(NA, n*k)
  for(i in 0:(k-1)) xdat.[(i*n+1):((i+1)*n)] <- p1[, (i+1)]
  ydat. <- rep(NA,n*k)
  for(i in 0:(k-1)) ydat.[(i*n+1):((i+1)*n)] <- 1*(y <= i)
  bw <- npregbw(ydat = ydat., xdat = xdat., ckertype = "epanechnikov")
  return(bw$bw)
}

#' @keywords internal
resid.ordi_quasi <- function(model, line_args, ...){
  y  <- as.numeric(model$model[,1]) - 1
  p1 <- t(apply(model$fitted.values, 1, cumsum))

  h       <- bandwidthord(y = y, p1 = p1)
  x.input <- seq(0, 1, length.out = 101)
  y_curve <- marginm(x = x.input, y = y, p1 = p1, h = h)

  xlim_vec <- c(0, 1)

  .resid_quasi_finalize(
    x         = x.input,
    y         = y_curve,
    main      = "Quasi, Ordinal",
    xlim      = xlim_vec,
    line_args = line_args,
    ...
  )
}

###########
## zpois ##
###########
marginzerop <- function(u, y, pzero, meanpoisson,h) {
  n <- length(y)
  inv1 <- ifelse(u < (pzero + (1 - pzero) * (ppois(0, lambda = meanpoisson))), 0,
                 qpois(pmax((u - pzero) / (1 - pzero), 0), lambda = meanpoisson)
  )
  inv1m <- ifelse(inv1 > 0, inv1 - 1, 0)
  n1 <- cbind(inv1, inv1m)
  p1 <- cbind(
    (pzero + (1 - pzero) * (ppois(inv1, meanpoisson))),
    (pzero + (1 - pzero) * (ppois(inv1m, meanpoisson)))
  )
  ind1 <- apply(abs(p1 - u), 1, which.min)
  wei <- 1 * ((p1[cbind(1:n, ind1)] - u)^2 < 5 * h^2) *
    (1 - ((p1[cbind(1:n, ind1)] - u)^2) / h^2 / 5)
  l <- sum(wei * 1 * (y <= n1[cbind(1:n, ind1)])) / sum(wei)
  l
}
marginzerop <- Vectorize(marginzerop, "u")

### Bandwidth selection
bandwidth0p <- function(y, pzero, meanpoisson) {
  n <- length(y)
  newout <- unlist(sapply(split(cbind(qpois(0.1, meanpoisson), qpois(0.9, meanpoisson)), 1:n), listvec))
  newlambda <- rep(meanpoisson, times = qpois(0.9, meanpoisson) - qpois(0.1, meanpoisson) + 1)
  newzero <- rep(pzero, times = qpois(0.9, meanpoisson) - qpois(0.1, meanpoisson) + 1)
  newx <- newzero + (1 - newzero) * ppois(newout, newlambda)
  newy <- 1 * (rep(y, times = qpois(0.9, meanpoisson) - qpois(0.1, meanpoisson) + 1) <= newout)
  return(npregbw(ydat = newy[which(newx <= 0.9 & newx >= 0.1)], xdat = newx[which(newx <= 0.9 & newx > 0.1)], ckertype = "epanechnikov")$bw)
}

#' @keywords internal
resid.zpois_quasi <- function(model, line_args, ...){
  y           <- model$model[,1]
  meanpoisson <- predict(model, type = "count")
  pzero       <- predict(model, type = "zero")

  h       <- bandwidth0p(y = y, pzero = pzero, meanpoisson = meanpoisson)
  x.input <- seq(0, 1, length.out = 101)

  y_curve <- marginzerop(
    u           = x.input,
    y           = y,
    pzero       = pzero,
    meanpoisson = meanpoisson,
    h           = h
  )

  xlim_vec <- c(
    min(pzero + (1 - pzero) * ppois(0, meanpoisson)),
    1
  )

  .resid_quasi_finalize(
    x         = x.input,
    y         = y_curve,
    main      = "Quasi, 0-Inflated Poisson",
    xlim      = xlim_vec,
    line_args = line_args,
    ...
  )
}


###########
##  znb  ##
###########
marginzerop.nb <- function(u, y, pzero, size1f, mu.hat, h) {
  n <- length(y)
  inv1 <- ifelse(u < (pzero + (1 - pzero) * (pnbinom(0, size = size1f, mu = mu.hat))), 0,
                 qnbinom(pmax((u - pzero) / (1 - pzero), 0), mu = mu.hat, size = size1f)
  )
  inv1m <- ifelse(inv1 > 0, inv1 - 1, 0)
  n1 <- cbind(inv1, inv1m)
  p1 <- cbind(
    (pzero + (1 - pzero) * (pnbinom(inv1, size=size1f, mu=mu.hat))),
    (pzero + (1 - pzero) * (pnbinom(inv1m, size=size1f, mu=mu.hat)))
  )
  ind1 <- apply(abs(p1 - u), 1, which.min)
  wei <- 1 * ((p1[cbind(1:n, ind1)] - u)^2 < 5 * h^2) *
    (1 - ((p1[cbind(1:n, ind1)] - u)^2) / h^2 / 5)
  l <- sum(wei * 1 * (y <= n1[cbind(1:n, ind1)])) / sum(wei)
  l
}

marginzerop.nb <- Vectorize(marginzerop.nb, "u")

### Bandwidth selection
bandwidth0p.nb <- function(y, pzero, mu.hat, size1f) {
  n <- length(y)
  newout <- unlist(sapply(split(cbind(qnbinom(0.1, mu=mu.hat, size=size1f), qnbinom(0.9, mu=mu.hat, size=size1f)), 1:n), listvec))
  newlambda <- rep(mu.hat, times = qnbinom(0.9, mu=mu.hat, size=size1f) - qnbinom(0.1, mu=mu.hat, size=size1f) + 1)
  newzero <- rep(pzero, times = qnbinom(0.9, mu=mu.hat, size=size1f) - qnbinom(0.1, mu=mu.hat, size=size1f) + 1)
  newx <- newzero + (1 - newzero) * ppois(newout, newlambda)
  newy <- 1 * (rep(y, times = qnbinom(0.9, mu=mu.hat, size=size1f) - qnbinom(0.1, mu=mu.hat, size=size1f) + 1) <= newout)

  return(npregbw(ydat = newy[which(newx <= 0.9 & newx >= 0.1)], xdat = newx[which(newx <= 0.9 & newx > 0.1)], ckertype = "epanechnikov")$bw)
}

#' @keywords internal
resid.znb_quasi <- function(model, line_args, ...){
  y      <- model$model[,1]
  mu.hat <- predict(model, type = "count")
  size1f <- model$theta
  pzero  <- predict(model, type = "zero")

  h       <- bandwidth0p.nb(y = y, pzero = pzero, mu.hat = mu.hat, size1f = size1f)
  x.input <- seq(0,1, length.out=101)

  y_curve <- marginzerop.nb(
    u      = x.input,
    y      = y,
    pzero  = pzero,
    size1f = size1f,
    mu.hat = mu.hat,
    h      = h
  )

  xlim_vec <- c(
    min(pzero + (1 - pzero) * pnbinom(0, size = size1f, mu = mu.hat)),
    1
  )

  .resid_quasi_finalize(
    x         = x.input,
    y         = y_curve,
    main      = "Quasi, 0-Inflated NB",
    xlim      = xlim_vec,
    line_args = line_args,
    ...
  )
}