
clear all
close all force

% Catheter Data
 load catheter.dat
% catheter=[...
% 16	0	0	2
% 9	    0	0	0
% 16	0	0	0
% 18	0	0	4
% 10	1	0	2
% 29	0	0	2
% 19	1	0	0
% 14	0	0	2
% 24	0	0	2
% 17	0	0	0
% 47	1	0	4
% 51	1	0	2
% 28	0	0	2
% 10	1	0	4
% 9  	0	0	2
% 34	0	0	2
% 17	1	0	2
% 27	0	0	4
% 40	0	0	4
% 26	0	0	2
% 27	0	1	2
% 14	0	1	2
% 22	1	1	2
% 16	0	1	0
% 20	1	1	2
% 14	1	1	2
% 23	0	1	4
% 9	    0	1	2
% 12	0	1	4
% 11	0	1	0
% 18	1	1	2
% 31	0	1	4
% 14	0	1	2
% 16	1	1	4
% 15	1	1	2
% 9	    0	1	0
% ];		 
% 
days=catheter(:,1);  censored=catheter(:,2);  
type=catheter(:,3);  french=catheter(:,4);
%  
X=[type french];

[b,logL,H,stats] = coxphfit(X,days,...
  'censoring',censored,'baseline',0) 

% 
% 
% b =
% 
%     0.7251
%    -0.4442
% 
% 
% logL =
% 
%   -61.9182
% 
% 
% H =
% 
%     9.0000         0
%     9.0000    0.0406
%     9.0000    0.0813
%     9.0000    0.1219
%     9.0000    0.1626
%    11.0000    0.2133
%    12.0000    0.2700
%    14.0000    0.3279
%    14.0000    0.3857
%    14.0000    0.4436
%    16.0000    0.5178
%    16.0000    0.5920
%    16.0000    0.6662
%    17.0000    0.7698
%    18.0000    0.8912
%    23.0000    1.1123
%    24.0000    1.3519
%    26.0000    1.6176
%    27.0000    1.9161
%    27.0000    2.2145
%    28.0000    2.6431
%    29.0000    3.1636
%    31.0000    3.8258
%    34.0000    4.6872
%    40.0000    6.0213
% 
% 
% stats = 
% 
%     covb: [2x2 double]
%     beta: [2x1 double]
%       se: [2x1 double]
%        z: [2x1 double]
%        p: [2x1 double]
% 
%   

      stats.covb 
      
%       ans =
% 
%     0.2077   -0.0144
%    -0.0144    0.0345

      stats.beta 
%     ans =
% 
%     0.7251
%    -0.4442
  
      
      stats.se 
%     ans =
% 
%     0.4557
%     0.1858
%      
      stats.z 
      
%   ans =
% 
%     1.5911
%    -2.3915
%     
      
      stats.p 
       
% ans =
% 
%     0.1116
%     0.0168
       
       
% H is a 2-column matrix that contains cummulative hazard estimate.
% The first column of H containins 'days' values and the second column contains
% the estimated baseline cumulative hazard evaluated at values of 'days'.

% We also selected two subjects from the study to
% calculate survival curves corresponding to their 
% covariates.

X(11,:) % 1 0 4 (censored group=0, french=14+4)
%  
X(33,:) %0 1 2 (observed, group=1, french=14+2)
%  
  
% First we find a cummulative hazards at zero values
% of predictors as well as for subjects #11 and #33
Hbaseline(:,2)   =  H(:,2); %base
Hsubj11(:,2) =  H(:,2) .* exp(X(11,:)*b); %subject #11  
Hsubj33(:,2) =  H(:,2) .* exp(X(33,:)*b); %subject #33

  
% b =
% 
%     0.7251
%    -0.4442

 
% 
 Sbaseline    = exp(-Hbaseline(:,2));  
 Ssubj11   = exp(-Hsubj11(:,2));
 Ssubj33  = exp(-Hsubj33(:,2)); 


%
figure(1)
stairs(H(:,1),Sbaseline,'k-','linewidth',2)
hold on
stairs(H(:,1),Ssubj11,'b-','linewidth',2)
stairs(H(:,1),Ssubj33,'r-','linewidth',2)
xlabel('$t$ (days)','Interpreter','LaTeX')
ylabel('$\hat S(t)$','Interpreter','LaTeX')
legend('baseline','subject #11', 'subject #33', 'location','SW')
axis tight
print -depsc 'cox1.eps'

figure(2)
% 
[pars]=mle(days,'distribution','lognormal','censoring',censored)
%  23.0815    2.1182

% 
stairs(H(:,1),exp(-H(:,2) .*exp(mean(X)*b)),'LineWidth',2)
xx = linspace(0,50);
line(xx,1-logncdf(xx, pars(1),pars(2)),'color','r','LineWidth',2)
xlim([0,50])
legend('Average Survivor Function','Lognormal Survivor Function')
print -depsc 'cox2.eps'

%%
addpath 'C:\STAT1\Survival\Survivalmat\'
ple( days(type==1), censored(type==1))
hold on
ple( days(type==0), censored(type==0))
              
 [p chi2] = manteltwo(days(type==0), censored(type==0), ...
     days(type==1), censored(type==1))
%  p =
% 
%     0.1922
% 
% 
% chi2 =
% 
%     1.7003

 %%
logrank([days(type==1), censored(type==1)],  ...
    [days(type==0), censored(type==0)] )
% LOG-RANK TEST FOR KAPLAN-MEIER SURVIVAL FUNCTIONS
%  
% --------------------------------------------------------------------------------
% UL			S.E.		z		p-value		alpha
% --------------------------------------------------------------------------------
% 2.72357			2.08870		1.06457		0.28707		0.050
% --------------------------------------------------------------------------------
% 		The survival functions are not statistically different



[b,logL,H,stats] = coxphfit(type,days,...
  'censoring',censored,'baseline',0) 

% stats = 
% 
%     covb: 0.1959
%     beta: 0.5556
%       se: 0.4426
%        z: 1.2552
%        p: 0.2094
% 

