Announcements

1. Announcements

2. Stat’s Seminars UpdateRest of the Fall • Nov 4: Competing Risks, Frailty, and Violations of Proportional Hazards Tom Greene • Nov 11: Alternatives to Cox regression Tom Greene & Nan Hu • Dec 2: Chasing Tail(s): Statistical Intuition for High-Throughput Studies Ken Rice, University of Washington • Dec 9: Case Control Study Design Carroll Sweeney • Dec 16: Analyses for case-control studies Greg Stoddard Statistical Methods in Progress Meetings: (also in WB – usually 238A or B) Charlie Casper: Nov 8, 10 AM Mary Slattery: Nov 22, 10 AM Andrew Redd: Dec 6, 10 AM Causal Modeling Workgroup: 3rd Friday of Each Month (3PM, Div. EPI conference room, WB)

3. Other Announcements • Slides sets for seminars to be archived at http://www.ccts.utah.edu/biostats/ • R User’s group forming • Andrew Redd

4. Competing Risks, Frailty, and Violations of Proportional Hazards Tom Greene

5. References • Aalen O, Borgan O, Gjessing H. Survival and Event History Analysis: A Process Point of View, Springer, 2008 • Martinussen T, Scheike T. Dynamic Regression Models for Survival Data. Springer, 2010

6. Issues for Today • Understanding what standard Kaplan-Meier curves and competing risk cumulative incidence curves really estimate in competing risk situations • Same issue for standard hazard ratios and competing risk hazard ratios for subdistributions • When to use standard methods or competing risk methods?

7. Issues for Today • What is the difference between individual-level hazard ratios and population hazard ratios? • What are the implications of variation in frailty for population hazard ratios in Cox regression ? • What are the implications of variation in frailty for standard survival analysis in clinical trials? • What if there are also competing risks?

8. Questions for Today • How can one check for non-proportional hazards in Cox regression? • How can Cox regression be modified to account for non-proportional hazards?

9. Quantities in Survival Analysis S(t) = Pr(T > t) for a survival time T T may be subject to right censoring time C: we observe X = min(T,C) The hazard rate λ(t) = event rate in time interval (t, t+ε) for the subgroup still at risk at time t The cumulative hazard H(t) is the area under the λ(u) vs. u curve up to time t: Relationship of survival and cumulative hazard:

10. Estimating S(t) and H(t) dj = # events at time tj; nj= # at risk at time tj Kaplan-Meier estimator: Nelson-Aalen estimator (discrete version):

11. Cox Model λ(t|Z) = λ0(t) × exp(βZ) • Z : Covariate(s) of interest • λ0 (t): Baseline hazard function when Z = 0 • β: Cox regression coefficients • exp(β): Hazard ratios under the Cox model • Two key assumptions: • T and C are independent given Z • Proportional hazards

12. Competing Risks Events • An event whose occurrence precludes the occurrence of an event of interest • Usually a clinical event, possibly having a direct or indirect relationship with the event being analyzed

13. Competing Risk Subdistribution • The subdistribution, or cumulative incidence probability, corresponding to the fraction of all patients with the event of interest, after accounting for “true censoring” • Estimated by • tj: Event times (all event types) • dj: Number of pts with event of interest at tj • nj:Number at risk at tj • : Kaplan-Meier Curve for the composite outcome T = min(T1,T2)

14. Traditional Kaplan Meier and Competing Risk Cumulative Incidence Pts > 55 Years of Age Example: Study of chronic kidney disease; Death and End Stage Renal Disease (ESRD) are competing Risks Traditional Kaplan-Meier, ESRD censoring death Cumulative Incidence of ESRD, Death as competing risk Traditional Kaplan-Meier, Death censoring ESRD Cumulative Incidence, Death as competing risk Months Follow-up

15. Interpretation of Kaplan Meier curves with competing risks • The Kaplan Meier cumulative incidence curve indicates the event probability for a “virtual population” in which the competing risk does not occur, but which is linked to the studied population such that the target event occurs after a competing event at the same rate as for patients not experiencing the competing event. • This virtual population is hard to fathom if the competing risk and target event are associated

16. Competing Risk Cumulative Incidence Curves for Time to ESRD or Death prior to ESRD ESRD or Death ESRD or Death Age < 55 years Age > 55 years ESRD ESRD Death Death Months Follow-up Months Follow-up

17. Competing Risk Formulation Decomposes Cohort Into 3 or More States Age > 55 years Alive without ESRD Alive without ESRD Age < 55 years Alive with ESRD prior to Death Alive with ESRD prior to Death Death prior to ESRD Death prior to ESRD Months Follow-up Months Follow-up

18. Conventional and Competing Risk Hazard Calculation Hazard Increment at time tj nj(at risk) Competing Risk Calculation nj(cmp.rsk) nj(event) nj(at risk) Conventional Risk Calculation nj(cmp.rsk) nj(event) tj tj- Time

19. Analysis of Subdistriubtion Hazard • Typically implemented using a proportional model (Grey’s model), analogous to proportional Hazards Cox regression • Implemented in R library cuminc • Can be viewed as assigning + infinity to survival times for those with competing risks

20. What to do in practice? Example: RCT of effect of BP drug on CV composite; non-CV death is competing risk. The trial’s focus is on the CV outcome, but the treatment may also affect non-CV death. • Useful to display competing risk cumulative incidence • The standard analysis (using log rank or Cox) gives the correct HR in the hypothetical “independent competing risk” universe • This remains true even if the treatment affects the competing risk • This is not true of the HR under Gray’s competing risk model • Treatment can appear beneficial by increasing the competing risk even if it has no effect on the target endpoint • Hence, the “standard approach” is typically used to compare hazard rates in RCTs

21. A Parable for Statisticians A man is walking home one dark night when he sees a sage searching for something frantically under a streetlight. “Teacher, what have you lost?” the man asks. “I am searching for my key,” the sage says worriedly. “I’ll help you look,” the man says and joins the search. Soon both are down on their knees under the streetlight, looking for the lost key. After some time, the man asks, “Tell me, teacher, do you remember where exactly you dropped the key?” The sage waves his arm back toward the darkness and says, “Over there, in my house. I lost the key inside my house…” Shocked and exasperated, the man shouts, “Then why are you searching for the key out here in the street?” “Because there is more light here than inside my house,” the teacher answers non-chalantly.

22. Population and Individual Level Hazards • The hazard in the Cox model λ(t|Z) = λ0 (t) × exp(βZ) represents the “population hazard” given Z. • The population hazard represents an average of individual level hazards among patients remaining at-risk • The variation in individual level hazards, conditional on Z, is captured by the concept of frailty.

23. Proportional Frailty Model • Assume the hazard function for an subject with frailty W is α(t|W) = W × α(t). • As before: • The population hazard, λ(t) is related to the standardized individual hazard α(t) by: Standardized individual hazard, when W = 1

24. Proportional Frailty Model • Special Case: Assume W is gamma distributed with mean 1 and variance δ • Then Relationship of λ(t) with α(t) if α(t) = t2 δ=0 δ=0.4 δ=1 δ=3 Time

25. Relationship of Population and Individual Hazard Ratios • Suppose the proportional frailty model holds in two groups, with gamma frailties with mean 1 and variance δ within each group. • If the individual hazards in 2 groups are W1rα(t) and W2α(t), so the individual HR is r, then Attenuates towards 1 as t increases

26. Attenuation in Population HR Due to Frailty Variance δ=0 Assumes r = 4, α(t) = 0.5 (i.e., ~40% with events in 1 yr) δ=0.40 δ=1.00 δ=3.00

27. Impact of Immediate Attenuation of Patient-Level HR to 1 Under Frailty Variation • Assumes • δ=1 • r(t)=0.5 for • first year, • r(t) = 1 after • 1 year • 3) α(t) = 1

28. Competing Risks and False ProtectivityDifferential Survival Bias WBαB(t) Y0, Y2, Y3 independent gammas with mean 1 and variances δ0, δ2, and δ3 Event B WB= Y0 + Y2 Alive without Event WC= Y0 + Y3 WCαC(t) Event C Then A treatment which has no effect on αB(t) but increases αC(t) will have lower λB, and thus appear protective for event B.

29. Practical Implications • Variation in population hazards ratio over time depends both on variation in individual hazards ratio and to frailty selection • It is often found that HRs attenuate towards 1 or reverse over time, or among “survivors” who reach an advance stage of a chronic disease (e.g., those with end stage renal disease)

30. Relative Mortality for Norwegian Men 1901-1905 compared to 1991 Horizontal line indicates relative risk of 1 From “Statistics Norway” as reproduced in Aalen O, Survival and Event History Analysis

31. Practical Implications • Prevailing practice is to avoid covariate adjustment for survival outcomes in RCTs. • But adjustment for strong prognostic factors can: • Reduce conservative bias in estimated treatment effect • Increase power • Reduce differential survival bias

32. Practical Implications • Three approaches to estimating effects on individual hazards: • Analysis of repeat event data • Joint analysis of longitudinal and time-to-event outcomes • Use background knowledge to posit assumptions for frailty distribution and/or on the mathematical form of α(t)

33. Evaluation of Proportional Hazards • Parametric models for change in hazard ratios over time • Non-parametric smooths of Schoenfeld residuals • Non-parametric models for multiplicative hazards

34. Parametric models for change in hazard ratios over time Example: • Cox Regression of Effects of Dose Group (Ktv_grp) and baseline serum albumin (Balb) in the HEMO Study • RCT with 871 deaths in 1871 patients; planned follow-up 1.5 to 7 years.

35. Parametric models for change in hazard ratios over time 1) Test for linear interactions of predictors with follow-up time proc phreg data=demsum01 ; model fu_yr * ev_d(0) = ktv_grpbalbktv_grptbalbt; ktv_grpt = ktv_grp*fu_yr; balbt = balb*fu_yr; Parameter Standard Parameter DF Estimate Error Chi-Square Pr > ChiSq KTV_GRP 1 -0.06933 0.12112 0.3277 0.5670 BALB 1 -1.50469 0.17279 75.8320 <.0001 Ktv_grpt 1 0.00686 0.04406 0.0243 0.8762 Balbt 1 0.16874 0.06378 7.0001 0.0082 HR for baseline albumin (in g/dL) is 0.22 at time 0, but attenuates by a factor of exp(0.1687) = 1.184 per year

36. Parametric models for change in hazard ratios over time 2) Test for interactions of predictors with time period (> 1 yr vs. < 1yr) proc phreg data=demsum01 ; model fu_yr * ev_d(0) = Ktv_grp1 Balb1 Ktv_grp2 Balb2; if fu_yr > 1 then period =1; if . < fu_yr <= 1 then period = 0; Ktv_grp1 = ktv_grp*(1-period); Balb1 = balb*(1-period); Ktv_grp2 = ktv_grp*period; Balb2 = balb*period; PropHazKtv: test Ktv_grp1=Ktv_grp2; PropHazBalb: test Balb1 = Balb2;

37. Parametric models for change in hazard ratios over time 2) Test for interactions of predictors with time period (> 1 yr vs. < 1yr) Parameter Standard Parameter DF Estimate Error Chi-Square Pr > ChiSq Ktv_grp1 1 -0.07569 0.13413 0.3185 0.5725 Balb1 1 -1.58168 0.18855 70.3661 <.0001 Ktv_grp2 1 -0.04680 0.07861 0.3545 0.5516 Balb2 1 -0.95959 0.11553 68.9860 <.0001 Wald Label Chi-Square DF Pr > ChiSq PropHazKtv 0.0345 1 0.8526 PropHazBalb 7.9139 1 0.0049

38. Plots of Schoenfeld Residuals R Code: require(survival) HEMOCox<-coxph(Surv(fu_yr,EV_D) ~ KTV_GRP+BALB,data=hemodat) HEMOPropchk<-cox.zph(HEMOCox) HEMOPropchk plot(HEMOPropchk,var="KTV_GRP") plot(HEMOPropchk,var="BALB") plot(HEMOPropchk,var=“KTV_GRP",resid=FALSE) plot(HEMOPropchk,var="BALB",resid=FALSE) rho chisq p KTV_GRP 0.0049 0.0209 0.88510 BALB 0.0968 8.2343 0.00411 GLOBAL NA 8.2570 0.01611

39. Schoenfeld Residual Plots with Cubic Spline Smooths: R Output For KTV_GRP For Baseline Albumin

40. Schoenfeld Residual Plots with Cubic Spline Smooths: R Output (Omitting the residuals) For KTV_GRP For Baseline Albumin

41. Multiplicative Hazards with Time Varying Coefficients • Standard Cox Proportional Hazards Model λ(t|Z) = λ0 (t) × exp(β Z) • Cox Proportional Hazards Model with Time-Dependent Covariates λ(t|Z) = λ0 (t) × exp(β Z(t)) • Multiplicative Hazards Model with Fixed Covariates and Time-Varying Coefficients λ(t|Z) = λ0 (t) × exp(β(t) Z) • Multiplicative Hazards Model with Time-Dependent Covariates and Time-Varying Coefficients λ(t|Z) = λ0 (t) × exp(β(t) Z(t))

42. Multiplicative Hazards Model with Fixed Covariates and Time-Varying Coefficients • Model λ(t|Z) = λ0 (t) × exp(β(t) Z) • Useful if proportional hazards assumption in doubt, and you don’t want to assume a particular parametric model for change in HR over time • Estimands are cumulative Cox regression coefficients • Can use timereg package in R (if you are careful to center predictor variables)

43. Multiplicative Hazards Model with Fixed Covariates and Time-Varying Coefficients summary(fit) Multiplicative Hazard Model Test for nonparametric terms Test for non-significant effects Supremum-test of significance p-value H_0: B(t)=0 cKTV_GRP 1.38 0.931 cBALB 11.30 0.000 Test for time invariant effects Kolmogorov-Smirnov test p-value H_0:constant effect cKTV_GRP 0.163 0.986 cBALB 0.833 0.038 Cramer von Mises test p-value H_0:constant effect cKTV_GRP 0.015 0.993 cBALB 1.040 0.026

44. Multiplicative Hazards Model with Fixed Covariates and Time-Varying Coefficients R Code: > cBALB<-BALB – mean(BALB) > cKTV_GRP <- KTV_GRP – mean(KTV_GRP) > require(timereg) > fit<-timecox(Surv(fu_yr,EV_D)~KTV_GRP+cBALB,max.time=5) > summary(fit) > par(mfrow(1,2) > plot(fit,c(2,3))

45. Multiplicative Hazards Model with Fixed Covariates and Time-Varying Coefficients

46. Next Week Alternatives to the Cox Model • Aalen’s Additive Hazards Model • Accelerated Failure Time Models • Proportional Odds Models for Survival Data