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]]