FOCUS AREAS

Categorical Data Analysis: Chapter 13

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

   data colds;
      input sex $ residence $ periods count @@;
      datalines;  
   female rural 0 45 female rural 1 64  female rural 2 71
   female urban 0 80 female urban 1 104 female urban 2 116
   male rural   0 84 male   rural 1 124 male   rural 2 82
   male urban   0 106 male  urban 1 117 male   urban 2 87
   ;
   run;

   proc catmod; 
      weight count;
      response means;
      model periods = sex residence sex*residence /freq prob;     
   run;

   proc catmod; 
      weight count;
      response means;
      model periods = sex residence;     
   run; 

   proc catmod;
      population sex residence; 
      weight count;
      response means;
      model periods = sex;     
   run;    

   data cpain;
      input dstatus $ invest $ treat $ status $ count @@;
      datalines;
   I   A active poor  3 I   A active  fair 2 I A active moderate 2  
   I   A active good  1 I   A active  excel 0 
   I   A placebo poor 7 I   A placebo fair 0 I A placebo moderate 1  
   I   A placebo good 1 I   A placebo excel 1 
   I   B active poor  1 I   B active  fair 6 I B active moderate 1  
   I   B active good  5 I   B active  excel 3 
   I   B placebo poor 5 I   B placebo fair 4 I B placebo moderate 2  
   I   B placebo good 3 I   B placebo excel 3 
   II  A active poor  1 II  A active  fair 0 II A active moderate 1  
   II  A active good  2 II  A active  excel 2 
   II  A placebo poor 1 II  A placebo fair 1 II A placebo moderate 0  
   II  A placebo good 1 II  A placebo excel 1 
   II  B active poor  0 II  B active  fair  1 II B active moderate 1  
   II  B active good  1 II  B active  excel 6 
   II  B placebo poor 3 II  B placebo fair 1 II B placebo moderate 1  
   II  B placebo good 5 II  B placebo excel 0 
   III A active poor  2 III A active  fair 0 III A active moderate 3  
   III A active good  3 III A active  excel 2 
   III A placebo poor 5 III A placebo fair 0 III A placebo moderate 0  
   III A placebo good 8 III A placebo excel 1 
   III B active poor  2 III B active  fair 4 III B active moderate 1  
   III B active good 10 III B active  excel 3 
   III B placebo poor 2 III B placebo fair 5 III B placebo moderate 1  
   III B placebo good 4 III B placebo excel 2 
   IV  A active poor  8 IV  A active  fair 1 IV A active moderate 3  
   IV  A active good  4 IV  A active  excel 0 
   IV  A placebo poor 5 IV  A placebo fair 0 IV A placebo moderate 3  
   IV  A placebo good 3 IV  A placebo excel 0 
   IV  B active poor  1 IV  B active  fair 5 IV B active moderate 2  
   IV  B active good  3 IV  B active  excel 1 
   IV  B placebo poor 3 IV  B placebo fair 4 IV B placebo moderate 3  
   IV  B placebo good 4 IV  B placebo excel 2 
   ;

   proc catmod order=data;
     weight count;
     response 1 2 3 4 5;
     model status=dstatus|invest|treat;
   run;   

   proc catmod order=data;
      weight count;
      response 1 2 3 4 5;
      model status=dstatus|invest|treat@2; 
   run;

   proc catmod order=data;
      weight count;
      response 1 2 3 4 5;
      model status=dstatus invest treat;
   run;

   proc catmod order=data;
      weight count;
      response 1 2 3 4 5;
      model status=dstatus invest treat; 
      contrast 'I versus II'    dstatus 1 -1  0;
      contrast 'I versus III'   dstatus 1  0 -1;
      contrast 'I versus IV'    dstatus 2  1  1; 
      contrast 'II versus III'  dstatus 0  1 -1;
      contrast 'II versus IV'   dstatus 1  2  1;
      contrast 'III versus IV'  dstatus 1  1  2; 
      contrast 'dstatus'  dstatus 1 0 0 ,
                          dstatus 0 1 0 ,
                          dstatus 0 0 1 ;
   run;

   data byss;
      input workplace $ em_years $ smoking $ status $ count @@;
      datalines ;
   dusty        <10   yes  yes 30  dusty     <10   yes no  203
   dusty        <10   no   yes  7  dusty     <10   no  no  119
   dusty       >=10   yes  yes 57  dusty     >=10  yes no  161
   dusty       >=10   no   yes 11  dusty     >=10  no  no   81
   notdusty     <10   yes  yes 14  notdusty  <10   yes no 1340
   notdusty     <10   no   yes 12  notdusty  <10   no  no 1004
   notdusty    >=10   yes  yes 24  notdusty  >=10  yes no 1360
   notdusty    >=10   no   yes 10  notdusty  >=10  no  no  986
   ; 
   run;

   proc catmod order=data;
      weight count;
      response marginals;
      model status = workplace|em_years|smoking;
   run; 

   proc catmod order=data;
       weight count;
       response marginals;
       model status = workplace|em_years workplace|smoking
                      em_years|smoking;
   run;  

    proc catmod order=data;
       weight count ;
       response marginals;
       model status = workplace|em_years workplace|smoking;
    run;     

   proc catmod order=data;
       weight count;
       response marginals;
       model status = workplace
                     em_years(workplace='dusty')
                     em_years(workplace='notdusty')                   
                     smoking(workplace='dusty')
                     smoking(workplace='notdusty') ;                   
   run;

   proc catmod order=data;
      weight count ;
      response marginals;
      model status = workplace
                     em_years(workplace='dusty')
                     smoking(workplace='dusty') / pred;
   run;

   data pain (drop=h0-h8);
      input center initial $ treat $ h0-h8;
      array hours h0-h8;
      do i=1 to 9;
         no_hours=i-1; count=hours(i); output;
      end;    
      datalines; 
   1 some placebo  1 0 3 0 2 2 4 4 2 
   1 some treat_a  2 1 0 2 1 2 4 5 1 
   1 some treat_b  0 0 0 1 0 3 7 6 2 
   1 some treat_ba 0 0 0 0 1 3 5 4 6 
   1 lot  placebo  6 1 2 2 2 3 7 3 0 
   1 lot  treat_a  6 3 1 2 4 4 7 1 0 
   1 lot  treat_b  3 1 0 4 2 3 11 4 0  
   1 lot  treat_ba 0 0 0 1 1 7 9 6 2 
   2 some placebo  2 0 2 1 3 1 2 5 4 
   2 some treat_a  0 0 0 1 1 1 8 1 7 
   2 some treat_b  0 2 0 1 0 1 4 6 6 
   2 some treat_ba 0 0 0 1 3 0 4 7 5 
   2 lot  placebo  7 2 3 2 3 2 3 2 2 
   2 lot  treat_a  3 1 0 0 3 2 9 7 1 
   2 lot  treat_b  0 0 0 1 1 5 8 7 4 
   2 lot  treat_ba 0 1 0 0 1 2 8 9 5 
   3 some placebo  5 0 0 1 3 1 4 4 5 
   3 some treat_a  1 0 0 1 3 5 3 3 6 
   3 some treat_b  3 0 1 1 0 0 3 7 11 
   3 some treat_ba 0 0 0 1 1 4 2 4 13 
   3 lot  placebo  6 0 2 2 2 6 1 2 1 
   3 lot  treat_a  4 2 1 5 1 1 3 2 3 
   3 lot  treat_b  5 0 2 3 1 0 2 6 7 
   3 lot  treat_ba 3 2 1 0 0 2 5 9 4 
   4 some placebo  1 0 1 1 4 1 1 0 10 
   4 some treat_a  0 0 0 1 0 2 2 1 13
   4 some treat_b  0 0 0 1 1 1 1 5 11 
   4 some treat_ba 1 0 0 0 0 2 2 2 14 
   4 lot  placebo  4 0 1 3 2 1 1 2 2 
   4 lot  treat_a  0 1 3 1 1 6 1 3 6 
   4 lot  treat_b  0 0 0 0 2 7 2 2 9  
   4 lot  treat_ba 1 0 3 0 1 2 3 4 8 
   ;
   proc print data=pain(obs=9);
   run;

   proc catmod;
      weight count; 
      response 0 .125 .25 .375 .5 .625 .75 .875 1;
      model no_hours =  center initial treat
                        treat*initial;
   run;

   proc catmod;
      weight count; 
      response 0 .125 .25 .375 .5 .625 .75 .875 1;
      model no_hours = center initial 
                       treat(initial);  
   run;

   proc catmod;
      weight count; 
      response 0 .125 .25 .375 .5 .625 .75 .875 1;
      model no_hours = center initial 
                       treat(initial);  
      contrast 'lot: a-placebo'  treat(initial) -1  1  0  0  0  0 ;
      contrast 'lot: b-placebo'  treat(initial) -1  0  1  0  0  0 ; 
      contrast 'lot: ba-placebo' treat(initial) -2 -1 -1  0  0  0 ; 
      contrast 'lot: ba-a'       treat(initial) -1 -2 -1  0  0  0 ;
      contrast 'lot: ba-b'       treat(initial) -1 -1 -2  0  0  0 ;
      contrast 'some:a-placebo'  treat(initial)  0  0  0 -1  1  0 ;
      contrast 'some:b-placebo'  treat(initial)  0  0  0 -1  0  1 ; 
      contrast 'some:ba-placebo' treat(initial)  0  0  0 -2 -1 -1 ; 
      contrast 'some:ba-a'       treat(initial)  0  0  0 -1 -2 -1 ;
      contrast 'some:ba-b'       treat(initial)  0  0  0 -1 -1 -2 ;
   run;

   proc catmod;
      weight count; 
      response 0 .125 .25 .375 .5 .625 .75 .875 1;
      model no_hours = center initial 
                       treat(initial);  
      contrast 'interact:a-placebo'  treat(initial) -1  1  0  1 -1  0 ;
      contrast 'interact:b-placebo'  treat(initial) -1  0  1  1  0 -1 ; 
      contrast 'interact:ba-placebo' treat(initial) -2 -1 -1  2  1  1 ; 
      contrast 'average:a-placebo'   treat(initial) -1  1  0 -1  1  0 ;
      contrast 'average:b-placebo'   treat(initial) -1  0  1 -1  0  1 ; 
      contrast 'average:ba-placebo'  treat(initial) -2 -1 -1 -2 -1 -1 ;
      contrast 'average:ba-a'        treat(initial) -1 -2 -1 -1 -2 -1 ;
      contrast 'average:ba-b'        treat(initial) -1 -1 -2 -1 -1 -2 ; 
      contrast 'interaction'         treat(initial) -1  1  0  1 -1  0 ,
                                     treat(initial) -1  0  1  1  0 -1 , 
                                     treat(initial) -2 -1 -1  2  1  1 ; 
      contrast 'treatment effect'    treat(initial) -1  1  0 -1  1  0 ,
                                     treat(initial) -1  0  1 -1  0  1 , 
                                     treat(initial) -2 -1 -1 -2 -1 -1 ; 
   run;

   data wbeing;
        input #1 b1-b5 _type_ $ _name_ $8. #2 b6-b10;
        datalines;
    7.24978   7.18991   7.35960   7.31937   7.55184 parms
    7.93726   7.92509   7.82815   7.73696   8.16791 
    0.01110   0.00101   0.00177  -0.00018  -0.00082 cov       b1
    0.00189  -0.00123   0.00434   0.00158  -0.00050 
    0.00101   0.02342   0.00144   0.00369   0.25300 cov       b2
    0.00118  -0.00629  -0.00059   0.00212  -0.00098 
    0.00177   0.00144   0.01060   0.00157   0.00226 cov       b3
    0.00140  -0.00088  -0.00055   0.00211   0.00239 
   -0.00018   0.00369   0.00157   0.02298   0.00918 cov       b4
   -0.00140  -0.00232   0.00023   0.00066  -0.00010 
   -0.00082   0.00253   0.00226   0.00918   0.01921 cov       b5
    0.00039   0.00034  -0.00013   0.00240   0.00213 
    0.00189   0.00118   0.00140  -0.00140   0.00039 cov       b6
    0.00739   0.00019   0.00146  -0.00082   0.00076 
   -0.00123  -0.00629  -0.00088  -0.00232   0.00034 cov       b7
    0.00019   0.01172   0.00183   0.00029   0.00083 
    0.00434  -0.00059  -0.00055   0.00023  -0.00013 cov       b8
    0.00146   0.00183   0.01050  -0.00173   0.00011 
    0.00158   0.00212   0.00211   0.00066   0.00240 cov       b9
   -0.00082   0.00029  -0.00173   0.01335   0.00140 
   -0.00050  -0.00098   0.00239  -0.00010   0.00213 cov      b10
    0.00076   0.00083   0.00011   0.00140   0.01430 
   ;

   proc catmod data=wbeing;
        response read b1-b10;
        factors sex $ 2, age $  5 / 
             _response_ = sex|age 
        profile = (female '25-34',
                   female '35-44',
                   female '45-54',
                   female '55-64',
                   female '65-74',
                   male   '25-34',
                   male   '35-44',
                   male   '45-54',
                   male   '55-64',
                   male   '65-74');
        model _f_ = _response_; 
   run;


   proc catmod data=wbeing;
      response read b1-b10;
      factors sex $ 2, age $ 5 / 
         _response_ = sex age  
         profile = (male '25-34' ,
                    male '35-44',
                    male '45-54' ,
                    male '55-64',
                    male '65-74' ,
                    female '25-34',
                    female '35-44',
                    female '45-54',
                    female '55-64' ,
                    female '65-74');
      model _f_ = _response_;
      contrast '25-34 vs. 35-44' all_parms  0 0 1 -1  0  0;
      contrast '25-34 vs. 45,54' all_parms  0 0 1  0 -1  0;
      contrast '25-34 vs. 55,64' all_parms  0 0 1  0  0 -1;
      contrast '25-34 vs. 65,74' all_parms  0 0 2  1  1  1;
   run;

   proc catmod data=wbeing;
      response read b1-b10;
      factors sex $ 2, age $ 5 / 
         _response_ = sex age  
         profile = (male '25-34' ,
                    male '35-44',
                    male '45-54' ,
                    male '55-64',
                    male '65-74' ,
                    female '25-34',
                    female '35-44',
                    female '45-54',
                    female '55-64' ,
                    female '65-74');
      model _f_ = _response_;
      contrast '25-64 the same'    all_parms  0 0 1 -1  0  0,
                                   all_parms  0 0 1  0 -1  0,
                                   all_parms  0 0 1  0  0 -1;
   run;

   proc catmod data=wbeing; 
      response read b1-b10;
      factors sex $ 2, age $  5 / 
         _response_ = sex age 
            profile = (female '25-34',
                       female '35-44',
                       female '45-54',
                       female '55-64',
                       female '65-74',
                       male   '25-34',
                       male   '35-44',
                       male   '45-54',
                       male   '55-64',
                       male   '65-74');
      model _f_ = ( 1 0 0 ,
                    1 0 0 ,
                    1 0 0 ,
                    1 0 0 ,
                    1 0 1 ,
                    1 1 0 ,
                    1 1 0 ,
                    1 1 0 ,
                    1 1 0 ,
                    1 1 1 ) (1='Intercept', 2='Sex', 3='65-74')
                            / pred;

   proc freq data=cpain; 
      weight count; 
      tables dstatus*invest*treat*status/ measures; 
   run; 

   data MannWhitney;
      input b1-b8 _type_ $ _name_ $8.;
      datalines;
   .6000  .6011  .6042 .8389 .5130 .5947 .5000 .4922 parms  
   .03091 .0000  .0000 .0000 .0000 .0000 .0000 .0000 cov b1
   .0000  .00918 .0000 .0000 .0000 .0000 .0000 .0000 cov b2
   .0000  .0000  .3280 .0000 .0000 .0000 .0000 .0000 cov b3
   .0000  .0000  .0000 .0084 .0000 .0000 .0000 .0000 cov b4
   .0000  .0000  .0000 .0000 .0129 .0000 .0000 .0000 cov b5
   .0000  .0000  .0000 .0000 .0000 .0093 .0000 .0000 cov b6 
   .0000  .0000  .0000 .0000 .0000 .0000 .0101 .0000 cov b7
   .0000  .0000  .0000 .0000 .0000 .0000 .0000 .0112 cov b8
   ;

   proc catmod data=MannWhitney;
        response read b1-b8;
        factors diagnosis $ 4 , invest $ 2 / 
           _response_ = diagnosis invest
        profile = (I   A,
                   I   B,
                   II  A,
                   II  B,
                   III A,
                   III B,
                   IV  A,
                   IV  B);
        model _f_ = _response_ / cov;
   run;

   proc catmod data=MannWhitney;
      response read b1-b8;
      factors diagnosis $ 4 , invest $ 2 / 
         _response_ = diagnosis invest
      profile = (I   A,
                 I   B,
                 II  A,
                 II  B,
                 III A,
                 III B,
                 IV  A,
                 IV  B);
      model _f_ = ( 1 ,
                    1 ,
                    1 ,
                    1 ,
                    1 ,
                    1 ,
                    1 ,
                    1 ) ;
   run;


Statistics and Operations Research