nll_frailty_correlated function

Negative log-likelihood function: correlated frailty model

Function calculating the negative log-likelihood (nll) for patterns of mortality in infected and uninfected treatments when there is unobserved variation in background mortality and mortality due to infection, where these two sources of variation are positively correlated.

nll_frailty_correlated(
  a1 = a1,
  b1 = b1,
  a2 = a2,
  b2 = b2,
  theta01 = theta01,
  theta02 = theta02,
  rho = rho,
  data = data,
  time = time,
  censor = censor,
  infected_treatment = infected_treatment,
  d1 = "Weibull",
  d2 = "Weibull"
)

Arguments

• a1, b1: location and scale parameters describing background mortality
• a2, b2: location and scale parameters describing mortality due to infection
• theta01: parameter describing variance of unobserved variation in background rate of mortality
• theta02: parameter describing variance of unobserved variation in rate of mortality due to infection
• rho: parameter for strength of positive correlation between theta01 and theta02; rho > 0
• data: name of data frame containing survival data
• time: name of data frame column identifying time of event; time > 0
• censor: name of data frame column idenifying if event was death (0) or right-censoring (1)
• infected_treatment: name of data frame column identifying if data are from an infected (1) or uninfected (0) treatment
• d1, d2: names of probability distributions chosen to describe background mortality and mortality due to infection, respectively; both default to the Weibull distribution

Details

This function assumes both sources of unobserved variation follow the Gamma distribution where both have a mean = 1.0, and variances 'theta01' and 'theta02', respectively. The nll function is calculated based on seven parameters; location and scale parameters for background mortality and mortality due to infection, respectively, the unobserved variance of in each source of mortality, and the strength of the positive correlation between them.

Examples

# step #1: parameterise nll function to be passed to 'mle2'
    m01_prep_function <- function(
      a1 = a1, b1 = b1, a2 = a2, b2 = b2,
      theta01 = theta01, theta02 = theta02, rho = rho){
        nll_frailty_correlated(
          a1 = a1, b1 = b1, a2 = a2, b2 = b2,
          theta01 = theta01, theta02 = theta02, rho = rho,
          data = data_lorenz,
          time = t,
          censor = censored,
          infected_treatment = g,
          d1 = "Gumbel",
          d2 = "Gumbel"
          )}

# step #2: send 'prep_function' to 'mle2' for maximum likelihood estimation
#          lower bounds of estimates set to 1e-6
  m01 <- mle2(m01_prep_function,
    start = list(a1 = 20, b1 = 5, a2 = 20, b2 = 4,
                 theta01 = 1, theta02 = 1, rho = 1),
    method = "L-BFGS-B",
    lower = list(a1 = 1e-6, b1 = 1e-6, a2 = 1e-6, b2 = 1e-6,
                 theta01 = 1e-6, theta02 = 1e-6, rho = 1e-6)
                 )

  summary(m01)

# NB no std. errors estimated and estimate of 'theta01' at lower boundary

# rerun model with theta01 set at lower boundary
    m02 <- mle2(m01_prep_function,
      start = list(a1 = 20, b1 = 5, a2 = 20, b2 = 4,
                   theta01 = 1, theta02 = 1, rho = 1),
      fixed = list(theta01 = 1e-6),
      method = "L-BFGS-B",
      lower = list(a1 = 1e-6, b1 = 1e-6, a2 = 1e-6,
                   b2 = 1e-6, theta02 = 1e-6, rho = 1e-6)
      )

summary(m02)

# NB std. error of 'rho' crosses lower boundary

# rerun model with theta01 and rho set at lower limits
    m03 <- mle2(m01_prep_function,
      start = list(a1 = 20, b1 = 5, a2 = 20, b2 = 4,
                   theta01 = 1, theta02 = 1, rho = 1),
      fixed = list(theta01 = 1e-6, rho = 1e-6),
      method = "L-BFGS-B",
      lower = list(a1 = 1e-6, b1 = 1e-6, a2 = 1e-6,
                   b2 = 1e-6, theta02 = 1e-6)
      )

    summary(m03)

# result of m03 corresponds with estimates from 'nll_frailty' model,
  # i.e., where it is assumed there is no unobserved variation in the
  # rate of background mortality and where the gamma distribution describes
  # the unobserved variation in virulence

# step #1: parameterise nll function to be passed to 'mle2'
    m04_prep_function <- function(a1, b1, a2, b2, theta){
      nll_frailty(a1, b1, a2, b2, theta,
        data = data_lorenz,
        time = t,
        censor = censored,
        infected_treatment = g,
        d1 = "Gumbel",
        d2 = "Gumbel",
        d3 = "Gamma"
        )}

# step #2: send 'prep_function' to 'mle2' for maximum likelihood estimation
    m04 <- mle2(m04_prep_function,
      start = list(a1 = 20, b1 = 5, a2 = 20, b2 = 5, theta = 2)
      )

    summary(m04)

# compare estimated values m03 vs. m04
    coef(m03) ; coef(m04)

See Also

nll_frailty