bernoulli_cusum function

Risk-adjusted Bernoulli CUSUM

Risk-adjusted Bernoulli CUSUM

This function can be used to construct a risk-adjusted Bernoulli CUSUM chart for survival data. It requires the specification of one of the following combinations of parameters as arguments to the function:

  • glmmod & theta
  • p0 & theta
  • p0 & p1
bernoulli_cusum(data, followup, glmmod, theta, p0, p1, h, stoptime, assist,
  twosided = FALSE)

Arguments

  • data: A data.frame containing the following named columns for each subject:

    • entrytime:: time of entry into study (numeric);
    • survtime:: time from entry until event (numeric);
    • censorid:: censoring indicator (0 = right censored, 1 = observed) (integer);

    and optionally additional covariates used for risk-adjustment.

  • followup: The value of the follow-up time to be used to determine event time. Event time will be equal to entrytime + followup for each subject.

  • glmmod: Generalized linear regression model used for risk-adjustment as produced by the function glm(). Suggested:

    glm(as.formula("(survtime <= followup) & (censorid == 1) ~covariates"), data = data).

    Alternatively, a list containing the following elements:

    • formula:: a formula() in the form ~ covariates;
    • coefficients:: a named vector specifying risk adjustment coefficients for covariates. Names must be the same as in formula and colnames of data.

  • theta: The θ\theta value used to specify the odds ratio eθe^\theta under the alternative hypothesis. If θ>=0\theta >= 0, the chart will try to detect an increase in hazard ratio (upper one-sided). If θ<0\theta < 0, the chart will look for a decrease in hazard ratio (lower one-sided). Note that

p1=p0eθ1−p0+p0eθ.p1=(p0∗eθ)/(1−p0+p0∗eθ). p_1 = \frac{p_0 e^\theta}{1-p_0 +p_0 e^\theta}.p1 = (p0 * e^\theta)/(1-p0+p0 * e^\theta).
  • p0: The baseline failure probability at entrytime + followup for individuals.
  • p1: The alternative hypothesis failure probability at entrytime + followup for individuals.
  • h: (optional): Control limit to be used for the procedure.
  • stoptime: (optional): Time after which the value of the chart should no longer be determined.
  • assist: (optional): Output of the function parameter_assist()
  • twosided: (optional): Should a two-sided Bernoulli CUSUM be constructed? Default is FALSE.

Returns

An object of class bercusum containing:

  • CUSUM: A data.frame containing the following named columns:

    • time:: times at which chart is constructed;
    • value:: value of the chart at corresponding times;
    • numobs:: number of observations at corresponding times.

  • call: the call used to obtain output;

  • glmmod: coefficients of the glm() used for risk-adjustment, if specified;

  • stopind: indicator for whether the chart was stopped by the control limit.

Details

The Bernoulli CUSUM chart is given by

Sn=max⁡(0,Sn−1+Wn),Sn=max(0,Sn−1+Wn), S_n = \max(0, S_{n-1} + W_n),S_n = max(0, S_{n-1} + W_n),

where

Wn=Xnln⁡(p1(1−p0)p0(1−p1))+ln⁡(1−p11−p0)Wn=Xnln((p1∗(1−p0))/(p0∗(1−p1)))+ln((1−p1)/(1−p0)) W_n = X_n \ln \left( \frac{p_1 (1-p_0)}{p_0(1-p_1)} \right) + \ln \left( \frac{1-p_1}{1-p_0} \right)W_n = X_n ln((p_1 * (1-p_0))/(p_0 * (1-p_1))) + ln((1-p_1)/(1-p_0))

and XnX_n is the outcome of the nn-th (chronological) subject in the data. In terms of the Odds Ratio:

Wn=Xnln⁡(eθ)+ln⁡(11−p0+eθp0)Wn=Xnln(exp(theta))+ln((1)/(1−p0−exp(theta)∗p0)) W_n = X_n \ln \left( e^\theta \right) + \ln \left( \frac{1}{1-p_0 + e^\theta p_0} \right)W_n = X_n ln(exp(theta)) + ln((1)/(1-p_0 - exp(theta) * p_0))

For a risk-adjusted procedure (when glmmod is specified), a patient specific baseline failure probability p(0i)p_(0i) is modelled using logistic regression first. Instead of the standard practice of displaying patient numbering on the x-axis, the time of outcome is displayed.

Examples

#We consider patient outcomes 100 days after their entry into the study.
followup <- 100
#Determine a risk-adjustment model using a generalized linear model.
#Outcome (failure within 100 days) is regressed on the available covariates:
exprfitber <- as.formula("(survtime <= followup) & (censorid == 1)~ age + sex + BMI")
glmmodber <- glm(exprfitber, data = surgerydat, family = binomial(link = "logit"))
#Construct the Bernoulli CUSUM on the 1st hospital in the data set.
bercus <- bernoulli_cusum(data = subset(surgerydat, unit == 1), glmmod = glmmodber,
 followup = followup, theta = log(2))
#Plot the Bernoulli CUSUM
plot(bercus)

See Also

plot.bercusum, runlength.bercusum

Author(s)

Daniel Gomon

success packageRead PDF manual