kmplot.m

function varargout=kmplot(varargin)
% KMPLOT Plot the Kaplan-Meier estimation of the survival function
% Survival times are data that measure follow-up time from a defined
% starting point to the occurrence of a given event, for example the time
% from the beginning to the end of a remission period or the time from the
% diagnosis of a disease to death. Standard statistical techniques cannot
% usually be applied because the underlying distribution is rarely Normal
% and the data are often "censored". A survival time is described as
% censored when there is a follow-up time but the event has not yet
% occurred or is not known to have occurred. For example, if remission time
% is being studied and the patient is still in remission at the end of the
% study, then that patient�s remission time would be censored. If a patient
% for some reason drops out of a study before the end of the study period,
% then that patient�s follow-up time would also be considered to be
% censored. The survival function S(t) is defined as the probability of
% surviving at least to time t. The graph of S(t) against t is called the
% survival curve. The Kaplan�Meier method can be used to estimate this
% curve from the observed survival times without the assumption of an
% underlying probability distribution.
%
% Syntax: 	kmplot(x,alpha,censflag)
%      
%     Inputs:
%           X (mandatory)- Nx2 data matrix:
%                          (X:,1) = survival time of the i-th subject
%                          (X:,2) = censored flag 
%                                   (0 if not censored; 1 if censored)
%           note that if X is a vector, all the flags of the second column
%           will be set to 0 (all data are not censored).
%           ALPHA (optional) - significance level (default 0.05) 
%           CENSFLAG (optional) - Censored Plot flag (default 0). If 0
%           censored data will be plotted spreaded on the horizontal
%           segment; if 1 they will be plotted at the given time of censoring.
%     Outputs:
%           Kaplan-Meier plot
%
%      Example: (+ indicate that patient is censored)
%  
%                   ---------------------
%                   Patient     Survival
%                               time       
%                   ---------------------
%                      1        7    
%                      2        12   
%                      3        7+    
%                      4        12+  
%                      5        11+  
%                      6        8    
%                      7        9    
%                      8        6
%                      9        7+
%                     10        2    
%                   ----------------------
%    X=[7 0; 12 0; 7 1; 12 1; 11 1; 8 0; 9 0; 6 0; 7 1; 2 0];
%
%    Calling on Matlab the function: kmplot(X) the function will plot the
%    Kaplan-Meier estimation of the survival function
%
%           Created by Giuseppe Cardillo
%           giuseppe.cardillo-edta@poste.it
%
% To cite this file, this would be an appropriate format:Curve
% Cardillo G. (2008). KMPLOT: Kaplan-Meier estimation of the survival
% function.
% http://www.mathworks.com/matlabcentral/fileexchange/22293

%Input Error handling
args=cell(varargin);
nu=numel(args);
if isempty(nu) 
    error('Warning: Data vectors are required')
elseif nu>3
    if nu>4
        error('Warning: Max two input data are required')
    end
end
default.values = {[7 0; 12 0; 7 1; 12 1; 11 1; 8 0; 9 0; 6 0; 7 1; 2 0],0.05,0,1};
default.values(1:nu) = args;
[x alpha cflag flag] = deal(default.values{:});
if ~all(isfinite(x(:))) || ~all(isnumeric(x(:)))
    error('Warning: all X values must be numeric and finite')
end
if isvector(x) 
    x(:,2)=0;
else
    if ~isequal(size(x,2),2)
        error('KMPLOT requires Nx2 matrix data.');
    end
    if ~all(x(:,2)==0 | x(:,2)==1)
        error('Warning: all X(:,2) values must be 0 or 1')
    end
end
if nu>1
    if isempty(alpha)
        alpha=0.05;
    else
        if ~isscalar(alpha) || ~isnumeric(alpha) || ~isfinite(alpha)
            error('Warning: it is required a numeric, finite and scalar ALPHA value.');
        end
        if alpha <= 0 || alpha >= 1 %check if alpha is between 0 and 1
            error('Warning: ALPHA must be comprised between 0 and 1.')
        end
    end
end
if nu==3
    if isempty(cflag)
        cflag=0;
    else
        if ~isscalar(cflag) || ~isnumeric(cflag) || ~isfinite(cflag)
            error('Warning: it is required a numeric, finite and scalar CENSFLAG value.');
        end
        if cflag~=0 && cflag~=1
            error('Warning: CENSFLAG value must be 0 or 1')
        end
    end
end    
clear args default nu
%string for LEGEND function
str1=[num2str((1-alpha)*100) '% confidence interval']; 

%sort data by survival time
x=sortrows(x,1);
%table of patients observed for each survival time
%the TABULATE function sets up this matrix:
%table1=[time count percent(on total)]
table1=[0 size(x,1) 1; tabulate(x(:,1))];
%if all observed time are integers remove not observed time added by
%TABULATE function
table1(table1(:,3)==0,:)=[];

%Table of censored data
table12=tabulate(x(x(:,2)==1));
if ~isempty(table12)
    % remove not observed time added by TABULATE function
    table12(table12(:,3)==0,:)=[];
    % setup the vector of the censored data
    [cens,loc]=ismember(table1(:,1),table12(:,1)); %find censored data
end    

%the percents stored in the the third column are unuseful;
%so, place in the third column how many subjects are still alive at the
%beginning of the i-th interval.
a1=[table1(1,2); -1.*table1(2:end,2)];
table1(:,3)=cumsum(a1); table1(2:end,3)=table1(1:end-1,3);
%number of deaths in the intervals (don't take in account the censored
%data)
if ~isempty(table12)
    table1(cens,2)=table1(cens,2)-table12(loc(cens),2);
end
%finally, delete the first row that is now useless
table1(1,:)=[];

t1=[0;table1(:,1)]; %this is the x variable (time);
%this is the y variable (survival function)
T1=[1;cumprod(1-(table1(:,2)./table1(:,3)))];
if flag %if this function was not called by LOGRANK function
    %compute the standard error of the survival function
    SE=[0;T1(2:end).*sqrt(cumsum(table1(:,2)./(table1(:,3).* ...
        (table1(:,3)-table1(:,2)))))];
end

%censored data plotting
if ~isempty(table12) 
    %if there are censored data after max(t1), add a new cell into the t1,
    %T1 and SE arrays
    if table12(end,1)>=t1(end,1)
        t1(end+1,1)=table12(end,1)+1;
        T1(end+1,1)=T1(end,1);
        if flag %if this function was not called by LOGRANK function
            SE(end+1,1)=SE(end,1);
        end
    end
    if ~cflag
        %vectors preallocation
        xcg=zeros(1,sum(table12(:,2))); ycg=xcg; J=1;
        %for each censored data into the i-th time interval...
        for I=1:size(table12,1)
            %compute how many position into the array they must occupy
            JJ=J+table12(I,2)-1;
            %find the correct time interval in which censored data must be
            %placed
            A=find(t1<=table12(I,1),1,'last');
            B=find(t1>table12(I,1),1,'first');
            %equally divide this interval
            int=linspace(table12(I,1),t1(B,1),table12(I,2)+2);
            %put all in the vectors of the plotting variables
            xcg(J:JJ)=int(2:end-1);
            ycg(J:JJ)=T1(A);
            %update the counter
            J=JJ+1;
        end
    else
        xcg=table1(table1(:,2)==0,1);
        ycg=T1(table1(:,2)==0);
    end
else
    if ~flag %if this function was called by LOGRANK function
        xcg=[]; ycg=[];
    end
end
%compute the hazard rate
c1=T1.*numel(x);
c2=-(diff(log(c1(1:end-1)))./diff(t1(1:end-1)));
lambda=mean(c2(c2~=0));

if flag %if this function was not called by LOGRANK function
    %compute the (1-alpha)*100% confidence interval curves
    cv=realsqrt(2)*erfcinv(alpha); %critical value
    %lower curve (remember that: the lower curve values can't be negative)
    lowc=max(0,T1-SE.*cv);
    %if the lower curve reaches the 0 earlier than survival function, trim the
    %data.
    if isequal(lowc(end-1:end),[0; 0])
        lowcend=find(lowc==0,1,'first');
    else
        lowcend=length(lowc);
    end
    %upper curve (remember that the upper curve values can't be >1)
    upc=min(1,T1+SE.*cv);
    %eventually, correct the data.
    if isequal(upc(end),1) 
        cupend=find(upc<1,1,'last');
        upc(cupend:end)=upc(cupend);
    end

    %compute the median survival time (if exist...)
    if isempty(T1(T1==0.5)) %if there is not a point where T=0.5...
        I=find(T1>0.5,1,'last'); %find the first point where T>0.5
        J=find(T1<0.5,1,'first'); %find the first point where T<0.5
        if isempty(J) %if all points are >0.5...
            mt=0; %...there is no median time
        else 
            %compute the median time by linear interpolation.
            p=polyfit([t1(I) t1(J)],[T1(I) T1(J)],1);
            mt=(0.5-p(2))/p(1);
            str2=['Median time ' num2str(mt)]; %string for LEGEND function
        end
    else
        mt=t1(T1==0.5);
        str2=['Median time ' num2str(mt)]; %string for LEGEND function
    end

    %plot all the data
    clf
    hold on
    S2=stairs(t1(1:lowcend),lowc(1:lowcend),'g--'); %lower confidence interval curve
    stairs(t1,upc,'g--'); %upper confidence interval curve
    S1=stairs(t1,T1,'b'); %Kaplan-Meier survival function
    if mt>0 %if exist a median time...
        S3=plot([0 mt mt],[0.5 0.5 0],'k:'); 
    end
    if ~isempty(table12) %if there are censored data...
        S4=plot(xcg,ycg,'r+');
    else
        S4=[];
    end
    hold off

    %set the axis properly
    xmax=max(t1)+1;
    axis([0 xmax 0 1.2]);
    axis square
    %add labels and legend
    txt=sprintf('Kaplan-Meier estimate of survival function (hazard rate: %0.4f)\n',lambda);
    title(txt,'FontName','Arial','FontSize',14,'FontWeight','Bold'); 
    ylabel('Estimated survival function','FontName','Arial','FontSize',14,'FontWeight','Bold'); 
    xlabel('Time','FontName','Arial','FontSize',14,'FontWeight','Bold'); 
    if mt
        if isempty(S4)
            legend([S1 S2 S3],'Data',str1,str2)
        else
            legend([S1 S2 S3 S4],'Data',str1,str2,'Censored')
        end
    else
        if isempty(S4)
            legend([S1 S2],'Data',str1)
        else
            legend([S1 S2 S4],'Data',str1,'Censored')
        end
    end
    disp('HAZARD RATE IS AN EXPERIMENTAL FUNCTION!!!!')
end
if nargout
    varargout(1)={table1};
    varargout(2)={table12};
    varargout(3)={t1};
    varargout(4)={T1};
    varargout(5)={xcg};
    varargout(6)={ycg};
    varargout(7)={lambda};
end