From: Peter Jepsen <PJ_at_dce.au.dk>

Date: Thu, 31 Jan 2008 00:21:00 +0100

***

-----Oprindelig meddelelse-----

Fra: r-help-bounces_at_r-project.org [mailto:r-help-bounces_at_r-project.org] På vegne af Thomas Lumley Sendt: 30. januar 2008 00:31

Til: Peter Jepsen

Cc: r-help

Emne: Re: [R] Direct adjusted survival?

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Wed 30 Jan 2008 - 23:24:15 GMT

Date: Thu, 31 Jan 2008 00:21:00 +0100

Thank you for your reply, Thomas. I'm not quite sure whether survexp() does that. It seems that the idea of survexp() is to take the ratetable from a mortality table or Cox model based on one dataset and apply it to another dataset. I'm trying to adjust for confounding, so I want to take the ratetable from a Cox model based on one dataset and apply it to the SAME dataset. Here's an example of how I try to achieve this:

require(survival)

data(pbc)

## compare data to Cox model

## fit to randomised patients in Mayo PBC data
coxph(Surv(time,status)~edtrt,data=pbc)

m<-coxph(Surv(time,status)~edtrt+log(bili),data=pbc)
m # log(bili) is a strong confounder

plot(survfit(Surv(time,status)~edtrt,data=pbc))
lines(survexp(~edtrt+ratetable(edtrt=edtrt,bili=bili),data=pbc,ratetable=m,cohort=TRUE),col="purple")

The lines that I hoped to be the survival probabilities for each edtrt-group adjusted for confounding by log(bili) are nearly identical to the KM-lines, and they certainly don't appear adjusted for the very strong confounding by log(bili). I'm not quite sure what they are, though.

Ghali et al. claim to have an S-plus implementation of the 'direct adjusted survival' method (Ghali WA, Quan H, Brant R, van Melle G, Norris CM, Faris PD, et al. Comparison of 2 methods for calculating adjusted survival curves from proportional hazards models. JAMA 2001;286:1494-1497). I have found the function here: http://stat.ubc.ca/~rollin/stats/S/surv.html. It is inserted below, but please note that I have made one modification.

I'm still very new to R, so I don't follow exactly what happens. It seems that avg.surv() wants edtrt as a factor that takes integer values?! I realize that this is changes the Cox model specification, but, anyway, this code produces a result that is much closer what I expected:

pbc$edtrt.fac<-as.factor(pbc$edtrt*2)

m2<-coxph(Surv(time,status)~edtrt.fac+log(bili),data=pbc)
fits<-avg.surv(m2, var.name="edtrt.fac", data=pbc, var.values=c("0","1","2"))
matlines(fits$time,fits$fits)

However, avg.surv() does not provide standard errors, hence my question regarding the Zhang paper. If anyone can help me sort out what is going on, I'd be very thankful.

Best regards,

Peter.

***

avg.surv <- function(cfit, var.name, var.values, data, weights)
{

if(missing(data)) { if(!is.null(cfit$model)) mframe <- cfit$model else mframe <- model.frame(cfit, sys.parent()) } else mframe <- model.frame(cfit, data) var.num <- match(var.name, names(mframe)) data.patterns <- apply(data.matrix(mframe[, - c(1, var.num)]), 1, paste, collapse = ",") data.patterns <- factor(data.patterns,levels=unique(data.patterns)) if(missing(weights)) weights <- table(data.patterns) else weights <- tapply(weights, data.patterns, sum) kp <- !duplicated(data.patterns) mframe <- mframe[kp,] obs.var <- mframe[,var.num] lps <- (cfit$linear.predictor)[kp] tframe <- mframe[rep(1,length(var.values)),] tframe[,var.num] <- var.values xmat <- model.matrix(cfit,data=tframe)[,-1] tlps <- as.vector(xmat%*%cfit$coef) names(tlps) <- var.values obs.off <- tlps[as.character(obs.var)] explp.off <- exp(outer(lps - obs.off ,tlps,"+")) bfit <- survfit(cfit, se.fit = F) # Changed from "survfit.coxph" to "survfit" fits <- outer(bfit$surv,explp.off,function(x,y) x^y) avg.fits <- apply(sweep(fits,2,weights,"*"),c(1,3),sum)/sum(weights) dimnames(avg.fits) <- list(NULL,var.values) list(time=bfit$time,fits=avg.fits)}

***

-----Oprindelig meddelelse-----

Fra: r-help-bounces_at_r-project.org [mailto:r-help-bounces_at_r-project.org] På vegne af Thomas Lumley Sendt: 30. januar 2008 00:31

Til: Peter Jepsen

Cc: r-help

Emne: Re: [R] Direct adjusted survival?

On Wed, 30 Jan 2008, Peter Jepsen wrote:

*>
*

> I am trying to find an R function to compute 'direct adjusted survival'

*> with standard errors. A SAS-macro to do this is presented in Zhang X,
**> Loberiza FR, Klein JP, Zhang MJ. A SAS macro for estimation of direct
**> adjusted survival curves based on a stratified Cox regression model.
**> Comput Methods Programs Biomed 2007;88:95-101. It appears that this
**> method is not implemented in R. Can anyone prove me wrong?
**>
*

This looks like what survexp() does. It's hard to be sure, since I can only find the abstract online.

-thomas

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Wed 30 Jan 2008 - 23:24:15 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Wed 30 Jan 2008 - 23:30:09 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*