# What function is used to estimate regression coefficients of a CoxPH model via coxph() in R

I am trying to understand the step by step process of how to estimate the beta coefficients in a CoxPH model. I first get the beta estimates using the coxph() function in r for a sample dataset of six records. Then I manually calculated the partial likelihood function, L(β) and partial log-likelihood function , l(β). If I use these functions in R to estimate beta coefficients using optim(), it gives the corrects estimated value; fun1() for L(β) and fun2() for l(β)

Here is my code:

##########coxph() in r to get beta estimate#############
test_cox_sbp <- coxph(Surv(Time, Status) ~ SBP , data = testd_sbp)
summary(test_cox_sbp)

########## L(β)#############
fun1 = function(beta1){
h1 = exp(120*beta1)
h2 = exp(130*beta1)
h3 = exp(140*beta1)
h4 = exp(150*beta1)
h5 = exp(160*beta1)
h6 = exp(145*beta1)
p1=h1/(h1+h2+h3+h4+h5+h6)
p2=h3/(h3+h4+h5+h6)
p3=h4/(h4+h5+h6)
return(p1*p2*p3)
}

result <- optim(par=0, fn = fun1, method = "L-BFGS-B",
control=list(fnscale = -1),
lower = -3, upper = 1)
########## l(β)#############
fun2 = function(beta1){
h1 = exp(120*beta1)
h2 = exp(130*beta1)
h3 = exp(140*beta1)
h4 = exp(150*beta1)
h5 = exp(160*beta1)
h6 = exp(145*beta1)
p1 = ln(h1)-ln(h1+h2+h3+h4+h5+h6)
p2 = ln(h3)-ln(h3+h4+h5+h6)
p3 = ln(h4)-ln(h4+h5+h6)
return(p1+p2+p3)
}

result <- optim(par=0, fn = fun2, method = "L-BFGS-B",
control=list(fnscale = -1),
lower = -3, upper = 1)
result$par  All the functions above including the built-in function in R return beta1 = -0.2206. But when I use the score function U(β) which is the first derivative of l(β); fun3 in the code below , it doesn't return the correct value ########## U(β)############# fun3 = function(beta1){ h1 = exp(120*beta1) h2 = exp(130*beta1) h3 = exp(140*beta1) h4 = exp(150*beta1) h5 = exp(160*beta1) h6 = exp(145*beta1) a1 = 120 * h1 a2 = 130 * h2 a3 = 140 * h3 a4 = 150 * h4 a5 = 160 * h5 a6 = 145 * h6 p1 = 120-((a1+a2+a3+a4+a5+a6)/(h1+h2+h3+h4+h5+h6)) p2 = 140-((a3+a4+a5+a6)/(h3+h4+h5+h6)) p3 = 150-((a4+a5+a6)/(h4+h5+h6)) return(p1+p2+p3) } result <- optim(par=0, fn = fun3, method = "L-BFGS-B", control=list(fnscale = -1), lower = -3, upper = 1) result$par


If this is not the correct way to estimate the beta coefficient but then how does the coxph() function in R return beta coefficient estimates. What is the main function behind the coxph() function in R

• Here is a gist that attempts to reproduce the optimization examples for fun1, fun2, and fun3 in Python. For fun1 and fun2 I also obtained -0.22061141, but for fun3 I obtained -3 which happens to be one of the bounds of parameter constraints. Oct 19, 2022 at 2:49
• Looks like CoxPHFitter in lifelines estimates the Cox proportional hazards model in Python. Oct 19, 2022 at 2:52

It looks like your code relating to fun3 is trying to “optimize” the score function. You need instead to solve for the value of $$\beta$$ at which the score function equals 0.
root <- uniroot(fun3,c(-3,1))

That's the value of $$\beta$$ you found by minimizing the (log) partial likelihood.