Well, I would suggest using the code already in place in the survival
package. Here is my code for your problem.
I'm using a copy of the larynx data as found from the web resources for
the Klein and Moeschberger book.
larynx <- read.table("larynx.dat", skip=12,
col.names=c("stage", "time", "age", "year", "death"))
larynx$stage <- factor(larynx$stage)
larynx$age.grp <- cut(larynx$age, c(0, 55, 70, 90),
labels=c("40-55", "55-70", "70-90"))
fit1 <- coxph(Surv(time, death) ~ age.grp + stage + I(year >=74),
larynx) #fit a Cox model
km <- survfit(Surv(time, death) ~ stage, data=larynx) # KM
plot(km)
direct <- survexp( ~stage, data=larynx, ratetable=fit1) # "direct
adjusted" survival
lines(direct, col=2)
---
A few further comments
1. It would help those of us who advise if you would simplify your
code. I didn't need to see 20+ lines of data set manipulation and curve
drawing. Give the essence of the problem.
2. The "direct" estimate above was first created for expected
survival curves based on population rate tables and publised by Edererer
in 1961. The same idea was rediscovered in the context of the Cox model
and renamed the "direct adjusted" estimate; I like to give credit to the
oridingal.
3. I did not try to debug your function.
Terry Therneau
--- begin included message ---
Hello R users,
I am trying to obtain a direct adjusted survival curve. I am sending my
whole code (see below). It's basically the larynx cancer data with Stage
1-4.I am using the cox model using coxph option, see the fit3 coxph.
When I use the avg.surv option on fit3, I get the following error:
"fits<-avg.surv(fit3, var.name="stage.fac", var.values=c(1,2,3,4),
data=larynx)
Error in `contrasts<-`(`*tmp*`, value = contrasts.arg[[nn]]) :
contrasts apply only to factors"
If try using var.values = c ("1","2","3","4") option, I get the
following error:
"Error in xmat %*% cfit$coef : non-conformable arguments
In addition: Warning message:
In model.matrix.default(Terms, mf, contrasts = contrast.arg) :
variable 'stage.fac' converted to a factor"
Please advise me on what I am doing wrong!
Regards,
Manish Dalwani
University of Colorado
library(survival)
library(ggplot2)
setwd("/Users/manishdalwani/surv_6646/")
#file.show("Larynx.txt")
#Stage of disease (1=stage 1, 2=stage2, 3=stage 3, 4=stage 4)
#Time to death or on-study time, months
#Age at diagnosis of larynx cancer
#Year of diagnosis of larynx cancer
#Death indicator (0=alive, 1=dead)
larynx <-read.table(file ="larynx.dat",col.names
=c("stage","time","age","year","death"),colClasses
=c("factor","numeric","numeric","numeric","numeric"))
#the death indicator needs to be numeric for the survfit function to work
# splitting age into three intervals: <55, 50-70, >70
age.larynx <- rep(0, length(larynx[,1]))
age.larynx [larynx$age < 55] <- 1
age.larynx [larynx$age >= 55 & larynx$age <= 70] <- 2
age.larynx [larynx$age >70] <- 3
larynx <- cbind(larynx, age.larynx)
larynx$age.larynx <- factor(larynx$age.larynx, labels=c("1", "2", "3"))
larynx$stage.fac<-as.factor(larynx$stage)
year.larynx <- rep(0, length(larynx[,1]))
year.larynx [larynx$year >= 74] <- 1
year.larynx [larynx$year < 74] <- 2
larynx <- cbind(larynx, year.larynx)
head(larynx)
summary(larynx)
kapmeier.fit <- survfit( Surv(time, death) ~ stage, data = larynx, type
= "kaplan-meier")
sumkap <- summary(kapmeier.fit)
attributes(kapmeier.fit)
attributes(sumkap)
plot(kapmeier.fit, lty=2:3, fun="cumhaz",
xlab="Months",ylab="Cumulative Hazard of death")
legend(4,.4,c("stage1","stage2","stage3","stage4"),lty=1:2)
plot(kapmeier.fit, lty=2:3,
xlab="Months",ylab="Survival Function")
#construct a data frame for the plots
plotdata <- data.frame(time = sumkap$time, surv = sumkap$surv, strata =
sumkap$strata)
fact.stage<-factor(larynx$stage)
fit1<-coxph(Surv(time, death) ~ fact.stage + age.larynx +
factor(year>=74), data=larynx, method=c("efron"))
summary(fit1)
avg.surv <- function(cfit, var.name, var.values, data, weights)
{
if(missing(data)){
if(!is.null(cfit$model))
mframe <- cfit$model
elsemframe <-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.coxph(cfit, se.fit = F)
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)
}
fits<-avg.surv(fit1, var.name="fact.stage", data=larynx,
var.values=c("0","1","2","3"))
matlines(fits$time,fits$fits)
fit2<-coxph(Surv(time, death) ~ stage + age.larynx + factor(year>=74),
data=larynx, method=c("efron"))
summary(fit2)
fit3<-coxph(Surv(time, death) ~ stage.fac + age.larynx + year.larynx),
data=larynx, method=c("efron"))
summary(fit3)
fit4<-coxph(Surv(time,death)~factor(stage)+factor(age>=74)+factor(year>=74),data=larynx,method=c("efron"))
summary(fit4)
lines(survexp(~stage+ratetable(stage,age.larynx,year>=74),data=larynx,ratetable=fit2,cohort=TRUE),col="purple")
fits<-avg.surv(fit3, var.name="stage.fac", var.values=c(1,2,3,4),
data=larynx)
matlines(fits$time,fits$fits)
______________________________________________
R-help@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.
