% HEART DATA

% This database contains 13 attributes (which have been extracted from
% a larger set of 75)       
%   
% 
% 
% Attribute Information:
% ------------------------
%       -- 1. age       
%       -- 2. sex       
%       -- 3. chest pain type  (4 values: 1-4)       
%       -- 4. resting blood pressure  
%       -- 5. serum cholestoral in mg/dl      
%       -- 6. fasting blood sugar > 120 mg/dl       
%       -- 7. resting electrocardiographic results  (values 0,1,2) 
%       -- 8. maximum heart rate achieved  
%       -- 9. exercise induced angina    
%       -- 10. oldpeak = ST depression induced by exercise relative to rest   
%       -- 11. the slope of the peak exercise ST segment     (1<2<3)
%       -- 12. number of major vessels (0-3) colored by flourosopy        
%       -- 13.  thal: 3 = normal; 6 = fixed defect; 7 = reversable defect     
% 
% Attributes types
% -----------------
% Real: 1, 4, 5, 8, 10, 12
% Ordered: 11 
% Binary: 2, 6, 9
% Nominal:3, 7, 13
% 
% Variable to be predicted
% ------------------------
% Absence (1) or presence (2) of heart disease
% 
% Cost Matrix
% 
% 	                        true absence  true presence
% predict absence	      0	                  5
% predict presence        1                   0
% 
% h
% No missing values; 270 observations
clear all
close all force
msize=7; lw=2;
load('../Logisticdat/heart.dat')
% save '../Logisticdat/heart.mat' heart -mat
x1 = heart(:,1); x2 = heart(:,2);  x3 = heart(:,3); 
x4 = heart(:,4); x5 = heart(:,5);  x6 = heart(:,6); 
x7 = heart(:,7); x8 = heart(:,8);  x9 = heart(:,9); 
x10 = heart(:,10); x11 = heart(:,11);  x12 = heart(:,12); 
x13 = heart(:,13); y = heart(:,14)-1;  

n=length(y)
for i = 1:n
    dux31(i)=0; dux32(i)=0; dux33(i)=0;
    if  (x3(i)==2) dux31(i)=1;
    elseif (x3(i)==3) dux32(i) =1;
    elseif (x3(i)==4) dux33(i)=1;
    end
    %
    %x7 nominal, values 0, 1, 2 (baseline 0)
    dux71(i)=0; dux72(i)=0;  
    if  (x7(i)==1) dux71(i)=1;
    elseif (x7(i)==2) dux72(i) =1;
    end
    %
    %x11 ordinal, values 1<2<3  -- baseline=1
    dux111(i)=0; dux112(i)=0;  
    if  (x11(i)>=2) dux111(i)=1;
    elseif (x11(i)>=3) dux112(i) =1;
    end
    %x13 nominal, values 3, 6, 7 (baseline 3)
    dux131(i)=0; dux132(i)=0;  
    if  (x13(i)==6) dux131(i)=1;
    elseif (x13(i)==7) dux132(i) =1;
    end 
end

X=[x1 x2 dux31' dux32' dux33'  x4 x5 x6 dux71' dux72' x8 x9 x10 dux111' dux112' ...
    x12 dux131' dux132'];


[b, dev, stats]=glmfit(X,y, 'binomial','link','logit')

Xdes =[ones(size(y)) X];
 lin = Xdes * b  %linear predictor, n x 1 vector
% 
figure(1)
plot(lin, y,'o','MarkerSize',msize,...
          'MarkerEdgeColor','k','MarkerFaceColor','g')          
 xx = -8:0.01:8;
 mp =  exp(xx)./(1 + exp(xx));
 hold on
 plot(xx, mp,'r-','LineWidth',lw)
 axis([-10 10   0 1])
 xlabel('Linear Predictor: X b','Interpreter','LaTeX')
 ylabel('Probability of Heart Disease','Interpreter','LaTeX') 
 legend('Observations','Logistic Fit','Position','NorthEast')
 %print -depsc 'C:\Springer\Logistic\Logisticeps\arrh1.eps'

%%
figure;
 pp=0.3:0.001:0.7;
 yhat=[];
for p =  0.3:0.001:0.7  
yhat = [yhat   sum((glmval(b,X,'logit') > p)==y)];
end
plot(pp, yhat/n,'-')
hold on
mdl =  stepwiseglm(X,y,...
    'constant','upper','linear',...
'PredictorVars' ,{'x1', 'x2', 'dux31', 'dux32',' dux33' , 'x4', ...
'x5', 'x6', 'dux71',  'dux72', 'x8', 'x9',  'x10', 'dux111', 'dux112', ...
    'x12', 'dux131', 'dux132'},...
'Distribution','Binomial')

% 1. Adding  dux33, Deviance = 299.702, Chi2Stat = 71.25725, PValue = 3.135573e-17
% 2. Adding dux132, Deviance = 253.7734, Chi2Stat = 45.92859, PValue = 1.226427e-11
% 3. Adding x12, Deviance = 219.3152, Chi2Stat = 34.45825, PValue = 4.354921e-09
% 4. Adding dux111, Deviance = 201.184, Chi2Stat = 18.13116, PValue = 2.061998e-05
% 5. Adding x2, Deviance = 193.1786, Chi2Stat = 8.005407, PValue = 0.004663788
% 6. Adding x4, Deviance = 187.099, Chi2Stat = 6.079615, PValue = 0.01367506
% 7. Adding x8, Deviance = 183.0411, Chi2Stat = 4.057905, PValue = 0.04396502
%%
X=[x2 dux31' dux32' dux33'  x4   x8   dux111' dux112' ...
    x12 dux131' dux132'];
[b, dev, stats]=glmfit(X,y, 'binomial','link','logit')
Xdes =[ones(size(y)) X];
%% 
 yhat=[];
for p =  0.3:0.001:0.7  
yhat = [yhat   sum((glmval(b,X,'logit') > p)==y)];
end
plot(pp, yhat/n,'r-')
%%

 






