\documentclass{article}[11pt] \usepackage{Sweave} \usepackage{amsmath} \addtolength{\textwidth}{1in} \addtolength{\oddsidemargin}{-.5in} \setlength{\evensidemargin}{\oddsidemargin} \SweaveOpts{keep.source=TRUE, fig=FALSE} % Ross Ihaka suggestions \DefineVerbatimEnvironment{Sinput}{Verbatim} {xleftmargin=2em} \DefineVerbatimEnvironment{Soutput}{Verbatim}{xleftmargin=2em} \DefineVerbatimEnvironment{Scode}{Verbatim}{xleftmargin=2em} \fvset{listparameters={\setlength{\topsep}{0pt}}} \renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}} \SweaveOpts{width=6,height=4} \setkeys{Gin}{width=\textwidth} %\VignetteIndexEntry{Adjusted Survival Curves} <>= options(continue=" ", width=60) options(SweaveHooks=list(fig=function() par(mar=c(4.1, 4.1, .3, 1.1)))) pdf.options(pointsize=8) #text in graph about the same as regular text if (!exists('coxph')) library(survival) @ \title{Adjusted Survival Curves} \author{Terry M Therneau, Cynthia S Crowson, Elizabeth J Atkinson} \date{} % no date \newcommand{\myfig}[1]{\includegraphics[height=!, width=\textwidth] {adjcurve-#1.pdf}} \begin{document} \maketitle \section{Introduction} Suppose we want to investigate to what extent some factor influences survival, as an example we might compare the experience of diabetic patients who are using metformin versus those on injected insulin as their primary treatment modality. There is some evidence that metformin has a positive influence, particularly in cancers, but the ascertainment is confounded by the fact that it is a first line therapy: the patients on metformin will on average be younger and have had a diabetes diagnosis for a shorter amount of time than those using insulin. ``Young people live longer'' is not a particularly novel observation. The ideal way to test this is with a controlled clinical trial. This is of course not always possible, and assessments using available data that includes and adjusts for such confounders is also needed. There is extensive literature --- and debate --- on this topic in the areas of modeling and testing. The subtopic of how to create honest survival curve estimates in the presence of confounders is less well known, and is the focus of this note. Assume that we have an effect of interest $x$ and a set of possible confounding variables $c$. There are two main approaches to adjustment. The first approach, sometimes known as \emph{marginal} analysis, modifies the data so that the confounders $c$ are balanced across the factor of interest $x$, we can then proceed with simple analysis of survival versus $x$ using the reformulated data, ignoring the confounders. The second approach seeks to understand and model the effect of each confounder, with this we can then correct for them. This is often called the \emph{conditional} approach since we are examining the conditional survivals given $x$ and $c$. From these conditional survivals average curves are created that balance on the confounder. As shown below, these differ essentially in the \emph{order} in which the two necessary operations are done, balancing and survival curve creation. \begin{center} \begin{tabular}{rccc} Marginal: & balance data on $c$ & $\longrightarrow$ & form survival curves by $x$ \\ Conditional: & form survival curves by $x$ and $c$& $\longrightarrow$ & average curves over $c$ to balance \end{tabular} \end{center} Many analyses use combinations of these, of course, balancing on some factors and modeling others. All analyses are marginal analyses with respect to important predictors that are unknown to us, and most are conditional with respect to at least a few variables. \begin{figure}[tb] \myfig{flc1} \caption{Survival of \Sexpr{nrow(flchain)} residents of Olmsted County, broken into three cohorts based on FLC value.} \label{flc1} \end{figure} \section{Free Light Chain} Our example data set for this comparison uses a particular assay of plasma immunoglobulins and is based on work of Dr Angela Dispenzieri and her colleagues at the Mayo Clinic \cite{Dispenzieri12}. In brief, plasma cells are responsible for the production of immunoglobulins; they comprise a small portion ($<1$\%) of the total blood and marrow hematapoetic cell population in normal patients. The normal human repertoire is estimated to contain over $10^{8}$ unique immunoglobulins, conferring a broad range of immune protection. In multiple myeloma, the most common form of plasma cell malignancy, almost all of the circulating antigen will be the product of a single malignant clone. An electrophoresis examination of circulating immunoglobulins will exhibit a ``spike'' corresponding to this unique molecule. This anomaly is used both as a diagnostic method and in monitoring the course of the disease under treatment. The presence of a similar, albeit much smaller, spike in normal patients has been a long term research interest of the Mayo Clinic hematology research group \cite{Kyle93}. In 1995 Dr Robert Kyle undertook a population based study of this, and collected serum samples on 19,261 of the 24,539 residents of Olmsted County, Minnesota, aged 50 years or more \cite{Kyle06}. In 2010 Dr. Dispenzieri assayed a sub fraction of the immunoglobulins, the free light chain (FLC), on 15,748 of these subjects who had sufficient remaining sera from the original sample collection. All studies took place under the oversight of the appropriate Institutional Review Boards, which ensure rigorous safety and ethical standards in research. \begin{table} \centering \begin{tabular}{r|cccc} & 50--59 & 60--69 & 70--79 & 80+ \\ \hline <>= group3 <- factor(1+ 1*(flchain$flc.grp >7) + 1*(flchain$flc.grp >9), levels=1:3, labels=c("FLC < 3.38", "3.38 - 4.71", "FLC > 4.71")) age1 <- cut(flchain$age, c(49,59,69,79, 110)) levels(age1) <- c(paste(c(50,60,70), c(59,69,79), sep='-'), '80+') temp1 <- table(group3, age1) temp2 <- round(100* temp1/rowSums(temp1)) pfun <- function(x,y) { paste(ifelse(x<1000, "\\phantom{0}", ""), x, " (", ifelse(y<10, "\\phantom{0}", ""), y, ") ", sep="") } cat(paste(c("FLC low", pfun(temp1[1,], temp2[1,])), collapse=" & "), "\\\\\n") cat(paste(c("FLC med", pfun(temp1[2,], temp2[2,])), collapse=" & "), "\\\\\n") cat(paste(c("FLC high", pfun(temp1[3,], temp2[3,])), collapse=" & "), "\n") @ \end{tabular} \caption{Comparison of the age distributions (percents) for each of the three groups.} \label{tflc1} \end{table} Our data set is an age and sex stratified sample of this second data set and contains 7874 of the subjects from the original study. It is included in the survival package as the data frame \texttt{flchain}. Figure \ref{flc1} shows the survival curves for three subgroups of the patients: those whose total free light chain (FLC) is in the upper 10\% of all values found in the full study, those in the 70--89th percentile, and the remainder. There is a clear survival effect. Average free light chain amounts rise with age, however, at least in part because it is eliminated through the kidneys and renal function declines with age. Table \ref{tflc1} shows the FLC by age distribution. In the highest decile of FLC (group 3) over half the subjects are age 70 or older compared to only 23\% in those below the 70th percentile. How much of the survival difference is truly associated with FLC and how much is simply an artifact of age? The data set contains 3 subjects whose blood sample was obtained on the day of their death. It is rather odd to think of a sample obtained on the final day as ``predicting'' death, or indeed for any results obtained during a patient's final mortality cascade. There are also a few patients with no follow-up beyond the clinic visit at which the assay occurred. Therefore, all the analyses below have excluded subjects with less than 7 days of follow-up. <>= fdata <- flchain[flchain$futime > 7,] fdata$group <- factor(1+ 1*(fdata$flc.grp >7) + 1*(fdata$flc.grp >9), levels=1:3, labels=c("FLC < 3.38", "3.38 - 4.71", "FLC > 4.71")) fdata$age2 <- cut(fdata$age, c(49,54, 59,64, 69,74,79, 89, 110), labels = c(paste(c(50,55,60,65,70,75,80), c(54,59,64,69,74,79,89), sep='-'), "90+")) sfit1 <- survfit(Surv(futime, death) ~ group, fdata) plot(sfit1, mark.time=F, col=c(1,2,4), lty=1, lwd=2, xscale=365.25, xlab="Years from Sample", ylab="Survival") text(c(11.1, 10.5, 7.5), c(.88, .57, .4), c("FLC < 3.38", "3.38 - 4.71", "FLC > 4.71"), col=c(1,2,4)) @ \section{Marginal methods} \subsection{Selection} \begin{figure}[tb] \myfig{flc2} \caption{Survival curves from a balanced subset are shown as solid lines, dashed lines are curves for the unweighted data set (as found in figure \ref{flc1}).} \label{flc2} \end{figure} <>= temp <- with(fdata, table(group, age2, sex)) size <- apply(temp, 2:3, min) set.seed(1978) select <- NULL dd <- dim(temp) for (i in 1:dd[1]) { for (j in 1:dd[2]) { for (k in 1:dd[3]) { indx <- which(as.numeric(fdata$group)==i & as.numeric(fdata$age2) ==j & as.numeric(fdata$sex) ==k) select <- c(select, sample(indx, size[j,k])) } } } data2 <- fdata[select,] sfit2 <- survfit(Surv(futime, death) ~ group, data2) plot(sfit2, mark.time=F, col=c(1,2,4), lty=1, lwd=2, xscale=365.25, xlab="Years from Sample", ylab="Survival") lines(sfit1, mark.time=F, col=c(1,2,4), lty=2, lwd=1, xscale=365.25) legend(2,.4, levels(fdata$group), lty=1, col=c(1,2,4), bty='n', lwd=2) @ \begin{table} \centering \begin{tabular}{crrrrrrrr} \multicolumn{3}{c}{Females} \\ & \multicolumn{8}{c}{Age} \\ FLC group & 50--54& 55--59& 60--64 & 65--69 & 70--74 & 75--79 & 80--89& 90+ \\ \hline <>= tab1 <- with(fdata, table(group, age2, sex)) cat("Low&", paste(tab1[1,,1], collapse=" &"), "\\\\\n") cat("Med&", paste(tab1[2,,1], collapse=" &"), "\\\\\n") cat("High&", paste(tab1[3,,1], collapse=" &"), "\\\\\n") @ \\ \multicolumn{3}{c}{Males} \\ % & 50--54& 55--59& 60--64 & 65--69 & 70--74 & 75--79 & 80--89& 90+ \\ \hline <>= cat("Low&", paste(tab1[1,,2], collapse=" &"), "\\\\\n") cat("Med&", paste(tab1[2,,2], collapse=" &"), "\\\\\n") cat("High&", paste(tab1[3,,2], collapse=" &"), "\n") @ \end{tabular} \caption{Detailed age and sex distribution for the study population} \label{tab2} \end{table} One way to balance the data set is to select a subset of subjects such that the subset is balanced. Table \ref{tab2} is shows an expanded age/sex distribution for the study. The balanced subset has all \Sexpr{tab1[3,1,1]} females aged 50--54 from FLC group 3, a random sample of \Sexpr{tab1[3,1,1]} out of the \Sexpr{tab1[1,1,1]} females in FLC group 1, and \Sexpr{tab1[3,1,1]} out of \Sexpr{tab1[2,1,1]} of the females in FLC group 2 for that age group. The same is done for all of the triplets in the table: select the largest number possible for which we can get perfect balance. \begin{table}[tb] \centering \begin{tabular}{ccccccc} &\multicolumn{2}{c}{FLC low} & \multicolumn{2}{c}{FLC med}& \multicolumn{2}{c}{FLC high} \\ & Total & Subset & Total & Subset & Total & Subset \\ \hline <>= tab3 <- with(fdata, table(age2, group)) tab3 <- round(100*scale(tab3, center=F, scale=colSums(tab3))) tab4 <- with(data2, table(age2, group)) tab4 <- round(100*scale(tab4, center=F, scale=colSums(tab4))) tab5 <- cbind(tab3[,1], tab4[,1], tab3[,2], tab4[,2], tab3[,3], tab4[,3]) pfun <- function(x) paste(ifelse(x<10, paste("\\phantom{0}", x), x), collapse=" &") dtemp <- dimnames(tab5)[[1]] for (j in 1:7) cat(dtemp[j], " &", pfun(tab5[j,]), "\\\\\n") cat(dtemp[8], " & ", pfun(tab5[8,]), "\n") @ \end{tabular} \caption{Age distributions (\%) of the original data set along with that of the subset, for the three FLC groups.} \label{tflc2} \end{table} The survival curves for the balanced subset are shown in figure \ref{flc2}. Two features stand out. First we see that adjustment for age and sex has reduced the apparent survival difference between the groups by about half, but a clinically significant effect for high FLC values remains. The second is that the curve for group 1 has moved a lot while that for group 3 has changed hardly at all. This is a consequence of an overall shift in the distribution of ages, as shown in table \ref{tflc2}. For group 3, the subsetting process has hardly shifted the age distribution, whereas for group 1 it has moved the average age from young to old. The subsetting approach is most often labeled as a ``case-control'' method, and is applied in the situation where one of the groups is small and precious (a rare disease say) and the other group (the controls) is much larger. Since there are excess controls the final age/sex distribution will match that of the cases, and the resulting survival curves have a ready interpretation. Interpretation of the curves for the FLC data is much less clear; we have gained an unbiased comparison, but for a very peculiar population of subjects. <>= # I can't seem to put this all into an Sexpr z1 <- with(fdata,table(age, sex, group)) z2<- apply(z1, 1:2, min) ztemp <- 3*sum(z2) z1b <- with(fdata, table(age>64, sex, group)) ztemp2 <- sum(apply(z1b, 1:2, min)) @ A second deficit of the matching approach is that the choice of our matching criteria is somewhat arbitrary. We want to make the subsets fine enough that subjects are biologically homogeneous within age group, yet large enough so as to limit the number of subjects excluded. For instance, if we had divided the data set only into age 50--64 versus $\ge$ 65, the upper group would still be quite heterogeneous with respect to expected survival (a 70 vs 90 year old for instance) and the distribution of ages within the upper group would remain unbalanced. The resulting data set would have been larger, however, %with \Sexpr{ztemp2} in each group rather than %\Sexpr{nrow(data2)/3}. due to less stringent balancing criteria. On the other hand, if we had done a very fine division, say by age in weeks, the resulting matched data set will be reduced to almost zero; we become hostage to the random fluctuations in sample size for these small groups. In this large data set we have the luxury of dividing finely enough to create acceptably homogeneous groups, but this is not common. @ One advantage of matched subsets is that standard variance calculations for the curves are correct; the values provided by the usual Kaplan-Meier program need no further processing. We can also use the usual statistical tests to check for differences between the curves. <<>>= survdiff(Surv(futime, death) ~ group, data=data2) @ \subsection{Reweighting} \label{sect:logistic} Another way to adjust the data is by reweighting. Let $\pi(a,s)$, $a$ = age group, $s$ = sex be a target population age/sex distribution for our graph, and $p(a,s,i)$ the observed probability of each group/age/sex combination in the data. Each of these will sum to 1. Then if each observation in the data set is given a case weight of \begin{equation} w_{asi} = \frac{\pi(a,s)}{p(a,s,i)} \label{wt1} \end{equation} the weighted age/sex distribution for each of the groups will equal the target distribution $\pi$. When results from multiple studies are to be compared, the target distribution $\pi$ will often be based on external criteria, e.g., the US 2000 population, and the result from the re weighted sample is called the \emph{direct adjustment} method. An obvious advantage of this approach, as compared to subset selection, is that the resulting curves represent a tangible and well defined group. As an example, we will adjust our curves to match the age/sex distribution of the 2000 US population, a common reference target in epidemiology studies. The \texttt{uspop2} data set is found in later releases of the survival package in R. It is an array of counts with dimensions of age, sex, and calendar year. We only want ages of 50 and over, and the population data set has collapsed ages of 100 and over into a single category. We create a table \texttt{tab100} of observed age/sex counts within group for our own data, using the same upper age threshold. New weights are the values $\pi/p$ = \texttt{pi.us/tab100}. <<>>= refpop <- uspop2[as.character(50:100),c("female", "male"), "2000"] pi.us <- refpop/sum(refpop) age100 <- factor(ifelse(fdata$age >100, 100, fdata$age), levels=50:100) tab100 <- with(fdata, table(age100, sex, group))/ nrow(fdata) us.wt <- rep(pi.us, 3)/ tab100 #new weights by age,sex, group range(us.wt) @ There are infinite weights! This is because the US population has coverage at all ages, but our data set does not have representatives in every age/sex/FLC group combination. (There are 32 zeros in tab100). Let us repeat the process, collapsing the US population from single years into the 8 age groups used previously in table \ref{tab2}. Merging the per age/sex/group weights found in the 3-dimensional array \texttt{us.wt} into the data set as per-subject weights uses matrix subscripts, a useful but less known feature of R. <<>>= temp <- as.numeric(cut(50:100, c(49, 54, 59, 64, 69, 74, 79, 89, 110)+.5)) pi.us<- tapply(refpop, list(temp[row(refpop)], col(refpop)), sum)/sum(refpop) tab2 <- with(fdata, table(age2, sex, group))/ nrow(fdata) us.wt <- rep(pi.us, 3)/ tab2 range(us.wt) index <- with(fdata, cbind(as.numeric(age2), as.numeric(sex), as.numeric(group))) fdata$uswt <- us.wt[index] sfit3a <-survfit(Surv(futime, death) ~ group, data=fdata, weight=uswt) @ \begin{figure}[tb] \myfig{flc3a} \caption{Population totals for the US reference (red) and for the observed data set (black).} \label{flc3a} \end{figure} Since the present data set is itself population based and has excellent coverage of the county, it is sensible to use the overall age/sex distribution of the sample itself as our target distribution $\pi$. If we compare the target distribution $\pi$ based on the US population and on the Olmsted County population they are quite similar as shown in figure \ref{flc3a}. Not surprisingly, the population adjusted survival curves based on these two reference populations nearly overlap. <>= tab1 <- with(fdata, table(age2, sex))/ nrow(fdata) matplot(1:8, cbind(pi.us, tab1), pch="fmfm", col=c(2,2,1,1), xlab="Age group", ylab="Fraction of population", xaxt='n') axis(1, 1:8, levels(fdata$age2)) tab2 <- with(fdata, table(age2, sex, group))/nrow(fdata) tab3 <- with(fdata, table(group)) / nrow(fdata) rwt <- rep(tab1,3)/tab2 round(rwt[,1,], 1) #show female data fdata$rwt <- rwt[index] # add per subject weights to the data set sfit3 <- survfit(Surv(futime, death) ~ group, data=fdata, weight=rwt) @ \begin{figure}[tb] \myfig{flc3} \caption{Survival curves for the three groups using reweighted data are shown with solid lines, the original unweighted analysis as dashed lines. The heavier solid line adjusts to the Olmsted population and the lighter one to the US population.} \label{flc3} \end{figure} <>= plot(sfit3, mark.time=F, col=c(1,2,4), lty=1, lwd=2, xscale=365.25, xlab="Years from Sample", ylab="Survival") lines(sfit3a, mark.time=F, col=c(1,2,4), lty=1, lwd=1, xscale=365.25) lines(sfit1, mark.time=F, col=c(1,2,4), lty=2, lwd=1, xscale=365.25) legend(2,.4, levels(fdata$group), lty=1, col=c(1,2,4), bty='n', lwd=2) @ We see that for the low FLC group there are somewhat larger weights for the older ages, whereas the high FLC group requires substantial weights for the youngest ages in order to achieve balance. The resulting survival curve is shown in figure \ref{flc3}. The distance between the adjusted curves is similar to the results from subset selection, which is as expected since both approaches are correcting for the same bias, but results are now for an overall population distribution that matches Olmsted County. The curves estimate what the results would have looked like, had each of the FLC groups contained the full distribution of ages. Estimation based on reweighted data is a common theme in survey sampling. Correct standard errors for the curves are readily computed using methods from that literature, and are available in some software packages. In R the \texttt{svykm} routine in the \texttt{survey} package handles both this simple case and more complex sampling schemes. Tests of the curves can be done using a weighted Cox model; the robust variance produced by \texttt{coxph} is identical to the standard Horvitz-Thompsen variance estimate used in survey sampling \cite{Binder92}. The score test from \texttt{coxph} is equivalent to the log-rank test that is often used to label curves. (In the example below the sykm function is only run if the survey package is already loaded, as the variance calculation is very slow for this large data set.) <<>>= id <- 1:nrow(fdata) cfit <- coxph(Surv(futime, death) ~ group + cluster(id), data=fdata, weight=rwt) summary(cfit)$robscore if (exists("svykm")) { #true if the survey package is loaded sdes <- svydesign(id = ~0, weights=~rwt, data=fdata) dfit <- svykm(Surv(futime, death) ~ group, design=sdes, se=TRUE) } @ Note: including the \texttt{cluster} term in the coxph call causes it to treat the weights as resampling values and thus use the proper survey sampling style variance. The default without that term would be to treat the case weights as replication counts. (This same alternate variance estimate is also called for when there are correlated observations; many users will be more familiar with the cluster statement in that context.) \paragraph{Inverse probability weighting} Notice that when using the overall population as the target distribution $\pi$ we can use Bayes rule to rewrite the weights as \begin{align*} \frac{1}{w_{asi}} &= \frac{{\rm Pr}({\rm age}=a, {\rm sex} =s, {\rm group}=i)} {{\rm Pr}({\rm age}=a, {\rm sex} =s)} \\ &= {\rm Pr}({\rm group}=i | {\rm age}=a, {\rm sex} =s) \end{align*} This last is precisely the probability estimated by a logistic regression model, leading to \emph{inverse probability weighting} as an alternate label for this approach. We can reproduce the weights calculated just above with three logistic regression models. <>= options(na.action="na.exclude") gg <- as.numeric(fdata$group) lfit1 <- glm(I(gg==1) ~ factor(age2) * sex, data=fdata, family="binomial") lfit2 <- glm(I(gg==2) ~ factor(age2) * sex, data=fdata, family="binomial") lfit3 <- glm(I(gg==3) ~ factor(age2) * sex, data=fdata, family="binomial") temp <- ifelse(gg==1, predict(lfit1, type='response'), ifelse(gg==2, predict(lfit2, type='response'), predict(lfit3, type='response'))) all.equal(1/temp, fdata$rwt) @ If there were only 2 groups then only a single regression model is needed since P(group 2) = 1 - P(group 1). Note the setting of na.action, which causes the predicted vector to have the same length as the original data even when there are missing values. This simplifies using the derived weights with the original data set. An advantage of the regression framework is that one can easily accommodate more variables by using a model with additive terms and only a few selected interactions, and the model can contain continuous as well as categorical predictors. The disadvantage is that such models are often used without the necessary work to check their validity. For instance models with \texttt{age + sex} could have been used above. This makes the assumption that the odds of being a member of group 1 is linear in age and with the same slope for males and females; ditto for the models for group 2 and group 3. How well does this work? Since the goal of reweighting is to standardize the ages, a reasonable check is to compute and plot the reweighted age distribution for each flc group. \begin{figure}[tb] \myfig{flc4} \caption{The re-weighted age distribution using logistic regression with continuous age, for females, FLC groups 1--3. The target distribution is shown as a ``+''. The original unadjusted distribution is shown as dashed lines.} \label{flc4} \end{figure} Figure \ref{flc4} shows the result. The reweighted age distribution is not perfectly balanced, i.e., the 123 symbols do not over plot one another, but in this case the simple linear model has done an excellent job. We emphasize that whenever the reweighting is based on a simplified model then such a check is obligatory. It is quite common that a simple model is not sufficient and the resulting weight adjustment is inadequate. Like a duct tape auto repair, proceeding forward as though the underlying problem has been addressed is then most unwise. <>= lfit1b <-glm(I(gg==1) ~ age + sex, data=fdata, family="binomial") lfit2b <- glm(I(gg==2) ~ age +sex, data=fdata, family="binomial") lfit3b <- glm(I(gg==3) ~ age + sex, data=fdata, family="binomial") # weights for each group using simple logistic twt <- ifelse(gg==1, 1/predict(lfit1b, type="response"), ifelse(gg==2, 1/predict(lfit2b, type="response"), 1/predict(lfit3b, type="response"))) tdata <- data.frame(fdata, lwt=twt) #grouped plot for the females temp <- tdata[tdata$sex=='F',] temp$gg <- as.numeric(temp$group) c1 <- with(temp[temp$gg==1,], tapply(lwt, age2, sum)) c2 <- with(temp[temp$gg==2,], tapply(lwt, age2, sum)) c3 <- with(temp[temp$gg==3,], tapply(lwt, age2, sum)) xtemp <- outer(1:8, c(-.1, 0, .1), "+") #avoid overplotting ytemp <- 100* cbind(c1/sum(c1), c2/sum(c2), c3/sum(c3)) matplot(xtemp, ytemp, col=c(1,2,4), xlab="Age group", ylab="Weighted frequency (%)", xaxt='n') ztab <- table(fdata$age2) points(1:8, 100*ztab/sum(ztab), pch='+', cex=1.5, lty=2) # Add the unadjusted temp <- tab2[,1,] temp <- scale(temp, center=F, scale=colSums(temp)) matlines(1:8, 100*temp, pch='o', col=c(1,2,4), lty=2) axis(1, 1:8, levels(fdata$age2)) @ \paragraph{Rescaled weights} As the weights were defined in equation \ref{wt1}, the sum of weights for each of the groups is \Sexpr{nrow(fdata)}, the number of observations in the data set. Since the number of subjects in group 3 is one seventh of that in group 1, the average weight in group 3 is much larger. An alternative is to define weights in terms of the \emph{within} group distribution rather than the overall distribution, leading to the rescaled weights $w^*$ \begin{align} w^* &= \frac{\pi(a,s)}{p(a,s|i)} \label{wt2} \\ &= \frac{{\rm P}({\rm group}=i)} {{\rm P}({\rm group}=i | {\rm age}=a, {\rm sex}=s)} \label{wt2b} \end{align} Each group's weights are rescaled by the overall prevalence of the group. In its simplest form, the weights in each group are scaled to add up to the number of subjects in the group. <<>>= # compute new weights wtscale <- table(fdata$group)/ tapply(fdata$rwt, fdata$group, sum) wt2 <- c(fdata$rwt * wtscale[fdata$group]) c("rescaled cv"= sd(wt2)/mean(wt2), "rwt cv"=sd(fdata$rwt)/mean(fdata$rwt)) cfit2a <- coxph(Surv(futime, death) ~ group + cluster(id), data=fdata, weight= rwt) cfit2b <- coxph(Surv(futime, death) ~ group + cluster(id), data=fdata, weight=wt2) round(c(cfit2a$rscore, cfit2b$rscore),1) @ The rescaling results in weights that are much less variable across groups. This operation has no impact on the individual survival curves or their standard errors, since within group we have multiplied all weights by a constant. When comparing curves across groups, however, the rescaled weights reduce the standard error of the test statistic. This results in increased power for the robust score test, although in this particular data set the improvement is not very large. \section{Conditional method} In the marginal approach we first balance the data set and then compute results on the adjusted data. In the conditional approach we first compute a predicted survival curve for each subject that accounts for flc group, age and sex, and then take a weighted average of the curves to get an overall estimate for each flc group. For both methods a central consideration is the population of interest, which drives the weights. Modeling has not removed the question of \emph{who} these curves should represent, it has simply changed the order of operation between the weighting step and the survival curves step. \subsection{Stratification} Our first approach is to subset the data into homogeneous age/sex strata, compute survival curves within each strata, and then combine results. We will use the same age/sex combinations as before. The interpretation of these groups is different, however. In the marginal approach it was important to find age/sex groups for which the probability of membership within each FLC group was constant within the strata (independent of age and sex, within strata), in this case it is important that the survival for each FLC group is constant in each age/sex stratum. Homogeneity of membership within each stratum and homogeneity of survival within each stratum may lead to different partitions for some data sets. Computing curves for all the combinations is easy. <>= allfit <- survfit(Surv(futime/365.25, death) ~ group + age2 + sex, fdata) temp <- summary(allfit)$table temp[1:6, c(1,4)] #abbrev printout to fit page @ The resultant survival object has 48 curves: 8 age groups * 2 sexes * 3 FLC groups. To get a single curve for the first FLC group we need to take a weighted average over the 16 age/sex combinations that apply to that group, and similarly for the second and third FLC subset. Combining the curves is a bit of a nuisance computationally because each of them is reported on a different set of time points. A solution is to use the \texttt{summary} function for survfit objects along with the \texttt{times} argument of that function. This feature was originally designed to allow printout of curves at selected time points (6 months, 1 year, \ldots), but can also be used to select a common set of time points for averaging. We will arbitrarily use 4 per year, which is sufficient to create a visually smooth plot over the time span of interest. By default \texttt{summary} does not return data for times beyond the end of a curve, i.e., when there are no subjects left at risk; the \texttt{extend} argument causes a full set of times to always be reported. As seen in the printout above, the computed curves are in sex within age within group order. The overall curve is a weighted average chosen to match the original age/sex distribution of the population. <>= xtime <- seq(0, 14, length=57) #four points/year for 14 years smat <- matrix(0, nrow=57, ncol=3) # survival curves serr <- smat #matrix of standard errors pi <- with(fdata, table(age2, sex))/nrow(fdata) #overall dist for (i in 1:3) { temp <- allfit[1:16 + (i-1)*16] #curves for group i for (j in 1:16) { stemp <- summary(temp[j], times=xtime, extend=T) smat[,i] <- smat[,i] + pi[j]*stemp$surv serr[,i] <- serr[,i] + pi[i]*stemp$std.err^2 } } serr <- sqrt(serr) matplot(xtime, smat, type='l', lwd=2, col=c(1,2,4), ylim=c(0,1), lty=1, xlab="Years from sample", ylab="Survival") lines(sfit1, mark.time=F, lty=2, col=c(1,2,4), xscale=365.25) @ \begin{figure}[tb] \myfig{flc5} \caption{Estimated curves from a stratified model, along with those from the uncorrected fit as dashed lines.} \label{flc5} \end{figure} Figure \ref{flc5} shows the resulting averaged curves. Overlaid are the curves for the unadjusted model. Very careful comparison of these curves with the weighted estimate shows that they have almost identical spread, with just a tiny amount of downward shift. There are two major disadvantages to the stratified curves. The first is that when the original data set is small or the number of confounders is large, it is not always feasible to stratify into a large enough set of groups that each will be homogeneous. The second is a technical aspect of the standard error estimate. Since the curves are formed from disjoint sets of observations they are independent and the variance of the weighted average is then a weighted sum of variances. However, when a Kaplan-Meier curve drops to zero the usual standard error estimate at that point involves 0/0 and becomes undefined, leading to the NaN (not a number) value in R. Thus the overall standard error becomes undefined if any of the component curves falls to zero. In the above example this happens at about the half way point of the graph. (Other software packages carry forward the se value from the last no-zero point on the curve, but the statistical validity of this is uncertain.) To test for overall difference between the curves we can use a stratified test statistic, which is a sum of the test statistics computed within each subgroup. The most common choice is the stratified log-rank statistic which is shown below. The score test from a stratified Cox model would give the same result. <<>>= survdiff(Surv(futime, death) ~ group + strata(age2, sex), fdata) @ \subsection{Modeling} The other approach for conditional estimation is to model the risks due to the confounders. Though we have left it till last, this is usually the first (and most often the only) approach used by most data analysts. Let's start with the very simplest method: a stratified Cox model. <>= cfit4a <- coxph(Surv(futime, death) ~ age + sex + strata(group), data=fdata) surv4a <- survfit(cfit4a) plot(surv4a, col=c(1,2,4), mark.time=F, xscale=365.25, xlab="Years post sample", ylab="Survival") @ This is a very fast and easy way to produce a set of curves, but it has three problems. First is the assumption that this simple model adequately accounts for the effects of age and sex on survival. That is, it assumes that the effect of age on mortality is linear, the sex difference is constant across all ages, and that the coefficients for both are identical for the three FLC groups. The second problem with this approach is that it produces the predicted curve for a single hypothetical subject of age \Sexpr{round(cfit4a[['means']][1], 1)} years and sex \Sexpr{round(cfit4a[['means']][2],2)}, the means of the covariates, under each of the 3 FLC scenarios. However, we are interested in the adjusted survival of a \emph{cohort} of subjects in each range of FLC, and the survival of an ``average'' subject is not the average survival of a cohort. The third and most serious issue is that it is not clear exactly what these ``adjusted'' curves represent --- exactly who \emph{is} this subject a sex of \Sexpr{round(cfit4a[['means']][2],2)}? Multiple authors have commented on this problem, see Thomsen et al \cite{Thomsen91}, Nieto and Coresh \cite{Nieto96} or chapter 10 of Therneau and Grambsh \cite{Therneau00} for examples. Even worse is a Cox model that treated the FLC group as a covariate, since that will impose a an additional constraint of proportional hazards across the 3 FLC groups. \begin{figure} \myfig{flc6} \caption{Curves for the three groups, adjusted for age and sex via a risk model. Dotted lines show the curves from marginal adjustment. Solid curves are for the simple risk model \texttt{cfit4a}.} \label{flc6} \end{figure} We can address this last problem problem by doing a proper average. A Cox model fit can produce the predicted curves for any age/sex combination. The key idea is to produce a predicted survival curve for every subject of some hypothetical population, and then take the average of these curves. The most straightforward approach is to retrieve the predicted individual curves for all \Sexpr{nrow(fdata)} subjects in the data set, assuming each of the three FLC strata one by one, and take a simple average for each strata. For this particular data set that is a bit slow since it involves \Sexpr{nrow(fdata)} curves. However there are only 98 unique age/sex pairs in the data, it is sufficient to obtain the 98 * 3 FLC groups unique curves and take a weighted average. We will make use of the survexp function, which is designed for just this purpose. Start by creating a data set which has one row for each age/sex combination along with its count. Then replicate it into 3 copies, assigning one copy to each of the three FLC strata. <>= tab4a <- with(fdata, table(age, sex)) uage <- as.numeric(dimnames(tab4a)[[1]]) tdata <- data.frame(age = uage[row(tab4a)], sex = c("F","M")[col(tab4a)], count= c(tab4a)) tdata3 <- tdata[rep(1:nrow(tdata), 3),] #three copies tdata3$group <- factor(rep(1:3, each=nrow(tdata)), labels=levels(fdata$group)) sfit4a <- survexp(~group, data=tdata3, weight = count, ratetable=cfit4a) plot(sfit4a, mark.time=F, col=c(1,2,4), lty=1, lwd=2, xscale=365.25, xlab="Years from Sample", ylab="Survival") lines(sfit3, mark.time=F, col=c(1,2,4), lty=2, lwd=1, xscale=365.25) legend(2,.4, c("FLC low", "FLC med", "FLC high"), lty=1, col=c(1,2,4), bty='n', lwd=2) @ Figure \ref{flc6} shows the result. Comparing this to the prior 3 adjustments shown in figures \ref{flc3}, \ref{flc4}, and \ref{flc5} we see that this result is different. Why? Part of the reason is due to the fact that $E[f(X)] \ne f(E[X])$ for any non-linear operation $f$, so that averages of survival curves and survival curves of averages will never be exactly the same. This may explain the small difference between the stratified and the marginal approaches of figures \ref{flc3} and \ref{flc5}, which were based on the same subsets. The Cox based result is systematically higher than the stratified one, however, so something more is indicated. Aside: An alternate computational approach is to create the individual survival curves using the \texttt{survfit} function and then take averages. <<>>= tfit <- survfit(cfit4a, newdata=tdata, se.fit=FALSE) curves <- vector('list', 3) twt <- c(tab4a)/sum(tab4a) for (i in 1:3) { temp <- tfit[i,] curves[[i]] <- list(time=temp$time, surv= c(temp$surv %*% twt)) } @ The above code is a bit sneaky. I know that the result from the survfit function contains a matrix \texttt{tfit\$surv} of 104 columns, one for each row in the tdata data frame, each column containing the curves for the three strata one after the other. Sub setting \texttt{tfit} results in the matrix for a single flc group. Outside of R an approach like the above may be needed, however. \begin{figure} \myfig{flc6b} \caption{Left panel: comparison of Cox model based adjustment (solid) with the curves based on marginal adjustment (dashed). The Cox model curves without (black) and with (red) an age*sex interaction term overlay. Right panel: plot of the predicted relative risks from a Cox model \texttt{crate} versus population values from the Minnesota rate table.} \label{flc6b} \end{figure} So why are the modeling results so different than either reweighting or stratification? Suspicion first falls on the use of a simple linear model for age and sex, so start by fitting a slightly more refined model that allows for a different slope for the two sexes, but is still linear in age. In this particular data set an external check on the fit is also available via the Minnesota death rate tables, which are included with the survival package as \texttt{survexp.mn}. This is an array that contains daily death rates by age, sex, and calender year. <>= par(mfrow=c(1,2)) cfit4b <- coxph(Surv(futime, death) ~ age*sex + strata(group), fdata) sfit4b <- survexp(~group, data=tdata3, ratetable=cfit4b, weights=count) plot(sfit4b, fun='event', xscale=365.25, xlab="Years from sample", ylab="Deaths") lines(sfit3, mark.time=FALSE, fun='event', xscale=365.25, lty=2) lines(sfit4a, fun='event', xscale=365.25, col=2) temp <- median(fdata$sample.yr) mrate <- survexp.mn[as.character(uage),, as.character(temp)] crate <- predict(cfit4b, newdata=tdata, reference='sample', type='lp') crate <- matrix(crate, ncol=2)[,2:1] # mrate has males then females, match it # crate contains estimated log(hazards) relative to a baseline, # and mrate absolute hazards, make both relative to a 70 year old for (i in 1:2) { mrate[,i] <- log(mrate[,i]/ mrate[21,2]) crate[,i] <- crate[,i] - crate[21,2] } matplot(mrate, crate, col=2:1, type='l') abline(0, 1, lty=2, col=4) @ The resulting curves are shown in the left panel of figure \ref{flc6b} and reveal that addition of an interaction term did not change the predictions, and that the Cox model result for the highest risk group is distinctly different predicted survival for the highest FLC group is distinctly different when using model based prediction. The right hand panel of the figure shows that though there are slight differences with the Minnesota values, linearity of the age effect is very well supported. So where exactly does the model go wrong? Since this is such a large data set we have the luxury of looking at subsets. This would be a very large number of curves to plot --- age by sex by FLC = 48 --- so an overlay of the observed and expected curves by group would be too confusing. Instead we will summarize each of the groups according to their observed and predicted number of events. <>= obs <- with(fdata, tapply(death, list(age2, sex, group), sum)) pred<- with(fdata, tapply(predict(cfit4b, type='expected'), list(age2, sex, group), sum)) excess <- matrix(obs/pred, nrow=8) #collapse 3 way array to 2 dimnames(excess) <- list(dimnames(obs)[[1]], c("low F", "low M", "med F", "med M", "high F", "high M")) round(excess, 1) @ The excess risks, defined as the observed/expected number of deaths, are mostly modest ranging from .8 to 1.2. The primary exception exception is the high FLC group for ages 50--59 which has values of 1.6 to 2.5; the Cox model fit has greatly overestimated the survival for the age 50--54 and 55--59 groups. Since this is also the age category with the highest count in the data set, this overestimation will have a large impact on the overall curve for high FLC subset, which is exactly where the the deviation in figure \ref{flc6b} is observed to lie. There is also mild evidence for a linear trend in age for the low FLC females, in the other direction. Altogether this suggests that the model might need to have a different age coefficient for each of the three FLC groups. <<>>= cfit5a <- coxph(Surv(futime, death) ~ group:age +sex + strata(group), fdata) cfit5b <- coxph(Surv(futime, death) ~ group:(age +sex) + strata(group), fdata) cfit5c <- coxph(Surv(futime, death) ~ group:(age *sex) + strata(group), fdata) options(show.signif.stars=FALSE) # see footnote anova(cfit4a, cfit5a, cfit5b, cfit5c) temp <- coef(cfit5a) names(temp) <- c("sex", "ageL", "ageM", "ageH") round(temp,3) @ The model with separate age coefficients for each FLC group gives a major improvement in goodness of fit, but adding separate sex coefficients per group or further interactions does not add importantly beyond that. \footnote{There are certain TV shows that make one dumber just by watching them; adding stars to the output has the same effect on statisticians.} \begin{figure} \myfig{flc7} \caption{Adjusted survival for the 3 FLC groups based on the improved Cox model fit. Dashed lines show the predictions from the marginal model.} \label{flc7} \end{figure} A recheck of the observed/expected values now shows a much more random pattern, though some excess remains in the upper right corner. The updated survival curves are shown in figure \ref{flc7} and now are in closer concordance with the marginal fit. <>= pred5a <- with(fdata, tapply(predict(cfit5a, type='expected'), list(age2, sex, group), sum)) excess5a <- matrix(obs/pred5a, nrow=8, dimnames=dimnames(excess)) round(excess5a, 1) sfit5 <- survexp(~group, data=tdata3, ratetable=cfit5a, weights=count) plot(sfit3, fun='event', xscale=365.25, mark.time=FALSE, lty=2, xlab="Years from sample", ylab="Deaths") lines(sfit5, fun='event', xscale=365.25) @ One problem with the model based estimate is that standard errors for the curves are complex. Standard errors of the individual curves for each age/sex/FLC combination are a standard output of the survfit function, but the collection of curves is correlated since they all depend on a common estimate of the model's coefficient vector $\beta$. Curves with disparate ages are anti-correlated (an increase in the age coefficient of the model would raise one and lower the other) whereas those for close ages are positively correlated. A proper variance for the unweighted average has been derived by Gail and Byar \cite{Gail86}, but this has not been implemented in any of the standard packages, nor extended to the weighted case. A bootstrap estimate would appear to be the most feasible. \section{Conclusions} When two populations need to be adjusted and one is much larger than the other, the balanced subset method has been popular. It is most often seen in the context of a case-control study, with cases as the rarer group and a set of matched controls selected from the larger one. This method has the advantage that the usual standard error estimates from a standard package are appropriate, so no further work is required. However, in the general situation it leads to a correct answer but for the wrong problem, i.e., not for a population in which we are interested. The population reweighted estimate is flexible, has a readily available variance in some statistical packages (but not all), and the result is directly interpretable. It is the method we recommend in general. The approach can be extended to a large number of balancing factors by using a regression model to derive the weights. Exploration and checking of said model for adequacy is an important step in this case. The biggest downside to the method arises when there is a subset which is rare in the data sample but frequent in the adjusting population. In this case subjects in that subset will be assigned large weights, and the resulting curves will have high variance. The stratified method is closely related to reweighting (not shown). It does not do well if the sample size is small, however. Risk set modeling is a very flexible method, but is also the one where it is easiest to go wrong by using an inadequate model, and variance estimation is also difficult. To the extent that the fitted model is relevant, it allows for interpolation and extrapolation to a reference population with a different distribution of covariates than the one in the training data. It may be applicable in cases such as rare subsets where population reweighting is problematic, with the understanding that one is depending heavily on extrapolation in this case, which is always dangerous. \bibliographystyle{plain} \bibliography{refer} \end{document}