%leukoplakia as logit

clear all
close all
load leukoplakia.dat
%
tab = array2table(leukoplakia, 'VariableNames',{'leukoplak','smoker','alcohol'})
% nullmodel
modelspec1 = 'leukoplak  ~ 1'
mdl1 = fitglm(tab,modelspec1,'distribution','binomial','link','comploglog')
dev1 = mdl1.Deviance

% Estimated Coefficients:
%                    Estimate       SE        tStat     pValue 
%                    ________    ________    _______    _______
% 
%     (Intercept)    0.038353    0.089332    0.42933    0.66768
% 
% 212 observations, 211 error degrees of freedom
% Dispersion: 1
% 
% dev1 =   275.4946
% 
%%
modelspec2 = 'leukoplak  ~ 1 + smoker'
mdl2 = fitglm(tab,modelspec2,'distribution','binomial','link','comploglog')
dev2 = mdl2.Deviance

dev1 - dev2

% Estimated Coefficients:
%                    Estimate      SE        tStat       pValue  
%                    ________    _______    _______    __________
% 
%     (Intercept)    -0.2295     0.12449    -1.8436      0.065243
%     smoker         0.65666     0.18182     3.6116    0.00030436
% 
% 212 observations, 210 error degrees of freedom
% Dispersion: 1
% Chi^2-statistic vs. constant model: 12.9, p-value = 0.000321,  
%  
% 
% dev2 =   262.5533
% 
% 
% ans =    12.9413
% 
%%
modelspec3 = 'leukoplak  ~ 1 + alcohol'
mdl3 = fitglm(tab,modelspec3,'distribution','binomial','link','comploglog')
dev3 = mdl3.Deviance

dev1 - dev3
%                    Estimate      SE        tStat      pValue 
%                    ________    _______    _______    ________
% 
%     (Intercept)    0.27073     0.14507     1.8662    0.062009
%     alcohol        -0.2264     0.11728    -1.9303    0.053565
% 
% 
% 212 observations, 210 error degrees of freedom
% Dispersion: 1
% Chi^2-statistic vs. constant model: 4.25, p-value = 0.0394
  1-chi2cdf(4.25, 1)  % 0.0393
% dev3 =   271.2495
%  ans =    4.2451
%%
modelspec4 = 'leukoplak  ~ 1 + smoker + alcohol'
mdl4 = fitglm(tab,modelspec4,'distribution','binomial','link','comploglog')

% Estimated Coefficients:
%                    Estimate       SE        tStat       pValue  
%                    _________    _______    ________    _________
% 
%     (Intercept)     -0.12897    0.20201    -0.63844      0.52319
%     smoker           0.60259    0.19651      3.0665    0.0021656
%     alcohol        -0.076527    0.12669    -0.60403      0.54582
% 
% 
% 212 observations, 209 error degrees of freedom
% Dispersion: 1
% Chi^2-statistic vs. constant model: 13.3, p-value = 0.00128
% 

dev4 = mdl4.Deviance %262.1663

dev1 - dev4  %13.3283
1-chi2cdf(13.3283, 2)  %0.0013
dev2 - dev4  %0.3870
dev3 - dev4  % 9.0832
 

   
%%
%saturated
modelspec5 = 'leukoplak  ~ 1 + smoker +alcohol + smoker:alcohol'
mdl5 = fitglm(tab,modelspec5,'distribution','binomial','link','comploglog')

dev5 = mdl5.Deviance
% Generalized linear regression model:
%     comploglog(leukoplak) ~ 1 + smoker*alcohol
%     Distribution = Binomial
% 
% Estimated Coefficients:
%                       Estimate      SE        tStat      pValue 
%                       ________    _______    _______    ________
% 
%     (Intercept)       0.075951     0.2326    0.32653     0.74403
%     smoker             0.19692    0.30385    0.64809     0.51693
%     alcohol             -0.242    0.16346    -1.4805     0.13875
%     smoker:alcohol     0.48203    0.26957     1.7882    0.073749
% 
% 
% 212 observations, 208 error degrees of freedom
% Dispersion: 1
% Chi^2-statistic vs. constant model: 16.7, p-value = 0.000808
% 
% dev5 =
% 
%   258.7783
% 
% 
 

   


dev4 - dev5   %3.3881

%%
%  of one variable does not depend on the value of the other. For example,
% the eect of desiring no more children is 2 in all four age groups. (This
% assumption must obviously be tested, and we shall see that it is not consistent
% with the data.)
% The deviance of the additive model is 16.8 on three d.f. With this value
% we can calculate three dierent tests of interest, all of which involve com-
% parisons between nested models.
%  As we move from model D to A + D the deviance decreases by 37.2
% while we lose three d.f. This statistic tests the hypothesis H0 : i = 0
% for all i, concerning the net eect of age after adjusting for desire for
% more children, and is highly signicant.
%  As we move from model A to A+D we reduce the deviance by 49.7 at
% the expense of one d.f. This chi-squared statistic tests the hypothesis
% H0 : 2 = 0 concerning the net eect of desire for more children after
% adjusting for age. This value is highly signicant, so we reject the
% hypothesis of no net eects.
%  Finally, the deviance of 16.8 on three d.f. is a measure of goodness of
% t of the additive model: it compares this model with the saturated
% model, which adds an interaction between the two factors. Since t
