R/prevost.R

In ellessenne/random: A Random Collection of R Things

#' @title Incidence of Interval Cancers Modelled as a Function of Mean Sojourn Time
#' @description Model 1 from Provost _et al_. (1998), Am J Epidemiol
#' @param par Parameter \eqn{\lambda} of the model
#' @param J Incidence of preclinical disease
#' @param data Dataset
#' @source \url{https://academic.oup.com/aje/article-lookup/doi/10.1093/oxfordjournals.aje.a009687}
#' @return A numeric value representing the negative log-likelihood of the model
#' @export
#'
#' @examples
#' data("crcscreen", package = "tag")
#' df <- na.omit(crcscreen[crcscreen$ageD == "55-64" & crcscreen$t > 0, ])
#' J <- crcscreen[crcscreen$ageD == "55-64" & crcscreen$group == "Controls", ]
#' J <- J$r / J$py
#' optimise(
#'   f = prevost_m1,
#'   interval = c(0, 100),
#'   data = df,
#'   J = J
#' )
prevost_m1 <- function(par, J, data) {
  lambda <- par[1]
  i <- (J * data$py) * (1 - exp(-lambda * (data$t - 0.5)))
  ll <- sum(data$r * log(i) - i - log(data$r))
  return(-ll)
}

#' @title Incidence of Interval Cancers and Screen-detected Cancers Modelled as a Function of Sensitivity and Mean Sojourn Time
#' @description Model 2 from Provost _et al_. (1998), Am J Epidemiol
#' @param par Parameter \eqn{\lambda} of the model
#' @param J Incidence of preclinical disease
#' @param data Dataset
#' @source \url{https://academic.oup.com/aje/article-lookup/doi/10.1093/oxfordjournals.aje.a009687}
#' @return A numeric value representing the negative log-likelihood of the model
#' @export
#'
#' @examples
#' data("crcscreen", package = "tag")
#' df <- na.omit(crcscreen[crcscreen$ageD == "55-64", ])
#' J <- crcscreen[crcscreen$ageD == "55-64" & crcscreen$group == "Controls", ]
#' J <- J$r / J$py
#' optim(
#'   par = c(0.5, 0.5),
#'   f = prevost_m2,
#'   data = df,
#'   J = J
#' )
prevost_m2 <- function(par, J, data) {
  lambda <- par[1]
  S <- par[2]
  screened <- data[data$t == 0, ]
  data <- data[data$t > 0, ]
  avgT <- mean(as.numeric(c(substr(screened$ageD, 4, 5), substr(screened$ageD, 1, 2))))
  n <- screened$py
  c <- screened$r
  P <- (n * S * J * (exp(-lambda * avgT) - exp(-J * avgT)) / (J - lambda)) / (exp(-J * avgT) + J * (exp(-lambda * avgT - exp(-J * avgT))) / (J - lambda))
  I <- (J * data$py) * (1 - exp(-lambda * (data$t - 0.5))) + c * (1 - S) / S * (exp(-lambda * (data$t - 1)) - exp(-lambda * data$t))
  ll <- c * log(P / n) + (n - c) * log(1 - P / n) + sum(data$r * log(I) - I - log(data$r))
  return(-ll)
}