 
clear all
close all
 disp('Binary Response: Bumpus Sparrow')
% Bumpus' Sparrows Data. After an unusually severe storm in February
% of 1898, a number of house sparrows, Passer domesticus, were brought
% to the Anatomical Laboratory of Brown University, Providence, Rhode
% Island. Seventy-two of these birds revived; sixty-four perished. This
% event is described by Hermon Carey Bumpus, the first PhD graduate
% of Clark University, whose paper (Bumpus, 1898) has served as an example
% of natural selection in action. The data set provided by Bumpus
% included several anatomic measurements on 136 birds (as data structure
% bumpus.mat) and had been analyzed since by many diverse researches.
% surv 1 if survived, 0 if perished
% sex 11 = male young; 12 = male adult; 2 = female
% lbt Length (mm) from tip of the beak to the tip of the tail
% ae Alar extent (mm) from tip to tip of the extended wings
% wei Weight (g)
% lbh Length of beak and head (mm), from tip of the beak to the occiput
% hum Length of Humerus [arm/wing bone] (in)
% fem Length of Femur [thigh bone] (in)
% tib Length of Tibiotarsus [leg bone linked to femur] (in)
% wos Width of Skull (in), from the postorbital bone of
% one side to the postorbital bone of the other
% kos Length of Keel of Sternum [an extension of breastbone] (in)
%-----------------figure defaults
lw = 2.5;  
set(0, 'DefaultAxesFontSize', 16);
fs = 15;
msize = 8;
randn('seed',3)



load 'C:\STAT\Logistic\Logisticdat\bumpus.mat'
 mind = bumpus.sex >= 11 %male sparrows
 Y = bumpus.surv(mind);
 Xlbt = bumpus.lbt(mind);
 Xae = bumpus.ae(mind);
 Xwei = bumpus.wei(mind);
  Xlbh = bumpus.lbh(mind);
 Xhum = bumpus.hum(mind);
 Xfem = bumpus.fem(mind);
  Xtib = bumpus.tib(mind);
 Xwos = bumpus.wos(mind);
 Xkos = bumpus.kos(mind);

 
 
X =[Xlbt  Xae Xwei  Xlbh Xhum  Xfem Xtib Xwos Xkos];

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 = -6:0.01:4;
 mp =  exp(xx)./(1 + exp(xx));
 hold on
 plot(xx, mp,'r-','LineWidth',lw)
 axis([-6 4 0 1])
 xlabel('Score','Interpreter','LaTeX')
 ylabel('Probability of Survival','Interpreter','LaTeX') 
 legend('Observations','Logistic Fit',2)
 %%
 stepwiseglm(X,Y,'constant','upper','linear','distr','binomial','link','logit')
%  Generalized Linear regression model:
%     logit(y) ~ 1 + x1 + x3 + x5 + x9
%     Distribution = Binomial
% 
% Estimated Coefficients:
%                    Estimate      SE        tStat       pValue  
%                    ________    _______    _______    __________
% 
%     (Intercept)       48.45      18.26     2.6533     0.0079711
%     x1             -0.63391    0.16509    -3.8398    0.00012316
%     x3             -0.75313    0.30162    -2.4969      0.012527
%     x5                69.25     20.182     3.4312    0.00060091
%     x9               26.306     11.839      2.222      0.026282
% 
% 
% 87 observations, 82 error degrees of freedom
% Dispersion: 1
%  
 
 
 %%
 %print -depsc 'C:\Springer\Logistic\Logisticeps\arrh1.eps'
%Xlbt  Xae Xwei  Xlbh Xhum  Xfem Xtib Xwos Xko
%    26.2380
%    -0.7488 lbt
%     0.0928 ae
%    -0.8576 wei
%     0.5247 lbh
%    34.3716 hum
%    23.8114 fem
%    -1.8177 tib
%    28.1994 wos
%    23.4389 kos
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 
  
  figure(2)
  range = [2:10];
  errorbar(range-1, b(range), 1.96*sb(range),'o','MarkerSize',msize,...
       'LineWidth',2,   'MarkerEdgeColor','k','MarkerFaceColor','g')  
  hold on
  plot([-1:10],0*[-1:10],'r-')
  xlabel('Index','Interpreter','LaTeX');  ylabel('$b \pm 1.96 s_b$','Interpreter','LaTeX')
  axis([0.5 10.5 -50.4 100.3])
  hold off
  
  figure(3)
  range = [2:5];
  errorbar(range-1, b(range), 1.96*sb(range),'o','MarkerSize',msize,...
       'LineWidth',2,   'MarkerEdgeColor','k','MarkerFaceColor','g')  
  hold on
  plot([-1:6],0*[-1:6],'r-')
  xlabel('Index','Interpreter','LaTeX'); ylabel('$b \pm 1.96 s_b$','Interpreter','LaTeX')
   axis([-0.5 5.5 -2 2])
  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])
%%
mdl=stepwiseglm(X,Y,'constant','distrib','binomial','link','probit',...
    'criterion','deviance','penter',0.05,'premove',0.10)
