clear all
close all

% Despite advances over the past decade, including the advent of 5-[HT_3] receptor antagonists, combination therapy,
% and multimodal strategies,  post-operative nausea and vomiting (PONV) remains  serious and frequent adverse event 
% associated with surgery and anesthesia. 
% PONV can be very distressing for patients, can lead to medical complications and impose economic burdens.  
% A meta analysis of several studies give rates of 37% for nausea and 20% for vomiting in patients undergoing general anaesthesia.  
% There are, however, considerable variations in the reported incidence of PONV, which can be attributed to a number of factors.
% Risk factors for PONV can be divided into patient risk factors,
% procedural risk factors, anaesthetic risk factors and post-operative risk
% factors. Main and well understood risk factors are gender, history of
% motion sickness/PONV, smoking status, and use of postoperative opioids.

 disp('POST OPERATIVE NAUSEA AND VOMITING')
 lw = 3; 
 set(0, 'DefaultAxesFontSize', 16);
 fs = 15;
 msize = 10;

[vars names] = xlsread('C:\Brani\WorkHSI\Lela\PONV.xls');
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 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); 
   
   
  mean(ApfelScore) %0.2934
  std(ApfelScore) %0.1202
  trimmean(ApfelScore,20) %0.2891
  iqr(ApfelScore) %0.1585
  hist(ApfelScore, 30)
   X=[PONV0to24 SinclairScore ApfelScore LelaScore];
corr(X)

%ans =
%
 %   1.0000    0.0366    0.1179    0.1928
%    0.0366    1.0000    0.5780    0.2962
 %   0.1179    0.5780    1.0000    0.6533
 %   0.1928    0.2962    0.6533    1.0000
 
    
%%%
X = ApfelScore(PONV0to24 == 0 );
Y = ApfelScore(PONV0to24 == 1 );

%X = LelaScore(PONV0to24 == 0 );
%Y = LelaScore(PONV0to24 == 1 );

n1 = length(X) %n1=623
n2 = length(Y) %n2=293
 Xbar = mean(X) %Xbar=13.1268
 Ybar = mean(Y) %Ybar = 14.9693
s1 =std(X) %s1=4.3246
s2=std(Y) %s2=4.4885
pval = 2 * fcdf( min(s1^2, s2^2)/max(s1^2, s2^2), n1-1, n2-1 ) %0.4496
sp = sqrt(    ((n1-1)*s1^2 + (n2-1)*s2^2)/(n1 +n2 - 2)  )  %sp=4.3776
t = (Xbar - Ybar)/(sp * sqrt(1/n1 + 1/n2) ) %t=-5.9414
pval = tcdf(t, n1 + n2 - 2)  %pval= 2.0068e-009

%%Sensitivity and Specificity//Risks

a = sum(LelaScore > 13 & PONV0to24 == 1)   % entry a in the table [a b;c d]
b = sum(LelaScore > 13 & PONV0to24 == 0)   % entry b in the table [a b;c d]
c = sum(LelaScore <= 13 & PONV0to24 == 1)  % entry c in the table [a b;c d]
d = sum(LelaScore <= 13 & PONV0to24 == 0)  % entry d in the table [a b;c d]
[se, sp, pre, ppv, npv, lr, ag] =  sesp(a,b,c,d)
[rd rdl rdu rr rrl rru or orl oru] = risk(a, b, c, d)
% 
%                                     PONV   no PONV
% ---------------------------------------------------------------
% LelaScore  >  13        185           296
% LelaScore <= 13        108           327
% ---------------------------------------------------------------
%   
%      Se          Sp          Pre         PPV       NPV       LRP        Ag   
%   -----------------------------------------------------------------------------------
%     0.6314   0.5249   0.3199   0.3846   0.7517  1.3289  0.5590
% 
% BAYES
%
preout = 0.4
ppvout = se * preout/(  se * preout + (1-sp)*(1-preout))
% ppvout = 0.4698


% 
% rd = 0.1363
% [ 0.0769,  0.1958]
% 
% rr = 1.5491
% [1.2699, 1.8898]
% 
% or = 1.8924
% [ 1.4235,  2.5157]
% 
%		     mean	sd	MC_error	val2.5pc	median	val97.5pc	start	sample
%	or	   1.909	0.2786	2.77E-4	   1.423	1.888	2.512	             1001	1000000
%	rd   	0.1357	0.03028	3.022E-5	0.07605	0.1358	0.1947     	1001	1000000
%	rr	   1.555	0.1582	1.568E-4	1.271	1.546	1.891	             1001	1000000


a12 = sum(LelaScore > 12 & PONV0to24 == 1)   %
b12 = sum(LelaScore > 12 & PONV0to24 == 0)   % 
c12 = sum(LelaScore <= 12 & PONV0to24 == 1)  % 
d12 = sum(LelaScore <= 12 & PONV0to24 == 0)  % 

[se12, sp12, pre12, ppv12, npv12, lr12, ag12] =  sesp(a12,b12,c12,d12)

 %    Se        Sp        Pre       PPV       NPV       LRP        Ag   
 %   0.6587    0.4880    0.3199    0.3770    0.7525    1.2864  0.5426


a14 = sum(LelaScore > 14 & PONV0to24 == 1)   %
b14 = sum(LelaScore > 14 & PONV0to24 == 0)   % 
c14 = sum(LelaScore <= 14 & PONV0to24 == 1)  % 
d14 = sum(LelaScore <= 14 & PONV0to24 == 0)  % 
[se14, sp14, pre14, ppv14, npv14, lr14, ag14] =  sesp(a14,b14,c14,d14)

 %    Se        Sp        Pre       PPV       NPV       LRP        Ag   
%    0.6314    0.5249    0.3199    0.3846    0.7517    1.3289  0.5590


figure;
plot( [0  1-sp14  1-sp 1-sp12, 1],[0, se14, se, se12, 1],'ko')
hold on
plot( [0  1-sp14  1-sp 1-sp12, 1],[0, se14, se, se12, 1],'k-')
hold off

aur( [0  1-sp14  1-sp 1-sp12, 1],[0, se14, se, se12, 1])
% 0.5786
%%
%clear all
close all
load neuronfires
whos
%  Name        Size            Bytes  Class     Attributes
%
 % Y         989x1              7912  double              

T  = diff(Y,1);
hist(T,30)

%lambfdahat = 1 

theo = 1 - expcdf(4, 1)
emp = sum( T > 4 )/length(T)
% theo = 0.0183  emp = 0.0182

