Compare Survival Models for Psychiatric Patient Data

Compare Kaplan-Meier and Gompertz-Makeham survival models for right-censored lifetime data

A group of 26 psychiatric patients were followed until death or the end of the study. Patients still alive at the end of the study period contribute right-censored lifetimes. Compare a nonparametric Kaplan-Meier estimate with a parametric Gompertz-Makeham fit to assess median survival time.

Retrieve a description of the dataset:

In[1]:=

ExampleData[{"Statistics","Psychiatric"},"LongDescription"]

Out[1]=

Survival data on 26 psychiatric patients admitted to the University of Iowa hospitals between 1935 and 1948.

Load the data:

In[2]:=

d=EventData@@(ExampleData[{"Statistics","Psychiatric"}]〚All,{3,4}〛)

Out[2]=

EventData

Data points: 26

Censoring: Right



Fit a nonparametric survival model using the default Kaplan-Meier estimator:

In[3]:=

=SurvivalModelFit[d]

Out[3]=

SurvivalModel

Time: 1 to 40

Data points: 26



Fit a Gompertz–Makeham distribution to the censored data by maximum likelihood. In this model, the hazard rate is

(

) =

c·t



, where 𝑡 is time, 𝑎 is a constant baseline risk, 𝑏 sets the size of the age-dependent contribution, and 𝑐 controls how quickly that contribution grows over time, making the model suitable for data in which mortality risk increases with age:

In[4]:=

gm=EstimatedDistribution[d,GompertzMakehamDistribution[a,b]]

Out[4]=

GompertzMakehamDistribution[0.0469025,0.18742]

Overlay the Kaplan–Meier step function with its 95% pointwise confidence band and the fitted Gompertz–Makeham survival curve. The close agreement across the observation window suggests the parametric model is a reasonable fit:

In[5]:=

Plot{[t],["PointwiseBands"][t],SurvivalFunction[gm,t]},{t,0,40},



Out[5]=

	Kaplan–Meier
	95% pointwise band
	Gompertz–Makeham

Retrieve the median survival time from each model. The Kaplan-Meier estimate is exact at observed event times; the parametric estimate is smoother but close:

In[6]:=

["MedianSurvival"]

Out[6]=

In[7]:=

Median[gm]

Out[7]=

32.9879

The Gompertz-Makeham hazard rate increases with time, reflecting an aging effect. Plot it to confirm the shape assumed by the model:

In[8]:=

PlotHazardFunction[gm,t],{t,0,40},



Out[8]=

Retrieve the lower and upper 95% pointwise confidence bounds at the median (t = 30) to quantify uncertainty around the Kaplan-Meier estimate:

In[9]:=

["PointwiseBands"][30]

Out[9]=

{0.367631,0.739082}

Publisher Information

Contributed by: Wolfram Staff

License Information

This work is licensed under a Creative Commons Attribution 4.0 International License

Wolfram Language Example Repository

Compare Survival Models for Psychiatric Patient Data

See Also

Publisher Information

License Information

Compare Survival Models for Psychiatric Patient Data

See Also

Related Symbols

Publisher Information

License Information