5 Modeling Survival Data with Parametric Regression

Transcription

1 5 Modeling Survival Data with Parametric Regression Models 5. The Accelerated Failure Time Model Before talking about parametric regression models for survival data, let us introduce the accelerated failure time (AFT) Model. Denote by S (t) ands 2 (t) the survival functions of two populations. The AFT models says that there is a constant c>0 such that S (t) =S 2 (ct) for all t 0. (5.) This model implies that the aging rate of population is c times as much as that of population 2. (For example, if S (t) is the survival function for the dog population and S 2 (t) isthesurvival function for the human population, then the conventional wisdom that a year for a dog is equivalent to 7 years for a human implies c =7,andS (t) =S 2 (7t). So the probability that a dog can survive 0 years or beyond is the same as the probability that a human subject can survive 70 years or beyond) Let µ i be the mean survival time for population i and let ϕ i be the population quantiles such that S i (t)(ϕ i )=θ for some θ (0, ). Then µ 2 = = c = c S 2 (t)dt S 2 (cu)du (t = cu) S (u)du = cµ and S 2 (ϕ 2 )=θ = S (ϕ )=S 2 (cϕ ). Assume that S 2 (t) is a strictly decreasing function. Then we have ϕ 2 = cϕ. PAGE 00

2 This simple argument tells us that under the accelerated failure time model (5.), the expected survival time, median survival time of population 2 all are c times as much as those of population. Suppose we have a sample of size n from a target population. For subject i (i =, 2,..., n), we have observed values of covariates z i,z i2,..., z ip and possibly censored survival time T i.the procedure Proc Lifereg in SAS fits models to data specified by the following equations log(t i )=β 0 + β z i β p z ip + σε i, (5.2) where β 0,..., β p are the regression coefficients of interest, σ is a scale parameter and ε i are the random disturbance terms, usually assumed to be independent and identically distributed with some density function f(ε). The reason why we take logarithm of T i is obvious considering the fact that the survival times are always positive (with probability ). Equation (5.2) is very similar to a linear regression model for the log-transformed response variable Y i =log(t i ). In a linear regression, the random error term e i is usually assumed to be i.i.d. from N(0,σ 2 )sothate i can be written as e i = σε i,whereε i are i.i.d. from N(0, ) in this case. At this moment, let us see how the regression coefficients in model (5.2) can be interpreted in general. We will investigate their interpretation more closely later when we consider more specific models (i.e., with different distributional assumptions for ε i ). For this purpose, let us consider β k (k =,,p). Holding other covariate values fixed, let us increase covariate z k by one unit from z k to z k + and denote by T and T 2 the corresponding survival times for the two populations with covariate values z k and z k + (with other covariate values fixed). Then T and T 2 canbeexpressedas T = e β 0+β z +...+β k z k +...+β pz p e σε = c e σε T 2 = e β 0+β z +...+β k (z k +)+...+β pz p e σε 2 = c 2 e σε 2 where c 2 and c are two constants related by c 2 = c e β k. The corresponding survival functions PAGE 0

3 are S (t) = P [T t] =P [c e σε t] =P [e σε c t], S 2 (t) = P [T 2 t] =P [c 2 e σε 2 t] =P [e σε 2 c 2 t] Since ε and ε 2 have the same distribution, and c 2 = c e β k,wehave S 2 (e β k t)=p [e σε 2 c 2 e β k t]=p [e σε 2 c e β k e β k t]=p [e σε 2 t] =P [e σε t] =S (t). Therefore, we have accelerated failure time model between populations (covariate value=z k ) and 2 (covariate value=z k +)withc =e β k. So if we increase the covariate value of zk by one unit while holding other covariate values unchanged, the corresponding average survival time µ 2 and µ will be related by µ 2 =e β k µ. If β k is small, then µ 2 µ µ =e β k β k. Similarly we have for the population quantiles ϕ i ϕ 2 ϕ ϕ =e β k β k. Therefore, when β k is small, it can interpreted as the percentage increase if β k > 0 or percentage decrease if β k < 0 in the average survival time and/or median survival time when we increase the covariate value of z k by one unit. Thus the greater value of the covariate with positive β k is more beneficial in improving survival time for the target population. This interpretation of β k is very similar to that in a linear regression model. 5.2 Some Popular AFT Models We can assume different distributions for the disturbance term ε i in model (5.2). For example, we can assume ε i i.i.d. N(0, ). This assumption is equivalent to assuming that T i has log-normal PAGE 02

4 distribution (of course, conditional on the covariates z s). In this section, we will introduce some popular parametric models for T i (equivalently for ε i ). The following table gives some of these distributions: Distribution of ε Distribution of T Syntax in Proc Lifereg extreme values (2 par.) Weibull dist = weibull extreme values ( par.) exponential dist = exponential log-gamma gamma dist = gamma logistic log-logistic dist = llogistic normal log-normal dist = lnormal In Proc Lifereg of SAS, all models are named for the distribution of T rather than the distribution of ε. Although these above models fitted by Proc Lifereg all are AFT models (so the regression coefficients have a unified interpretation), different distributions assume different shapes for the hazard function. The exponential model The simplest model is the exponential model where T at z = 0 (usually referred to as the baseline) has exponential distribution with constant hazard exp( β 0 ). This is equivalent to assuming that σ = and ε has a standard extreme value distribution f(ε) =e ε eε, which has the density function shown in Figure 5.. (So e ε has the standard exponential distribution with constant hazard.) From this specification, it is easy to see that the distribution of T at any covariate vector z is exponential with constant hazard (independent of t) λ(t z) =e β 0 β z β pz p. PAGE 03

5 Figure 5.: The density function of the standard extreme value distribution f(x) x So automatically, we get proportional hazards models. For a given set of covariates (z,z 2,..., z p ), the corresponding survival function is S(t z) =e λ(t z)t, where λ(t z) =e β 0 β z β pz p.letβ j = β j.thenequivalently λ(t z) =e β 0 +β z + +β p zp. Therefore, if we increase the value of covariate z k (k =,,p) by on unit from z k to z k + while holding other covariate values fixed, then the ratio of the corresponding hazards is equal to λ(t z k +) λ(t z k ) =e β k. Thus e β k can be interpreted as the hazard ratio corresponding to one unit increase in the covariate z k,orequivalently,β k can be interpreted as the increase in log-hazard as the value of covariate z k increases by one unit (while other covariate values being held constant). Note: Another SAS procedure Proc Phreg fits a proportional hazards model to the data and outputs the regression coefficient estimates in log-hazard form (i.e., in β k). Therefore, if an PAGE 04

6 exponential model fits the data well (Proc Phreg will also fit the data well in this case.) then the regression coefficient estimates in outputs from Proc Lifereg using dist=exponential and Proc Phreg should be just opposite to each other (opposite sign but almost the same absolute value). We should be able to shift back and forth between these two models. Example (Autologous and Allogeneic Bone Marrow Transplants for Hodgkin s and Non- Hodgkin s Lymphoma, pages -2 of the textbook): Data on 43 bone marrow transplant patients were collected. Patients had either Hodgkin s disease or Non-Hodgkin s Lymphoma, and were given either an allogeneic (Allo) transplant (from a HLA match sibling donor) or autogeneic (Auto) transplant (their own marrow was cleansed and returned to them after a high dose of chemotherapy). Other covariates are Karnofsky score (a subjective measure of how well the patient is doing, ranging from 0-00) and waiting time (in months) from diagnosis to transplant. It is of substantial interest to see the difference in leukemia-free survival (in days) between those patients given an Allo or Auto transplant, after adjusting for patients disease status, Karnofsky score and waiting time. The data were given in Table.5 of the textbook. We used the following SAS program to fit an exponential model to the data title "Exponential fit"; proc lifereg data=bone; model time*status(0) = allo hodgkins kscore wtime / dist=exponential; run; where allo= for allogeneic transplant and allo=0 for autologous transplant, hodgkins= for Hodgkin s disease and hodgkins=0 for Non-Hodgkin s Lymphoma, kscore is the Karnofsky score and wtime is the waiting time. We got the following output: Exponential fit 7:35 Tuesday, March, 2005 The LIFEREG Procedure Model Information Data Set WORK.BONE Dependent Variable Log(time) Censoring Variable status Censoring Value(s) 0 PAGE 05

8 The Weibull model The only difference between the Weibull model and the exponential model is that the scale parameter σ is estimated rather than being set to be one. In this case, the distribution of σε is an extreme value distribution with scale parameter σ. The survival function of T at covariate value z =(,z,,z p ) T can be shown to be { S(t z) =exp [ te ] } zt β σ, where β =(β 0,,β p ) T hazard function is the vector of regression coefficients. Equivalently, in terms of (log)- logλ(t z) = ( σ ) logt logσ z T (β/σ). Let α =/σ, β0 = logσ β 0 /σ, andβj = β j /σ for j =,,p(again, pay close attention to the negative sign). Then we have logλ(t z) =(α )logt + β0 + z β + + z p βp. Thus we also get a proportional hazards model and the coefficient βk (k =,,p)alsohasthe interpretation that it is the increase in log-hazard when the value of covariate z k increases by one unit while other covariate values being held unchanged. The function is the baseline hazard (i.e., whenz =0). λ 0 (t) =t α e β 0 = αt α e β 0/σ = αt α e αβ 0 Note: If the Weibull model is a reasonable model for your data and you use Proc Lifereg and Proc Phreg to fit the data, then the regression coefficient estimates not only have opposite signs (except possibly for the intercept) but also have different magnitude (depending on whether σ>orσ<). Since β k = β k /σ for k =,,p.testingh 0 : β k = 0 is equivalent to testing H 0 : β k =0. If we are interested in calculating standard error for the estimate of β k and constructing a confidence interval for β k, we can use delta method for this purpose. PAGE 07

10 Scale Weibull Shape If we think that the Weibull model is a reasonable one, then the likelihood ratio test statistic is 2*( ( )) = 2.56 and the p-value = 0.096, not a strong evidence against the exponential model. From this model, we see similar results in the transplant methods. The log-normal model The log-normal model simply assumes that ε N(0, ). Let λ 0 (t) be the hazard function of T when β =0(β 0 = β = = β p = 0). Then it can be shown that λ 0 (t) has the following functional form where φ(x) = x2 /2 2π e is λ 0 (t) = φ ( ) log(t) σ [ ( )] Φ log(t), σ σt the probability density function and Φ(x) = x 2π e u2 /2 du is the cumulative distribution function of the standard normal distribution. Then the log-hazard function of T at any covariate value z canbeexpressedas logλ(t z) =logλ 0 (te zt β ) z T β. (5.3) Obviously we no longer have a proportional hazards model. Note: The function λ 0 (t) is not the baseline hazard function. If such a function is desired, it can be obtained from equation (5.3) by setting z =0. Some typical patterns that the hazard function λ 0 (t) assumes are presented in Figure 5.2. The inverted U-shaped of the log-normal hazard if often appropriate for repeated events such as a residential move (i.e, the interest is the time to next move). Immediately after a move, the hazard for another move is likely to be low, then increases with time, and eventually begins to decline since people tend to not move as they get older. The survival function S(t z) at any covariate value z can be expressed as Φ [S(t z)] = β 0 + β z + + β pz p αlog(t), (5.4) PAGE 09

11 Figure 5.2: Typical hazard functions for a log-normal model Hazard functions sigma=0.5 sigma= sigma= time or equivalently S(t z) =Φ[β 0 + β z + + β pz p αlog(t)], where α =/σ and β j = β j /σ for j =0,,,p. Thisisaprobit regression model with intercept depending on t. Note: Equation (5.4) indicates that the coefficients β j can be estimated using Proc Logistic or Proc Genmod by specifying probit link function. Specifically, pick a time point of interest, say t 0. Then dichotomize each subject based on his/her survival status at t 0 (in this case αlog(t 0 ) is absorbed into the intercept). Of cause, there are some limitations using this approach. First, there should not be censoring prior to time t 0. Second, the scale parameter σ is not estimable. Third, we will lose efficiency since we did not use all the information on the exact timing of the events. Since normal distribution and the logistic distribution we will introduce soon behave similarly to each other, the parameters βk here have similar interpretation as the parameters in log-logistic model. Example (Bone marrow transplant data revisited) If we want to fit a log-normal for the bone marrow transplant data, we use the following tt SAS program: title "Log-normal fit"; PAGE 0

13 The Log-Logistic Model The log-logistic model assumes that the disturbance term ε has a standard logistic distribution f(ε) = e ɛ ( + e ɛ ) 2. The density is plotted in Figure 5.3. Graphically, it looks like the standard normal density except that the standard normal density puts more mass around its mean value (0). The hazard function of T at any covariate value z has a closed form: where α =/σ. λ(t z) = αtα e zt β/σ +t α e zt β/σ, Figure 5.3: Density function of standard logistic distribution f(x) x Since the logistic distribution looks similar to the normal distribution, it is expected that the hazard function of log-logistic distribution would also look like that of log-normal distribution, i.e., would have an inverted U-shaped hazard. However, this is the case only for σ<. The hazard function of T when β = 0 for some values of σ s is presented in Figure 5.4. PAGE 2

14 Figure 5.4: Hazard functions for the log-logistic distribution Hazard functions sigma=0.5 sigma= sigma= time The random variable T has a very simple survival function at covariate value z S(t z) = +(te zt β ) /σ. Some simple algebra then shows that [ ] S(t z) log = β0 + β S(t z) z + + βpz p αlog(t), (5.5) where β j = β j /σ for j =0,,,p. This is nothing but a logistic regression model with the intercept depending on t. Since S(t z) is the probability of surviving to time t for any given time t, the ratio S(t z)/( S(t z)) is often called the odds of surviving to time t. Therefore, with one unit increase in z k while other covariates being held fixed, the odds ratio is given by S(t z k +)/( S(t z k +)) S(t z k )/( S(t z k )) =e β k for all t 0, which is a constant over time. Therefore, we have a proportional odds models. Hence βk can be interpreted as the log odds ratio (for surviving) with one unit increase in z k and β is the log odds ratio of dying before time t with one unit increase in z k. At the times when the event of failure is rare (such as the early phase of a study), βk can also be approximately interpreted as the log relative risk of dying. The log-logistic model is the only one that is both an AFT model and a proportional odds model. PAGE 3

16 Type III Analysis of Effects Wald Effect DF Chi-Square Pr > ChiSq allo hodgkins kscore <.000 wtime Analysis of Parameter Estimates Standard 95% Confidence Chi- Parameter DF Estimate Error Limits Square Pr > ChiSq Intercept allo hodgkins kscore <.000 wtime Scale We got the same conclusion that two transplants are not significantly different after adjusting for other covariates. The effects of other covariates are similar too. The gamma model The procedure Proc Lifereg in SAS actually fits a generalized gamma model (not a standard gamma model) to the data by assuming T 0 =e ε to have the following density function f(t) = δ (t δ /δ 2 ) /δ2 exp( t δ /δ 2 )/(tγ(/δ 2 )), where δ is called shape with label Gamma shape parameter in the output of Proc Lifereg when we specify dist=gamma. The hazard function of this gamma distribution does not have a closed form and is presented graphically in Figure 5.5 for some δ s. Note: Theyarenot the hazard functions of the standard gamma distribution. Clearly from this plot, the hazard function is an inverted U-shaped function of time when δ< and takes a U-shape when δ>. This feature makes gamma model very appropriate to model the survival times, especially for human. In practice, the hazard function is determined jointly by the scale parameter σ and the shape parameter δ and we need to examine the resulting PAGE 5

19 order to show this: Γ(x) 2πx x 2 e x as x. 3. δ = σ T z has the standard gamma distribution with the following density: f(t z) = tk e t γ γ K Γ(K), where K =/δ 2 is the shape parameter and γ = δ 2 e zt β is the scale parameter in the standard gamma distribution. However, Proc Lifereg will not fit a standard gamma distribution to data. We have to use grid search to fit a standard gamma distribution. Specifically, use the output from a generalized gamma distribution to get an idea about the true value of σ = δ. Then form a grid. For each grid point of σ (or δ), say, σ =.2, fit a gamma distribution with the specification dist=gamma noshape shape=.2 noscale scale=.2;. Select the model that gives the largest log-likelihood. Categorical Variables and Class statement in Proc Lifereg Ifacovariatez k is categorical, then we can use class statement in Proc Lifereg to tell the procedure that z k is categorical and just enter z k in the model statement in a usual way. Then the output of Proc Lifereg will provide a χ 2 statistic and a p-value to test the null hypothesis that the survival time is not associated with z k. Then it reports the estimates, standard errors, and p-values, etc., for the contrasts between each level and the highest level of the category. However, we need to create appropriate variables for the interaction. 5.3 Goodness-of-fit Using Likelihood Ratio Test The following fact on the models we described in this chapter allows us to perform likelihood ratio tests in order to pick a right model:. generalized gamma (δ, σ) standard gamma (δ = σ) exponential distribution (δ = σ = ). PAGE 8

20 2. generalized gamma (δ, σ) Weibull (δ =,σ) exponential distribution (δ = σ =). 3. generalized gamma (δ, σ) log-normal (δ =0,σ). We can use the above nested models to conduct likelihood ratio test for the bone marrow transplant data: Maximum likelihood values for different models Model Maximum Likelihood Gamma Log-logistic -6.5 Log-normal Weibull -6.2 Exponential Assuming the Gamma model is a reasonable model for the data, the LRT indicates that the log-normal model is equally good and the Weibull model is also acceptable. Since the log-normal model and the Weibull model have the same number of parameters, we might want to take the log-normal model as the final model based on the larger maximum log-likelihood value. PAGE 9

SUMAN DUVVURU STAT 567 PROJECT REPORT SURVIVAL ANALYSIS OF HEROIN ADDICTS Background and introduction: Current illicit drug use among teens is continuing to increase in many countries around the world.

Chapter 7 Survival Models Our final chapter concerns models for the analysis of data which have three main characteristics: (1) the dependent variable or response is the waiting time until the occurrence

Ordinal Regression Chapter 4 Many variables of interest are ordinal. That is, you can rank the values, but the real distance between categories is unknown. Diseases are graded on scales from least severe

Chapter 3 Logit Models for Binary Data We now turn our attention to regression models for dichotomous data, including logistic regression and probit analysis. These models are appropriate when the response

Nominal and ordinal logistic regression April 26 Nominal and ordinal logistic regression Our goal for today is to briefly go over ways to extend the logistic regression model to the case where the outcome

Chapter 4 Poisson Models for Count Data In this chapter we study log-linear models for count data under the assumption of a Poisson error structure. These models have many applications, not only to the

Chapter 95 Introduction Most statisticians have a set of probability tables that they refer to in doing their statistical wor. This procedure provides you with a set of electronic statistical tables that

Chapter 50 Two Correlated Proportions (Mcemar Test) Introduction This procedure computes confidence intervals and hypothesis tests for the comparison of the marginal frequencies of two factors (each with

Chapter 552 Gamma Distribution Fitting Introduction This module fits the gamma probability distributions to a complete or censored set of individual or grouped data values. It outputs various statistics

Introduction to General and Generalized Linear Models General Linear Models - part I Henrik Madsen Poul Thyregod Informatics and Mathematical Modelling Technical University of Denmark DK-2800 Kgs. Lyngby

Chapter 3 RANDOM VARIATE GENERATION In order to do a Monte Carlo simulation either by hand or by computer, techniques must be developed for generating values of random variables having known distributions.

Introduction to Event History Analysis DUSTIN BROWN POPULATION RESEARCH CENTER Objectives Introduce event history analysis Describe some common survival (hazard) distributions Introduce some useful Stata

Introduction to Hypothesis Testing Point estimation and confidence intervals are useful statistical inference procedures. Another type of inference is used frequently used concerns tests of hypotheses.

Regression Analysis Prof. Soumen Maity Department of Mathematics Indian Institute of Technology, Kharagpur Lecture - 2 Simple Linear Regression Hi, this is my second lecture in module one and on simple

Week 1 Survival Distributions, Hazard Functions, Cumulative Hazards 1.1 Definitions: The goals of this unit are to introduce notation, discuss ways of probabilistically describing the distribution of a

Chapter 408 Confidence Intervals for Exponential Reliability Introduction This routine calculates the number of events needed to obtain a specified width of a confidence interval for the reliability (proportion

A Tutorial on Logistic Regression Ying So, SAS Institute Inc., Cary, NC ABSTRACT Many procedures in SAS/STAT can be used to perform logistic regression analysis: CATMOD, GENMOD,LOGISTIC, and PROBIT. Each

Logistic Regression Introduction The simple and multiple linear regression methods we studied in Chapters 10 and 11 are used to model the relationship between a quantitative response variable and one or

Simple Linear Regression Inference 1 Inference requirements The Normality assumption of the stochastic term e is needed for inference even if it is not a OLS requirement. Therefore we have: Interpretation

Section 6.1 Joint Distribution Functions We often care about more than one random variable at a time. DEFINITION: For any two random variables X and Y the joint cumulative probability distribution function

Generalized Linear Models We have previously worked with regression models where the response variable is quantitative and normally distributed. Now we turn our attention to two types of models where the

Multiple Regression - Selecting the Best Equation When fitting a multiple linear regression model, a researcher will likely include independent variables that are not important in predicting the dependent

Survival Analysis Mark Lunt Arthritis Research UK Centre for Excellence in Epidemiology University of Manchester 01/12/2015 Survival Analysis is concerned with the length of time before an event occurs.

Statistics in Retail Finance 1 Overview > So far we have focussed mainly on application scorecards. In this chapter we shall look at behavioural models. We shall cover the following topics:- Behavioural

One-Degree-of-Freedom Tests Test for group occasion interactions has (number of groups 1) number of occasions 1) degrees of freedom. This can dilute the significance of a departure from the null hypothesis.

Final Mathematics 51, Section 1, Fall 24 Instructor: D.A. Levin Name YOU MUST SHOW YOUR WORK TO RECEIVE CREDIT. A CORRECT ANSWER WITHOUT SHOWING YOUR REASONING WILL NOT RECEIVE CREDIT. Problem Points Possible

Survival Analysis Using SPSS By Hui Bian Office for Faculty Excellence Survival analysis What is survival analysis Event history analysis Time series analysis When use survival analysis Research interest

Modeling the Claim Duration of Income Protection Insurance Policyholders Using Parametric Mixture Models Abstract This paper considers the modeling of claim durations for existing claimants under income

Unit 29 Chi-Square Goodness-of-Fit Test Objectives: To perform the chi-square hypothesis test concerning proportions corresponding to more than two categories of a qualitative variable To perform the Bonferroni

Numerical Summarization of Data OPRE 6301 Motivation... In the previous session, we used graphical techniques to describe data. For example: While this histogram provides useful insight, other interesting

Logit and Probit Brad 1 1 Department of Political Science University of California, Davis April 21, 2009 Logit, redux Logit resolves the functional form problem (in terms of the response function in the

A Primer on Mathematical Statistics and Univariate Distributions; The Normal Distribution; The GLM with the Normal Distribution PSYC 943 (930): Fundamentals of Multivariate Modeling Lecture 4: September

Checking proportionality for Cox s regression model by Hui Hong Zhang Thesis for the degree of Master of Science (Master i Modellering og dataanalyse) Department of Mathematics Faculty of Mathematics and

Chapter 450 on-inferiority Tests for Two Means using Differences Introduction This procedure computes power and sample size for non-inferiority tests in two-sample designs in which the outcome is a continuous

125: Chi-Square Goodness of Fit Tests CD12-1 125: CHI-SQUARE GOODNESS OF FIT TESTS In this section, the χ 2 distribution is used for testing the goodness of fit of a set of data to a specific probability

136 Poisson Regression Analysis 13. Poisson Regression Analysis We have so far considered situations where the outcome variable is numeric and Normally distributed, or binary. In clinical work one often