clc
STATS = AsthmaData1;
[m,n] = size(STATS);
SD = STATS(:,6); 
TH = STATS(:,8); 
ST = STATS(:,4);
GP = STATS(:,7); 
CS = STATS(:,10); 
OB = STATS(:,13); 
CS = CS.*OB;
% Separate population (A) into Patients (P) and Controls (C)
% Program could be simplified if we dealt with columns
% of STATS directly, but would be hard to read
p = 0; q = 0;
for j = 1:m
   G1(j) = 2-ST(j);
   q = q + G1(j);
   G2(j) = ST(j)-1;
   p = p + G2(j);
end
ii = 0; jj = 0; kk =0;
for j = 1:m
   if GP(j) ~= 0.0
      if OB(j) <= 1.5
         kk = kk+1;
         SDA(kk) = SD(j);
         THA(kk) = TH(j);
         GPA(kk) = GP(j);
         CSA(kk) = CS(j);
         if G1(j) > 0
            ii = ii+1;
            SDC(ii) = SD(j);
            THC(ii) = TH(j);
            GPC(ii) = GP(j);
            CSC(ii) = CS(j);
         else
            jj = jj+1;
            SDP(jj) = SD(j);
            THP(jj) = TH(j);
            GPP(jj) = GP(j);
            CSP(jj) = CS(j);
         end
      end
   end
end
sprintf('total population %3.0f ', kk)
sprintf('controls  %3.0f ', ii)
sprintf('patients  %3.0f ', jj)
% Limits and continuous variables for given ranges
%GPA Gap variable: t axis (t0 <= g0 <= t <= g1 <= t1; t is used for continuous form a+b*t+b*t^2
g0 = min(GPA) ; g1 = max(GPA);
t0 = floor(g0); t1 = ceil(g1);
tp = linspace(t0,t1,100);
%SDA Small Diameter: z axis (z0 <= d0 <= z <= d1 <= z1);z is used for continuous form a+b*z+c*z^2
d0 = min(SDA) ; d1 = max(SDA);
z0 = floor(d0); z1 = ceil(d1);
zp = linspace(z0,z1,100);
% THA Thickness: y axis, discrete variable is h
h0 = min(THA) ; h1 = max(THA);
y0 = floor(h0); y1 = ceil(h1);
% CSA Cross section: u axis, discrete variable is c
c0 = min(CSA) ; c1 = max(CSA);
u0 = floor(c0); u1 = ceil(c1);
% residuals and correlations
% (A) For all data
% (1) y = a1 + b1*z + c1*z^2
%[ya,ya1,ya2, res1, Sres1] = parabola2(SDA,THA,zp);
%[a1,b1,c1,res1,Sres1] = parabola(SDA,THA);
% (3) y = a3 + b3*t + c3*t^2
%[ya,ya1,ya2,res3,Sres3] = parabola2(GPA,THA,tp);
% (P) For the Patient population
% (5) y = a5 + b5*z + c5*z^2
%[yp,yp1,yp2,res5,Sres5] = parabola2(SDP,THP,zp);
% (7) y = a7 + b7*t + c7*t^2
%[yp,yp1,yp2,res7,Sres7] = parabola2(GPP,THP,tp);
% (C) For the Control population
% (9) y = a9 + b9*z + c9*z^2
%[yc,yc1,yc2,res9,Sres9] = parabola2(SDC,THC,zp);
%(11) y = a11 + b11*z + c11*z^2
%[yc,yc1,yc2,res11,Sres11] = parabola2(GPC,THC,tp);
% (16) u = a16 + b16*z + c16*z^2
%[ua,ua1,ua2,res16,Sres16] = parabola2(SDA,CSA,zp);
% (17) u = a17 + b17*z + c17*z^2
%[uc,uc1,uc2,res17,Sres17] = parabola2(SDC,CSC,zp);
% (18) u = a18 + b18*z + c18*z^2
%[up,up1,up2,res18,Sres18] = parabola2(SDP,CSP,zp);
% (19) u = a19 + b19*t + c19*t^2
[ua,ua1,ua2,res19,Sres19] = parabola2(GPA,CSA,tp);
% (20) u = a20 + b20*t + c20*t^2
[uc,uc1,uc2,res20,Sres20] = parabola2(GPC,CSC,tp);
% (21) u = a21 + b21*t + c21*t^2
[up,up1,up2,res21,Sres21] = parabola2(GPP,CSP,tp);
%---------PLOTS, pvalues and confidence intervals
% 1: for thickness vs. small diameter
%figure
%subplot(2,2,1)
%siga = Sres1/(kk-3); sigp = Sres5/(jj-3); sigc = Sres9/(ii-3);
%pdt = ((kk-6)/3)*(Sres1/(Sres5+Sres9) -1);
%pv  = pvalue(pdt,kk-6);
%plot(zp,yp1,'r--',zp,yp2,'r--',zp,yp,'r--',SDP,THP,'rx')%,...
%hold on
%plot(zp,ya,'g-')
%plot(zp,yc2,'b-.',zp,yc1,'b-.',zp,yc,'b-.',SDC,THC,'bo')
%axis([z0 z1 1. 2.5])
%title(sprintf('Thknss vs. out_diam (.95 confidence); o/-.-: controls (%3.0f), x/--: asthmatics (%3.0f). F-stat = %4.2f, p-value <= %4.3f',ii,jj,pdt,pv))
%xlabel(sprintf(' Outer Diameter (mm) '))
%ylabel(sprintf(' Wall Thickness (mm) '))
%hold off
%%--------------------------------
%% 2: for thickness vs. gap
%figure
%subplot(2,2,2)
%siga = Sres3/(kk-3); sigp = Sres7/(jj-3); sigc = Sres11/(ii-3);
%pgt = ((kk-6)/3)*(Sres3/(Sres7+Sres11) -1)
%pv  = pvalue(pgt,kk-6);
%plot(tp,ya,'g-')
%hold on
%plot(tp,yp2,'r--',tp,yp1,'r--',tp,yp,'r--',GPP,THP,'rx')
%plot(tp,yc2,'b-.',tp,yc1,'b-.',tp,yc,'b-.',GPC,THC,'bo')
%axis([t0 t1 1. 2.5])
%title(sprintf('Thcknss vs. inner diam(.95 conf.); o/-.-: controls (%3.0f), x/--: asthmatics (%3.0f). F-stat = %4.2f, p-value <= %4.3f',ii,jj,pgt,pv))
%xlabel(sprintf(' Inner Diameter (mm) '))
%ylabel(sprintf(' Wall Thickness (mm) '))
%hold off
%--------------------------------
%% 3: for small diameter vs. cross section
%figure
%subplot(2,2,3)
%pdc = ((kk-6)/3)*(Sres16/(Sres17+Sres18) -1)
%pv  = pvalue(pdc,kk-6);
%plot(zp,ua,'g-')
%hold on
%plot(zp,up2,'r--',zp,up1,'r--',SDP,CSP,'rx', zp,up,'r--')
%plot(zp,uc2,'b-.',zp,uc1,'b-.',SDC,CSC,'bo', zp,uc,'b--')
%hold off
%axis([z0 z1 u0 1.5])
%title(sprintf('CrsSect vs. outdiam (.95 conf.); o/-.-: controls (%3.0f), x/-- :asthmatics (%3.0f). F-stat = %4.2f, p-value <= %4.3f',ii,jj,pdc,pv))
%xlabel(sprintf(' Outer Diameter (mm)  [ Obliquity <= 1.5 ]'))
%ylabel(sprintf(' Cross Section (cm^2) '))
%--------------------------------
%% 4: for gap vs. cross section
%figure
%subplot(2,2,4)
pgc = ((kk-6)/3)*(Sres19/(Sres20+Sres21) -1)
pv  = pvalue(pgc,kk-6);
plot(tp,ua,'g-')
hold on
plot(tp,up1,'r--',tp,up2,'r--',GPP,CSP,'rx', tp,up,'r--')
plot(tp,uc1,'b-.',tp,uc2,'b-.',GPC,CSC,'bo', tp,uc,'b-.')
hold off
axis([t0 t1 u0 1.5])
pv = .1
title(sprintf('CrsSect vs. inndiam (.95 conf.); o/-.-: controls (%3.0f), x/-- :asthmatics (%3.0f). F-stat = %4.2f, p-value > %4.3f',ii,jj,pgc,pv))
xlabel(sprintf(' Inner Diameter (mm) [ Obliquity <= 1.5 ]'))
ylabel(sprintf(' Cross Section (cm^2) '))
%--------------------------------
