%---------------------------------------MFILE LOGISTIC REGRESSION-------------------------------------------
%--------BLISS DATA: USING MATLAB's glmfit
%--------and glmval

clear all
close all
 disp('Logistic Regression 2: Bliss Beetle Data')
 lw = 2.5; 
 set(0, 'DefaultAxesFontSize', 16);
 fs = 15;
 msize = 10;

%Insects were exposed to gaseous carbon disulphide for a period of 5 hours.
%Eight experiments were run with different concentrations of carbon disulphide.
%Bliss (1935). The calculation of the dosage-mortality curve. Annals of Applied Biology.
%
%log(concentration)    #killed     #exposed
% beetle=[...
% 3.893 6 59;...
% 3.970 13 60;...
% 4.042 18 62;...
% 4.108 28 56;...
% 4.171 52 63;...
% 4.230 53 59;...
% 4.285 61 62;...
% 4.338 60 60];
%log10(concentration)    #killed     #exposed
beetle=[...
1.6907 6 59;...
1.7242 13 60;...
1.7552 18 62;...
1.7842 28 56;...
1.8113 52 63;...
1.8369 53 59;...
1.8610 61 62;...
1.8839 60 60];

%%%%%%%%%%%%%%%%%
xi = beetle(:,1);  yi=beetle(:,2); ni=beetle(:,3);
figure(1)
[b, dev, stats]        = glmfit(xi, [yi ni],'binomial','link','logit');
[b1, dev1, stats1] = glmfit(xi, [yi ni],'binomial','link','probit');
[b2, dev2, stats2] = glmfit(xi, [yi ni],'binomial','link','comploglog');

[dev, dev1, dev2]

xs = 1.5:0.01:2.0;
ys   = glmval(b, xs, 'logit');
y1s = glmval(b1, xs, 'probit');
y2s = glmval(b2, xs, 'comploglog');
% Plot 
plot(xs,ys,'r-','LineWidth',lw)
hold on
plot(xs,y1s,'k--','LineWidth',lw)
plot(xs,y2s,'-.','LineWidth',lw)
plot(xi, yi./ni, 'o','MarkerSize',msize,...
      'MarkerEdgeColor','k','MarkerFaceColor','g')
axis([1.5 2.0 0 1])
grid on
xlabel('Log concentration')
ylabel('Proportion killed')
legend('Logit','Probit','Clog-log','Obs. proportions','Location','NorthWest')
%%
%print -depsc 'C:\BESTAT\Logistic\Logisticeps\bliss1.eps'
%% 
phati = glmval(b,xi,'logit');
yhati = ni .* phati;
devires =sign(yi - yhati) .* sqrt(  2 * abs( yi .* log(yhati./yi ) + ...
                 (ni - yi) .* log( (ni - yhati)./(ni-yi+eps)  )))
devh =  sum(devires.^2) %
df = 8 -  2;
pval = 1 - chi2cdf(devh, 8 - 2 )  %small p-value critical for the model, too big deviance
% pval = 0.0811
figure(2)
ys = sort(abs(devires));
qs = ((1:8)+8 -1/8)/(2 * 8 + 1/2)
xs = norminv(qs) % $\Phi^{-1}\left( \frac{i+n-1/8}{2 n + 1/2}\right).$
plot(xs, ys, '*','MarkerSize',msize)

% 
%%
phati = glmval(b2,xi,'comploglog');
yhati = ni .* phati;
devires =sign(yi - yhati) .* sqrt(  2 * abs( yi .* log(yhati./yi ) + ...
                 (ni - yi) .* log( (ni - yhati)./(ni-yi+eps)  )))
devh =  sum(devires.^2) %  0.7511
df = 8 -  2;
pval = 1 - chi2cdf(devh, 8 - 2 )  %small p-value critical for the model, too big deviance
% 
figure(3)
ys = sort(abs(devires));
qs = ((1:8)+8 -1/8)/(2 * 8 + 1/2)
xs = norminv(qs) % $\Phi^{-1}\left( \frac{i+n-1/8}{2 n + 1/2}\right).$
plot(xs, ys, '*','MarkerSize',msize)


