Adjusted survival curves (Thanks to sample code: https://rpubs.com/daspringate/survival ) Thanks to Moderator/Admin's Great Work! For a successful solution I used advice that could be understood: 1. Peter Dalgaard: The code does not work, because the covariates are not factors. 2. Jeff Newmiller: "Change the columns into factors before you give them to the coxph function, e.g. df$treatment <- as.factor(df$treatment)" And I study David Winsemius's instructions.THANKS!!! Code works: library(survival) library(survminer) df<-read.csv("F:/R/data/edgr-orig.csv", header = TRUE, sep = ";") head(df) # "age" means the age groups ID start stop censor sex age stage treatment 1 0 66 0 2 1 3 1 2 0 18 0 1 2 4 2 3 0 43 1 2 3 3 1 4 0 47 1 2 3 NA 2 5 0 26 0 1 4 3 NA # Change continuous var. as factors df$sex<-as.factor(df$sex) df$age<-as.factor(df$age) df$stage<-as.factor(df$stage) df$treatment<-as.factor(df$treatment) S <- Surv( ??time = df$start, ??time2 = df$stop, ??event = df$censor) model <- coxph( S ~ treatment + age + sex + stage, data = df) plot(survfit(model), ?????las=1, ?????xscale = 1.00, ?????xlab = "Months after diagnosis", ?????ylab = "Proportion survived", ?????main = "Baseline Hazard Curve") treat <- with(df, ??????????????data.frame( ??????????????treatment = levels(treatment), ??????????????age = rep(levels(age)[1], 2), ??????????????sex = rep(levels(sex)[1], 2), ??????????????stage = rep(levels(stage)[1], 2))) plot(survfit(model, newdata = treat), ?????las=1, ?????xscale = 1.00, ?????conf.int = TRUE, ?????xlab = "Months after diagnosis", ?????ylab = "Proportion survived", ?????col = c("red", "green")) legend(8, 0.9, ???????legend = c("Beta blockers", ??????????????????"No beta blockers"), ???????lty = 1, ???????col = c("green", "red")) [[alternative HTML version deleted]]