%   Duchenne Muscular Dystrophy Dataset. 
%   Duchenne Muscular Dystrophy, DMD, is a genetically transmitted disease, passed from a mother to her children. Affected female offspring usually suffer no apparent symptoms and may unknowingly carry the disease. Male offspring with the disease die at a young age. Not all cases of the disease come from an affected mother. A fraction, perhaps one third, of the cases arise spontaneously, to be genetically-transmitted by an affected female. This is the most widely held view at present. The incidence of DMD is about 1 in 10,000 male births. The population risk that a woman is a DMD carrier is about 3 in 10,000.
% 
% This dataset is modified data from Percy et al (1981) (entries with missing data excluded).
%  It contains 194 observations corresponding to blood samples collected in a project to develop a screening program for female relatives of boys with Duchenne Muscular Dystrophy   (DMD). The program
%   was implemented in Canada and its goal was to inform a woman of her chances of being a carrier based on serum markers as well as her family pedigree. Another question of interest was whether age should be taken into account. Enzyme levels were measured in known carriers (67 samples) and in a group of non-carriers (127 samples).
% 
% The first two serum markers, creatine kinase and hemopexin 
% (ck,h), are inexpensive to obtain, while the last two,
% pyruvate kinase and lactate dehydroginase ( pk,ld), are expensive.
% 
% %It was of interest to measure how much pk and ld add toward
% %predicting the carrier status.
% The importance of age is also of interest.
% The columns in the data set  dmd.csv|dat|mat|xlsx   are
%  
%     Column    Variable Name    Description  
% ==============================
%      	1   age	   age of a woman in the study 
% 	    2   ck     creatine kinase level  
%         3   h 	   hemopexin  
%         4   pk     pyruvate kinase	 
%         5   ld     lactate dehydroginase 
%         6   carrier   indicator if a woman is DMD carrier 
% =============================== 	
% 
% Maire E. Percy, David F. Andrews, Margaret W. Thompson, John M. Opitz (1981).
% Duchenne muscular dystrophy carrier detection using logistic discrimination: Serum creatine kinase and hemopexin in combination.
%  American Journal of Medical Genetics, 8, 4, 397 -- 409. 
% 

clear all
close all force
 disp('Binary Response: DMD')

%-----------------figure defaults
lw = 2.5;  
set(0, 'DefaultAxesFontSize', 16);
fs = 15;
msize = 8;
randn('seed',3)
%-------------------------------------------------------------

load 'dmd.dat'
y = dmd(:,6);
X = dmd(:,1:5);   %Design matrix n x (p-1) without 
                           %vector 1 (intercept)
 Xdes =[ones(size(y)) X];  %with the intercept: n x p
 [n p] = size(Xdes); 
 alpha = 0.05;  %alpha for CIs

[b, dev, stats]=glmfit(X,y, 'binomial','link','logit')
 lin = Xdes * b  %linear predictor, n x 1 vector

figure(1)
plot(lin, y,'o','MarkerSize',msize,...
          'MarkerEdgeColor','k','MarkerFaceColor','g')          
 xx = -10:0.01:15;
 mp =  exp(xx)./(1 + exp(xx));
 hold on
 plot(xx, mp,'r-','LineWidth',lw)
 axis([-10 15 0 1])
 xlabel('Linear Predictor: $b_0 + b_1 X_1 + ... + b_{5}   X_{5}$','Interpreter','LaTeX')
 ylabel('Probability of DMD','Interpreter','LaTeX') 
 legend('Observations','Logistic Fit','Location','NorthWest')
 %print -depsc 'C:\Springer\Logistic\Logisticeps\arrh1.eps'

 
%% 
  phat = exp(lin)./(1 + exp(lin));
  V = diag( phat .* (1 - phat) ); 
  % sqrtV = diag( sqrt(phat .* (1 - phat) ))
  sb = sqrt(   diag( inv( Xdes' * V * Xdes ) ) )
  % inv( Xdes' * V * Xdes ) is stats.covb
  % Wald tests for parameters beta
  z = b./sb  %tests for beta_i = 0, i=0,...,p-1
  pvals = 2 * normcdf(-abs(z))   %p-values 
%   pvals =
% 
%     0.0000
%     0.0010
%     0.0017
%     0.0024
%     0.0158
%     0.0294
  %%
%   figure(2)
%   range = [2:7];
%   errorbar(range-1, b(range), 1.96*sb(range),'o','MarkerSize',msize,...
%        'LineWidth',2,   'MarkerEdgeColor','k','MarkerFaceColor','g')  
%   hold on
%   plot([-1:12],0*[-1:12],'r-')
%   xlabel('Index','Interpreter','LaTeX');  ylabel('$b \pm 1.96 s_b$','Interpreter','LaTeX')
%   axis([0.5 6.5 -0.4 0.3])
%   hold off
  
  figure(3)
  range = [2:6];
  errorbar(range-1, b(range), 1.96*sb(range),'o','MarkerSize',msize,...
       'LineWidth',2,   'MarkerEdgeColor','k','MarkerFaceColor','g')  
  hold on
  plot([-0:6],0*[-0:6],'r-')
  xlabel('Index','Interpreter','LaTeX'); ylabel('$b \pm 1.96 s_b$','Interpreter','LaTeX')
   axis([-0 6.0 -0.1 0.3])
  hold off
  %print -depsc 'C:\Springer\Logistic\Logisticeps\arrhrb.eps'
 %%  
%  %(1-alpha)*100% CI for betas
   CIs = [b - norminv(1 - alpha/2) * sb , b + norminv(1-alpha/2) * sb]
   %(1-alpha)*100% CIs for odds ratios
   CIsOR = exp([b-norminv(1-alpha/2)*sb , b+norminv(1-alpha/2)*sb])
  
%    CIs =
% 
%   -26.0284  -13.0693
%     0.0601    0.2355
%     0.0176    0.0758
%     0.0292    0.1360
%     0.0219    0.2108
%     0.0012    0.0238
% 
% 
% CIsOR =
% 
%     0.0000    0.0000
%     1.0619    1.2655
%     1.0178    1.0787
%     1.0296    1.1457
%     1.0221    1.2347
%     1.0012    1.0240
% 
%    
   
%% 
 
 mdl=stepwiseglm(X,y,'constant','distrib','binomial','link','logit',...
   'upper','linear', 'criterion','deviance','penter',0.05,'premove',0.1)
% 1. Adding x2, Deviance = 153.6368, Chi2Stat = 96.4398, PValue = 9.199919e-23
% 2. Adding x3, Deviance = 116.3913, Chi2Stat = 37.24551, PValue = 1.041544e-09
% 3. Adding x1, Deviance = 99.0159, Chi2Stat = 17.3754, PValue = 3.06764e-05
% 4. Adding x4, Deviance = 88.8645, Chi2Stat = 10.1514, PValue = 0.0014419
% 5. Adding x5, Deviance = 84.0384, Chi2Stat = 4.82603, PValue = 0.0280332
% 
% mdl = 
% 
% 
% Generalized linear regression model:
%     logit(y) ~ 1 + x1 + x2 + x3 + x4 + x5
%     Distribution = Binomial
% 
% Estimated Coefficients:
%                    Estimate       SE         tStat       pValue  
%                    ________    _________    _______    __________
% 
%     (Intercept)     -19.549       3.3059    -5.9133    3.3541e-09
%     x1              0.14779     0.044749     3.3025    0.00095821
%     x2             0.046698     0.014845     3.1457     0.0016571
%     x3             0.082587     0.027239     3.0319     0.0024301
%     x4              0.11634     0.048199     2.4137       0.01579
%     x5             0.012505    0.0057431     2.1774      0.029449
% 
% 
% 194 observations, 188 error degrees of freedom
% Dispersion: 1
% Chi^2-statistic vs. constant model: 166, p-value = 5.11e-34

%% 
   %Log-likelihood
   loglik = sum( y  .* lin - log( 1 + exp(lin) ))  %-42.0192
   %fitting null model. 
   [b0, dev0, stats0] = glmfit(zeros(size(y)),y,'binomial','link','logit')
      %b0= -0.6395,  dev0=250.0766, stats=... (structure)
   loglik0 = sum( y  .* b0(1) - log(1 + exp(b0(1))) ) % -125.0383
%   
    G = -2 * (loglik0 -  loglik)  %  166.0382 
    dev0 - dev %  166.0382, the same as G, difference of deviances
%   
    %model deviance
    devi = -2 * sum( y  .* log( phat + eps) + (1-y ).*log(1 - phat + eps) )  
                       %  84.0384,  directly 
    dev                %  84.0384,  glmfit output
    -2 * loglik        %  84.0384,  as a link between loglik and deviance
%   
    
    %McFadden Pseudo R^2
    mcfadden = 1 - loglik/loglik0 % 0.6639
    1 - dev/dev0 % 0.6639
    G/dev0 % 0.6639
%   
    mcfaddenadj = 1 - (loglik - (p-1))/loglik0 
    %     0.6240 counterpart to R2adj, p=5 parameters, p-1
    % predictors
%   
    coxsnell  = 1 - (exp(loglik0)/exp(loglik) )^(2/n) % 0.5751
%   
    nagelkerke = (1-(exp(loglik0)/exp(loglik))^(2/n))/( 1-exp(loglik0)^(2/n))
   %  0.7938
%   
    effron = 1 -  sum((y  - phat).^2)/sum((y  - sum(y)/n ).^2) 
    % 0.7019
%   
%%   
%   %%%RESIDUAL ANALYSIS
%   
   ro = y - phat;  %ordinary residuals, generally noninformative,
   % the distribution of ordinary residuals variable and unknown.
   rpea = ro./sqrt( phat .*(1-phat) ); %Pearson Residuals 
   %sum of rpea^2 is Chi^2 statistics, when data in Binomial shape, it is Chi2 with df =
   %c-p, when in Bernoulli shape sum distribution is not chi2.
%   
%Studentized Pearson: Pearson residuals are divided with their estimated
%standard deviation. Estimated hat-matrix H=W^(1/2) X (X' W X)^(-1) X'
%W^(1/2), W = diag( phat_i (1- phat_i) ), Y = HY.
% 
%    H = sqrtV * Xdes * inv(Xdes' * V * Xdes) * Xdes' * sqrtV;
%    hii = diag(H);
%    rstpea = rpea./sqrt(1 - hii);
   
   
%Deviance Residuals
   rdev = sign(y  - phat) .* sqrt( - 2 * y  .* log(phat+eps) - 2*(1 - y) .* log(1 - phat+eps));
%Anscombe Residuals    
   ransc = (betainc(y , 2/3, 2/3) - betainc(phat, 2/3, 2/3) ) .* ( phat  .* (1-phat) + eps).^(1/6);
% 
% Deviance as a sum of squared rev. residuals
   sum(rdev.^2)  %78.2515
% or as an output from glmfit
   dev %78.2515
%   
%   
  figure(4)
  subplot(2,2,1)
  plot( phat, ro, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')
                                          hold on;  plot( [0,1],[0,0],'k-')  
                                          xlabel('$\hat p$','Interpreter','LaTeX')
                                          ylabel('Ordinary Res.','Interpreter','LaTeX')
  subplot(2,2,2)
  plot( phat, rpea, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')  
                                       hold on;  plot( [0,1],[0,0],'k-')  
                                       xlabel('$\hat p$','Interpreter','LaTeX')
                                       ylabel('Pearson Res.','Interpreter','LaTeX')
  subplot(2,2,3)
  plot( phat, rdev, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')  
                                       hold on;  plot( [0,1],[0,0],'k-')
                                      xlabel('$\hat p$','Interpreter','LaTeX')
                                      ylabel('Deviance Res.','Interpreter','LaTeX')
  subplot(2,2,4)
  plot( phat, ransc, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')
   hold on;  plot( [0,1],[0,0],'k-')
      xlabel('$\hat p$','Interpreter','LaTeX')
      ylabel('Anscombe Res.','Interpreter','LaTeX')
                                      
    %%                                   
%  If the model is adequate, the loess smoothing
%  of residuals (Pearson) should result in a function close to 0
      g = loess(phat,rpea,phat,1,1,1);
     g1 = loess(phat,rdev,phat,1,1,1);
%    loess.m -- not part of MATLAB, MatlabCENTRAL, author Richard Johnson
%    part of Data Visualization Toolbox
%    ====================================================================
%    apply loess curve fit -- nonparametric regression
%    x,y  data points
%    newx,g  fitted points
%    alpha  smoothing  typically 0.25 to 1.0
%    lambda  polynomial order 1 or 2
%    if robustFlag is present, use bisquare                                    
  figure(5)
  plot(phat, rpea, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')
  hold on
  plot(phat, g,'o','MarkerSize',3, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','r')
                                        
                                       
  plot( [0,1],[0,0],'k-')
  xlabel('$\hat p$','Interpreter','LaTeX')
  h=legend('Pearson Residuals','Loess Smoothed','Location','NorthEast');set(h,'Interpreter','LaTeX')

  
   
  figure(6)
  plot(phat, rdev, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')
  hold on
  plot(phat, g1,'o','MarkerSize',3, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','r')
  plot( [0,1],[0,0],'k-') 
  xlabel('$\hat p$','Interpreter','LaTeX')
  legend('Deviance Residuals','Loess-Smoothed','Location','NorthEast');
  set(h,'Interpreter','LaTeX')
  
  
  %%
%  print -depsc 'C:\Springer\Logistic\Logisticeps\arrh6.eps'
%   
%   
%   % INFLUENTIAL  OBSERVATIONS
%  Half normal Probability Plot of the deviance residuals are
%  used for checking adequacy of the linear part of the model and 
%  detection of outlying  residuals

% Half normal plot in GLM setting is an extension of Atkinsons (1985) half
% normal plot in regular linear regression models. In the plot the absolute
% values of the deviance residuals rDi are ordered and the kth largest absolute
% residual is plotted against z( (k+n - 1/8)/(2n+1/2) ). Here,
% (k+n - 1/8)/(2n+1/2) is suggested by Blum.  
% z(q) is the q-percentile of the standard normal distribution. 
% As in the linear regression model, the half
% normal plot in GLM setting can also be used to detect influential points.
% 


  figure(7)
  k = 1:n;
  q = norminv((k + n - 1/8)./(2 * n + 1/2));
  plot( q, sort(abs(rdev)),  'k-','LineWidth',1.5);
 
  % Simulated Envelope
  rand('state',1)
  env =[];
for i = 1:1119
  surrogate = binornd(1, phat);
  rdevsu = sign(surrogate - phat).*sqrt(- 2*surrogate .* log(phat+eps) ...
      - 2*(1 - surrogate) .* log(1 - phat+eps) );
  env = [env   sort(abs(rdevsu))];
end
envel=[min(env'); mean(env' ); max(env' )]'; 
hold on
plot( q , envel(:,1), 'r-');
plot( q , envel(:,2), 'g-');
plot( q , envel(:,3), 'r-');
xlabel('Half-Normal Quantiles','Interpreter','LaTeX'); 
ylabel('Abs. Dev. Residuals','Interpreter','LaTeX')
h=legend('Abs. Residuals','Simul. $95\%$ CI','Simul. Mean','Location','NorthWest');set(h,'Interpreter','LaTeX')
axis tight
 


%
% PROBABILITY OF Y=1 FOR A NEW OBSERVATION
Xh =[1  26 100  80  30 100]' ;   
% responses for a new person
pXh = exp(Xh' *  b)/(1 + exp(Xh' * b) ) %      0.5771
%(1-alpha) * 100% CI
ppXh = Xh' * b  %    0.3111
s2pXp = Xh' * inv( Xdes' * V * Xdes ) * Xh  %     1.6654
spXh = sqrt(s2pXp)   %     1.2905
% confidence interval on the linear part
li = [ppXh - norminv(1-alpha/2) * spXh   ppXh + norminv(1-alpha/2) * spXh]
%     -2.2183    2.8404
% transfromation to the CI for the mean response
exp(li)./(1 + exp(li)) %    0.0981    0.9448



% PREDICTING SINGLE FUTURE OBSERVATION

cutoff = sum(y)/n %0.3454
% often taken as 0.5
% 
pXh > cutoff  %Ynew = 1

