data cov;
	tment=0;
run;

proc phreg data=pbc3;
	model days*status(0)=tment/rl;
	baseline out=breslow cumhaz=breslow covariates=cov;
run;

data breslow; 
	set breslow; 
	daysyears = days/365.25; 
run; 


proc gplot data=breslow;
	plot breslow*daysyears=tment/haxis=axis1 vaxis=axis2;
	axis1 order=0 to 6 by 1 minor=none label=('Years');
	axis2 order=0 to 0.7 by 0.1 minor=none label=(a=90 'Cumulative baseline hazard');
	symbol1  v=none i=stepjl c=blue;
run;



#------------------------------------------------------------------#

# Fit a Cox model using the pbc3 data with treatment as a covariate
coxfit <- coxph(Surv(days, status != 0) ~ tment, data = pbc3, method = "breslow")

# Summary of fit
summary(coxfit)

# Extract cumulative baseline hazard
coxcumhaz <- basehaz(coxfit, centered = FALSE)

# Collect data for plot
coxdata <- data.frame(cumhaz = coxcumhaz$hazard, 
                      time = coxcumhaz$time, 
                      tment = rep("0", nrow(coxcumhaz)), 
                      type = rep("Breslow estimate", nrow(coxcumhaz)))

# Create Figure
fig <- ggplot(aes(x = time / 365.25, y = cumhaz, linetype = tment), data = coxdata) + 
  geom_step(size = 1) + 
  xlab("Time since randomization (years)") + 
  ylab("Cumulative baseline hazard") + 
  scale_linetype_discrete("Treatment", labels = c("Placebo")) + 
  scale_x_continuous(expand = expansion(mult = c(0.001, 0.05)),
                     limits = c(0, 6),
                     breaks = seq(0, 6, 1)) + 
  scale_y_continuous(expand = expansion(mult = c(0.005, 0.05))) +
  theme_general

fig