function [T1 T2]=logrank(varargin)
% LOGRANK Comparing survival curves of two groups using the log rank test
% Comparison of two survival curves can be done using a statistical
% hypothesis test called the log rank test. It is used to test the null
% hypothesis that there is no difference between the population survival
% curves (i.e. the probability of an event occurring at any time point is
% the same for each population). This function use the Kaplan-Meier
% procedure to estimate the survival function, so it is mandatory to download
% KMPLOT (http://www.mathworks.com/matlabcentral/fileexchange/22293).
%
% Syntax: 	logrank(x1,x2,alpha)
%      
%     Inputs:
%           X1 and X2 (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) 
%     Outputs:
%           Kaplan-Meier plot
%           Log-rank statistics
%
%      Example: 
%           load logrankdata x1 x2
%           logrank(x1,x2)
%
%LOG-RANK TEST FOR KAPLAN-MEIER SURVIVAL FUNCTIONS
%
%--------------------------------------------------------------------------------
%UL				S.E.			z				p-value			alpha
%--------------------------------------------------------------------------------
%6.57226		2.80788			2.16258			0.01529			0.050
%--------------------------------------------------------------------------------
%		The survival functions are statistically different
%
%           Created by Giuseppe Cardillo
%           giuseppe.cardillo-edta@poste.it
%
% To cite this file, this would be an appropriate format:
% Cardillo G. (2008). LogRank: Comparing survival curves of two groups
% using the log rank test
% http://www.mathworks.com/matlabcentral/fileexchange/22317

%Input Error handling
args=cell(varargin);
nu=numel(args);
if isempty(nu) || nu<2
    error('Warning: Data vectors are required')
elseif nu>3
    error('Warning: Max three input data are required')
end
default.values = {[],[],0.05};
default.values(1:nu) = args;
[x1 x2 alpha] = deal(default.values{:});
if ~all(isfinite(x1(:))) || ~all(isnumeric(x1(:))) ...
        || ~all(isfinite(x2(:))) || ~all(isnumeric(x2(:)))
    error('Warning: all X1 and X2 values must be numeric and finite')
end
if isvector(x1) 
    x1(:,2)=0;
else
    if ~isequal(size(x1,2),2)
        error('LOGRANK requires Nx2 matrix data.');
    end
    if ~all(x1(:,2)==0 | x1(:,2)==1)
        error('Warning: all X1(:,2) values must be 0 or 1')
    end
end
if isvector(x2) 
    x2(:,2)=0;
else
    if ~isequal(size(x2,2),2)
        error('LOGRANK requires Nx2 matrix data.');
    end
    if ~all(x2(:,2)==0 | x2(:,2)==1)
        error('Warning: all X2(:,2) values must be 0 or 1')
    end
end
if nu>2
    if ~isscalar(alpha) || ~isnumeric(alpha) || ~isfinite(alpha) || isempty(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
clear args default nu

%recall KMPLOT function to construct tables of data (table1 and table2),
%tables of censored data (table12 and table 22), Kaplan-Meier variables
%(t1, t2, T1 and T2) and Kaplan-Meier graphical data for censored data 
%(xcg and ycg).
try
    [table1 table12 t1 T1 xcg1 ycg1]=kmplot(x1,0.05,0);
    [table2 table22 t2 T2 xcg2 ycg2]=kmplot(x2,0.05,0);
catch ME
    disp('Download KMPLOT: http://www.mathworks.com/matlabcentral/fileexchange/22293')
    rethrow(ME);
end

%plot both Kaplan-Meier curves
clf
hold on
S1=stairs(t1,T1,'b'); %Kaplan-Meier curve for treatment 1
if ~isempty(table12)
    S3=plot(xcg1,ycg1,'k+'); %Censored data for treatment 1 (if there are)
else
    S3=[];
end

S2=stairs(t2,T2,'r'); %Kaplan-Meier curve for treatment 2
if ~isempty(table22)
    S3=plot(xcg2,ycg2,'k+'); %Censored data for treatment 2 (if there are)
end
hold off
%set the axis properly
xmax=max([t1;t2])+1;
axis([0 xmax 0 1.2]);
axis square
%add labels and legend
title('Kaplan-Meier estimate of survival functions')
ylabel('Estimated survival functions')
xlabel('Time')
if isempty(S3)
    legend([S1 S2],'Treatment 1','Treatment 2')
else
    legend([S1 S2 S3],'Treatment 1','Treatment 2','Censored')
end
clear S1 S2 S3 xmax T1 T2 t1 t2 xcg1 ycg1 xcg2 ycg2

%Full-blown LOGRANK procedure
%Merge the first columns of Table1 and Table2 (time intervals)
%and pick-up unique values
A=unique([table1(:,1);table2(:,1)]);
table=zeros(length(A),9); %matrix preallocation
%Out in the first column the time intervals
table(:,1)=A; 
%Put in the columns 2 and 3 and in the proper rows the deaths and alive
%taken from table1 columns 2 and 3
[c ia ib]=intersect(table1(:,1),A);
table(ib,2:3)=table1(ia,2:3);
%Put in the columns 4 and 5 and in the proper rows the deaths and alive
%taken from table2 columns 2 and 3
[c ia ib]=intersect(table2(:,1),A);
table(ib,4:5)=table2(ia,2:3);
%remove the rows where there arent't deaths in both treatments
table((table(:,2)==0 & table(:,4)==0),:)=[];
clear A c ia ib table1 table2
%fill the "pigeon-holes"
c=find(table(:,3)==0); %find the "pigeon-holes" of treatment 1
for I=1:length(c)
    if c(I)~=1
        %find the first interval time before the hole where there is almost 1
        %death
        J=find(table(1:c(I)-1,3)>0,1,'last');
        table(c(I),3)=table(J,3)-table(J,2);
        if ~isempty(table12)
        %find eventually censored data
            K=find((table12(:,1)<table(c(I),1) & table12(:,1)>=table(J,1)),1,'last');
            %Put in the hole how many subject were alive before the interval time
            %of the hole
            if ~isempty(K)
                table(c(I),3)=table(c(I),3)-sum(table12(K,2));
            end
        end
    else
        table(1,3)=length(x1);
    end
end
%Do the same for tratment 2
c=find(table(:,5)==0);
for I=1:length(c)
    if c(I)~=1
        J=find(table(1:c(I)-1,5)>0,1,'last');
        table(c(I),5)=table(J,5)-table(J,4);
        if ~isempty(table22)
            K=find((table22(:,1)<table(c(I),1) & table22(:,1)>=table(J,1)),1,'last');
            if ~isempty(K)
                table(c(I),5)=table(c(I),5)-sum(table22(K,2));
            end
        end
    else
        table(1,5)=length(x2);
    end
end
clear c I J K table12 table22

%Fill the table and compute the statistic variable
%Compute the total deaths and alive before the i-th time interval
table(:,6:7)=[sum(table(:,[2 4]),2) sum(table(:,[3 5]),2)];
%Compute the difference between observed deaths for treatment 1 and
%expected deaths in the hyphthesis that the treatments are similar
table(:,8)=table(:,2)-table(:,3).*table(:,6)./table(:,7);
%Log-rank statistic is the sum of column 8 values
UL=abs(sum(table(:,8)));
%Compute the contribute to the standard error
table(:,9)=prod(table(:,[3 5 6]),2).*(table(:,7)-table(:,6)) ...
    ./(table(:,7).^2.*(table(:,7)-ones(size(table,1),1)));
%find if there is some NaN (i.e. 0/0)
loc=isnan(table(:,9));
if any(loc)
    table(loc,9)=0;
end
SUL=sqrt(sum(table(:,9))); %Compute the totale standard error
z=abs((UL-0.5)/SUL); %normalized UL with Yates'es correction
p=2*(1-0.5*erfc(-z/realsqrt(2))); %p-value

%display results
disp('LOG-RANK TEST FOR KAPLAN-MEIER SURVIVAL FUNCTIONS')
disp(' ')
tr=repmat('-',1,80);
disp(tr)
fprintf('UL\t\t\tS.E.\t\tz\t\tp-value\t\talpha\n')
disp(tr)
fprintf('%0.5f\t\t\t%0.5f\t\t%0.5f\t\t%0.5f\t\t%0.3f\n',UL,SUL,z,p,alpha)
disp(tr)
if p<alpha
    fprintf('\t\tThe survival functions are statistically different\n')
 else
    fprintf('\t\tThe survival functions are not statistically different\n')
end

function [O1 O2 O3 O4 O5 O6]=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)
%      
%     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) 
%     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:
% 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
    error('Warning: Max two input data are required')
end
default.values = {[],0.05,1};
default.values(1:nu) = args;
[x alpha 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 ~isscalar(alpha) || ~isnumeric(alpha) || ~isfinite(alpha) || isempty(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
clear args default nu
%string for LEGEND function
str1=[num2str((1-alpha)*100) '% CI']; 

%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
    %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
    if ~flag %if this function was called by LOGRANK function
        xcg=[]; ycg=[];
    end
end

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=['Med T ' num2str(mt)]; %string for LEGEND function
        end
    else
        mt=t1(T1==0.5);
        str2=['Med T ' 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
    %title('Kaplan-Meier estimate of survival function')
    ylabel('$\hat S(t)$','Interpreter','LaTeX')
    xlabel('Time')
    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
else
    O1=table1;
    O2=table12;
    O3=t1;
    O4=T1;
    O5=xcg;
    O6=ycg;
end


