Nothing
Cox regression for right-truncated data"
In coxrt: Cox Proportional Hazards Regression for Right-Truncated Data
The coxrt package accompanies the paper of Vakulenko-Lagun, Mandel and Betensky (2019) and is designed for analysis of right-truncated data.
What are right-truncated data?
When we are interested in the analysis of lifetimes that start at an initiating event and end up when an event of interest occurs, but we cannot obtain a population of interest to follow them prospectively, retrospective ascertainment of cases could be a good alternative. However, this might result in right-truncated survival data.
Right truncation often happens when we cannot identify the population of interest since their initiating events are latent (like the times of infections), and the population becomes "visible" only when another event, that can be more easily identified, happens. Here we assume that for those who are sampled, the times of initiating events are known or can be reliably estimated.
One example of such type of data is the data set on AIDS cases resulted from contaminated blood transfusion [@Kalbfeisch]. Here, the time of interest is the incubation period of AIDS, i.e., the time from HIV infection to development of AIDS. Figure 1 shows schematically how these data were collected. The individuals who developed AIDS by June 30, 1986 (the date when the AIDS cases were retrospectively ascertained) were selected, satisfying by this the condition $T_i< R_i$, where $T_i$ is an incubation period for the $i^{th}$ person and $R_i$ is his/her truncation time, defined as the time from HIV infection to June 30, 1986.
In the two functions of the package, coxph.RT and coxph.RT.a0, we assume that $T_i$ are observed exactly and there is no censoring. The analysis is based only on sampled events. Those who belong to the population but have not developed AIDS by the time of collection of data remain unobserved.
knitr::include_graphics('slide1_116.png')
The sampling condition $T_i< R_i$ means that we tend to sample shorter lifetimes than longer ones. And an important complication of such sampling scheme is that the follow-up is limited: on Figure 1, $r_$, defined as the time from the earliest HIV infection (January 1983) to June 30 1986, is the maximum follow-up time in this data set. This is an important statistic, since it allows us to realize whether our sampling covers the whole range of possible incubation periods or not. If the follow-up is too short, most likely we do not capture longer incubation periods, and they will not only be underrepresented in our sample, but even worse, they are likely to be non-represented at all. We shall call the condition $P(T>r_)=0$ a positivity assumption.
For example in the AIDS data, positivity does not hold since the incubation period of AIDS can be longer than the maximum observed follow-up time in this data set that equals 3.5 years.
Under such short retrospective follow-up there was no chance to observe longer incubation times, which means that $P(T>r_)>0$, or equivalently, the support of observed $T^$ is shorter than the support of the original $T$ in the population.
Right truncation causes distortion of the original distribution
Even if positivity holds, truncation anyway results in distorted distribution, and it is important to account for this bias.
For example, on Figure 2 we illustrate the problem of distortion of the original distribution using simulated data. Here, the population lifetime $T$ follows Weibull distribution as rweibull(20000, shape=3, scale=1).
And the population truncation time $R$ was generated (independently of $T$) from rweibull(20000, shape=2.5, scale=1).
The sampled lifetimes $T^$ comprise all observations from the population that satisfy $T_iR)=0.5$.
Figure 2 shows densities of $T$, $R$ and $T^ $. It is clear that distributions of $T$ and $T^*$ are different.
Importantly, in this example positivity holds since the supports of $T$ and $R$ are the same, so that every possible lifetime has a chance to be sampled, and this kind of bias can be easily corrected.
knitr::include_graphics('densTpTsRp.png')
Why is it important to adjust for right truncation in regression?
Even if we do not aim to estimate the whole distribution and are interested only in estimation of the effects of some covariates on the incubation period, we still need to adjust for right truncation, since, as we illustrate below, the right truncation blurs the difference between subgroups.
Suppose we have one binary covariate $Z$. Figure 3 shows distributions in the two subpopulations defined by $Z$. The lifetime $T$ is generated from the proportional hazard regression model $h(t;z)=h_0(t)\exp(\beta z)$ with $\beta=2$:
library(ggplot2)
library(survival)
set.seed(17)
N <- 20000
Z <- as.numeric(runif(N,0,1) < 0.5) # binary explanatory variable
X <-rweibull(N, shape=3, scale=1*exp(-(2*Z)/3)) # shape a=3, scale b=1, beta=2
m <- max(X)
d.tr <- data.frame(time=c(X), variable=c(Z))
d.tr$variable=as.factor(d.tr$variable)
levels(d.tr$variable) <- c( "T | Z=0", "T | Z=1", " ")
p <- ggplot(d.tr) +
geom_density(aes(x=time, color=variable, fill=variable), alpha=0.5) +
scale_color_brewer(palette = "Set1")+
scale_fill_brewer(palette = "Set1")+
theme_classic() +
labs(x = "time")+
ggtitle("In the population")+
theme(axis.text.x = element_text(face="bold"),
axis.text.y = element_text(face="bold"),
plot.title = element_text(hjust = 0.5) )+ xlim(0, m)
p
Figure 4 shows sample distributions under right truncation with positivity. For this illustration we use the same population as in Figure 3, but truncate it using uniformly distributed $R$:
set.seed(0)
r <- runif(N, 0, m)
ind <- X<r
cat("P(T<R) = selection probability = ", mean(ind), "\n")
d.tr <- data.frame(time=c(X[ind],r), variable=c(Z[ind], rep(3, N)))
d.tr$variable=as.factor(d.tr$variable)
levels(d.tr$variable) <- c( "T* | Z*=0", "T* | Z*=1", "R")
p <- ggplot(d.tr) +
geom_density(aes(x=time, color=variable, fill=variable), alpha=0.5) +
scale_color_brewer(palette = "Set1")+
scale_fill_brewer(palette = "Set1")+
theme_classic() +
labs(x = "time")+
ggtitle(paste("In the sample with R~U[0,", round(m,1), "]", sep=""))+
theme(axis.text.x = element_text(face="bold"),
axis.text.y = element_text(face="bold"),
plot.title = element_text(hjust = 0.5) ) + xlim(0, m) #+ ylim(0,6.3)
p
Despite the 30% chance of truncation, both Figures 3 and 4 look relatively similar to each other.
On Figure 5, we shorten the follow-up time by using $r_*=1.1$. This leads to violation of positivity, since longer lifetimes have no chance to be sampled, and the sub-groups become more similar to each other. In such case, without any adjustment, the estimators of $\beta$ will be biased and it is likely we will not able to discover the difference between subgroups.
set.seed(0)
r <- runif(N, 0, m/2)
ind <- X<r
d.tr <- data.frame(time=c(X[ind],r), variable=c(Z[ind], rep(3, N)))
d.tr$variable=as.factor(d.tr$variable)
levels(d.tr$variable) <- c( "T* | Z*=0", "T* | Z*=1", "R")
p <- ggplot(d.tr) +
geom_density(aes(x=time, color=variable, fill=variable), alpha=0.5) +
scale_color_brewer(palette = "Set1")+
scale_fill_brewer(palette = "Set1")+
theme_classic() +
labs(x = "time")+
ggtitle(paste("In the sample with R~U[0,", round(m/2,1), "]", sep=""))+
theme(axis.text.x = element_text(face="bold"),
axis.text.y = element_text(face="bold"),
plot.title = element_text(hjust = 0.5) ) + xlim(0, m) #+ ylim(0,6.3)
p
On Figure 6, severe violation of positivity with $R\sim U(0,0.4)$ causes the groups to be almost identical:
set.seed(0)
r <- runif(N, 0, m*4/20)
ind <- X<r
d.tr <- data.frame(time=c(X[ind],r), variable=c(Z[ind], rep(3, N)))
d.tr$variable=as.factor(d.tr$variable)
levels(d.tr$variable) <- c( "T* | Z*=0", "T* | Z*=1", "R")
p <- ggplot(d.tr) +
geom_density(aes(x=time, color=variable, fill=variable), alpha=0.5) +
scale_color_brewer(palette = "Set1")+
scale_fill_brewer(palette = "Set1")+
theme_classic() +
labs(x = "time")+
ggtitle(paste("In the sample with R~U[0,", round(m*4/20,1), "]", sep=""))+
theme(axis.text.x = element_text(face="bold"),
axis.text.y = element_text(face="bold"),
plot.title = element_text(hjust = 0.5) ) + xlim(0, m) #+ ylim(0,6.3)
p
And the situation becomes more complicated when there are several dependent covariates, as described in @RT.
What does coxrt suggest for right-truncated data?
The coxrt package allows to estimate regression coefficients $\beta$ in the proportional hazards model $h(t;z)=h_0(t)\exp(\beta z)$.
Under positivity use coxph.RT
Function coxph.RT implements estimating equations denoted as IPW-S in @RT. This function assumes positivity (as, e.g., on Figure 4).
For example, for AIDS data we can fit proportional hazards regression of incubation time on age:
library(gss)
library(coxrt)
data(aids)
all <- data.frame(age=aids$age, ageg=as.numeric(aids$age<=59), T=aids$incu, R=aids$infe, hiv.mon =102-aids$infe)
all$T[all$T==0] <- 0.5 # as in Kalbfeisch and Lawless (1989)
s <- all[all$hiv.mon>60,] # select those who were infected in 1983 or later
# analysis assuming positivity
# we request bootstrap estimate of Asymptotic Standard Error (ASE) as well:
sol <- try(coxph.RT(T~ageg, right=R, data=s, bs=TRUE, nbs.rep=500) )
knitr::kable(sol$summary) # print the summary of fit based on the analytic ASE
knitr::kable(sol$bs$summary) # print the summary of fit based on the bootstrap sample distribution
Technically, for point estimation coxph.RT uses the standard coxph function from the survival package. However, the estimator of the standard error (SE) obtained from coxph even with robust option might underestimate for smaller samples (of $<1000$ observations).
Although coxph.RT includes both analytic and bootstrap estimates of a SE, the inference based on the analytic SE estimate, implemented in coxph.RT, is more reliable than based on its bootstrap counterpart, since bootstrap resampling does not always capture the variability in data that comes from near-zero sampling probabilities $P(R>T_i)$.
Since the Inverse-Probability-Weighting (IPW) approach, implemented in coxrt, relies on "blowing up" the sampled lifetimes $T_i$ in order to obtain the unbiased pseudo-population, which should recreate the same probability law as in the original population, coxph.RT in general will give biased estimators in case there is an unobservable right tail of the distribution. Hence the function coxph.RT.a0 aims to adjust for this violation of positivity by using the structural assumption of the Cox model in the population and the information on $a_0=P(T>r_*\mid Z=0)$.
Without positivity use coxph.RT.a0
For the situations without positivity (as in Figures 5 and 6), if $a_0$ can be hypothesized, based on the subject matter knowledge, or if $a_0$ is unknown but its interval can be guessed, this information can be used in coxph.RT.a0. This function implements estimating equations denoted as IPW-SA in the original paper. IPW-SA uses stabilized weights. Technically, coxph.RT.a0 solves non-linear equations using BBsolve function from BB package using default parameter settings.
For example, in order to estimate $\beta$, the effect of age at HIV on the incubation period, given that $a_0=0.2$, we run:
# analysis using adjusted estimating equations for a0=0.2
sol.02 <- try(coxph.RT.a0(T~ageg, right=R, data=s, a0=0.2, bs=TRUE, nbs.rep = 500))
knitr::kable(round(sol.02$bs$summary,2))
We can do a sensitivity analysis for $\hat{\beta}$ for different values of a sensitivity parameter $a_0$:
# for a0=0:
sol <- try(coxph.RT(T~ageg, right=R, data=s, bs=FALSE) )
# senstivity analysis for different values of a0
a_ <- seq(0.05, 0.55, by=0.05)
est <- CI.L <- CI.U <- NULL
start_time <- Sys.time()
for(q in 1:length(a_))
{
sol.a <- try(coxph.RT.a0(T~ageg, right=R, data=s, a0=a_[q], bs=TRUE))
if (sol.a$convergence!=0)
{
cat("a0=", a_[q], ". Error occurred in BBsolve.\n")
} else
{
est <- c(est, sol.a$coef)
CI.L <- c(CI.L, sol.a$bs$summary$CI.L)
CI.U <- c(CI.U, sol.a$bs$summary$CI.U)
}
}
end_time <- Sys.time()
end_time - start_time # It takes some time to solve nonlinear estimating equations and to run bootstrap for each value of a0 ...
require(ggplot2)
res.d <- data.frame(a0=c(0, a_), beta=c(sol$coef, est),
L.95=c(sol$summary$CI.L, CI.L), U.95=c(sol$summary$CI.U, CI.U))
knitr::kable(round(res.d, 2))
p <- ggplot(res.d, aes(x=a0, y=beta)) +
geom_ribbon(aes(ymin=L.95, ymax=U.95), alpha=0.2) +
geom_line() + geom_point() +
geom_hline(yintercept=0)
p + xlab(expression( paste(a[0], "=P(T>", r['*']," | z=0)" , sep="")) )+
ylab(expression( paste(hat(beta), "(", a[0], ")" , sep="")) ) +
scale_x_continuous(breaks=res.d$a0, labels=res.d$a0) +
theme(axis.text.x = element_text(face="bold", angle=45),
axis.text.y = element_text(face="bold"))
In coxph.RT.a0, only the bootstrap estimator of SE is available.
Some remarks
-
Positivity cannot be tested. This assumption is based on the subject matter knowledge.
-
The approach assumes independence between $T$ and ($R$, $Z$) in the original population. If there is dependence between $T$ and $R$, it can be discovered by Kendall's tau test which can be done using packages tranSurv or SurvTrunc, or by more powerful permutation tests which can be done using package permDep.
-
If a non-terminal event of interest cannot be observed because death preceded it, that is the event of interest is truncated by death, the approach implemented here can be used only if two lifetimes, time to an intermediate event (=$T$) and time to death (=$R$), are independent, otherwise the results might be biased.
Some limitations of the coxrt package
-
Although in practice, it is possible to have additional left or interval censoring in right-truncated data, the package assumes that the exact lifetimes, without any censoring, are observed.
-
The current implementation allows for only time-independent covariates.
-
The package estimates only covariate effects $\beta$, and does not provide estimates for the baseline hazard $h_0(t)$ in the Cox model.
References
Try the coxrt package in your browser
Any scripts or data that you put into this service are public.
coxrt documentation built on Jan. 11, 2020, 9:39 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
The coxrt package accompanies the paper of Vakulenko-Lagun, Mandel and Betensky (2019) and is designed for analysis of right-truncated data.
What are right-truncated data?
When we are interested in the analysis of lifetimes that start at an initiating event and end up when an event of interest occurs, but we cannot obtain a population of interest to follow them prospectively, retrospective ascertainment of cases could be a good alternative. However, this might result in right-truncated survival data.
Right truncation often happens when we cannot identify the population of interest since their initiating events are latent (like the times of infections), and the population becomes "visible" only when another event, that can be more easily identified, happens. Here we assume that for those who are sampled, the times of initiating events are known or can be reliably estimated.
One example of such type of data is the data set on AIDS cases resulted from contaminated blood transfusion [@Kalbfeisch]. Here, the time of interest is the incubation period of AIDS, i.e., the time from HIV infection to development of AIDS. Figure 1 shows schematically how these data were collected. The individuals who developed AIDS by June 30, 1986 (the date when the AIDS cases were retrospectively ascertained) were selected, satisfying by this the condition $T_i< R_i$, where $T_i$ is an incubation period for the $i^{th}$ person and $R_i$ is his/her truncation time, defined as the time from HIV infection to June 30, 1986.
In the two functions of the package, coxph.RT and coxph.RT.a0, we assume that $T_i$ are observed exactly and there is no censoring. The analysis is based only on sampled events. Those who belong to the population but have not developed AIDS by the time of collection of data remain unobserved.
knitr::include_graphics('slide1_116.png')
The sampling condition $T_i< R_i$ means that we tend to sample shorter lifetimes than longer ones. And an important complication of such sampling scheme is that the follow-up is limited: on Figure 1, $r_$, defined as the time from the earliest HIV infection (January 1983) to June 30 1986, is the maximum follow-up time in this data set. This is an important statistic, since it allows us to realize whether our sampling covers the whole range of possible incubation periods or not. If the follow-up is too short, most likely we do not capture longer incubation periods, and they will not only be underrepresented in our sample, but even worse, they are likely to be non-represented at all. We shall call the condition $P(T>r_)=0$ a positivity assumption.
For example in the AIDS data, positivity does not hold since the incubation period of AIDS can be longer than the maximum observed follow-up time in this data set that equals 3.5 years. Under such short retrospective follow-up there was no chance to observe longer incubation times, which means that $P(T>r_)>0$, or equivalently, the support of observed $T^$ is shorter than the support of the original $T$ in the population.
Right truncation causes distortion of the original distribution
Even if positivity holds, truncation anyway results in distorted distribution, and it is important to account for this bias.
For example, on Figure 2 we illustrate the problem of distortion of the original distribution using simulated data. Here, the population lifetime $T$ follows Weibull distribution as rweibull(20000, shape=3, scale=1).
And the population truncation time $R$ was generated (independently of $T$) from rweibull(20000, shape=2.5, scale=1).
The sampled lifetimes $T^$ comprise all observations from the population that satisfy $T_i
knitr::include_graphics('densTpTsRp.png')
Why is it important to adjust for right truncation in regression?
Even if we do not aim to estimate the whole distribution and are interested only in estimation of the effects of some covariates on the incubation period, we still need to adjust for right truncation, since, as we illustrate below, the right truncation blurs the difference between subgroups.
Suppose we have one binary covariate $Z$. Figure 3 shows distributions in the two subpopulations defined by $Z$. The lifetime $T$ is generated from the proportional hazard regression model $h(t;z)=h_0(t)\exp(\beta z)$ with $\beta=2$:
library(ggplot2) library(survival) set.seed(17) N <- 20000 Z <- as.numeric(runif(N,0,1) < 0.5) # binary explanatory variable X <-rweibull(N, shape=3, scale=1*exp(-(2*Z)/3)) # shape a=3, scale b=1, beta=2 m <- max(X) d.tr <- data.frame(time=c(X), variable=c(Z)) d.tr$variable=as.factor(d.tr$variable) levels(d.tr$variable) <- c( "T | Z=0", "T | Z=1", " ") p <- ggplot(d.tr) + geom_density(aes(x=time, color=variable, fill=variable), alpha=0.5) + scale_color_brewer(palette = "Set1")+ scale_fill_brewer(palette = "Set1")+ theme_classic() + labs(x = "time")+ ggtitle("In the population")+ theme(axis.text.x = element_text(face="bold"), axis.text.y = element_text(face="bold"), plot.title = element_text(hjust = 0.5) )+ xlim(0, m) p
Figure 4 shows sample distributions under right truncation with positivity. For this illustration we use the same population as in Figure 3, but truncate it using uniformly distributed $R$:
set.seed(0) r <- runif(N, 0, m) ind <- X<r cat("P(T<R) = selection probability = ", mean(ind), "\n") d.tr <- data.frame(time=c(X[ind],r), variable=c(Z[ind], rep(3, N))) d.tr$variable=as.factor(d.tr$variable) levels(d.tr$variable) <- c( "T* | Z*=0", "T* | Z*=1", "R")
p <- ggplot(d.tr) + geom_density(aes(x=time, color=variable, fill=variable), alpha=0.5) + scale_color_brewer(palette = "Set1")+ scale_fill_brewer(palette = "Set1")+ theme_classic() + labs(x = "time")+ ggtitle(paste("In the sample with R~U[0,", round(m,1), "]", sep=""))+ theme(axis.text.x = element_text(face="bold"), axis.text.y = element_text(face="bold"), plot.title = element_text(hjust = 0.5) ) + xlim(0, m) #+ ylim(0,6.3) p
Despite the 30% chance of truncation, both Figures 3 and 4 look relatively similar to each other.
On Figure 5, we shorten the follow-up time by using $r_*=1.1$. This leads to violation of positivity, since longer lifetimes have no chance to be sampled, and the sub-groups become more similar to each other. In such case, without any adjustment, the estimators of $\beta$ will be biased and it is likely we will not able to discover the difference between subgroups.
set.seed(0) r <- runif(N, 0, m/2) ind <- X<r d.tr <- data.frame(time=c(X[ind],r), variable=c(Z[ind], rep(3, N))) d.tr$variable=as.factor(d.tr$variable) levels(d.tr$variable) <- c( "T* | Z*=0", "T* | Z*=1", "R") p <- ggplot(d.tr) + geom_density(aes(x=time, color=variable, fill=variable), alpha=0.5) + scale_color_brewer(palette = "Set1")+ scale_fill_brewer(palette = "Set1")+ theme_classic() + labs(x = "time")+ ggtitle(paste("In the sample with R~U[0,", round(m/2,1), "]", sep=""))+ theme(axis.text.x = element_text(face="bold"), axis.text.y = element_text(face="bold"), plot.title = element_text(hjust = 0.5) ) + xlim(0, m) #+ ylim(0,6.3) p
On Figure 6, severe violation of positivity with $R\sim U(0,0.4)$ causes the groups to be almost identical:
set.seed(0) r <- runif(N, 0, m*4/20) ind <- X<r d.tr <- data.frame(time=c(X[ind],r), variable=c(Z[ind], rep(3, N))) d.tr$variable=as.factor(d.tr$variable) levels(d.tr$variable) <- c( "T* | Z*=0", "T* | Z*=1", "R") p <- ggplot(d.tr) + geom_density(aes(x=time, color=variable, fill=variable), alpha=0.5) + scale_color_brewer(palette = "Set1")+ scale_fill_brewer(palette = "Set1")+ theme_classic() + labs(x = "time")+ ggtitle(paste("In the sample with R~U[0,", round(m*4/20,1), "]", sep=""))+ theme(axis.text.x = element_text(face="bold"), axis.text.y = element_text(face="bold"), plot.title = element_text(hjust = 0.5) ) + xlim(0, m) #+ ylim(0,6.3) p
And the situation becomes more complicated when there are several dependent covariates, as described in @RT.
What does coxrt suggest for right-truncated data?
The coxrt package allows to estimate regression coefficients $\beta$ in the proportional hazards model $h(t;z)=h_0(t)\exp(\beta z)$.
Under positivity use coxph.RT
Function coxph.RT implements estimating equations denoted as IPW-S in @RT. This function assumes positivity (as, e.g., on Figure 4).
For example, for AIDS data we can fit proportional hazards regression of incubation time on age:
library(gss) library(coxrt) data(aids) all <- data.frame(age=aids$age, ageg=as.numeric(aids$age<=59), T=aids$incu, R=aids$infe, hiv.mon =102-aids$infe) all$T[all$T==0] <- 0.5 # as in Kalbfeisch and Lawless (1989) s <- all[all$hiv.mon>60,] # select those who were infected in 1983 or later # analysis assuming positivity # we request bootstrap estimate of Asymptotic Standard Error (ASE) as well: sol <- try(coxph.RT(T~ageg, right=R, data=s, bs=TRUE, nbs.rep=500) ) knitr::kable(sol$summary) # print the summary of fit based on the analytic ASE knitr::kable(sol$bs$summary) # print the summary of fit based on the bootstrap sample distribution
Technically, for point estimation coxph.RT uses the standard coxph function from the survival package. However, the estimator of the standard error (SE) obtained from coxph even with robust option might underestimate for smaller samples (of $<1000$ observations).
Although coxph.RT includes both analytic and bootstrap estimates of a SE, the inference based on the analytic SE estimate, implemented in coxph.RT, is more reliable than based on its bootstrap counterpart, since bootstrap resampling does not always capture the variability in data that comes from near-zero sampling probabilities $P(R>T_i)$.
Since the Inverse-Probability-Weighting (IPW) approach, implemented in coxrt, relies on "blowing up" the sampled lifetimes $T_i$ in order to obtain the unbiased pseudo-population, which should recreate the same probability law as in the original population, coxph.RT in general will give biased estimators in case there is an unobservable right tail of the distribution. Hence the function coxph.RT.a0 aims to adjust for this violation of positivity by using the structural assumption of the Cox model in the population and the information on $a_0=P(T>r_*\mid Z=0)$.
Without positivity use coxph.RT.a0
For the situations without positivity (as in Figures 5 and 6), if $a_0$ can be hypothesized, based on the subject matter knowledge, or if $a_0$ is unknown but its interval can be guessed, this information can be used in coxph.RT.a0. This function implements estimating equations denoted as IPW-SA in the original paper. IPW-SA uses stabilized weights. Technically, coxph.RT.a0 solves non-linear equations using BBsolve function from BB package using default parameter settings.
For example, in order to estimate $\beta$, the effect of age at HIV on the incubation period, given that $a_0=0.2$, we run:
# analysis using adjusted estimating equations for a0=0.2 sol.02 <- try(coxph.RT.a0(T~ageg, right=R, data=s, a0=0.2, bs=TRUE, nbs.rep = 500)) knitr::kable(round(sol.02$bs$summary,2))
We can do a sensitivity analysis for $\hat{\beta}$ for different values of a sensitivity parameter $a_0$:
# for a0=0: sol <- try(coxph.RT(T~ageg, right=R, data=s, bs=FALSE) ) # senstivity analysis for different values of a0 a_ <- seq(0.05, 0.55, by=0.05) est <- CI.L <- CI.U <- NULL start_time <- Sys.time() for(q in 1:length(a_)) { sol.a <- try(coxph.RT.a0(T~ageg, right=R, data=s, a0=a_[q], bs=TRUE)) if (sol.a$convergence!=0) { cat("a0=", a_[q], ". Error occurred in BBsolve.\n") } else { est <- c(est, sol.a$coef) CI.L <- c(CI.L, sol.a$bs$summary$CI.L) CI.U <- c(CI.U, sol.a$bs$summary$CI.U) } } end_time <- Sys.time() end_time - start_time # It takes some time to solve nonlinear estimating equations and to run bootstrap for each value of a0 ...
require(ggplot2) res.d <- data.frame(a0=c(0, a_), beta=c(sol$coef, est), L.95=c(sol$summary$CI.L, CI.L), U.95=c(sol$summary$CI.U, CI.U)) knitr::kable(round(res.d, 2)) p <- ggplot(res.d, aes(x=a0, y=beta)) + geom_ribbon(aes(ymin=L.95, ymax=U.95), alpha=0.2) + geom_line() + geom_point() + geom_hline(yintercept=0) p + xlab(expression( paste(a[0], "=P(T>", r['*']," | z=0)" , sep="")) )+ ylab(expression( paste(hat(beta), "(", a[0], ")" , sep="")) ) + scale_x_continuous(breaks=res.d$a0, labels=res.d$a0) + theme(axis.text.x = element_text(face="bold", angle=45), axis.text.y = element_text(face="bold"))
In coxph.RT.a0, only the bootstrap estimator of SE is available.
Some remarks
-
Positivity cannot be tested. This assumption is based on the subject matter knowledge.
-
The approach assumes independence between $T$ and ($R$, $Z$) in the original population. If there is dependence between $T$ and $R$, it can be discovered by Kendall's tau test which can be done using packages tranSurv or SurvTrunc, or by more powerful permutation tests which can be done using package permDep.
-
If a non-terminal event of interest cannot be observed because death preceded it, that is the event of interest is truncated by death, the approach implemented here can be used only if two lifetimes, time to an intermediate event (=$T$) and time to death (=$R$), are independent, otherwise the results might be biased.
Some limitations of the coxrt package
-
Although in practice, it is possible to have additional left or interval censoring in right-truncated data, the package assumes that the exact lifetimes, without any censoring, are observed.
-
The current implementation allows for only time-independent covariates.
-
The package estimates only covariate effects $\beta$, and does not provide estimates for the baseline hazard $h_0(t)$ in the Cox model.
References
Try the coxrt package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.