The goal is to plot separate survival curves for each level of a categorical variable in R, using a coxph model including the variable and a second categorical variable with >2 levels.
If the second categorical variable had had only 2 levels, it is possible to give the proportion of one level as the reference value. This is frowned upon by the R community, but done by default in some other statistical software.
An example using the mgus dataset:
library(survival) mgustest <- mgus mgustest$agecat <- cut(mgustest$age, breaks=quantile(mgustest$age)) mgustest$albcat <- cut(mgustest$alb, breaks=quantile(mgustest$alb, na.rm=T)) mfit.refs <- data.frame("sex"=mean(mgustest$sex == "female"), "albcat"=c('(1.8,3]','(3,3.2]','(3.2,3.5]','(3.5,5.1]')) mfit <- coxph(Surv(futime, death) ~ sex + albcat, data=mgustest) plot(survfit(mfit, newdata=mfit.refs), col=c("black","red","blue","green"))
The newdata data frame looks like this:
sex albcat 1 0.4190871 (1.8,3] 2 0.4190871 (3,3.2] 3 0.4190871 (3.2,3.5] 4 0.4190871 (3.5,5.1]
Assuming the model:
mfit2 <- coxph(Surv(futime, death) ~ agecat + albcat, data=mgustest)
How should I build the newdata data frame to be able to fit separate survival curves for albcat levels with survfit (I only want four curves representing the albact levels)?
Thoughts: How do other statistical software produce these survival curves without specifying reference values for agecat? For example, how does SPSS do it?
Disclaimer: I agree that it makes little sense to assume that every individual in a sample has the sex 0.42, when the only occurring values are 0 and 1. Still, this assumption is common in other statistical software, and it is useful to be able to replicate models resulting from such assumptions in R. I have googled for a couple of hours about the matter, but could not find an answer about how to deal with categorical variables with >2 levels.
Edit: I want to clarify based on the first comment. I do not want 4*4=16 survival curves like in this example:
mfit.refs2 <- data.frame( "agecat"=c(rep('(34,55]', 8),rep('(55,64]',8),rep('(64,72]', 8),rep('(72,90]', 8)), "albcat"=rep(c('(1.8,3]','(3,3.2]','(3.2,3.5]','(3.5,5.1]'), 4)) mfit2 <- coxph(Surv(futime, death) ~ agecat + albcat, data=mgustest) plot(survfit(mfit2, newdata=mfit.refs2), col=c("black","red","blue","green"))