Assess Stage Effects on Larynx Cancer Survival

Test for an increasing mortality trend across larynx cancer stages and estimate the relative risk between selected stages

This example uses a larynx cancer survival dataset to study how mortality risk changes with disease stage. We first summarize the observed outcomes across stages and test whether risk increases as stage advances. We then fit a Cox proportional hazards model to estimate the relative risk between stage 4 and stage 2 patients and compare fitted survival curves for patients of the same age.

Load the larynx cancer dataset and extract the stage, follow-up time, age, and event information used in the survival analysis:

In[1]:=

larynx=Normal@ResourceData["Sample Data: Larynx Cancer"];

In[2]:=

stage=larynx〚All,1〛;time=QuantityMagnitude[larynx〚All,2〛,"Months"];age=QuantityMagnitude[larynx〚All,3〛];status=larynx〚All,5〛;cens=Replace[larynx〚All,5〛,{"dead"1,"alive"0},{1}];e=EventData[time,cens]

Out[2]=

EventData

Data points: 90

Censoring: Right



Summarize the number of patients, deaths and censored observations in each stage:

In[3]:=

Dataset@Map["Stage"#〚1,1〛,"Patients"Length[#],"Deaths"Count[#〚All,5〛,"dead"],"Censored"Count[#〚All,5〛,"alive"]&,GatherBy[larynx,First]]

Out[3]=

Stage	Patients	Deaths	Censored
1	33	15	18
2	17	7	10
3	27	17	10
4	13	11	2

Fit a stratified Cox model so that each stage is treated as a separate ordered stratum for trend testing:

In[4]:=

stratifiedModel=CoxModelFit[{stage,e},{},s,StrataVariabless]

Out[4]=

CoxModel

Time: 0.1 to 10.7

Variables: 0



Use

LogRankTest

to obtain the score vector and covariance matrix needed for a test of increasing mortality trend across the ordered stages:

In[5]:=

{z,Σ}=LogRankTest[#["Response"]&/@stratifiedModel["StrataModels"],"Equal",{"TestStatisticList","CovarianceMatrix"}]

Out[5]=

{{-2.496,2.89383,-1.80297,1.40514},{{9.26374,-3.35847,-5.62602,-0.27925},{-3.35847,5.48866,-2.00997,-0.120225},{-5.62602,-2.00997,7.79326,-0.157267},{-0.27925,-0.120225,-0.157267,0.556743}}}

Compute the one-sided test statistic for ordered trends across the four stage levels:

In[6]:=

p=Length[stratifiedModel["StrataModels"]]

Out[6]=

In[7]:=

testStat=Sum[(j-1)z〚j〛,{j,p}]/Sqrt[Sum[(j-1)(g-1)Σ〚j,g〛,{j,p},{g,p}]]

Out[7]=

0.628969

Convert the test statistic to a one-sided p-value for testing whether mortality increases with stage:

In[8]:=

testPValue=NProbability[stat>testStat,statNormalDistribution[]]

Out[8]=

0.264685

Fit a Cox proportional hazards model with stage and age as predictors, treating stage as a categorical variable:

In[9]:=

model=CoxModelFit[{Transpose[{stage,age}],e},{s,a},{s,a},NominalVariabless]

Out[9]=

CoxModel

Time: 0.1 to 10.7

Variables: 4



Extract the fitted coefficients for stage 2 and stage 4:

In[10]:=

{stg2,stg4}=model["BestFitParameters"]〚{1,3}〛

Out[10]=

{0.346576,1.41189}

Estimate the relative risk of death for stage 4 patients compared with stage 2 patients as a hazard ratio:

In[11]:=

hazardRatio=Exp[stg4-stg2]

Out[11]=

2.90175

Use the covariance matrix of the fitted parameters to compute the standard error of the log-relative-risk estimate:

In[12]:=

cov=model["CovarianceMatrix"]

Out[12]=

{{0.174499,0.0638567,0.0708709,0.00193338},{0.0638567,0.163934,0.060948,0.000923715},{0.0708709,0.060948,0.64308,0.00133903},{0.00193338,0.000923715,0.00133903,0.000324239}}

In[13]:=

se=Sqrt[cov〚3,3〛+cov〚1,1〛-2cov〚1,3〛]

Out[13]=

0.822093

Construct a 95% confidence interval for the relative risk between stage 4 and stage 2 patients:

In[14]:=

zcrit=Quantile[NormalDistribution[],.975]

Out[14]=

1.95996

In[15]:=

ci=Exp[{stg4-stg2-zcritse,stg4-stg2+zcritse}]

Out[15]=

{0.579286,14.5354}

Since the entire interval lies above 1, the fitted model indicates a higher risk of death for stage 4 patients than for stage 2 patients.

Plot the fitted survival curves for 50-year-old patients at each cancer stage to compare stage-specific survival patterns visually:

In[16]:=

Plot[Evaluate@Table[model["SF"][{i,50}][t],{i,4}],{t,0,12},ExclusionsNone,PlotLegends(Row[{"Stage ",#}]&/@Range[4]),PlotLabel"Survival Time by Cancer Stage"]

Out[16]=