In umich-biostatistics/corrsurv: Two-Sample Tests for Correlated Survival Endpoints

options(digits = 3)

Introduction

This vignette presents a step-by-step guide to the analysis of correlated survival data using the methods of Tayob and Murray. By the end of the tutorial, you should be familiar with the functions available to users of corrsurv and in creating summaries, plots, and extracting other results.

Summary of methods

This package implements methods from the following three papers by Tayob and Murray:

Method 1: Murray, 2000

Description

Perform two-sample tests for treatment effects with paired censored survival data.

Reference

Murray, Susan. Nonparametric Rank-Based Methods for Group Sequential Monitoring of Paired Censored Survival Data. 2000. Biometrics, 56, pp. 984-990.

Jump to the pairtest() section for the methods from this paper.

Method 2: Tayob and Murray, 2014

Description

Perform the Tayob and Murray two-sample recurrent events test.

Reference

Tayob, N. and Murray, S., 2014. Nonparametric tests of treatment effect based on combined endpoints for mortality and recurrent events. Biostatistics, 16(1), pp.73-83.

Jump to the TM() section for the methods from this paper.

Method 3: Tayob and Murray, 2016

Description

Estimate the tau-restricted mean survival across multiple follow-up intervals.

Reference

Tayob, N. and Murray, S., 2016. Nonparametric restricted mean analysis across multiple follow-up intervals. Statistics & probability letters, 109, pp.152-158.

Jump to the TM2() section for the methods from this paper.

Install

First, install and load the corrsurv package from Github using the following code:

If you have not installed devtools, run the following:

install.packages("devtools")

If devtools is sucessfully installed, run:

if(!require(corrsurv)) {
  library(devtools)
  install_github('umich-biostatistics/corrsurv') 
}

Load the package:

library(corrsurv)

Results

pairtest()

Example 1: Analysis of package data set pairdata.

Enrolled 3711 patients with mild-to-severe nonproliferative or early proliferative diabetic retinopathy in both eyes from April 1980 to July 1985.

• One eye per patient was randomized to early photocoagulation
• The other eye deferred photocoagulation until detection of high-risk proliferative retinopathy
• Survival endpoint: Time to severe visual loss (defined as visual acuity < 5/200 at two consecutive visits) subject to censoring
• Eventually detected a benefit with early photocoagulation
• After 9 years of follow-up, only 5.5% of patients had experienced the event (lots of censoring)

Important Feature of ETDRS Data Set (pairdata): Paired Event Times

Paired Design Attractive:

• Historically used to minimize variability in treatment comparisons
• With positive correlation in pairs that is taken into account in an analysis:
  • higher power for fixed sample size than independent designs
  • smaller required sample size for fixed power than independent designs
  • At the time the ETDRS study was published, the investigators did not have access to these paired tests for censored survival data, but they recognized the statistical issues relating to the correlated structure of the data. Some exploratory analysis on their part suggested “that not taking pairing into account led to conservative tests (ETDRS Research Group, 1991a, p. 749).”

Here is a preview of the data structure:

head(pairdata)

Documentation for pairtest:

?pairtest

The function pairtest requires the time-to-event data for each pair (x1 and x2) and the event indicator for each of the times (delta1 for x1, and delta2 for x2). The event indicators (delta1, delta2) equal 1 if the event occurred and 0 if the patient was censored.

The call to pairtest produces the results for Murray's paired test:

eyeresults = pairtest(x1 = pairdata$x1, delta1 = pairdata$delta1, 
                      x2 = pairdata$x2, delta2 = pairdata$delta2, n = 3711)

Print a summary of the results:

summary(eyeresults)

If you need to access any result not printed by pairtest, use the $ operator to search available objects to extract. The test statistic and p-value for each test is returned from the call to pairtest in a list.

Here are the printed and stored results extracted from pairtest:

# Upper limit of integration
(integration_upper_limit = eyeresults$upperlim)

# Logrank statistic
(logrank_statistic = eyeresults$lgrk.stat.p)
(logrank_statistic_p = 2*(pnorm(logrank_statistic, lower.tail = T))) # since negative statistic

# Gehan statistic
(gehan_statistic = eyeresults$gehan.stat.p)
(gehan_statistic_p = 2*(pnorm(gehan_statistic, lower.tail = T))) # since negative statistic

# Years of Life statistic
(yls_statistic = eyeresults$yls.stat.p)
(yls_statistic_p = 2*(pnorm(yls_statistic, lower.tail = F))) # since positive statistic

# Pepe Fleming statistic
(pepe_flem_statistic = eyeresults$pf.stat.p)
(pepe_flem_statistic_p = 2*(pnorm(pepe_flem_statistic, lower.tail = F))) # since positive statistic

# Logrank statistic assuming independence
(logrank_assuming_indep = eyeresults$lgrk.nopair.stat.p)
(logrank_assuming_indep_p = 2*(pnorm(logrank_assuming_indep, lower.tail = T)) ) # since negative statistic

# Gehan statistic assuming independence
(gehan_assuming_indep = eyeresults$gehan.nopair.stat.p)
(gehan_assuming_indep_p = 2*(pnorm(gehan_assuming_indep, lower.tail = T))) # since negative statistic

# Years of Life assuming independence
(yls_assuming_indep = eyeresults$yls.nopair.stat.p)
(yls_assuming_indep_p = 2*(pnorm(yls_assuming_indep, lower.tail = F))) # since positive statistic

# Pepe Fleming assuming independence
(pf_assuming_indep = eyeresults$pf.nopair.stat.p)
(pf_assuming_indep_p = 2*(pnorm(pf_assuming_indep, lower.tail = F))) # since positive statistic

# group size adjustment
nstuff = ((3711*3711) / (3711 + 3711))^(.5)

# estimate area between survival curves and 95% conf. int.
(yls.diff = eyeresults$yls.num / nstuff)

(yls.upper = yls.diff + 1.96*sqrt(eyeresults$yls.type.var) / nstuff)
(yls.lower = yls.diff - 1.96*sqrt(eyeresults$yls.type.var) / nstuff)

# 95% CI assuming independence
(yls.upper2 = yls.diff + 1.96*sqrt(eyeresults$yls.nopair.type.var) / nstuff)
(yls.lower2 = yls.diff - 1.96*sqrt(eyeresults$yls.nopair.type.var) / nstuff)

Other options in summary and plot methods

The result of summary() of a pairtest object has a digits option to specify more or less precision than the default:

summary(eyeresults, digits = 6)

TM()

Repeated time-to-event data with two treatment groups.

The example data set contains 200 observations of 10 variables, where each row contains a subject ID, one covariate, censoring indicator, group indicator, and Z (a set of times to each follow up).

data("TMdata")

Preview the data set:

head(TMdata)

Documentation for TM:

?TM

This method requires some pre-processing of the data. Follow our example:

N = nrow(TMdata) # Number of subjects
X = TMdata$X
delta = TMdata$delta
table(delta)
Z = TMdata[,c('Z1', 'Z2', 'Z3', 'Z4', 'Z5', 'Z6')]
A = max(X) # length of the study
Treatment = as.numeric(TMdata$Group == 1) #1 if Case and 0 if Control
table(Treatment, delta)

The call to TM produces the results for the Tayob and Murray two-sample recurrent events test.

The two available methods for assessing treatment differences include "average" for mean difference test or "area" for area between the RMRL curves. The method = "average" default is assumed if not explicitly specified.

srec.average = TM(X = X, delta = delta, Z = Z, Group = Treatment, Tau = A/2, 
                  t = seq(from = 0, to = A-A/2, by = A/4))

Print a summary of the results:

summary(srec.average)

Now with method = "area":

srec.area = TM(X = X, delta = delta, Z = Z, Group = Treatment, Tau = A/2, 
               t = seq(from = 0, to = A-A/2, by = A/4), method = "area")

summary(srec.area)

If the area method is used, a plot option is also available:

plot(srec.area)

Other options in summary and plot methods

The result of summary() of a TM object has a digits option to specify more or less precision than the default:

summary(srec.area, digits = 6)

TM2()

The function TM2() allows the user to compute restricted mean survival across multiple follow-up intervals.

The reference paper is Tayob, N. and Murray, S., 2016. Nonparametric restricted mean analysis across multiple follow-up intervals. Statistics & probability letters, 109, pp.152-158.

This method allows for flexibility in choosing variance estimate for the restricted mean survival. The options are proposed, sandwich, and independence. See the reference paper for details.

The data set contains follow-up times over multiple intervals (X) and the censoring indicator (delta; 1 if subject died, 0 if censored). The data set is in the "long" format where multiple events are each recorded as a new row rather than a new column.

The data set contains a total of 100 follow-up times.

data("TM2data")
#?TM2data

The maximum recorded time in the data set is:

max(TM2data$X)

Similar to the other methods in this package, TM2 is the work-horse estimation function and calling summary() on the resultant fit object produces a printed summary of the output, and calling plot() produces a plot of the estimate and 95% bounds.

Documentation for TM2:

?TM2

An example of analysis with tau = 12 and all three variance estimates:

output = TM2(X = TM2data$X, delta = TM2data$delta, Tau = 12, 
             t = seq(from = 0, to = 24, by = 6), var_output = "all")
summary(output)
plot(output)

An example of analysis with tau = 12 and proposed variance estimate only (default):

output = TM2(X = TM2data$X, delta = TM2data$delta, Tau = 12, 
             t = seq(from = 0, to = 24, by = 6))
summary(output)
plot(output)

An example of analysis with tau = 12 and independent variance estimate only:

output = TM2(X = TM2data$X, delta = TM2data$delta, Tau = 12, 
             t = seq(from = 0, to = 24, by = 6), var_output = "independence")
summary(output)
plot(output)

An example of analysis with tau = 12 and sandwich variance estimate only:

output = TM2(X = TM2data$X, delta = TM2data$delta, Tau = 12, 
             t = seq(from = 0, to = 24, by = 6), var_output = "sandwich")

summary(output)
plot(output)

Other options in summary and plot methods

Use ? to explore other options to summary and plot:

?summary.TM2

As you can see from the Usage section of the document, the digits = argument can optionally be specified if you prefer more or less precision.

?plot.TM2

Optional arguments for the plot method include alpha = 0.05, conservative_index = 25, k = 500, n.sim = 1000.

Here is a summary from the documentation:

• alpha is the significance level for the confidence bands
• conservative_index is the minimum number of events after the start of the last interval
• k is the number of points for integration within follow-up window [0,Tau]
• n.sim is the number of samples simulated to calculate confidence bands

These arguments can be adjusted to suit the user's plotting needs.

umich-biostatistics/corrsurv documentation built on Jan. 11, 2020, 2:03 a.m.