/*  file: gts_alr3level.sas                      */
/*  old name: ~/gts/JQW1/cluster/alrtab6.sas                  */
/*  date: 07/8/2001                                           */
/*  name: John Preisser                                       */
/*  project: GTS (Green Tobacco Sickness)                     */
/*     "Detecting Patterns of Occupational Illness Clustering */
/*      with Alternating Logistic Regressions Applied to      */
/*      Longitudinal Data" by Preisser, Arcury and Quandt     */
/*      Amercian Journal of Epidemiology, 2003; 57:550-2.     */
/*  request: Apply ALR and produce results for Table 3 and 4  */
/*                                                            */
/*  input data files: ~/gts/datasets/daily2.ssd01             */
/*  output data files:                                        */

libname dat "/home/projects/gts/datasets";
libname formats "/home/projects/gts/sas_formats";

options ls=80 ps=58; 
options fmtsearch=(work.format, formats.formats, fmt.format);

proc format;

   value workty
   1 = 'prime        '
   2 = 'prime/barn'
   3 = 'barn'
   4 = 'top'
   5 = 'other'
       ;

   value hrc
   1 = 'O hours'
   2 = '.25-2.00 hours'
   3 = '2.25-4.00 hours'
   4 = '> 4 hours'
      ;
run;


data daily2;
    set dat.daily2;
    if (1 le cum_day le 34) then season=3;
    else if (35 le cum_day le 55) then season=2;
    else if cum_day ge 56 then season=1;
* 10 degree F;
    temp80=(temp-80)/10;
    if experien=1 then firstyr=1;
    else firstyr=0;
    if (experien=1 or experien=2) then lowexp=1;
    else lowexp=0;
    

run;
proc print data=daily2;
    where (experien=. or tobuser=. or workty=. or hrwork=. or hrwet=. or temp=. or humid=. or season=.);
    var siteid newpid experien tobuser workty hrwork hrwet temp humid season gts;
run;

proc sort data=daily2;
    by siteid newpid period;
run;

data wid;
   set daily2;
   by siteid newpid period;
   if first.siteid then wid=0;
   if first.newpid then wid+1;

   if hrwet > 0 then workwet=1;
   else workwet=0;


   if hrwet =0 then hrwet4=1;
   else if (0 lt hrwet le 2.00) then hrwet4=2;
   else if (2.00 lt hrwet le 4.00) then hrwet4=3;
   else if hrwet > 4 then hrwet4=4;

/*  Comment - workty already of daily2.ssd
   if worktype =1 then workty = 3;
   else if worktype=2 then workty = 1;
   else if worktype=3 then workty = 4;
   else if worktype=4 then workty = 2;
   else if worktype=5 then workty = 5;
*/
   
* combine barn and other into one category - make topping the reference category;

  if workty=1 then prime=1; else prime=0;
  if workty=2 then prbarn=1; else prbarn=0;
  if (workty=4 or workty=5) then othbarn=1; else othbarn=0;
   retain wid;
run;
/*
proc freq;
    table (prime prbarn othbarn)*gts;
run;
*/

* two ORs, same worker, or different worker;

ods output GEEEmpPEst = parestA;
ods trace on;
title 'first OR is within persons, 2nd OR is between persons';
title2 'Model A: adjusts for variables affecting exposure';
proc genmod data=wid descending;           
  class siteid newpid;
  model gts= prime prbarn othbarn lowexp
          /dist=bin noscale scale=1;
  repeated subject=siteid/logOr=NEST1 SUBCLUST=newpid;
run;
ods trace off;
ods listing;

data oddsA;
   set parestA;
   ORest = exp(Estimate);
   Uodds = exp(lowerCL);
   Lodds = exp(upperCL);
run;
proc print data=oddsA;
  title 'parameter estimates and empirical SEs for ALR-model A';
  var Parm Estimate stderr z probZ ORest Uodds Lodds;
run;



ods output GEEEmpPEst = parestB;
ods trace on;
title 'first OR is within persons, 2nd OR is between persons';
title2 'Model B: adjusts for exposure and adsorption variables';
proc genmod data=wid descending;           
  class siteid newpid;
  model gts= prime prbarn othbarn lowexp
             workwet temp80 tobuser
          /dist=bin noscale scale=1;
  repeated subject=siteid/logOr=NEST1 SUBCLUST=newpid;
run;
ods trace off;
ods listing;

data oddsB;
   set parestB;
   ORest = exp(Estimate);
   Uodds = exp(lowerCL);
   Lodds = exp(upperCL);
run;
proc print data=oddsB;
  title 'parameter estimates and empirical SEs for ALR-model B';
  var Parm Estimate stderr z probZ ORest Uodds Lodds;
run;