FOCUS AREAS

Categorical Data Analysis: Chapter 11

  

   options nodate nonumber ps=200 ls=80 formdlim=' ';

   data bacteria;
      input dose status $ count @@;
      ldose=log(dose);
      sq_ldose=ldose*ldose;
      datalines;
   1200      dead   0    1200      alive 5
   12000     dead   0    12000     alive 5
   120000    dead   2    120000    alive 3
   1200000   dead   4    1200000   alive 2
   12000000  dead   5    12000000  alive 1
   120000000 dead   5    120000000 alive 0
   ;
   proc print;
   run;

   proc logistic data=bacteria descending;
      freq count;
      model status = ldose sq_ldose / scale=none aggregate
            selection=forward start=1 details covb;
   run;

   data assay;
      input drug $ dose status $ count;
      int_n=(drug='n');
      int_s=(drug='s');
      ldose=log(dose);
      ldose_n=int_n*ldose;
      ldose_s=int_s*ldose;
      sqldose_n=int_n*ldose*ldose;
      sqldose_s=int_s*ldose*ldose;
      datalines;
   n 0.01   dead   0
   n 0.01   alive 30
   n  .03   dead   1
   n  .03   alive 29
   n  .10   dead   1
   n  .10   alive  9
   n  .30   dead   1
   n  .30   alive  9
   n 1.00   dead   4
   n 1.00   alive  6
   n 3.00   dead   4
   n 3.00   alive  6
   n 10.00  dead   5
   n 10.00  alive  5
   n 30.00  dead   7
   n 30.00  alive  3
   s   .30  dead   0
   s   .30  alive 10
   s  1.00  dead   0
   s  1.00  alive 10
   s 3.00   dead   1
   s 3.00   alive  9
   s 10.00  dead   4
   s 10.00  alive  6
   s 30.00  dead   5
   s 30.00  alive  5
   s 100.00 dead   8
   s 100.00 alive  2
   ;

   proc logistic data=assay descending;
      freq count;
      model status = int_n int_s ldose_n ldose_s
                     sqldose_n sqldose_s  / noint
                     scale=none aggregate
                     start=4 selection=forward details;
      eq_slope: test ldose_n=ldose_s;
   run;

   proc logistic data=assay descending outest=estimate
                      (drop= intercept _link_ _lnlike_) covout;
      freq count;
      model status = int_n int_s ldose /
                     noint scale=none aggregate covb;
   run;

   proc iml;
     use estimate;
     start fieller;
     title 'Confidence Intervals';
     use estimate;
     read all into beta where (_type_='PARMS');
     beta=beta`;
     read all into cov where (_type_='COV');
     ratio=(k`*beta) / (h`*beta);
     a=(h`*beta)**2-(3.84)*(h`*cov*h);
     b=2*(3.84*(k`*cov*h)-(k`*beta)*(h`*beta));
     c=(k`*beta)**2 -(3.84)*(k`*cov*k);
     disc=((b**2)-4*a*c);
     if (disc<=0 | a<=0) then do;
     print "confidence interval can't be computed", ratio;
     stop; end;
     sroot=sqrt(disc);
     l_b=((-b)-sroot)/(2*a);
     u_b=((-b)+sroot)/(2*a);
     interval=l_b||u_b;
     lname={"l_bound", "u_bound"};
     print "95 % ci for ratio based on fieller", ratio interval[colname=lname];
     finish fieller;
     k={ 1 -1 0 }`;
     h={ 0  0 1 }`;
     run fieller;
     k={-1 0 0 }`;
     h={ 0 0 1 }`;
     run fieller;
     k={ 0 -1 0 }`;
     h={ 0 0 1 }`;
     run fieller;

   data adverse;
      input diagnos $ dose status $ count @@;
      i_diagII=(diagnos='II');
      i_diagI= (diagnos='I');
      doseI=i_diagI*dose;
      doseII=i_diagII*dose;
      diagdose=i_diagII*dose;
      if doseI > 0 then ldoseI=log(doseI); else ldoseI=0;
      if doseII > 0 then ldoseII=log(doseII); else ldoseII=0;
      datalines;
   I    1   adverse  3 I    1 no 26
   I    5   adverse  7 I    5 no 26
   I    10  adverse 10 I   10 no 22
   I    12  adverse 14 I   12 no 18
   I    15  adverse 18 I   15 no 14
   II   1   adverse  6 II   1 no 26
   II   5   adverse 20 II   5 no 12
   II   10  adverse 26 II  10 no  6
   II   12  adverse 28 II  12 no  4
   II   15  adverse 31 II  15 no  1
   ;

   proc freq data=adverse;
      weight count;
      tables dose*status diagnos*status diagnos*dose*status /
             nopct nocol cmh;
   run;

   proc logistic data=adverse outest=estimate
                      (drop= intercept _link_ _lnlike_) covout;
      freq count;
      model status = i_diagI i_diagII doseI doseII  /
                     noint scale=none aggregate;
      eq_slope: test doseI=doseII;
   run;

   proc logistic data=adverse;
      freq count;
      model status = i_diagI i_diagII ldoseI ldoseII  /
                     noint scale=none aggregate;
      eq_slope: test ldoseI=ldoseII;
   run;

   proc iml;
      start fieller;
      title 'Confidence Intervals';
      use estimate;
      read all into beta where (_type_='PARMS');
      beta=beta`;
      read all into cov where (_type_='COV');
      ratio=(k`*beta) / (h`*beta);
      a=(h`*beta)**2-(3.84)*(h`*cov*h);
      b=2*(3.84*(k`*cov*h)-(k`*beta)*(h`*beta));
      c=(k`*beta)**2 -(3.84)*(k`*cov*k);
      disc=((b**2)-4*a*c);
      if (disc<=0 | a<=0) then do;
      print "confidence interval can't be computed", ratio;
      stop; end;
      sroot=sqrt(disc);
      l_b=((-b)-sroot)/(2*a);
      u_b=((-b)+sroot)/(2*a);
      interval=l_b||u_b;
      lname={"l_bound", "u_bound"};
      print "95 % ci for ratio based on fieller", ratio interval[colname=lname];
      finish fieller;
      k={ -1 0 0 0}`;
      h={  0 0 1 0}`;
      run fieller;
      k={ 0 -1 0 0}`;
      h={ 0  0 0 1}`;
      run fieller;


Statistics and Operations Research