Search code examples
rsurvival-analysis

Creating a grouping variable for a survival model in R


I am trying to plot a KM curve in R but first I need to fit the survival object. I have a dataset that consists of 100 rows where each row corresponds to a patient who is either in group A or group B. What I would like to do is to be able to plot (on the same plot) KM curves for group A versus group B versus group A+B (so everyone). The trouble that I am having is figuring out how to construct the group variable. I assume you can't do it in a single variable and so this is what I am trying, although it doesn't seem to be working correctly (I don't get everyone in group A and B).

set.seed(4)

n = 100
x = runif(n,0,1500)
y = runif(n,0,5)
survival = runif(n,1,1000)
censor = rbinom(n,1,.5)

dat = data.frame(x=x,y=y,survival=survival,censor=censor)

### Create a group indicator variable
# 1: Group A
# 2: Group B
# 3: Everyone else
group = rep(3,nrow(dat))
group[which(dat$x < 730.5)] = 1
group[which(dat$y >= 2)] = 2


### Kaplan Meier curves
# Need new group indicator variables
A = ifelse(group == 1,1,0)
B = ifelse(group == 2,1,0)
AB = ifelse(group != 3,1,0)


### Overall survival
os = survfit(Surv(dat$survival,dat$censor)~A + B + AB,data=dat) 

So if you run the example and type os you will see that the sample size in AB = 27, when what I want it to be is 17+56=73.


Solution

  • One simple way is to create a new column indicating the group (A or B) the row belongs to, and bind this with the whole population (A+B). Then simply run the model against the group.

    # Create a new variable to indicate the group and drop the group you don't need.
    dat$group = "C"
    dat$group = ifelse( dat$x < 730.5, "A", dat$group )
    dat$group = ifelse( dat$y >= 2, "B", dat$group )
    dat = subset( dat, dat$group != "C" )
    
    # Bind the sample with the population
    dat2 = dat
    dat2$group = "AB"
    data = rbind( dat2, dat )
    
    table( data$group )
    # A AB  B 
    # 17 73 56 
    
    # Plot 
    plot( survfit(Surv(data$survival,data$censor)~group,data=data) )