options ls=64;
%include 'gamfrail.sas';
/**********************************************************/
/* Skin graft data from                                  */
/* http://www.biostat.mcw.edu/homepgs/klein/13.5.1.html  */
/* Variables are patient, time, status, and quality of   */
/* match.  I'll use K and M's gamma frailty macro,       */
/* http://www.biostat.mcw.edu/software/gamfrail          */
/* Macro needs variables in this order: Time, status,    */
/* matching variable, covarites, and so I I move some    */
/* variables around.  Arguments are dat set name,        */
/**********************************************************/
data skin; infile 'skin.dat';
   input Patient Tm st match; 
   gp=patient; mtch=match; drop patient match; run;
title1 'Gamma frailty model applied to skin graft data';
proc IML symsize=800; option={y,n,n,n,n,n} ;
   %gamfrail(skin,{mtch},{0},.01,.0001,.1,option,o1,o2,o3);
/**********************************************************/
/* VA bladder tumor data from                            */
/* http://lib.stat.cmu.edu/datasets/Andrews/T45.1, after */
/* some reformatting.  Variables are treatment group,    */
/* followuptime, initial number of tumors, initial size  */
/* of tumors, and up to four recurrence times.           */
/**********************************************************/
data tumor ; infile 'tumor.dat' ;
   input group futime number size rt1 rt2 rt3 rt4; run;
title1 'Gamma frailty model applied to tumor Recurrence';
data tumor ; set tumor;
   time=0; status=1; person=_n_; gp=group;
   tl=futime; tn=0; 
   if rt1=. then; else 
      do; tn=tn+1; time=rt1; tl=tl-time; output; end;
   if rt2=. then; else 
      do; tn=tn+1; time=rt2-rt1; tl=tl-time; output; end;
   if rt3=. then; else 
      do; tn=tn+1; time=rt3-rt2; tl=tl-time; output; end;
   if rt4=. then; else 
      do; tn=tn+1; time=rt4-rt3; tl=tl-time; output; end;
   if tl>0 then do;status=0; tn=tn+1; time=tl; output; end;
   tn1=0; if tn=1 then tn1=1; tn2=0; if tn=2 then tn2=1;
   drop rt1 rt2 rt3 rt4 tl futime group number size;
run;
data tumor;  set tumor ; tn1=0; if tn=1 then tn1=1; tn2=0; 
   if tn=2 then tn2=1; if tn>3 then delete ; drop tn; run;
proc IML symsize=200; 
   cn={gp tn1 tn2}; init={0 0 0};
   option={y,n,n,n,n,n} ;
   %gamfrail(tumor,cn,init,.01,.00001,.1,option,o1,o2,o3);
/**********************************************************/
/* Data from http://lib.stat.cmu.edu/datasets/Arsenic     */
/* reformatted into ASCII on the course home page.  Data  */
/* reflect arsenic levels in toenail clippings; covar-    */
/* iates include age, sex, categorical measures of quant- */
/* ities used for drinking and cooking, arsenic in the    */
/* water, and arsenic in the nails.                       */
/**********************************************************/
title1 'Tobit model for arsenic in water';
data arsenic; infile 'arsenic.dat'; 
   input age sex drink cook arswater arsnails; status=1; 
   if arswater <.0001 then do; arswater=.0001; status=0; end;
   lower=.; if status=1 then lower=arswater;
   lnail=log(arsnails); lwater=log(arswater); run;
proc reg data=arsenic; model lwater=lnail; run;
proc lifereg data=arsenic;
   model (lower,arswater)=lnail/d=lognormal; run;
/**********************************************************/
/* Danger on injuries from lifting, from                  */
/* www2.tltc.ttu.edu/westfall/images/5349/alr_data.htm    */
/* Model the number of days lost as a latent normal var-  */
/* iable, rounded to zero if negative.  Important vari-   */
/* ables are alr (an index of lifting) and work days lost.*/
/**********************************************************/
title1 'Tobit model for Injuries from lifing';
data lifting;
   infile 'alr_data.dat';
   input subject hours status alr sex age weight armst 
      backst dynamend abdep shoulht knuckht dayslost;
   if dayslost=0 then lower = .; else lower = dayslost;
run;
/* Below dayslost is days of work lost to an injury, and 
ALR is an index of how much lifting was involved. */
proc lifereg;
   model (lower,dayslost) = alr/dist=normal;
run;
/**********************************************************/
/* Ulcer data from Collette, on web site.  Response is    */
/* time until ulcer diagnosis, many of which were made    */
/* at scheduled visits.  Variables are age, duration      */
/* (1 for short, 2 for long), treatment (A or B), time    */
/* (in months), and status.  Regular visits took place at */
/* 6 month intervals.                                     */
/* We will write out a data set with one of the three     */
/* following pattern of lines:                            */
/* R period                                               */
/* 1 1                                                    */
/* ----                                                   */
/* 0 1                                                    */
/* 1 2                                                    */
/* ----                                                   */
/* 0 1                                                    */
/* 0 2                                                    */
/* 1 3                                                    */
/* The final line corresponds to pi[i,3], which is by def-*/
/* inition 1; hence this line can be (and is) omitted.    */
/**********************************************************/
title1 'Ulcer recurrence data';
title2 'Incorrect answer ignoring interval censoring';
data ulcer; infile 'ulcerrecurrence.dat';
   input age duration ttt $ time result;
   trt=0; if ttt="B" then trt=1; run;
proc phreg data=ulcer; model time*result(1)=trt; run;
data ulcer; set ulcer;
   R=result-1; if time<6 then time=6; 
   if time>6 and time<12 then time=12;
   if time>=12 then do ; period=2; output; R=0; end;
   if time>=6 then do ; period=1; output; end;
run;
title2 'Correct answer accounting for interval censoring';
proc genmod data=ulcer descending; class period ; 
   model r=trt period/dist=binomial link=cloglog; run;