clear all
close all

 % PONV.M
clear all
close all
 disp('Binary Response:PONV')

%-----------------figure defaults
lw = 2.5;  
set(0, 'DefaultAxesFontSize', 16);
fs = 15;
msize = 8;
randn('seed',3)

%[vars names] = xlsread('C:\Springer\Logistic\Logisticdat\PONV.xls'); 
%  if PONV.mat is used, simply  load PONV
load 'C:\STAT\Logistic\Logisticdat\PONV.mat'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 1)patient	2) Ageover50	3)Age	4)gender	5)diagnosis	6)surgery	7)anaesthesiaduration	8)anaesthesiaoverhour	
    % 9) weight	10) BMI	11) smoking	12) Kinethosishist 13)	PONVhist	14) GLIKOKORTIKOID	15)  METOKLOPRAMIDop	
    % 16) SEROTONIN 17)	CORMACKLIHEENE	18) PONV0to2 19)	PONV2to24	20) PONV0to24	21) Headache0to2
    % 22)	Headache2to24	23) Headache0to24 24)	SinclairScore	25) ApfelScore	26) LelaScore
   patient = vars( :, 1 ); Ageover50=vars(:,2);	age=vars(:,3);	gender=vars(:,4); 	diagnosis=vars(:,5);
   surgery = vars(:,6); 	anaesthesiaduration=vars(:,7); 	anaesthesiaoverhour	= vars(:,8);
   weight=vars(:, 9); 	 BMI =vars(:, 10); 	 smoking=vars(:,11); 	 Kinethosishist=vars(:, 12); 
   PONVhist	=vars(:, 13);  GLIKOKORTIKOID=vars(:, 14);  METOKLOPRAMID =vars(:, 15); 	
   SEROTONIN =vars(:, 16); 	CORMACKLIHEENE	=vars(:, 17);  PONV0to2 =vars(:, 18); 
   PONV2to24=vars(:, 19); 	 PONV0to24=vars(:, 20); 	 Headache0to2=vars(:, 21); 
   Headache2to24=vars(:, 22); 	 Headache0to24 =vars(:, 23); 	SinclairScore=vars(:, 24); 	 ApfelScore=vars(:, 25); 
   LelaScore=vars(:, 26); 
   
%    [b, dev, stats] = glmfit(LelaScore, PONV0to2, 'binomial','link','logit')
%    varNames = {'age','gender','anaesthesiaoverhour','Kinethosishist','PONVhist'};
%    X=[age,gender,anaesthesiaoverhour,Kinethosishist,PONVhist]
%    [b, dev, stats]=glmfit(X,PONV0to24, 'binomial','link','logit')
%     lin = b(1) + b(2)*age + b(3)*gender + b(4)*anaesthesiaoverhour + ...
%                                                     b(5) * Kinethosishist + b(6) * PONVhist;
%     plot(lin, PONV0to24,'o','MarkerSize',msize,...
%                                           'MarkerEdgeColor','k',...
%                                            'MarkerFaceColor','g')
pred =  SinclairScore  + LelaScore/20;
[b, dev, stats] = glmfit(pred, PONV0to24, 'binomial','link','logit') 
plot(pred, PONV0to24,'o','MarkerSize',msize,...
                                           'MarkerEdgeColor','k',...
                                            'MarkerFaceColor','g')
xx = -1:0.01:3;
mp =  exp(b(1) + b(2)*xx)./(1 + exp(b(1) + b(2)*xx));
hold on
plot(xx, mp,'r-','LineWidth',lw)
xlabel('Adjusted Sinclair Score')
ylabel('Probability of PONV in 0-24 hours')
legend('PONV','Logistic Fit',2)
%print -depsc 'C:\Springer\Logistic\Logisticeps\PONV.eps'
%print -djpeg 'C:\Springer\Logistic\Logisticeps\PONV.jpg'
%% Diagnostics Aitkinson's Plot (half-normal Plot)
phati = glmval(b,pred,'logit');
ni=length(pred);
yi=PONV0to24;
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;
ys = sort(abs(devires));
qs = ((1:ni)+ni -1/8)/(2 * ni + 1/2);
xs = norminv(qs) % $\Phi^{-1}\left( \frac{i+n-1/8}{2 n + 1/2}\right).$
plot(xs, ys, '*','MarkerSize',msize)
axis([0 0.45 15 44])

   