0
$\begingroup$

Suppose I've a set of information,

id t  s   x          H0         h0        H0i         Hi Di       Eu          Elu
 1 2  1  53  0.01067033 0.01067033 0.01371871 0.01371871  2 1.039715  0.029300175
 2 4  0  60  0.02148908 0.01081875 0.02856059 0.04227930  1 1.019138  0.009121439
 3 5  0  45  0.04344224 0.02195316 0.05377424 0.09605354  2 1.038006  0.027655283
 4 5  0  50  0.04344224 0.02195316 0.05506429 0.15111783  2 1.036866  0.026556713
 5 6  0  52  0.05472922 0.01128698 0.07003181 0.22114964  2 1.035420  0.025161271
 6 7  1  51  0.07763398 0.02290476 0.09887098 0.32002062  2 1.033386  0.023194494

If I want to find the estimate of $\beta$ from this,

$$ Q(\beta)=\sum_{i=1}^{G}\sum_{j=1}^{n_{i}}[\delta_{ij}\{\ln h_{o}(t_{ij})+\beta'x_{ij}+E(\log U_{i})\}-H_{o}(t_{ij})e^{\beta'x_{ij}}E(U_{i})] $$ where, $t_{ij}$ is the $j^\mathrm{th}$ time point of the $i\mathrm{th}$ individual, $\delta_{ij}$ is the status, $x_{ij}$ is the vector of covariates.

Let initial $\beta$ be 0.02. How can we find the estimate of $\beta$ using R code?

$\endgroup$
  • $\begingroup$ How is $\hat{\beta}$ defined? Is it defined by $Q(\hat{\beta}) = 0$? Or by $Q'(\hat{\beta}) = 0$? By the way, does this have something to do with frailty models? $\endgroup$ – ocram Jan 6 '13 at 18:58
  • $\begingroup$ $Q'(\hat{\beta})=0$. Yes, it, is one part of M-step in shared frailty model. I am trying to do it in a way, but facing some problem for double summation and do not get the exact value. $\endgroup$ – Dihan Jan 6 '13 at 19:19
  • $\begingroup$ I have added the 'frailty' tag. $\endgroup$ – ocram Jan 6 '13 at 19:42
2
$\begingroup$

Some background

Using standard notations, the Cox model is defined as $$h(t) = h_0(t) \exp(\mathbf{x}^\prime \mathbf{\beta})$$ and the associated likelihood is \begin{align*} L(\mathbf{\beta}, h_0; \textrm{data}) & = \prod_{j=1}^n h_j(y_j)^{\delta_j} S_j(y_j) \\ & = \prod_{j=1}^n \left(h_0(t) \exp(\mathbf{x}^\prime_j \mathbf{\beta})\right)^{\delta_j} \exp \left( -H_0(y_j) \exp(\mathbf{x}^\prime_j \mathbf{\beta}) \right). \end{align*} The log-likelihood is thus given by $$ \ell(\mathbf{\beta}, h_0; \textrm{data}) = \sum_{j=1}^n \delta_j \left[\log(h_0(y_j)) + \mathbf{x}^\prime_j \mathbf{\beta} \right] - H_0(y_j) \exp(\mathbf{x}^\prime_j \mathbf{\beta}). $$

If $h_0$ is known (up to some parameters), this can be maximised in R using, e.g., 'survreg()'. On the other hand, if $h_0$ is left unspecified, then $\ell(\mathbf{\beta}, h_0; \textrm{data})$ has first to profiled into a partial log-likelihood, and the latter can be maximised in R using, e.g., coxph().


Your $Q(\mathbf{\beta})$ is basically the same as $\ell(\mathbf{\beta}, h_0; \textrm{data})$, except that the $\textrm{E}(\log(U_i))$'s enter as fixed offset terms. You can therefore rely on existing software since it is possible to include an offset in a model formula.

$\endgroup$
  • $\begingroup$ Thanks for your comment. But, I want to do it in step by step R function, not built in R code. So, I need to optimize the $Q(\beta)$. $\endgroup$ – Dihan Jan 6 '13 at 19:46
  • $\begingroup$ In that case, you need to enter minus the likelihood into R (Q <- function(beta, data) {... }) and then to rely on an optimisation routine like 'nml' or 'optim'. $\endgroup$ – ocram Jan 6 '13 at 19:50
  • $\begingroup$ Thanks. I am doing it by using 'optim'. Where some method ('BFGS., 'CG', 'SANN' etc) gives the result with warning but some does not. But do not understand the reason. $\endgroup$ – Dihan Jan 6 '13 at 19:57
  • $\begingroup$ There is not enough info to help you. However, I am afraid this becomes more suited to stackoverflow. $\endgroup$ – ocram Jan 6 '13 at 20:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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