Transcription

1 Using Time Dependent Covariates and Time Dependent Coefficients in the Cox Model Terry Therneau Cindy Crowson Mayo Clinic July 1, Introduction This vignette covers 3 different but interrelated concepts: ˆ An introduction to time dependent covariates, along with some of the most common mistakes. ˆ Tools for creating time-dependent covariates, or rather the data sets used to encode them. ˆ Time dependent coefficients. 2 Time dependent covariates One of the strengths of the Cox model is its ability to encompass covariates that change over time. The practical reason that time-dependent covariates work is based on the underlying way in which the Cox model works: at each event time the program compares the current covariate values of the subject who had the event to the current values of all others who were at risk at that time. One can think of it as a lottery model, where at each death time there is a drawing to decide which subject wins the event. Each subject s risk score exp(xβ) determines how likely they are to win, e.g., how many tickets they have purchased for the drawing. The model tries to assign a risk score to each subject that best predicts the outcome of each drawing based on ˆ The risk set: which subjects are present for each event; the set of those able to win the prize. ˆ The covariate values of each subject just prior to the event time. The model has a theoretical foundation in martingale theory, a mathematical construct which arose out of the study of games of chance. A key underlying condition for a martingale like game is that present actions depend only on the past. The decision of whether to play (is one in the risk set or not) and the size of a bet (covariates) can depend in any way on prior bets and patterns of won/lost, but cannot look into the future. If this holds then multiple properties can be proven about the resulting process. 1

2 The key rule for time dependent covariates in a Cox model is simple and essentially the same as that for gambling: you cannot look the future. A covariate may change in any way based on past data or outcomes, but it may not reach forward in time. A good example of this is found in a recent analysis from the Mayo Clinic study of aging (MCSA), a study which enrolled a stratified random sample from the population of Olmsted County and then has followed them forward in time. The occurence of mild cognitive impairment (MCI), dementia, and death are all of interest. The paper starts out with a table comparing baseline covariates for those who never progress to MCI versus those who ever did, there is also a table of baseline covariates versus survival. Both of these are fine: if you think in terms of an R formula they could be written as future past. A model that predicts survival as a function of ever versus never MCI is not correct, however; that is a model with a future occurence on both sides of the equation. One of the more well known examples of this error is analysis by treatment response: at the end of a trial a survival curve is made comparing those who had an early response to treatment (shrinkage of tumor, lowering of cholesterol, or whatever) to those who did not, and it discovered that responders have a better curve. A Cox model fit to the same data will demonstrate a strong significant effect. The problem arises because any early deaths, those that occur before response can be assessed, will all be assigned to the non-responder group, even deaths that have nothing to do with the condition under study. Below is a simple example based on the advanced lung cancer data set. Assume that subjects came in every monthly for 12 cycles of treatment, and randomly declare a response for 5% of the subjects at each visit. > set.seed(1953) # a good year > nvisit <- floor(pmin(lung$time/30.5, 12)) > response <- rbinom(nrow(lung), nvisit,.05) > 0 > badfit <- survfit(surv(time/365.25, status) ~ response, data=lung) > plot(badfit, mark.time=false, lty=1:2, xlab="years post diagnosis", ylab="survival") > legend(1.5,.85, c("responders", "Non-responders"), lty=2:1, bty='n') 2

3 Survival Responders Non responders Years post diagnosis What is most surprising about this error is the size of the false effect that is produced. A Cox model using the above data reports a hazard ratio of 1.9 fold with a p-value of less than 1 in The alarm about this incorrect approach has been sounded often [1, 2, 7] but the analysis is routinely re-discovered. A slightly subtler form of the error is discussed in Redmond et al [6]. Breast cancer chemotherapy patients were divided into three groups based on whether the patient eventually received > 85%, 65 85% or < 65% of the dose planned at the start of their treatment. The chemotherapy regiment spans 12 weeks of treatment and the early deaths, not surprisingly, do not get all their dose. If response is instead coded as a time-dependent covariate whose values depend only on the past, then the problem disappears. For treatment response this will be a variable that starts at 0 for all subjects and is recoded to 1 only when the response occurs. For dose it would measure cumulative dose to date. There are many variations on the error: interpolation of the values of a laboratory test linearly between observation times, removing subjects who do not finish the treatment plan, imputing the date of an adverse event as midway between observation times, etc. Using future data will often generate large positive or negative bias in the coefficients, but sometimes it generates little bias at all. It is nearly impossible to predict a priori which of these will occur in any given data set. Using such a covariate is similar to jogging across the Los Angeles freeway: disaster is not guaranteed but it is likely. The most common way to encode time-dependent covariates is to use the (start, stop] form of the model. 3

4 > fit <- coxph(surv(time1, time2, status) ~ age + creatinine, data=mydata) In data set mydata a patient might have the following observations subject time1 time2 status age creatinine In this case the variable age = age at entry to the study stays the same from line to line, while the value of creatinine varies and is treated as 1.3 over the interval (0, 15], 1.5 over (15, 46], etc. The intervals are open on the left and closed on the right, which means that the creatinine is taken to be 1.3 on day 15. The status variable describes whether or not each interval ends in an event. One common question with this data setup is whether we need to worry about correlated data, since a given subject has multiple observations. The answer is no, we do not. The reason is that this representation is simply a programming trick. The likelihood equations at any time point use only one copy of any subject, the program picks out the correct row of data at each time. There two exceptions to this rule: ˆ When subjects have multiple events, then the rows for events are correlated and a cluster variance is needed. ˆ When a subject appears in overlapping intervals. This however is almost always a data error, since it corresponds to two copies of the subject being present in the same strata at the same time, e.g., they could meet themselves on the sidewalk. A subject can be at risk in multiple strata at the same time, however. This corresponds to being simultaneously at risk for two distinct outcomes. 3 Building time-dependent sets with tmerge 3.1 The function A useful function for building data sets is tmerge, which is part of the survival library. The motivating case for tmerge came from a particular problem: the Rochester Epidemiology Project has tracked all subjects living in Olmsted County, Minnesota, from 1965 to the present. For an investigation of cumulative comorbidity we had three data sets ˆ base: demographic data such as sex and birth date ˆ timeline: one or more rows for each subject containing age intervals during which they were a resident of the county. The important variables are id, age1 and age2; each (age1, age2) pair marks an interval of residence. ˆ outcome: one row for each age/outcome pair of interest. The outcomes were 20 comorbid conditions as defined by NIH. 4

5 The structure for building the data is shown below. (The data for this example unfortunately cannot be included in the survival library so the code is shown but not executed.) > newd <- tmerge(data1=base, data2=timeline, id=repid, tstart=age1, tstop=age2, options(id="repid")) > newd <- tmerge(newd, outcome, id=repid, mtype = cumevent(age)) > newd <- with(subset(outcome, event='diabetes'), tmerge(newd, id=repid, diabetes= tdc(age))) > newd <- with(subset(outcome, event='arthritis'), tmerge(newd, id=repid, event =tdc(age))) The first call to tmerge adds the timeline for each observation to the baseline data. The tstart and tstop arguments refer to the starting and ending times for each subject and are taken from data set 2 (data2). The options argument tells the routine that the identifier variable in data set 1 is called repid, and will cause the identifier variable in the final output data set newd to have that name. By default, the names of the three key variables are id, tstart, and tstop ; these uniquely identify each row of the final data set. Each subsequent call adds a new variable to the data set. The second line creates an event variable which is a cumulative count of the number of comorbidities thus far, for each subject. The third line creates a time dependent covariate (tdc) which will be 0 until the age of diabetes and is 1 thereafter, the fourth line creates a time dependent variable for the presence of arthritis. This is the basic working approach for tmerge: first establish baseline covariates and the range over which the subject is at risk, and then add events and/or time dependent covariates to the data set one by one. These additions will often increase the number of rows in the data set. Say at some stage subject Smith has time intervals of (0,15), (15,40) and (40,100) and we add a time dependent covariate sbp (systolic blood pressure) which is evaluated at months 6, 12, and 24 with values of 134, 126, and 140, respectively. In the resulting data set Smith will have intervals of (0,6) (6,12) (12,15) (15,24) (24,40) (40,100) The value over the interval (0,6) will be the value of the variable sbp in data set 1 if the variable existed there it is assumed to contain the baseline value or NA otherwise CGD data set Chronic granulomatous disease (CGD) is a heterogeneous group of uncommon inherited disorders characterized by recurrent pyogenic infections that usually begin early in life and may lead to death in childhood. In 1986, Genentech, Inc. conducted a randomized, double-blind, placebocontrolled trial in 128 CGD patients who received Genentech s humanized interferon gamma (rifn-g) or placebo three times daily for a year. Data were collected on all serious infections until the end of followup, which occurred before day 400 for most patients. One patient was taken off on the day of his last infection; all others have some followup after their last episode. Below are the first 10 observations, see the help page for cgd0 for the full list of variable names. The last few columns contain the duration of follow-up for the subject followed by infection times. Subject 1 was followed for 414 days and had 2 infections on days 219 and 373, subject 2 had 7 infections, and subject 3 had none. 5

6 The data set above is included as cgd0 in the survival library. We want to turn this into a data set that has survival in a counting process form. ˆ Each row of the resulting data set represents a time interval (time1, time2] which is open on the left and closed on the right. Covariate values for that row are the covariate values that apply over that interval. ˆ The event variable for each row i is 1 if that time interval ends with an event and 0 otherwise. We don t want the variables etime1 etime7 in the final data set, so they are left out of the data1 argument in the first call. > newcgd <- tmerge(cgd0[, 1:13], cgd0, id=id, tstop=futime) > newcgd <- tmerge(newcgd, cgd0, id=id, infect = event(etime1)) > newcgd <- with(cgd0, tmerge(newcgd, id=id, infect = event(etime2))) > newcgd <- tmerge(newcgd, cgd0, id=id, infect = event(etime3)) > newcgd <- tmerge(newcgd, cgd0, id=id, infect = event(etime4), infect= event(etime5), infect=event(etime6), infect= event(etime7)) > attr(newcgd, "tcount") early late gap within boundary leading trailing tied infect infect infect infect > newcgd <- tmerge(newcgd, newcgd, id, enum=cumtdc(tstart)) > all.equal(newcgd[, c("id", "tstart", "tstop", "infect")], cgd [, c("id", "tstart", "tstop", "status")], check.attributes=false) [1] TRUE ˆ A standard way to build data sets is one addition at a time, as shown by the event(etime1) and event(etime2) lines above. When multiple additions are done using the same d the same data set, however, they can also be done using as single call as shown by the line that adds etime4 etime7. When there are multiple arguments they are processed sequentially. 6

7 ˆ Additions with a missing time value are skipped. ˆ The result of tmerge is a data frame with a few extra attributes. One of these, tcount, is designed to help visualize the process, and was printed out after the last step above. (Printing after every step may often be useful.) Assume that a subject already had 3 intervals of (2,5), (5,10) and (14,40). A new event added at time 1 would be early while one at time 50 is after any interval and would be recorded as late. An event at time 3 is within an interval, one at 5 is on the border of two intervals, one at 14 is at the leading edge of an interval and one at time 10 in on the trailing edge. In this data set all new additions fell strictly within prior intervals. We also see that etime6 and etime7 each added only a single event to the data. ˆ If two observations in data2 for a single person share exactly the same time, the created value will be the sum of the contributions. The tied column tells how often this happened; in some data sets this behavior might not be desired and one would need to break the ties before calling tmerge. ˆ Sometimes the where form of the call may be more convenient than using the data2 argument. An example is shown in the addition of etime2. ˆ The last addition above, after printing the tcount attribute, adds a simple time-dependent variable enum which is a running observation count for each subject. This can often be a useful variable in later models or processing, e.g. enum==1 select off the first row for each subject. ˆ The extra attributes of the data frame are ephemeral: they will be lost as soon as any further manipulation is done. This is intentional. The last line above shows that the created data set is identical to cgd, a (start, stop] version of the CGD data, also part of the survival library, which had been created by hand several years earlier. 3.2 Stanford heart transplant The jasa data set contains information from the Stanford heart transplant study, in the form that it appeared in the paper of Crowley and Hu [3]. Each row also contain also contains the age, time to transplant and time to last follow-up calculated from these dates. Patients were on medical treatment from their entry to the study until a matching heart became available, at which time they transferred to surgical treatment. As is often the case with real data, this data set contains a few anomalies that need to be dealt with when setting up an analysis data set. 1. The coefficients in table 6.1 of the definitive analysis found in Kalbfliesch and Prentice [4] will only be obtained if covariates are defined in precisely the same way. For age this is (age in days)/ years, and for year of enrollment it is the number of years since the start of the study: (entry date /10/1)/ One subject died on the day of entry. However (0,0) is an illegal time interval for the program. It suffices to have them die on day

8 3. A subject transplanted on day 10 is considered to have been on medical treatment for days 1 10 and as transplanted starting on day 11. That is, except for patient 38 who died during the procedure on day 5. They should be treated as a transplant death; the problem is resolved by moving the transplant day to 4.5. Since time is in days the fractional time of 0.5 could be any chosen value t with 0 < t < 1, it will not affect the results. > tdata <- jasa[, -(1:4)] #leave off the dates, temporary data set > tdata$futime <- pmax(.5, tdata$futime) # the death on day 0 > indx <- with(tdata, which(wait.time == futime)) > tdata$wait.time[indx] <- tdata$wait.time[indx] -.5 #the tied transplant > sdata <- tmerge(tdata, tdata, id=1:nrow(tdata), death = event(futime, fustat), trans = tdc(wait.time)) > attr(sdata, "tcount") early late gap within boundary leading trailing tied death trans > coxph(surv(tstart, tstop, death) ~ age + trans, sdata) Call: coxph(formula = Surv(tstart, tstop, death) ~ age + trans, data = sdata) coef exp(coef) se(coef) z p age trans Likelihood ratio test=5.17 on 2 df, p= n= 170, number of events= 75 This example shows one special case for the tmerge function that is moderately common: when the data1 and data2 arguments are the same, and the first created variable is an event code, then the range for each subject is inferred to be from 0 to the event time: an explicit tstop argument is not required. It also makes use of a two argument form of event. Each of the event, cumevent, tdc and cumtdc functions may have a second argument, which will be used as the value or increment to the event code or time dependent covariate. If not present a value of 1 is used. Also note that if the variable being created is already a part of data1, then our updates make changes to that variable. Be careful of this. This feature is what allowed for the infection indicator to be build up incrementally in the cgd example given earlier, but quite surprising results can occur when you think a new variable is being created de novo but its name already is in use. (As an example change trans to transplant in the code just above). For a variable that is not in data1, the starting point is either a vector of NA values or a vector of zeros; the first is used for a tdc (or cumtdc) call that has two arguments, and the zero vector for all other cases. 8

9 The tcount table for the above fit shows all the deaths at the trailing edge of their interval, not surprising since last follow-up was used to define the interval of risk. Three of the transplants happened on day 0, one of which we moved to 0.5, the other 2 listed as occurring on the leading edge of the follow-up interval. The other 67 transplants were strictly within the (0, last follow up) interval of each subject. As a further example of time dependent covariates consider the PBC data. The pbc data set contains baseline data and follow-up status for a set of subjects with primary biliary cirrhosis, while the pbcseq data set contains repeated laboratory values. The first data set contains data on 312 subjects in a clinical trial plus 106 that agreed to be followed off protocol, the second data set has data only on the protocol subjects. > temp <- subset(pbc, id <= 312, select=c(id:sex, stage)) > pbc2 <- tmerge(temp, temp, id=id, status = event(time, status)) > pbc2 <- tmerge(pbc2, pbcseq, id=id, ascites = tdc(day, ascites), bili = tdc(day, bili), albumin = tdc(day, albumin), protime = tdc(day, protime), alkphos = tdc(day, alk.phos)) > coef(coxph(surv(time, status==2) ~ log(bili) + log(protime), pbc)) log(bili) log(protime) > coef(coxph(surv(tstart, tstop, status==2) ~ log(bili) + log(protime), pbc2)) log(bili) log(protime) The coefficients of bilirubin and prothrombin time are somewhat larger in the time-dependent analysis than using only baseline values. In this autoimmune disease there is steady progression of liver damage, along with a steady rise in these two markers of dysfunction. The baseline analysis captures patients disease status at the start, the time-dependent also is able to account for those who progress more quickly. > attr(pbc2, "tcount") early late gap within boundary leading trailing ascites bili albumin protime alkphos tied ascites 0 bili 0 albumin 0 protime 0 alkphos 0 9

10 The tcount results are interesting. For the first addition of ascites we have 312 observations on a leading edge of follow up (all of the lab values at time 0) and 1495 within the subjects follow-up interval. This causes a new break point to be added at each laboratory date, and for subsequent additions these 1495 are all on a boundary of two intervals. Another 138 lab values occurred after the last follow-up date of the pbc data set and are ignored. The data for the pbcseq data set was gathered at a later calendar time. Since having lab tests done is a certain marker of survival, would a better analysis have first used this information to extend the last follow-up date for these 138 subjects? Not necessarily. Odd things happen in survival analysis when risk sets are extended piecemeal. A basic tenet of the Cox model (or a survival curve) is that if someone is marked as being at risk over some interval (s, t), this means that if they had had an event over that interval, we would have recorded it. Say someone ended their initial follow-up time at 3000 days and then had a lab test at 3350 days (subjects returned about once a year). If we only extend the time of those who had a test, then saying that this subject was at risk during the interval (3000, 3350) is false: if they had died in that interval, they would not have had the lab test and not obtained the extension. We need to extend the followup time of all subjects. In the study all subjects were actively followed up and the results of that endeavor are in the futime and status variables of the pbcseq data set. Extension needs to use those variables and not just the presence of another laboratory result. 4 Time dependent coefficients Time dependent covariates and time dependent coefficients are two different extensions of a Cox model, as shown in the two equations below. λ(t) = λ 0 (t)e βx(t) (1) λ(t) = λ 0 (t)e β(t)x (2) Equation (1) is a time dependent covariate, a commonly used and well understood usage. Equation (2) has a time dependent coefficient. These models are much less common, but represent one way to deal with non-proportional hazards the proportional hazard assumption is precisely that the coefficient does not change over time: β(t) = c. The cox.zph function will plot an estimate of β(t) for a study and is used to diagnose and understand non-proportional hazards. Here for example is a test case using the veterans cancer data. > options(show.signif.stars = FALSE) # display intelligence > vfit <- coxph(surv(time, status) ~ trt + prior + karno, veteran) > vfit Call: coxph(formula = Surv(time, status) ~ trt + prior + karno, data = veteran) coef exp(coef) se(coef) z p trt prior

11 karno e-11 Likelihood ratio test=43 on 3 df, p=2.41e-09 n= 137, number of events= 128 > quantile(veteran$karno) 0% 25% 50% 75% 100% > zp <- cox.zph(vfit, transform= function(time) log(time +20)) > zp rho chisq p trt prior karno GLOBAL NA > plot(zp[3]) > abline(0,0, col=2) Time Beta(t) for karno Karnofsky score is a very important predictor, but its effect is not constant over time as shown by both the test and the plot. Early on it has a large negative effect: the risk of someone 11

12 at the first quartile is approximately exp(35*.03377) = 3.2 fold times that of someone at the third quartile, but by 200 days this has waned. The cox.zph function does not produce a formal fit, however, only a graph and a linear test of whether the graph is flat. What if we want to actually fit the model? If β(t) is assumed to have a simple functional form we can fool an ordinary Cox model program in to doing the fit. The particular form β(t) = a + b log(t) has for instance often been assumed. Then β(t)x = ax + b log(t)x = ax + bz for the special time dependent covariate z = log(t)x. The time scale for the plot above of log(t + 20) was chosen to make the first part of the plot roughly linear. The simple linear model does not fit over the entire range, but we will forge ahead. (After all, most who fit the log(t) form have not bothered to even look at a plot.) An obvious but incorrect approach is > vfit2 <- coxph(surv(time, status) ~ trt + prior + karno + I(karno * log(time + 20)), data=veteran) This mistake has been made often enough the the coxph routine has been updated to print an error message for such attempts. The issue is that the above code does not actually create a time dependent covariate, rather it creates a time-static value for each subject based on their value for the covariate time; no differently than if we had constructed the variable outside of a coxph call. This variable most definitely breaks the rule about not looking into the future, and one would quickly find the circularity: large values of time predict long survival, because long survival leads to large values for time; the resulting model coefficient is large. A true time-dependent covariate can be constructed using the time-transform functionality of coxph. > vfit3 <- coxph(surv(time, status) ~ trt + prior + karno + tt(karno), data=veteran, tt = function(x, t,...) x * log(t+20)) > vfit3 Call: coxph(formula = Surv(time, status) ~ trt + prior + karno + tt(karno), data = veteran, tt = function(x, t,...) x * log(t + 20)) coef exp(coef) se(coef) z p trt prior karno e-05 tt(karno) Likelihood ratio test=53.8 on 4 df, p=5.7e-11 n= 137, number of events= 128 The time dependent coefficient is estimated to be β(t) = * log(t + 20). We can add said line to the above plot via abline(coef(vfit3)[3:4]). Not surprisingly, the result is rather too high for time >

13 (The same dichotomy exists in SAS phreg, by the way. A new covariate based on time will be fixed, while a programming statement within phreg that uses time as a variable will generate time-dependent objects. The error is less likely there because phreg s model statement has no equivalent to the I() function, i.e., you cannot create new variables on-the-fly.) 5 Predictable time-dependent covariates Occasionally one has a time-dependent covariate whose values in the future are predictable. The most obvious of these is patient age, occasionally this may also be true for the cumulative dose of a drug. If age is entered as a linear term in the model, then the effect of changing age can be ignored in a Cox model, due to the structure of the partial likelihood. Assume that subject i has an event at time t i, with other subject j R i at risk at that time, with a denoting age. The partial likelihood term is e β ai eβ (ai+ti) = β aj j R i e β (aj+ti) j R i e We see that using time-dependent age (the right hand version) or age at baseline (left hand), the partial likelihood term is identical since exp(βt i ) cancels out of the fraction. However, if the effect of age on risk is non-linear, this cancellation does not occur. Since age changes continuously, we would in theory need a very large data set to completely capture the effect, an interval per day to match the usual resolution for death times. In practice this level of resolution is not necessary; though we all grow older, risk does not increase so rapidly that we need to know our age to the day! One method to create a time-changing covariate is to use the time-transform feature of coxph. Below is an example using the pbc data set. The longest follow-up time in that data set is over 13 years, follow-up time is in days, and we might worry that the intermediate data set would be huge. The program only needs the value of the time dependent covariate(s) for each subject at the times of events, however, so the maximum number of rows in the intermediate data set is the number of subjects times the number of unique event times. > pfit1 <- coxph(surv(time, status==2) ~ log(bili) + ascites + age, pbc) > pfit2 <- coxph(surv(time, status==2) ~ log(bili) + ascites + tt(age), data=pbc, tt=function(x, t,...) { age <- x + t/ cbind(age=age, age2= (age-50)^2, age3= (age-50)^3) }) > pfit2 Call: coxph(formula = Surv(time, status == 2) ~ log(bili) + ascites + tt(age), data = pbc, tt = function(x, t,...) { age <- x + t/ cbind(age = age, age2 = (age - 50)^2, age3 = (age - 50)^3) }) 13

14 coef exp(coef) se(coef) z p log(bili) 1.05e e e < 2e-16 ascites 1.29e e e e-07 tt(age)age 6.66e e e e-06 tt(age)age2 3.04e e e tt(age)age3-9.17e e e Likelihood ratio test=179 on 5 df, p=0 n= 312, number of events= 125 (106 observations deleted due to missingness) > anova(pfit2) Analysis of Deviance Table Cox model: response is Surv(time, status == 2) Terms added sequentially (first to last) loglik Chisq Df Pr(> Chi ) NULL log(bili) ascites e-09 tt(age) < 2.2e-16 > # anova(pfit1, pfit2) #this fails > 2*(pfit2$loglik - pfit1$loglik)[2] [1] Since initial age is in years and time is in days, it was important to scale within the pspline function. The likelihood ratio of 10.8 on 2 degrees of freedom shows that the additional terms are mildly significant. When there are one or more terms on the right hand side of the equation marked with the tt() operator, the program will pre-compute the values of that variable for each unique event time. A user-defined function is called with arguments of ˆ the covariate: whatever is inside the tt() call ˆ the event time ˆ the event number: if there are multiple strata and the same event time occurs in two of them, they can be treated separately ˆ the weight for the observation, if the call used weights There is a single call to the function with a large x vector, it contains an element for each subject at risk at each event time. If there are multiple tt() terms in the formula, then the tt argument should be a list of functions with the requisite number of elements. 14

15 There are other interesting uses for the time-transform capability. One example is O Brien s logit-rank test procedure [5]. He proposed replacing the covariate at each event time with a logit transform of its ranks. This removes the influence of any outliers in the predictor x. For this case we ignore the event time argument and concentrate on the groupings. > function(x, t, riskset, weights){ obrien <- function(x) { r <- rank(x) (r-.5)/(.5+length(r)-r) } unlist(tapply(x, riskset, obrien)) } function(x, t, riskset, weights){ obrien <- function(x) { r <- rank(x) (r-.5)/(.5+length(r)-r) } unlist(tapply(x, riskset, obrien)) } This relies on the fact that the input arguments to tt() are ordered by the event number or riskset. This function is used as a default if no tt argument is present in the coxph call, but there are tt terms in the model formula. (Doing so allowed me to depreciate the survobrien function). Another interesting usage is to replace the data by simple ranks, not rescaled to 0 1. > function(x, t, riskset, weights) unlist(tapply(x, riskset, rank)) function(x, t, riskset, weights) unlist(tapply(x, riskset, rank)) The score statistic for this model is (C D)/2, where C and D are the number of concordant and discordant pairs, see the survconcordance function. The score statistic from this fit is then a test for significance of the concordance statistics, and is in fact the basis for the standard error reported by survconcordance. The O Brien test can be viewed as concordance statistic that gives equal weight to each event time, whereas the standard concordance weights each event proportionally to the size of the risk set. (The Cox score statistic depends on the mean x at each event time; since ranks go from 1 to number at risk the mean also scales.) Although handy, the computational impact of the tt argument should be considered before using it. The Cox model requires computation of a weighted mean and variance of the covariates at each event time, a process that is inherently O(ndp 2 ) where n = the sample size, d = the number of events and p= the number of covariates. Much of the algorithmic effort in coxph() is to use updating methods for the mean and variance matrices, reducing the compute time to O((n + d)p 2 ). When a tt term appears updating is not possible; for even moderate size data sets the impact of nd versus n + d can be surprising. The time-transform is a new addition and still has some rough edges. At this moment the x = T RU E argument is needed to get proper residuals and predicted values, and termplot 15

Tips for surviving the analysis of survival data Philip Twumasi-Ankrah, PhD Big picture In medical research and many other areas of research, we often confront continuous, ordinal or dichotomous outcomes

Survey, Statistics and Psychometrics Core Research Facility University of Nebraska-Lincoln Log-Rank Test for More Than Two Groups Prepared by Harlan Sayles (SRAM) Revised by Julia Soulakova (Statistics)

1 Multivariate Logistic Regression As in univariate logistic regression, let π(x) represent the probability of an event that depends on p covariates or independent variables. Then, using an inv.logit formulation

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

Randomized trials versus observational studies The case of postmenopausal hormone therapy and heart disease Miguel Hernán Harvard School of Public Health www.hsph.harvard.edu/causal Joint work with James

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.

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

STATISTICS 8, FINAL EXAM NAME: KEY Seat Number: Last six digits of Student ID#: Circle your Discussion Section: 1 2 3 4 Make sure you have 8 pages. You will be provided with a table as well, as a separate

Research methods - II 3 2. Simple Linear Regression Simple linear regression is a technique in parametric statistics that is commonly used for analyzing mean response of a variable Y which changes according

PholC60 September 001 DATA INTERPRETATION AND STATISTICS Books A easy and systematic introductory text is Essentials of Medical Statistics by Betty Kirkwood, published by Blackwell at about 14. DESCRIPTIVE

CORRELATION AND REGRESSION / 47 CHAPTER EIGHT CORRELATION AND REGRESSION Correlation and regression are statistical methods that are commonly used in the medical literature to compare two or more variables.

Cirrhosis and HCV Jonathan Israel M.D. Outline Relationship of fibrosis and cirrhosisprevalence and epidemiology. Sequelae of cirrhosis Diagnosis of cirrhosis Effect of cirrhosis on efficacy of treatment

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

MINITAB ASSISTANT WHITE PAPER This paper explains the research conducted by Minitab statisticians to develop the methods and data checks used in the Assistant in Minitab 17 Statistical Software. One-Way

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

Models Where the Fate of Every Individual is Known This class of models is important because they provide a theory for estimation of survival probability and other parameters from radio-tagged animals.

Systematic Reviews and Meta-analyses Introduction A systematic review (also called an overview) attempts to summarize the scientific evidence related to treatment, causation, diagnosis, or prognosis of

Data Mining: An Overview of Methods and Technologies for Increasing Profits in Direct Marketing C. Olivia Rud, VP, Fleet Bank ABSTRACT Data Mining is a new term for the common practice of searching through

Size of a study 15.1 Introduction It is important to ensure at the design stage that the proposed number of subjects to be recruited into any study will be appropriate to answer the main objective(s) of

Appendix D Basic Measurement And Statistics The following information was developed by Steven Rothke, PhD, Department of Psychology, Rehabilitation Institute of Chicago (RIC) and expanded by Mary F. Schmidt,

PEER REVIEW HISTORY BMJ Open publishes all reviews undertaken for accepted manuscripts. Reviewers are asked to complete a checklist review form (http://bmjopen.bmj.com/site/about/resources/checklist.pdf)

Comparing return to work outcomes between vocational rehabilitation providers after adjusting for case mix using statistical models Prepared by Jim Gaetjens Presented to the Institute of Actuaries of Australia

General Remarks This template of a data extraction form is intended to help you to start developing your own data extraction form, it certainly has to be adapted to your specific question. Delete unnecessary

Introduction to Statistics and Quantitative Research Methods Purpose of Presentation To aid in the understanding of basic statistics, including terminology, common terms, and common statistical methods.

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.

Goal of a Formal Sensitivity Analysis To replace a general qualitative statement that applies in all observational studies the association we observe between treatment and outcome does not imply causation

Richard L. Scheaffer University of Florida The reference material and many examples for this section are based on Chapter 8, Analyzing Association Between Categorical Variables, from Statistical Methods

COMMON CORE STATE STANDARDS FOR Mathematics (CCSSM) High School Statistics and Probability Mathematics High School Statistics and Probability Decisions or predictions are often based on data numbers in

CHAPTER 5. A POPULATION MEAN, CONFIDENCE INTERVALS AND HYPOTHESIS TESTING 5.1 Concepts When a number of animals or plots are exposed to a certain treatment, we usually estimate the effect of the treatment

S03-2008 The Difference Between Predictive Modeling and Regression Patricia B. Cerrito, University of Louisville, Louisville, KY ABSTRACT Predictive modeling includes regression, both logistic and linear,

Statistical analysis using Microsoft Excel Microsoft Excel spreadsheets have become somewhat of a standard for data storage, at least for smaller data sets. This, along with the program often being packaged

Cairo University Faculty of Economics and Political Science Statistics Department English Section Students' Opinion about Universities: The Faculty of Economics and Political Science (Case Study) Prepared

Case Study in Data Analysis Does a drug prevent cardiomegaly in heart failure? Harvey Motulsky hmotulsky@graphpad.com This is the first case in what I expect will be a series of case studies. While I mention

Simple Linear Regression in SPSS STAT 314 1. Ten Corvettes between 1 and 6 years old were randomly selected from last year s sales records in Virginia Beach, Virginia. The following data were obtained,

Genetic Discoveries and the Role of Health Economics Collaboration for Outcome Research and Evaluation (CORE) Faculty of Pharmaceutical Sciences University of British Columbia February 02, 2010 Content

CHAPTER 40 When Does it Make Sense to Perform a Meta-Analysis? Introduction Are the studies similar enough to combine? Can I combine studies with different designs? How many studies are enough to carry

In this presentation, you will be introduced to data mining and the relationship with meaningful use. Data mining refers to the art and science of intelligent data analysis. It is the application of machine

Parametric and Nonparametric: Demystifying the Terms By Tanya Hoskin, a statistician in the Mayo Clinic Department of Health Sciences Research who provides consultations through the Mayo Clinic CTSA BERD

Chapter 7 Summary and general discussion Summary and general discussion In this thesis, treatment of vitamin K antagonist-associated bleed with prothrombin complex concentrate was addressed. In this we

VI: 2 ELEMENTS FOR A PUBLIC SUMMARY Bicalutamide (CASODEX 1 ) is a hormonal therapy anticancer agent, used for the treatment of prostate cancer. Hormones are chemical messengers that help to control the

Normality Testing in Excel By Mark Harmon Copyright 2011 Mark Harmon No part of this publication may be reproduced or distributed without the express permission of the author. mark@excelmasterseries.com