 % DMD: Duchenne Muscular Dystrophy
% %---------------------------------------------------------
 close all
 clear all
% 
% 
   disp('DMD Example')
   lw = 3; 
   set(0, 'DefaultAxesFontSize', 16);
   fs = 15;
  msize = 10;
%  
% 
% 
alpha = 0.05;
load 'C:\Springer\Logistic\Logisticdat\dmd.dat'

	age = dmd(:,1); % age of a woman in the study 
    ck = dmd(:,2);   %creatine kinase level
     h = dmd(:,3);    % hemopexin
     pk = dmd(:,4);   % pyruvate kinase
     ld = dmd(:,5);     % lactate dehydroginase
     carrier = dmd(:,6);  % indicator if a woman is DMD carrier

     n = length(carrier)
%   
   [b, dev, stats] = glmfit([age ck h pk ld],carrier,'binomial','link','logit')
   Xdes = [ones(size(age)) age ck h pk ld];
   lin = b(1) + b(2)*age + b(3)*ck + b(4)*h + b(5)*pk + b(6) * ld;
   phat = exp(lin)./(1 + exp(lin));
   V = diag( phat .* (1 - phat) ); 
   sqrtV = diag( sqrt(phat .* (1 - phat) ));
   sb = sqrt(   diag( inv( Xdes' * V * Xdes ) ) )
   
%    Wald
   z = b./sb  %tests
   1 - normcdf( z  )    %pvals
%   
%   %CI
   [b - norminv(1 - alpha/2) * sb , b + norminv(1-alpha/2) * sb]
   % CIs for odds ratios
   exp([b - norminv(1 - alpha/2) * sb , b + norminv(1-alpha/2) * sb])
%   
   loglik = sum( carrier .* lin - log( 1 + exp(lin) ))
   [b0, dev0, stats0] = glmfit(zeros(size(carrier)),carrier,'binomial','link','logit')
   loglik0 = sum( carrier .* b0(1) - log(1 + exp(b0(1))) )
%   
    G = -2 * (loglik0 -  loglik)  % G, chisquare in SPSS, LR test, nested model chi2
    dev0 - dev % the same as G, difference of deviances
%   
    %Next 3 are the same -- model deviance
    devi = -2 * sum( carrier .* log( phat + eps) + (1-carrier).*log(1 - phat + eps) )  %by formula
    dev                %MATLAB
    -2 * loglik        %link between loglik and deviance
%   
    %Next 3 are the sam -- Pseudo R^2
    pseudor2 = -2 * (loglik0 -  loglik)/(- 2* loglik0)
    mcfadden = 1 - loglik/loglik0
    1 - dev/dev0
    
    G/dev0
%   
    mcfaddenadj = 1 - (loglik - 6)/loglik0  %counterpart to R2adj, p=5 parameters
%   
    coxsnell = 1 - (exp(loglik0)/exp(loglik) )^(2/n)
    coxsnell1 = 1 - (dev/dev0)^(2/n)
%   
    nagelkerke = (1 - (exp(loglik0)/exp(loglik) )^(2/n) ) /  ( 1 - exp(loglik0)^(2/n)  )
%   
    effron = 1 -  sum((carrier - phat).^2)/sum((carrier - sum(carrier)/n ).^2)
%   
%   
%   %%%RESIDUAL ANALYSIS
%   
   ro = carrier - 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.
   H = sqrtV * Xdes * inv(Xdes' * V * Xdes) * Xdes' * sqrtV;
   hii = diag(H)
   rstpea = rpea./sqrt(1 - hii)
%Deviance Residuals
   rdev = sign(carrier - phat) .* sqrt( - 2 * carrier .* log(phat+eps) - 2*(1 - carrier) .* log(1 - phat+eps) )
   ransc = (betainc(carrier, 2/3, 2/3) - betainc(phat, 2/3, 2/3) ) .* ( phat  .* (1-phat) + eps).^(1/6)
   ransc1 = 3 * (carrier - phat.^(2/3))./(2 * phat.^(1/6)+ eps);
%   
   sum(rdev.^2)
   dev
%   
%   
  figure(1)
  subplot(2,2,1)
  plot( phat, ro, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')
                                          hold on;  plot( [0,1],[0,0],'k-')
  subplot(2,2,2)
  plot( phat, rpea, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')  
                                       hold on;  plot( [0,1],[0,0],'k-')                                     
  subplot(2,2,3)
  plot( phat, rstpea, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')  
                                       hold on;  plot( [0,1],[0,0],'k-')
  subplot(2,2,4)
  plot( phat, rdev, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')
   hold on;  plot( [0,1],[0,0],'k-')
                                       
                                       
%  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)
%    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(2)
  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-')
  
%   
%   
%   % 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
  figure(3)
  k = 1:n;
  q = norminv((k + n - 1/8)./(2 * n + 1/2));
  plot( q, sort(abs(rdev)),  'k-');
 
  % Simulated Envelope
  rand('state',1)
  env =[];
for i = 1:19
  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-');

% PREDICTING NEW MEAN OBSERVATION
Xh =[1   42    111    80    23.6   218]' ;     
pXh = exp(Xh' *  b)/(1 + exp(Xh' * b) ) %prediction
%CI
ppXh = Xh' * b; %linear part
s2pXp = Xh' * inv( Xdes' * V * Xdes ) * Xh; %covar matrix of lin part X beta
spXh = sqrt(s2pXp);
% confidence interval on the linear part
li = [ppXh - norminv(1-alpha/2) * spXh   ppXh    ppXh + norminv(1-alpha/2) * spXh]
% transfromation to the mean response
exp(li)./(1 + exp(li))


% PREDICTING SINGLE FUTURE OBSERVATION

cutoff = sum(carrier)/n
% could be cutoff = 0.5
% 
pXh > sum(carrier)/n

        

