# Parameter estimates from extended cox model using the coxph() function from the survival package in R

I want to understand how the parameters from a Cox Model using the coxph() function from the survival package in R are estimated. I am following a book by Rizopoulos [1]. In the book the partial log likelihood for an extended Cox model is given below:

$l(\gamma,\alpha) = \sum_{i=1}^n \int_0^\infty \{R_i(t) \exp(\gamma^T w_i + \alpha y_i (t)) -\log[\sum_j R_j(t) \exp(\gamma^T w_j + \alpha y_j(t))] \} dN_i(t)$.

This is counting process integral notation, which I'm not really familiar with but I take to mean:

$l(\gamma,\alpha) = \left. \sum_{i=1}^n \exp(\gamma^T w_i + \alpha y_i (t)) \right|_{t = \text{time of death(i)}} - \left. \sum_{i=1}^n \log[\sum_j R_j(t) \exp(\gamma^T w_j + \alpha y_j(t))] \right|_{t = \text{time of death(i)}}$,

i.e. LHS is the covariate readings at time of death and RHS is the sum of covariate readings at time of death i over all j. Where $w_i$ is a factor and $y_i(t)$ is a covariate reading. I understand that there is no closed form solution for the MLE parameter estimates. I'd just like to know the form of the score equations so that I can input the parameter estimates and return zero. That would make me feel like I understand what I'm doing.

require(JM)
tdCox.pro<-coxph(Surv(start,stop,event)~pro+treat,data=prothro)

#I make some variables which may be of use in returning the score equation
deaths<-which(prothro$event==1) pred<-as.numeric(prothro$treat) - 1

LHS<-sum(exp(pred[deaths]*tdCox.pro$coefficients[2] +prothro$pro[deaths]*tdCox.pro$coefficients[1])) LHSscore<-sum(t(pred)[deaths]%*%exp(pred[deaths]*tdCox.pro$coefficients[2] + prothro$pro[deaths]*tdCox.pro$coefficients[1]))

x<-prothro$Time[deaths] y<-list() for(i in 1:292) {y[[i]]<-which(prothro$start < x[i] & x[i] <= prothro$stop)} RHS<-rep(0,292) for(i in 1:292) {RHS[i]<-sum(exp(pred[y[[i]]]*tdCox.pro$coefficients[2] + prothro$pro[y[[i]]]*tdCox.pro$coefficients[1]))}

RHSscore<-rep(0,292)
for(i in 1:292)
{RHSscore[i]<-sum(t(pred)[y[[i]]]%*%exp(pred[y[[i]]]*tdCox.pro$coefficients[2] + prothro$pro[y[[i]]]*tdCox.pro\$coefficients[1]))}

I'm sorry if my explanation of what I want was a bit awkward. Please tell me how to return zero from my score equations. Cheers!

[1] Rizopoulos, D. Joint models for longitudinal and time-to-event data: With applications in R CRC Press, 2012, pg47