discretised_logit_hazard.stan

Go to the documentation of this file.

 
real loglogistic_lcdf (real y, real alpha, real beta) {
  return -log1p((y / alpha) ^-beta);
}
 
vector lcdf_discretised(real mu, real sigma, int u, int dist) {
  vector[u] integer_lcdf;
  if (dist == 1) {
    real emu = exp(-mu);
    for (i in 1:u) {
      integer_lcdf[i] = exponential_lcdf(i | emu);
    }
  } else if (dist == 2) {
    for (i in 1:u) {
      integer_lcdf[i] = lognormal_lcdf(i | mu, sigma);
    }
  } else if (dist == 3) {
    real emu = exp(mu);
    for (i in 1:u) {
      integer_lcdf[i] = gamma_lcdf(i | emu, sigma);
    }
  } else if (dist == 4) {
    real emu = exp(mu);
    for (i in 1:u) {
      integer_lcdf[i] = loglogistic_lcdf(i | emu, sigma);
    }
  } else {
    reject("Unknown distribution function provided.");
  }
  return(integer_lcdf);
}
 
vector normalise_lcdf_as_uniform_double_censored(vector lcdf, int u,
  int max_strat) {
  vector[u] adjusted_lcdf;
  if (max_strat == 0) {
    // ignore (i.e. sum of probabilities does not add to 1)
    // NOTE: This cannot be used when using lcdf_to_logit_hazard as
    // it will return the same result as using strategy 2.
    adjusted_lcdf = lcdf;
  } else if (max_strat == 1) {
    // add to upper bound: cdf_u = 1 (i.e. lcdf_u = 0)
    adjusted_lcdf = lcdf;
    adjusted_lcdf[u] = 0;
    // normalise to account for double censoring
    adjusted_lcdf = adjusted_lcdf - log(2);
  } else if (max_strat == 2) {
    // normalize: cdf = cdf / (cdf[u] + cdf[u-1])
    // both u and u-1 are in the denominator to account for the 2 day width of
    // the (overlapping) intervals from the uniform interval approximation:
    // For a delay d, the interval is [d-1, d+1], e.g. the interval for delay 0 
    // goes from -1 to 1, and the interval for delay 1 goes from 0 to 2 in the 
    // LCDF, respectively. For a maximum delay of D days, we therefore evaluate 
    // the LCDF up to the integer u = D+1. The total mass of this LCDF is thus:
    // sum_{d=0}^D (F_{d+1} - F_{d-1}) = sum_{i=1}^u (F_{i} - F_{i-2}) = F_{u} + F_{u-1}
    // If u were infinite this would be equivalent to dividing by 2.
    if (u == 1) {
      adjusted_lcdf = lcdf - lcdf[1];
    }
    if (u > 1) {
      adjusted_lcdf = lcdf - log_sum_exp(lcdf[u], lcdf[u-1]);
    }
  } else {
    reject("Unknown strategy to handle probability mass beyond the maximum value.");
  }
  return(adjusted_lcdf);
}
 
vector lcdf_to_uniform_double_censored_log_prob(vector lcdf, int u) {
  vector[u] lpmf; 
  // p_0 = cdf_1 - cdf_0
  lpmf[1] = lcdf[1];
  if (u > 1) {
    // p_1 = cdf_2 - cdf_0
    lpmf[2] = lcdf[2];
    if (u > 2) {
      // p_u = cdf_u+1 - cdf_u-1
      lpmf[3:u] = log_diff_exp(lcdf[3:u], lcdf[1:(u-2)]);
    }
  }
  return(lpmf);
}
 
vector lprob_to_uniform_double_censored_log_hazard(vector lprob, vector lcdf,
   int u) {
  vector[u] lhaz;
  // h_0 = F_1
  lhaz[1] = lcdf[1];
  // h_d = p_d / (1 - sum^{d-1}_{i=0} p_i)
  // h_d = (F_d+1 - F_d-1) / (1 - sum^{d-1}_{i=0} F_i+1 - F_i-1)
  // h_d = (F_d+1 - F_d-1) / (1 - (F_d + F_d-1))
  // h_d = (F_d+1 - F_d-1) / ((1 - F_d) - F_d-1)
  // log(h_d) = log(F_d+1 - F_d-1) - log((1 - F_d) - F_d-1)
  if (u > 1) {
    vector[u-2] lccdf;
    // ccdf_i = 1 - cdf_i
    lccdf = log1m_exp(lcdf[1:(u-2)]);
    lhaz[2] = lprob[2] - lccdf[1];
    if (u > 2) {
      // Note that lprob[1] = log(p_0), lcdf[1] = log(F_1), lccdf[1] = log(1 - F_1)
      // e.g. lhaz[3] = lprob[3] - log_diff_exp(lccdf[2], lcdf[1])
      // means log(h_2) = log(p_2) - log(1 - F_2, F_1)
      lhaz[3:(u-1)] = lprob[3:(u-1)] - log_diff_exp(lccdf[2:(u-2)], lcdf[1:(u-3)]);
    }
  }
  return(lhaz);
}
 
vector log_hazard_to_logit_hazard(vector lhaz) {
  int n = num_elements(lhaz);
  vector[n] logit_haz;
  // Logit transformation
  logit_haz[1:(n-1)] = lhaz[1:(n-1)] - log1m_exp(lhaz[1:(n-1)]);
  // Set last logit transformed hazard to Inf (i.e h[n] = 1)
  logit_haz[n] = positive_infinity();
  return(logit_haz);
}
 
vector discretised_logit_hazard(real mu, real sigma, int dmax, int dist, 
                                int max_strat, int ref_as_p) {
  vector[dmax] lcdf;
  vector[dmax] lprob;
  vector[dmax] logit_haz; 
  lcdf = lcdf_discretised(mu, sigma, dmax, dist);
  lcdf = normalise_lcdf_as_uniform_double_censored(lcdf, dmax, max_strat);
  lprob = lcdf_to_uniform_double_censored_log_prob(lcdf, dmax);
  if (ref_as_p == 1) {
    // In the mode where there are no hazard effects downstream functions
    // make use of the log probability directly so we return it here without
    // converting to the logit hazard.
    logit_haz = lprob;
  }else{
    vector[dmax] lhaz;
    lhaz = lprob_to_uniform_double_censored_log_hazard(lprob, lcdf, dmax);
    logit_haz = log_hazard_to_logit_hazard(lhaz);
  }
  return(logit_haz);
}