% DISEASE OUTBREAK DATA .
%---------------------------------------------------------
close all
clear all


  disp('DISEASE OUTBREAK Example')
  lw = 3; 
  set(0, 'DefaultAxesFontSize', 16);
  fs = 15;
  msize = 10;
 


alpha = 0.05;
load C:\STAT1\Logistic\Logisticdat\outbreak.mat
%
X1  = outbreak(:,2);  %AGE
X0  = outbreak(:,3); % socioeconomic status 1 upper,2 middle ,3 lower
X4  = outbreak(:,4); %sector of the city 1 or 2 
Y    =  outbreak(:,5);  % disease status YES (1), NO (0)
   % X2, X3 dummy variables for socioeconomic status
   % 0 0 - upper;   1 0 - middle;  0 1 - lower
X5  = outbreak(:,6);  %has savings account 1 yes, 0 no 
% ntotal = 196
n = 98;
%--------------------------------------------------------------------------
              %constructing dummies X2 and X3
              X2 = zeros(size(X0)); X3 = zeros(size(X1));
              X2( X0 == 2 ) = 1;
              X3(X0 == 3) = 1;
 %-------------------------------------------------------------------------
  X1c  = X1(1:98);  X1p  = X1(99:end);
  X2c  = X2(1:98);  X2p  = X2(99:end);  
  X3c  = X3(1:98);  X3p  = X3(99:end);
  X4c  = X4(1:98);  X4p  = X4(99:end);
  Yc  = Y(1:98);  Yp  =Y(99:end);
  %------------------------------------------------------------------------
  %
  [b, dev, stats] = glmfit([X1c X2c X3c X4c],Yc,'binomial','link','logit')
  Xdes = [ones(size(X1c)) X1c X2c X3c X4c]
  lin = b(1) + b(2)*X1c + b(3)*X2c + b(4)*X3c + b(5)*X4c;
  pihat = exp(lin)./(1 + exp(lin));
  V = diag( pihat .* (1 - pihat) ); 
  sqrtV = diag( sqrt(pihat .* (1 - pihat) ))
  sb = sqrt(   diag( inv( Xdes' * V * Xdes ) ) )
  % inv( Xdes' * V * Xdes ) is the same as stats.covb
  % 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( Yc .* lin - log( 1 + exp(lin) ))
  [b0, dev0, stats0] = glmfit(zeros(size(Yc)),Yc,'binomial','link','logit')
  loglik0 = sum( Yc .* 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( Yc .* log( pihat ) + (1-Yc).*log(1 - pihat) )  %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
  G/dev0
  
  mcfaddenadj = 1 - (loglik - 5)/loglik0  %counterpart to R2adj, p=5 parameters
  
  coxsnell = 1 - (exp(loglik0)/exp(loglik) )^(2/n)
  
  nagelkerke = (1 - (exp(loglik0)/exp(loglik) )^(2/n) ) /  ( 1 - exp(loglik0)^(2/n)  )
  
  effron = 1 -  sum((Yc - pihat).^2)/sum((Yc - sum(Yc)/n ).^2)
  
  
  %%%RESIDUAL ANALYSIS
  
  ro = Yc - pihat;  %ordinary residuals, generally noninformative
  rpea = ro./sqrt( pihat .*(1-pihat) ) %sum of rpea^2 is Chi^2 statistics
  
  H = sqrtV * Xdes * inv(Xdes' * V * Xdes) * Xdes' * sqrtV;
  hii = diag(H)
  rstpea = rpea./sqrt(1 - hii)
  rdev = sign(Yc - pihat) .* sqrt( - 2 * Yc .* log(pihat) - 2*(1 - Yc) .* log(1 - pihat) )
  
  sum(rdev.^2)
  dev
  
  
  figure(1)
  subplot(2,2,1)
  plot( pihat, ro, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')
                                          hold on;  plot( [0,1],[0,0],'k-')
  subplot(2,2,2)
  plot( pihat, rpea, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')  
                                       hold on;  plot( [0,1],[0,0],'k-')                                     
  subplot(2,2,3)
  plot( pihat, rstpea, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')  
                                       hold on;  plot( [0,1],[0,0],'k-')
  subplot(2,2,4)
  plot( pihat, rdev, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')
   hold on;  plot( [0,1],[0,0],'k-')
                                       
                                       
                                       
     g = loess(pihat,rpea,pihat,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(pihat, rpea, 'o','MarkerSize',8, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','g')
  
  hold on
  plot(pihat, g,'o','MarkerSize',3, 'MarkerEdgeColor','k',...
                                           'MarkerFaceColor','r')
  plot( [0,1],[0,0],'k-')
  
  
  
  % INFLUENTIAL  OBSERVATIONS
  figure(3)
  k = 1:n;
  q = norminv((k + n - 1/8)./(2 * n + 1/2));
  plot( q, sort(abs(rdev)),  'k-');
 
  rand('state',0)
  env =[];
for i = 1:19
  surrogate = binornd(1, pihat);
  rdevsu = sign(surrogate - pihat) .* sqrt( - 2 * surrogate .* log(pihat) - 2*(1 - surrogate) .* log(1 - pihat) );
  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  30 0 0 2]' ;
pXh = exp(Xh' *  b)/(1 + exp(Xh' * b) )
s2pXp = Xh' * inv( Xdes' * V * Xdes ) * Xh;
spXh = sqrt(s2pXp);
ppXh = Xh' * b;
li = [ppXh - norminv(1-alpha/2) * spXh   ppXh    ppXh + norminv(1-alpha/2) * spXh]
exp(li)./(1 + exp(li))


% PREDICTING SINGLE FUTURE OBSERVATION
pXh > sum(Yc)/n