1
$\begingroup$

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

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Galen
    Oct 19, 2022 at 2:49
  • $\begingroup$ Looks like CoxPHFitter in lifelines estimates the Cox proportional hazards model in Python. $\endgroup$
    – Galen
    Oct 19, 2022 at 2:52

1 Answer 1

0
$\begingroup$

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))
root$root
# [1] -0.2206112

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.