R/logitexp.R

In trobinj/trtools: Miscellaneous Tools for Teaching Statistics

#' Logistic-exposure link function. 
#' 
#' This is a link function proposed by Schaffer (2004) for discrete time survival modeling. Let \eqn{T} be the number of discrete time units until an event (e.g., days), and assume that for a given unit the hazard function \eqn{P(T > t | T \ge t) = \pi} for all \eqn{t} (i.e., the probability of surviving time unit \eqn{t} is \eqn{\pi}). Then the probability of survived \eqn{e} time units is \eqn{\pi^e}. Assume a logistic regression model for survival on a given day so that \eqn{\pi = 1/(1 + exp(-\eta))}. Let \eqn{Y} be the number of observational units out of \eqn{m} that survive an exposure of \eqn{e} time units. Then the inverse link function is \eqn{E(Y/m) = 1/(1 + exp(-\eta))^e}, and the link function is \eqn{[\log(E(Y/m)^{1/e}/(1 - E(Y/m)^{1/e})] = \eta}.
#' 
#' @details The argument \code{exposure} specifies the number of time units. To be consistent with the motivation given above this should be a positive integer. 
#' 
#' @seealso \code{\link[trtools]{logitcomp}}
#' 
#' @note The code for this link function is also given as an example in the documentation for \code{family}. 
#' 
#' @source Schaffer, T. L. (2004). A unified approach to analyzing nest success. \emph{The Auk}, \emph{121(2)}, 526-540.
#' 
#' @examples
#' library(dplyr)
#' 
#' set.seed(123)
#' 
#' d <- data.frame(x = seq(0, 3, length = 100)) %>%
#'  mutate(days = sample(1:10, n(), replace = TRUE)) %>% 
#'  mutate(y = rbinom(n(), 1, plogis(x)^days))
#' 
#' m <- glm(y ~ x, data = d, 
#'  family = binomial(link = logitexp(exposure = d$days)))
#'  
#' # The logitcomp link function can be used to model the probability
#' # of not surviving. The following produces the same estimates except
#' # for a reversal of sign. 
#' 
#' m <- glm(1 - y ~ x, data = d, 
#'  family = binomial(link = logitcomp(m = d$days)))
#' @export
logitexp <- function(exposure = 1)
{
  linkfun <- function(mu) {
    qlogis(mu^(1/exposure))
  }
  linkinv <- function(eta) {
    plogis(eta)^exposure
  }
  mu.eta  <- function(eta) {
      exposure * plogis(eta)^(exposure-1) * binomial()$mu.eta(eta)
  }
  valideta <- function(eta) TRUE
  link <- paste0("logexp(", exposure, ")")
  structure(list(linkfun = linkfun, linkinv = linkinv,
    mu.eta = mu.eta, valideta = valideta, name = link), class = "link-glm")
}