MATLAB Answers

coxphfit survival numerical issue

3 views (last 30 days)
leo nidas
leo nidas on 15 Dec 2019
Commented: leo nidas on 16 Dec 2019
Below there is a code that works fine and it generates two correlated variables based on bivariate normality. I then want to naively use the Cox model to model one based on the other. When I ask for the survival curve of the one (say the response) at a given value of the other (say the covariate) fixed at a value close to its mean I get a behavior of the survival "cutting down" to 0 in spite of the fact that the cumulative hazard is looking good at the same range of support. The problem is that the survival (being exp(-H)) is getting really small due to large values of the H so MATLAB outputs that as 0 while it is not. Is there a way around it so that I can visualize the full tail of the survival curve at the third panel of the generated plot?
I know this is not the way to generate values from the Cox model. However I do want to generate bivariate normal data and then pretend one is the response and the other is the covariate. More specifically I would like to see how this process works for high correlation values, but the higher the correlation the stronger the problem. If I set "rhoreal=0" (which is the correlation) then it works fine (probably because of the smaller values of H).
Any ideas?
close all
clear all
clc
aa=3
rng(aa,'twister');
na=500;
nb=500;
rhoreal=0.8;
m1a=6;
m2a=6;
s1a=1;s2a=1;
mycov=rhoreal;
Wa=mvnrnd([m1a m2a]',[s1a^2 mycov;mycov s2a^2],na);
x1a=Wa(:,1);x2a=Wa(:,2);
ca=cov([x1a x2a]);cova=ca(1,2);
statusA=zeros(na,1); % Have this as a column of zeros equal to the length of the scores
%I have generated two correlated varianbles x1a x2a and I want to
%naively apply the Cox model with one as a response and the other
%as a covariate. I KNOW THAT THIS IS NOT THE WAY TO SIMULATE FROM A COX
%MODEL. The Cox model follows:
[b,logL,H0,st] = coxphfit(x1a,x2a,'Baseline',0);
S0=exp(-H0(:,2));
subplot(1,3,1)
stairs(H0(:,1),H0(:,2),'r')
title('Baseline Cumulative Hazard')
subplot(1,3,2)
stairs(H0(:,1),S0,'r')
title('Baseline Survival')
subplot(1,3,3)
Sz= S0.^(exp(b.*7))
stairs(H0(:,1),Sz,'r')
title('survival curve given the value of the covariate =7')
%Is there a way of completing the tail of the survival curve in the
%last panel. It "cuts down" to zero but this a numerical issue due
%the survival begin exp( -large number).

Answers (1)

Thiago Henrique Gomes Lobato
In your S0 variable for the 0.8 case the last non-zero value is 3.2114e-321, which is way lower than the maximum numerical accuracy you could get with doubles (eps; 2.2204e-16). Even if you somehow hack your way to make matlab "go below" there's absolutely no guarantee that your results would be right, in fact, if you have to deal with anything below eps it is almost certain that you have wrong results due to numerical accuracy. I would not trust your results if a numerical difference of 2e-317 (the difference of your two last non-zero S0) is suppose to produce meanigful results, the model is too ill-conditioned with those parameters, and even though some partial resutls may appear reasonable that is absolutely no guarantee that it actually is.
  1 Comment
leo nidas
leo nidas on 16 Dec 2019
Thank you for your answer.
There is an easy "hack" after all that provides a guarantee that the results are right because I can bypass the need of calculating such small quantities. I only needed to realize that instead of projecting S through the linear term of the Cox model directly, I could go through the cumulative hazard first and then use S=exp(-H). Namely, add in the previous code the following part to visualize the tail of the truncated survival curve
hold on
Hz= H0(:,2).*(exp(b.*7))
Sz=exp(-Hz)
stairs(H0(:,1),Sz,'g')

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!