clc
STATS = AsthmaData1;
[m,n] = size(STATS);
SD = STATS(:,6); ii(1) = 6;
TH = STATS(:,8); ii(4) = 8;
ST = STATS(:,4);
RR = STATS(:,12);ii(3) = 12;
GP = STATS(:,7); ii(2) = 7;
CS = STATS(:,10); ii(5) = 10;
OB = STATS(:,13); ii(6) = 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
STA = [SD,TH,GP,RR,CS];
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);
         RRA(kk) = RR(j);
         CSA(kk) = CS(j);
         if G1(j) > 0 
            ii = ii+1;
            SDC(ii) = SD(j);
            THC(ii) = TH(j);
            GPC(ii) = GP(j);
            RRC(ii) = RR(j);
            CSC(ii) = CS(j);
         else
            jj = jj+1;
            SDP(jj) = SD(j);
            THP(jj) = TH(j);
            GPP(jj) = GP(j);
            RRP(jj) = RR(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);
t = 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);
z = linspace(z0,z1,100);
% RRA Ratio: x-axis, discrete variable r
r0 = min(RRA) ; r1 = max(RRA);
x0 = floor(r0); x1 = ceil(r1);
% 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
[a1,b1,c1,res1,Sres1] = parabola(SDA,THA);
% (2) x = a2 + b2*z + c2*z^2
[a2,b2,c2,res2,Sres2] = parabola(SDA,RRA);
% (3) y = a3 + b3*t + c3*t^2
[a3,b3,c3,res3,Sres3] = parabola(GPA,THA);
% (4) x = a4 + b4*z + c4*z^2
[a4,b4,c4,res4,Sres4] = parabola(GPA,RRA);
% (P) For the Patient population
% (5) y = a5 + b5*z + c5*z^2
[a5,b5,c5,res5,Sres5] = parabola(SDP,THP);
% (6) x = a6 + b6*z + c6*z^2
[a6,b6,c6,res6,Sres6] = parabola(SDP,RRP);
% (7) y = a7 + b7*t + c7*t^2
[a7,b7,c7,res7,Sres7] = parabola(GPP,THP);
% (8) x = a8 + b8*z + c8*z^2
[a8,b8,c8,res8,Sres8] = parabola(GPP,RRP);
% (C) For the Control population
% (9) y = a9 + b9*z + c9*z^2
[a9,b9,c9,res9,Sres9] = parabola(SDC,THC);
%(10) y = a10 + b10*z + c10*z^2
[a10,b10,c10,res10,Sres10] = parabola(SDC,RRC);
%(11) y = a11 + b11*z + c11*z^2
[a11,b11,c11,res11,Sres11] = parabola(GPC,THC);
%(12) y = a12 + b12*z + c12*z^2
[a12,b12,c12,res12,Sres12] = parabola(GPC,RRC);
%(13) z = a13 + b13*t + c13*t^2
[a13,b13,c13,res13,Sres13] = parabola(GPA,SDA);
%(14) z = a14 + b14*t + c14*t^2
[a14,b14,c14,res14,Sres14] = parabola(GPP,SDP);
%(15) z = a15 + b15*t + c15*t^2
[a15,b15,c15,res15,Sres15] = parabola(GPC,SDC);
% (16) u = a16 + b16*z + c16*z^2
[a16,b16,c16,res16,Sres16] = parabola(SDA,CSA);
% (17) u = a17 + b17*z + c17*z^2
[a17,b17,c17,res17,Sres17] = parabola(SDC,CSC);
% (18) u = a18 + b18*z + c18*z^2
[a18,b18,c18,res18,Sres18] = parabola(SDP,CSP);
% (19) u = a19 + b19*t + c19*t^2
[a19,b19,c19,res19,Sres19] = parabola(GPA,CSA);
% (20) u = a20 + b20*t + c20*t^2
[a20,b20,c20,res20,Sres20] = parabola(GPC,CSC);
% (21) u = a21 + b21*t + c21*t^2
[a21,b21,c21,res21,Sres21] = parabola(GPP,CSP);
%---------PLOTS, pvalues and confidence intervals
% 1: for thickness vs. small diameter
%figure
%pdt = ((kk-6)/3)*(Sres1/(Sres5+Sres9) -1);
%pv  = pvalue(pdt,kk-6);
%plot(z,a1+b1*z+c1*z.*z,'g-',z,a5+b5*z+c5*z.*z,'r--',SDP,THP,'rx',...
%   z,a9+b9*z+c9*z.*z,'b-.',SDC,THC,'bo')
%axis([z0 z1 1. 2.5])
%title(sprintf('Wall thickness vs. outer diameter; 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) '))
%%--------------------------------
%% 2: for thickness vs. gap
%figure
%pgt = ((kk-6)/3)*(Sres3/(Sres7+Sres11) -1)
%pv  = pvalue(pgt,kk-6);
%plot(t,a3+b3*t+c3*t.*t,'g-',t,a7+b7*t+c7*t.*t,'r--',GPP,THP,'rx',...
%   t,a11+b11*t+c11*t.*t,'b-.',GPC,THC,'bo')
%axis([t0 t1 1. 2.5])
%title(sprintf('Wall thickness vs. inner diameter; 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) '))
%--------------------------------
%% 3: for small diameter vs. cross section
figure
pdc = ((kk-6)/3)*(Sres16/(Sres17+Sres18) -1)
pv  = pvalue(pdc,kk-6);
plot(z,a16+b16*z+c16*z.*z,'g-',z,a17+b17*z+c17*z.*z,'r--',SDP,CSP,'rx',...
   z,a18+b18*z+c18*z.*z,'b-.',SDC,CSC,'bo')
axis([z0 z1 u0 1.5])
title(sprintf('Cross Section vs. outer diameter; 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
pgc = ((kk-6)/3)*(Sres19/(Sres20+Sres21) -1)
pv  = pvalue(pgc,kk-6);
plot(t,a19+b19*t+c19*t.*t,'g-',t,a20+b20*t+c20*t.*t,'r--',GPP,CSP,'rx',...
   t,a21+b21*t+c21*t.*t,'b-.',GPC,CSC,'bo')
axis([t0 t1 u0 1.5])
pv = .1
title(sprintf('Cross Section vs. inner diameter; 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) '))
%--------------------------------
%
%
%% plot(1,1): thickness vs. small diameter: regressions for all, patients, controls
%% plot(2,1):       "         gap         :              "
%% plot(1,2): ratio     vs. small diameter:
%% plot(2,2): ratio     vs.   gap
%figure
%plot(SDP,res5,'x',SDC,res9,'o')
%
%figure
%plot(SDA,res1,'+')
%%axis([z0 z1 y0 y1])
%%subplot(2,2,1);
%%subplot(2,2,2); 
%%subplot(2,2,3);
%figure
%plot(z,a2+b2*z+c2*z.*z,z,a6+b6*z+c6*z.*z,SDP,RRP,'x',...
%   z,a10+b10*z+c10*z.*z,SDC,RRC,'o')
%axis([z0 z1 x0 x1])
%%subplot(2,2,4);
%figure
%plot(t,a4+b4*t+c4*t.*t,t,a8+b8*t+c8*t.*t,GPP,RRP,'x',...
%   t,a12+b12*t+c12*t.*t,GPC,RRC,'o')
%axis([t0 t1 x0 x1])
%%subplot(2,2,5);
%figure
%plot(t,a13+b13*t+c13*t.*t,t,a14+b14*t+c14*t.*t,GPP,SDP,'x',...
%   t,a15+b15*t+c12*t.*t,GPC,SDC,'o')
%axis([t0 t1 z0 z1])
%
%% various quantities:
%% compute F-statistic
%% 3: for ratio vs. sm.diam.
%pdr = ((kk-6)/3)*(Sres2/(Sres6+Sres10) -1)
%% 4: for ratio vs. gap
%pgr = ((kk-6)/3)*(Sres4/(Sres8+Sres12) -1)
%% 5: for small diameter vs. gap
%pgd = ((kk-6)/3)*(Sres13/(Sres14+Sres15) -1)
%%
