# Estimate of parameter $\beta$

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?

• 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? – ocram Jan 6 '13 at 18:58
• $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. – Dihan Jan 6 '13 at 19:19
• I have added the 'frailty' tag. – ocram Jan 6 '13 at 19:42

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.

• 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)$. – Dihan Jan 6 '13 at 19:46
• 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'. – ocram Jan 6 '13 at 19:50
• 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. – Dihan Jan 6 '13 at 19:57
• There is not enough info to help you. However, I am afraid this becomes more suited to stackoverflow. – ocram Jan 6 '13 at 20:00