How to compare clinical trial data to a natural history control

Question

Context: I am interested in understanding how to analyze data from clinical trials with a natural history control, meaning, a clinical trial where a group of people (say, those with genetic risk for a disease) are recruited, all of them are treated with a drug, and their survival (say, age at disease onset, or death due to this disease) is compared to the survival of untreated historical controls with the same genetic risk.

Before I start, I just want to acknowledge that I've read FDA Guidance documents (e.g. Rare Diseases: Common Issues in Drug Development) touching on this topic and I'm well aware that there are many types of bias that can creep into natural history comparisons and that, for this reason, FDA only rarely accepts such trials as evidence of a drug's efficacy, instead usually considering randomized trials as far stronger evidence. The question of how to determine whether a natural history cohort can fairly be compared to a prospective cohort treated with a drug is an important and complex question. But that is not the topic of my question today. Instead, I am finding that I am stumped on a yet much simpler question, the question of if a natural history cohort is considered comparable and unbiased and so on, then how, statistically speaking, would one do the comparison?

So here's a scenario.

• For simplicity, say my disease has a single, known, genetic cause, and people with this genotype are perfectly healthy until some age, and then they suddenly become very sick. The age of disease onset is highly variable, and hazard varies as a function of age. Everyone in every dataset described below has this disease-causing genotype.
• In dataset A, I have data on people in a hypothetical clinical trial. They enrolled at one age (starting_age), were treated with a drug, and were followed a variable amount of time until a second age (last_age) at which they either became sick (event == 1) or withdrew from the trial (event == 0).
• In an "ideal ideal" scenario, perhaps the natural history cohort would be a large number of people with this genotype, followed from birth until they either died of this disease or of an unrelated cause. That sort of dataset doesn't exist. Instead, consider two possible options for what the natural history cohort could be (B and C)...
• In dataset B, I have data on people who enrolled in a prospective study at some random age (starting_age), were not treated with any drug, and were simply followed for some variable amount of time until a second age (last_age) at which they either became sick (event == 1) or withdrew from the study (event == 0). The distribution of starting ages and the number of years of followup cannot be assumed to follow the same distributions as in dataset A.
• In dataset C, I have purely retrospective data; here there was no prospective followup. We simply observed each person once. Some people were observed to become sick at some particular age (age, and event == 1) while others were seen to still be healthy at some particular age (age, and event == 0).

If I understand the terminology correctly, the data in A and B would be considered left-truncated and right-censored; the data in C would be considered right-censored only. According to terms used in Cain 2011, I believe that dataset B is a "prevalent cohort" while dataset C is an "incident cohort".

My questions are:

1. How (if at all) can one construct a survival function (or a hazard function, etc) for dataset A or B, given that the data are left-truncated at different ages for each individual?
2. What statistical test (or bootstrapping method) would be used to test whether people in dataset A survive significantly longer or have significantly lower hazard than those in B or C?
3. If all you had was dataset B or C, and you wanted to estimate the statistical power of a given clinical trial with N individuals for a given hazard ratio, how would you go about it?
4. (Bonus question) if people also had, say, different genetic mutations associated with different hazards as a function of age, how would you incorporate that variable into the answers above?

Some things I've tried so far:

• I considered Flora's Z statistic [Flora 1978] which some have applied in similar situations, but I am concerned that it doesn't account for sampling variance in the natural history cohort being used as a reference.
• I Googled a few different combinations of terms such as survival analysis with different left truncation times and found some pages that discuss similar problems, such as Survival Analysis: Left-Truncated Data, Surviving Left Truncation using PROC PHREG, but these did not explain the underlying math of their solution.
• I finally came upon Cain 2011 which discusses the issue in detail and, helpfully, has R code for handling left-truncation in the supplement. They implemented their own function for an MLE incorporating left-truncation, but the claim is that left-truncation can already be handled in a Cox proportional hazards model using functionality built into the R survival library, for instance: coxph(Surv(time=agestart,time2=x,event=cx,type='counting')~bmi). Here a Surv object is created using time as the starting age, and time2 as the age where event either did or didn't happen. That sounds sensible, though the help file for Surv doesn't explain what it is actually doing here — for example it does not explain what the counting model entails. I turned to the survival manual and, as a reference for counting, got pointed to Andersen & Gill 1982 which in turn was a bit over my head and does not seem to discuss left truncation, or at least not by that name. One peculiarity is the "counting" model seems built to fit cases where an individual can have multiple events in their lifetime, which is not the case in my example, but perhaps that doesn't matter.
• Based on the above, it seemed like a Cox "counting" model might be the right way to compare dataset A to B, and I was able to do such a comparison in R (see code below) though I am still struggling to understand whether I am doing the right thing here, and I am not sure how (if at all) such an approach could be applied to compare dataset A to C.
• Finally, I searched Cross Validated for survival left truncation and found a large number of instances where people had asked similar questions as I'm asking (1, 2, 3, 4, 4, 5, 6, 7) but most were unanswered; one pointed to the Cain reference above, and another pointed to Klein & Moeschberger 2003 which was helpful (see e.g. p. 123 and p. 312) and seems to support the notion that A and B can be compared using a fairly simple Cox approach, though it doesn't address the comparability of A to C, nor the power calculation question (though, if the statistical test question is answered, I could presumably get to the power with some bootstrapping).

Below are some hypothetical data in R illustrating this scenario. I've included 20 rows for each dataset, though in case it matters, in the real life scenario I am imagining, datasets A and B could be perhaps on the order of 50 or 100 patients, and dataset C might be on the order of 500 or 1,000 patients.

# A) hypothetical data from clinical trial
data_a = read.table(sep='|',header=T,textConnection("
indiv_id|starting_age|last_age|event
1|33|42|0
2|45|49|0
3|47|52|1
4|30|34|0
5|37|44|0
6|34|37|0
7|29|34|0
8|58|66|0
9|58|60|0
10|66|75|0
11|37|41|0
12|37|46|0
13|58|62|0
14|44|48|1
15|45|50|0
16|56|65|0
17|54|63|0
18|36|41|0
19|47|55|1
20|45|55|0
"))

# B) hypothetical data from a prospective natural history study
data_b = read.table(sep='|',header=T,textConnection("
indiv_id|starting_age|last_age|event
101|19|28|0
102|39|52|0
103|38|41|1
104|18|27|0
105|20|24|0
106|16|20|0
107|39|41|0
108|48|50|0
109|40|50|0
110|38|41|1
111|40|43|1
112|26|29|0
113|37|39|0
114|21|30|0
115|36|41|1
116|46|48|0
117|27|32|0
118|26|29|0
119|29|38|0
120|47|58|0
"))


# C) hypothetical data from a retrospective natural history study
data_c = read.table(sep='|',header=T,textConnection("
indiv_id|age|event
201|43|1
202|53|1
203|64|1
204|45|1
205|88|1
206|70|1
207|66|1
208|55|1
209|51|1
210|48|1
211|63|1
212|36|0
213|61|0
214|63|1
215|63|1
216|57|1
217|74|0
218|63|1
219|59|1
220|57|1
"))


# one possible approach to compare A and B using Cox counting model
data_a$drug = TRUE
data_b$drug = FALSE
nh_compare = rbind(data_a, data_b)
m = survfit(Surv(time=starting_age,time2=last_age,event=event,type='counting')~drug, data=nh_compare)
summary(m)
coxph(Surv(time=starting_age,time2=last_age,event=event,type='counting')~drug, data=nh_compare)

Final note: because a reputation of 10 is required to post >2 links on Cross Validated, I have removed all of the (many) hyperlinks that should appearabove. A version of this post with links is available on my blog and I will update this post to include the links if or when I am permitted to do so. Update: thanks to everyone who upvoted! I now have reputation >10 and have updated this post to include links.

Answer 1

Years later, I've arrived at a satisfactory answer. This did indeed turn out to be to use a Cox proportional hazards counting model, which allows you to account for different left-truncation times (ages at which you started following the individuals) in addition to different right-censoring times. As noted in the question, this is implemented in R in the survival package coxph function, where time is the left-truncation time, time2 is the right-censoring time, event is what happened at the right-hand time, and you use type='counting' to specify the Cox counting model.

The answers to specific questions raised in my post are:

1. Survival can indeed be computed with survfit — the model accounts for how the number of at-risk individuals can both grow and shrink over time as people enter and exit the age ranges where they were followed. An example for the toy dataset posted above would be plot(survfit(Surv(time=starting_age,time2=last_age,event=event,type='counting')~1, data=data_a))

2. The model can compare the left-truncated prospective data and the non-left-truncated retrospective data, if you simply assume that the retrospective data are equivalent to following people from birth. This assumption may not be perfect but this is an inherent limitation of the dataset that no model is going to get around. Example code for the toy dataset above would be:

data_c$starting_age = 0
data_c$last_age = data_c$age
data_c$drug = FALSE
nh_compare = rbind(data_a, data_c[,c('indiv_id','starting_age','last_age','event','drug')])
coxph(Surv(time=starting_age,time2=last_age,event=event,type='counting')~drug, data=nh_compare)

3. There appears to be no closed-form power calculation, instead we did this by bootstrapping. Our code to do this for our specific dataset appears here.

4. coxph allows for covariates, so for example in our code we used coxph(Surv(time=ascertainment_age,time2=surv_age,event=surv_status,type='counting')~asc+family_mutation,data=prore) where family_mutation is a covariate.

We have published a paper where used this approach to calculate power for preventive clinical trials in genetic prion disease. You can read the details on bioRxiv and our R code is all in a public GitHub repo:

https://github.com/ericminikel/prnp_onset/

Citation:

Minikel EV, Vallabh SM, Orseth MC, Brandel JP, Haïk S, Laplanche JL, Zerr I, Parchi P, Capellari S, Safar J, Kenny J, Fong JC, Takada LT, Ponto C, Hermann P, Knipper T, Stehmann C, Kitamoto T, Ae R, Hamaguchi T, Sanjo N, Tsukamoto T, Mizusawa H, Collins SJ, Chiesa R, Roiter I, de Pedro-Cuesta J, Calero M, Geschwind MD, Yamada M, Nakamura Y, Mead S. Age at onset in genetic prion disease and the design of preventive clinical trials. Neurology. 2019 Jun 6. pii: 10.1212/WNL.0000000000007745. doi: 10.1212/WNL.0000000000007745. PubMed PMID: 31171647.

Answer 2

Eric, broadly speaking, your problem sounds severe enough that a search for off-the-shelf solutions seems misguided. Rather, you almost surely need recourse to bespoke modeling to exploit your special domain knowledge about the pathophysiology of the disease. Unless you use a modeling approach that enables you to bring such knowledge to bear, you might not stand a chance against the formidible 'opponent' you are facing!

Your best first step might be figuring out just what 'special domain knowledge' you actually possess. Can you simulate the process that generated your data (i.e., the data-generating process or DGP), including the (left-truncation) process governing the entry of individuals into your data set? Once you can simulate the DGP, Bayesian methods should enable you to 'challenge' your simulation model with data—e.g., to estimate your model's parameters. Notwithstanding Odd Aalen's fin de siècle skepticism about Bayesian methods for survival analysis [1], I note that there is now at least one text on such approaches [2].

If I were faced with such a problem, I would be inclined first to explore it through simulation and Bayesian inference. Perhaps I would learn enough in that process to formulate simpler process models that might yield to more traditional frequentist estimation approaches. The interplay between simpler models and richer simulations might indeed yield its own valuable forms of insight and understanding.

I hope you will eventually update us all here on what approach you eventually adopt, and how it works out!

1. Aalen OO. Medical statistics - no time for complacency. Stat Methods Med Res. 2000;9(1):31-40. doi:10.1177/096228020000900105.

2. Ibrahim JG, Chen M-H, Sinha D. Bayesian Survival Analysis. New York: Springer; 2010.