Re: Re: Mathematica code for survival functions
- To: mathgroup at smc.vnet.net
- Subject: [mg78433] Re: [mg78316] Re: [mg78212] Mathematica code for survival functions
- From: DrMajorBob <drmajorbob at bigfoot.com>
- Date: Sun, 1 Jul 2007 07:48:23 -0400 (EDT)
- References: <200706260809.EAA03911@smc.vnet.net> <3327004.1183023178098.JavaMail.root@m35>
- Reply-to: drmajorbob at bigfoot.com
In the absence of censoring (removing patients from a study), we should
have
deaths = Most@survivors - Rest@survivors
padded with a 0 to give deaths and survivors the same length.
Darren's example obviously was censored, since that relationship doesn't
hold.
His code confused me, so I did a test to determine what it was doing.
(Maybe this will be of help to others; maybe not!)
Clear[t,s,d]
times=Array[t,5];
survivors=Array[s,5];
deaths=Array[d,5];
KMFit[times,survivors,deaths,x]
\[Piecewise]1 0<=x<t[1]
(-d[1]+s[1])/s[1] t[1]<=x<t[2]
((-d[1]+s[1]) (-d[2]+s[2]))/(s[1] s[2]) t[2]<=x<t[3]
((-d[1]+s[1]) (-d[2]+s[2]) (-d[3]+s[3]))/(s[1] s[2] s[3]) t[3]<=x<t[4]
((-d[1]+s[1]) (-d[2]+s[2]) (-d[3]+s[3]) (-d[4]+s[4]))/(s[1] s[2] s[3]
s[4]) t[4]<=x<t[5]
((-d[1]+s[1]) (-d[2]+s[2]) (-d[3]+s[3]) (-d[4]+s[4]) (-d[5]+s[5]))/(s[1]
s[2] s[3] s[4] s[5]) t[5]<=x<\[Infinity]
Uncensored, the result is much simpler:
Clear[t, s, d]
times = Array[t, 5];
atRisk = Array[s, 5];
deaths = Array[d, 5];
KMFit[times, atRisk, deaths, x] //. {d[k_] :> s[k] - s[k + 1],
s[6] -> s[5]}
\[Piecewise] {
{1, 0 <= x < t[1]},
{s[2]/s[1], t[1] <= x < t[2]},
{s[3]/s[1], t[2] <= x < t[3]},
{s[4]/s[1], t[3] <= x < t[4]},
{s[5]/s[1], t[4] <= x < t[5] || t[5] <= x < \[Infinity]}
}
or
Clear[t,s,d]
times=Array[t,5];
survivors=Array[s,5];
deaths=Join[Most@survivors-Rest@survivors,{0}]
KMFit[times,survivors,deaths,x]
{s[1]-s[2],s[2]-s[3],s[3]-s[4],s[4]-s[5],0}
\[Piecewise]1 0<=x<t[1]
s[2]/s[1] t[1]<=x<t[2]
s[3]/s[1] t[2]<=x<t[3]
s[4]/s[1] t[3]<=x<t[4]
s[5]/s[1] t[4]<=x<t[5]||t[5]<=x<\[Infinity]
Here's a code that, to me, seems less confusing. (Probably because I wrote
it!)
Clear[kmFit]
kmFit[t_, s_, d_, x_] := With[{times = Join[{0}, t, {Infinity}]},
Piecewise@
Transpose@{FoldList[Times, 1, (s - d)/s],
LessEqual @@@
Transpose@{Most@times, ConstantArray[x, Length@times - 1],
Rest@times}}
]
And the following is more efficient, as it doesn't reassemble the parts of
Piecewise every time x changes.
kmFit[t_, s_, d_] :=
kmFit[t, s, d] =
Function[x, Evaluate@With[{times = Join[{0}, t, {Infinity}]},
Piecewise[
Transpose@{FoldList[Times, 1, (s - d)/s],
LessEqual @@@
Transpose@{Most@times, ConstantArray[x, Length@times - 1],
Rest@times}}]
]]
kmFit[times,survivors,deaths]
Function[x$,\[Piecewise]1 0<=x$<=t[1]
s[2]/s[1] t[1]<=x$<=t[2]
s[3]/s[1] t[2]<=x$<=t[3]
s[4]/s[1] t[3]<=x$<=t[4]
s[5]/s[1] t[4]<=x$<=t[5]||t[5]<=x$<=\[Infinity]
]
I'm not sure how best to get rid of that redundant Or in the last piece.
Bobby
On Thu, 28 Jun 2007 03:24:38 -0500, Darren Glosemeyer
<darreng at wolfram.com> wrote:
> Coleman, Mark wrote:
>> Greetings,
>>
>> Is anyone on MathGroup aware of Mathematica code that generates either
>> Kaplan-Meier survival curves or life-table survival curves? Thanks.
>>
>> Best regards,
>> Mark
>>
>>
>
> Here is a function that will give the Kaplan Meier fit as a Piecewise
> function. The first three arguments are lists of equal length for the
> times, survivor counts, and death counts. Your data may require some
> processing to get it into that form.
>
> In[1]:= KMFit[times_, survivors_, deaths_, var_] := Piecewise[
> MapThread[{#1, #2 <= var < #3} &,
> {Join[{1}, With[{list = (survivors - deaths)/survivors},
> FoldList[Times, First[list], Drop[list, 1]]]], Join[{0},
> times],
> Join[times, {Infinity}]}]]
>
> Here is an example based on data from Table 1.4 in The Statistical
> Analysis of Failure Time Data by Kalbfleisch and Prentice.
>
> In[2]:= t = {143, 164, 188, 190, 192, 206, 209, 213, 216, 220, 227, 230,
> 234,
> 246, 265, 304};
>
> In[3]:= s = {19, 18, 17, 15, 14, 13, 12, 11, 10, 8, 7, 6, 5, 3, 2, 1};
>
> In[4]:= d = {1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
>
> In[5]:= KMFit[t, s, d, x]
>
> 18 17
> Out[5]= Piecewise[{{1, 0 <= x < 143}, {--, 143 <= x < 164}, {--, 164 <=
> x < 188},
> 19 19
> 15 14 13
> > {--, 188 <= x < 190}, {--, 190 <= x < 192}, {--, 192 <= x < 206},
> 19 19 19
> 12 11 10
> > {--, 206 <= x < 209}, {--, 209 <= x < 213}, {--, 213 <= x < 216},
> 19 19 19
> 9 63 27
> > {--, 216 <= x < 220}, {---, 220 <= x < 227}, {--, 227 <= x < 230},
> 19 152 76
> 45 9 3
> > {---, 230 <= x < 234}, {--, 234 <= x < 246}, {--, 246 <= x < 265},
> 152 38 19
> 3
> > {--, 265 <= x < 304}}]
> 38
>
>
> Darren Glosemeyer
> Wolfram Research
>
>
--
DrMajorBob at bigfoot.com