Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Cox partial likelihood function, how to solve?

• To: mathgroup at smc.vnet.net
• Subject: [mg38517] Re: [mg38490] Cox partial likelihood function, how to solve?
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Fri, 20 Dec 2002 23:40:25 -0500 (EST)
• References: <200212200924.EAA20409@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

Søren Merser wrote:
> 
> Hi there
> I need a litle help here please
> 
> I've have the Mathematica expression below which is a Cox partial likelihood
> function f  (for some survival data)
> but I can't get Mathematica to solve the following problems:
> 
> Calculation of coefficient estimates of exp (-0.130) and age (0.611).
> Results known from a statistical program (R)
>     U=D[f, exp, age]
>     NSolve[U==0, exp, age]
> 
> I'd then would like to calculate the information matrix which is the
> negative second derivative of the log-likelihood function f
> I've tried
> 
>     I = -D[f, {exp,2},{age,2}]
> 
> but the result isn't right  when i replace the exp and age like
> se = Sqrt[Info]-1 /. exp-> -.130 /. age-> .611
> (se should be .920 and .997)
> 
> the methods above works OK when there is only one (exposure) variable - but
> not with additional covars
> 
> regards soren
> 
> likelihood function:
> 
> \!\(Log[\[ExponentialE]\^\(5\ age + 4\ exp\)\/\(\((1 +
> \[ExponentialE]\^\(age \
> + exp\))\)\ \((1 + 2\ \[ExponentialE]\^\(age + exp\))\)\ \((3 + 2\ \
> \[ExponentialE]\^age + \[ExponentialE]\^exp + 4\ \[ExponentialE]\^\(age +
> exp\
> \))\)\^5\)]\)


This works fine for me. Did you forget to take a matrix inverse?

Lets set this up usning simple Mathematica code.

f = Log[E^(5*age + 4*exp)/((1 + E^(age + exp))*(1 + 2*E^(age + exp))*
   (3 + 2*E^age + E^exp + 4*E^(age + exp))^5)];
vars = {exp, age};
grad = Map[D[f,#]&, vars];

We get the parameter values in Mathematica as follows.

vals = FindRoot[Thread[grad==0], {exp,0}, {age,0}]
Out[58]= {exp -> -0.130119, age -> 0.61097}

Now form the information matrix, giving it a valid Mathematica name (not
"I").

infomat = -Map[D[grad,#]&, vars];

Next we find the estimated covariance matrix.

covmat = Inverse[infomat /. vals];

The estimated standard errors are then given by the square roots of the
main diagonal and may be found as below.

standarderrors = Flatten[MapIndexed[Sqrt[#1[[#2]]]&, covmat]]
Out[68]= {0.920345, 0.997452}

I'm only guessing from your naming of the result as "se" that this is
what you want. You really did not make this very clear. Also from your
notation it appears that you are subtracting a matrix of ones from a
square root, or maybe an identity matrix. But presumably the actual
formula uses the "-1" in a way to denote matrix inversion, before taking
any square roots.


Daniel Lichtblau
Wolfram Research

• References:
  • Cox partial likelihood function, how to solve?
    • From: Søren Merser <merser@image.dk>