% HEART1 DATA

% This database contains 10 attributes (which have been extracted from
% a larger set and further modified to recode some categorical variables
% as binary.     
%   
% 
% 
% Predictor Information:
% ------------------------
%       Column
%       -- 1. age       
%       -- 2. sex       
%       -- 3. resting blood pressure  
%       -- 4. serum cholestoral in mg/dl      
%       -- 5. fasting blood sugar > 120 mg/dl    1 yes         
%       -- 6. maximum heart rate achieved  
%       -- 7. exercise induced angina 1 yes 
%       -- 8. oldpeak = ST depression induced by exercise relative to rest   
%       -- 9. number of major vessels (0-3) colored by flourosopy        
%       -- 10.  thal: 0 = normal; 1 =  defective     
% Variable to be predicted
% ------------------------
%       --11.  Absence (0) or presence (1) of heart disease
% 
% No missing values; 270 observations
clear all
close all force
msize=7; lw=2;
load('../Logisticdat/heart1.dat')
  %%
x1 = heart1(:,1); 
x2 = heart1(:,2);  
x3 = heart1(:,3); 
x4 = heart1(:,4); 
x5 = heart1(:,5);  
x6 = heart1(:,6); 
x7 = heart1(:,7); 
x8 = heart1(:,8);  
x9 = heart1(:,9); 
x10 = heart1(:,10); 
y = heart1(:,11);  
n= length(y);

 
X=[x1 x2 x3 x4 x5 x6 x7 x8 x9 x10];


[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', 'x3',  'x4', ...
'x5', 'x6', 'x7', 'x8', 'x9',  'x10'},...
'Distribution','Binomial')

% Estimated Coefficients:
%                    Estimate        SE         tStat       pValue  
%                    _________    _________    _______    __________
% 
%     (Intercept)      -1.7931        1.678    -1.0686       0.28524
%     x2                1.2614      0.47779       2.64     0.0082909
%     x4             0.0081298    0.0037763     2.1529      0.031329
%     x6             -0.023797    0.0087828    -2.7095     0.0067394
%     x7                1.1397      0.39589     2.8788     0.0039925
%     x8               0.47914      0.17645     2.7153     0.0066206
%     x9                1.0092      0.21862     4.6163    3.9057e-06
%     x10                 1.45      0.37753     3.8408    0.00012265


 %%
X=[x2  x4   x8   x6 x7 x8 x9 x10];
[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-')
%%

 






