xlsSeerComputeSurvivalFuntion(xlsfilename)

%By Dr. Rex Cheung -- Mainline, PA, USA. Fall of 2012. 
%This program is a SCOPE sub-routine to plot the survival data stored in
%an excel file
%This download should include this code, a randomly generated excel data
%file for testing purposes, and the example output.
%This program expects two excel columns 1. follow up time and 2. failure
%indicator (1=failure has occured and 0=no failure has occurred.
%A randomly generated excel data file example is included for testing.

%Type this at the command line to run: [a b]= xlsSeerComputeSurvivalFuntion('survivalTestData.xls')

function [f, x]= xlsSeerComputeSurvivalFuntion(xlsfilename)
%This program expects a copy of excel file that could be read and
%written to. So do not use your master file.
%__________________________________________________________________________

%Open the excel file and read the follow up time (futime) column from xls file
uiwait(msgbox('select the column of f/u time data')); %wait for the user to interact
[a b] =xlsread(xlsfilename,-1); %[number text]=xlsread()
columnname=a(1); %extract the name of the column, interchange a and b position in the line above to capture
                  %the information
tempinputcol=a(1:end,1); %same as a(2:end) fro string, use a(1:end,1) for numbers,

futime=tempinputcol;
%__________________________________________________________________________
%Open the excel file and read the failure indicator (failureIndicator) column from xls file
uiwait(msgbox('select the column of failure indicator data')); %wait for the user to interact
[c d] =xlsread(xlsfilename,-1); %[number text]=xlsread()
columnname=c(1); %extract the name of the column, interchange a and b position in the line above to capture
                  %the information
tempinputcol2=c(1:end,1); %same as a(2:end) fro string, use a(1:end,1) for numbers,

failureIndicator=tempinputcol2;

%compute the censored and censoredfutime
censored = 1-failureIndicator;
%concatenate the column of futime and censored into temp
temp= horzcat(futime, censored);
len=length(futime);
for i=1:len
    display(i);
    if (temp(i,2)==1)
        censoredfutime(i)=temp(i,1);
    else
        %do nothing
    end
end
censoredfutime(censoredfutime==0)=[]; %remove the zeros

% Calculate and plot the empirical cdf and confidence bounds obtained using Matlab ecdf()
t=futime;
[f,x,flo,fup] = ecdf(t,'censoring',censored);
survivalCurve = 1-f; %f is the computed failure curve taking into account censoring                  
flo_survivalCurve = 1-flo; %flo and fup are the lower and upper 95% confidence
fup_survivalCurve = 1-fup; %interval bounds for failure, this converts into survival bounds.

%Note x is the time corresponds to where the failure or censoring occured.
stairs(x,survivalCurve,'LineWidth',2);
hold on;
stairs(x,flo_survivalCurve,'r:','LineWidth',2);
stairs(x,fup_survivalCurve,'r:','LineWidth',2);

%save the results into an excel file. User could copy and paste as needed
   zz = num2cell(survivalCurve);
   yy={'survivalCurve',zz{:}};
   xlswrite(xlsfilename, yy',1,'E1'); 
   
   z = num2cell(x);
   y={'failure time',z{:}};
   xlswrite(xlsfilename, y',1,'D1'); 


% Find and plot the censored points on the survival curve with  +  markers

% check if the censoring occurs by the first failure, if so assign fc=1
for i=1:length(censoredfutime)
    if (censoredfutime(i)<=x(1))
            fc(i)=1;
    else
            %do nothing
    end
end
% check if the censoring occurs after the first failure, 
% and by the last filure, if so assign fc depending on where censoring
% occurs relative to the failure times
for j=1:length(censoredfutime)
    for m=1:(length(x)-1)
        if ((censoredfutime(j)>x(m))&&(censoredfutime(j)<=x(m+1)))
            fc(j)=f(m);
        else
            %do nothing
        end
    end
end
% check if the censoring occurs after the last failure, 
%if so assign fc as the f of last failure
max=length(x);
for k=1:length(censoredfutime)
    if (censoredfutime(k)>x(max))
            fc(k)=f(length(x));
    end
end

%plot the censored data as + markers
xc = censoredfutime;
plot(xc,1-fc,'r+');
hold off;