options ls=80 nodate pageno=1 NOFMTERR; 
 
 options MPRINT; * checks macro values;

  %include "X:\wthomas On My Mac\ PC SAS examples\SAS Macros\mackerel.sas";
  %include 
"C:\Documents and Settings\Administrator\My Documents\SAS Macros\bootstrapmacros.sas";

 "C:\Documents and Settings\Administrator\Desktop\";

 *  scatterplot matrix with histograms down diagonal;

 ods graphics on;
  proc corr data=A nomiss plots=matrix(histogram);
    var start_pay start_month start_age experience education;
  run;
  ods graphics off;

   ODS graphics on;
 proc univariate noprint data=a; *stacked histograms;
   class age_group;
   var PCE;
   histogram PCE  / nrows=3;
 run; 
 ODS graphics off;
ODS graphics on;
 proc univariate noprint data=pubh.mult_myeloma;
   var hemoglobin;
   histogram   hemoglobin / RTINCLUDE height=3 endpoints = 4 to 16 by 1; *height=axis labels cex;
    inset q1 median mean q3 / position=NW noframe;
 ODS graphics off;


  Proc SGplot data=  ;
   reg  y = glu   x= bodyfat /group = IR MARKERATTRS=(symbol="circlefilled" size=4);
   refline  17.6 37.5 63.9 78.32 / axis=x;

 Proc SGplot data=  ;
   loess  y = glu  x= percent-bodyfat /group = IR ;
 Proc SGplot data=  ; * spaghetti plot;
   series x=year y=income / group =family_id LineAttrs= (pattern=1 color="black");

 Proc SGscatter data=baseline;
   matrix  sysbp  bmi heart_rate totchol glucose /MARKERATTRS=(size=4); *diagonal=(histogram) ;

 Proc SGpanel data=baseline;
   panelby male;
   loess  y = sysbp  x= age 
    /degree=2 smooth=.25 MARKERATTRS=(size=6) lineattrs=(color=red thickness=4);
              
 ******************************************** 
  misc statements
  array g[5] gfr1-gfr5;
  array m[r,c] row11 - row1c  row21-row2c  etc;

   * to hard code a date: start = MDY(1,1,2003)
    or   cutoff = '01JAN2003'D;

 *************************************
 data options

 RENAME = ( oldname1=newname1 oldname2=newname2 ) ;

  data a;
    set b(RENAME = ( oldname1=newname1 oldname2=newname2 );

  proc print data=pubh.nhanes_all20 (obs=5);
  * option (firstobs=K obs=J) prints obs from K to J, so need J > K; 
 ******************************************** ;


  *To write data To_R to a CSV;

 proc export data = To_R
  outfile = "C:\Documents and Settings\Administrator\Desktop\To_R.csv"
  DBMS=CSV replace;

  * To write  data To_Investigator to an Excel spreadsheet;

  * Character MRN: in data step, FORMAT mrn :$10. ;
  
 proc export data=To_Investigator
   outfile="C:\Documents and Settings\Administrator\Desktop\To_Investigator.xls"
   dbms=excel replace;

 proc import datafile="X:\wthomas On My Mac\6470\Grades\HW3 grades.xls"
   out=midterm_grades
   dbms=excel replace;

 PROC IMPORT OUT= WORK.normal_ranges 
            DATAFILE= "C:\normal_ranges.xls" 
            DBMS=EXCEL REPLACE;
    * SHEET="normal_ranges$"; 
     GETNAMES=YES;
     MIXED=YES;  
     SCANTEXT=YES;   
     USEDATE=YES;    
     SCANTIME=YES;   

  data cfrd.IFG00add; * 59 additional OGTT results, 3 Aug 09;
   infile "X:\wthomas On My Mac\ PC Projects\CFRD1\CF database\OGTT_96_05_AdditionalTests.csv"
      DSD dlm="," firstobs=2 missover lrecl=300;
   input NewID :$10. OGTT_Date :MMDDYY10.  min000 min030 min060 min090 min120   
    Sex$ BirthDate :MMDDYY10. exercise_time :time10. ;
   format BirthDate OGTT_Date date9. ; *easier to read;
   format  NewID :$10. ;
   format exercise_time time10. ;


 Convert character data to numeric or date: (LSB 10.8, aer00.sas)

     alb = INPUT(UMALI_mg_L, 12.);
     PFTdate = INPUT(PFT_date, MMDDYY10.); * this was format of char variable;
     format PFTdate  MMDDYY10. ;

  Convert numeric to character:   newvar = PUT(oldvar, format);

    predose = put(Dayordinal, 4.); * converts numeric var to character.;
    prescribed = input(dose, 4.); * convert char to numeric; 

********************************************;
   BMI = body weight in kilograms/height in meters squared
********************************************;
 *  Rank only one variable per PROC.  With > 1 var, ranks assigned only to all non-missing obs;
  Proc Rank data=C out=D groups=3;
     var factor1 ;
     ranks Qfactor1 ;

  proc tabulate data=A4;
   class high_F2 bmi_class trt;
   var delta_f2 log_f20;
   table high_F2* bmi_class, trt*(delta_f2*mean log_f20*mean);
   table high_F2* bmi_class, trt*(delta_f2*n  log_f20*n );

 proc logistic descending data=b; 
  class age bact_status / param=GLM ; * this option gives 0/1 indicator variables;
  freq kids;
  model death = age bact_status age*bact_status / CLodds = PL  lackfit
    aggregate scale=N ;  * calculates scale to check overdispersion, but does not adjust;
  oddsratio bact_status / at (age=15 25 35 45);


********************************************;
 proc power;
   onecorr test=pearson npartialvars=2 power=. ntotal=116 corr=.27 ;
   twosamplemeans test=diff meandiff = . stddev = 1.89 npergroup = 6  power = .80; 
   twosamplemeans test=diff meandiff = . stddev= 5.8 2.8 4.4 1.8 226 0.11 1.0
        npergroup = 250  power = .80 alpha=.025;
   pairedmeans test=diff meandiff=25 PAIREDSTDDEVS= (30 38)  corr=0.6 power=.80  npairs=. ; *SD1, SD2;
   pairedmeans test=diff meandiff=10 stddev=7.7 9.1  corr=0.6 power=.80  npairs=. ;
   TWOSAMPLEFREQ test=fisher power=. groupns= (20 22) GROUPPROPORTIONS= .75 | .225 .25 .3 .375 .4 .50;
   ONESAMPLEMEANS test=t dist=normal mean=.6 stddev=1 power=. ntotal=45 ; 
   run;

********************************************;
 proc SQL;
   create table B as
   select distinct newid,  count(GTT_date) as IGT_events
   from A
   group by newid;

  proc sql;
    create table htx2 as
    select * from aer.MA_comparison, htx1 
    where htx1.MRN=MA_comparison.newid and (htx1.startdate le MA_comparison.collectdate le htx1.enddate);

   ********************************************;

  * Back transforming Proc GLM LSmeans from log;
 proc GLM data=a;
    class dose;
    model log_response = dose;
    lsmeans dose / CL;
    ods output LSMeanCL=logCI;

 proc print data=logCI;
*
Obs   Effect    Dependent     dose        LowerCL         LSMean        UpperCL

 1     dose    log_response   10         0.273384       1.495460       2.717537
 2     dose    log_response   25         1.328628       2.550704       3.772780
 3     dose    log_response   50         2.185752       3.407828       4.629905
 4     dose    log_response   100        5.464444       6.686520       7.908596
 5     dose    log_response   250       10.036386      11.258462      12.480539
 ;
 data b;
   set logCI;
   geometric_mean = exp(lsmean);
   left_CI=exp(LowerCL);
   right_CI = exp(UpperCL);
 proc print data=m;
   var Dependent dose geometric_mean  left_CI   right_CI;
********************************************;
  * Back transforming Proc Mixed  LSmeans from log;
 proc mixed data=&data;
    title3 "Data = &data, Outcome = &endpoint, &hour hours";
   class fiber;
   model &endpoint = fiber;
   lsmeans fiber / diff CL;
   ODS output lsmeans=M;
   ODS select tests3;

 data means;
   set M;
   geometric_mean = exp(estimate);
   left_CI=exp(Lower);
   right_CI = exp(Upper);
   keep fiber geometric_mean left_CI right_CI;

******************************************************;
 Proc Transpose data=long out=wide prefix =income_ ;
   ID year; * values become names of variables in output data (wide), combined with prefix;
   VAR income; * repeated measurements (will be transposed);
   BY family_id; * subject identifier: one observation for each level of BY variable(s);

proc transpose data=sig out=tab;
    var <values in column, to be put in rows>;
    id <characterstrings in column, to become variable names>;
  * see pre_cvh02.sas ;

  * counting obs with transpose, instead of SQL;
 * TRANSPOSE from ogtt7.sas (Frohnert);
  proc sort data=short; * unique values of BY variables, one obs per ogtt_date;
   by newid pair  ifg_status;

 proc transpose data= short out=sideways;
   var ogtt_date;
   by newid pair  ifg_status;

 data rates2;
   set sideways;
   ogtt_count1=N(OF col1-col9); * count ogtt_dates per patient;

 run;
********************************************;
  ods output CovParms=covariance lsmeans=xbar  Diffs=delta;
  proc mixed noclprint data=wiat1.two_groups_longit;
    class patid week treatment depression_meds transplant_organ diabetes delay;
   model comp_itt = 
        treatment depression_meds transplant_organ diabetes delay 
        treatment*week / ddfm=kenwardroger;
    random patid / subject=patid v vcorr; 
    lsmeans treatment*week  /diff  ; 
********************************************;
 * compares paired dm - controls over regional FA and posterior TBSS measurements ;

 proc mixed data=kodl.regional covtest;
  class proc pair dm region;
  model FA= dm region dm*region / ddfm=kenwardroger;
  random pair / subject=proc;
  repeated / subject=proc(pair) type=UN r rcorr ; 
  lsmeans dm;
  lsmeans dm*region / slice=region;


  %include 
 "C:\Documents and Settings\Administrator\My Documents\SAS Macros\GlobTest.sas";


 ******************************************** ;
   %Gpdf(figname=dose_response_resid.pdf,height=4, cex= 1)
 proc sort data=g;
   by dose;
 proc BOXPLOT data=g;
    plot residual*dose /boxstyle=skeletal noserifs boxconnect=mean vref=0;

 %points(pdf_gif=pdf ,figname=infmort3S.pdf , cex=1 )
 symbol3 font=marker value=C h=.6;
 symbol4  value=square h=1; 
 proc Gplot data=two;
  plot studentized_resid * fitted = region/ vref=0;

 ******************************************** 
 Proc REG plots -- macro doesn't work

  Goptions reset=all vsize=5in ftext=simplex device=pdf lfactor=2
   noborder gsfname=graphout gsfmode= replace;
   filename graphout "C:\Documents and Settings\William Thomas\Desktop\Florida05.pdf";
   symbol1 value=circle  height=1;
 proc reg data=pred;
   model log_buchanan =  log_pop log_votes p_bush p_nader white_percent logpblack  
           logphispanic income percent_college percent_65;* / selection=cp ;
   plot student. * predicted.;
*******************************************************

*  scatterplot with smooth;
 Goptions reset=all vsize=4in ftext=simplex device=pdf lfactor=2 
   noborder gsfname=graphout gsfmode= replace; 
   filename graphout "C: ... minnows.pdf"; 
  symbol1 value=circle h=1.5; 
 Proc Gplot data=pubh.minnows ; 
  plot wt_g * length_cm / haxis = 4.5 to 7.5 by 0.5; 
  symbol1 interpol = sm60s width=2; 

 * interaction plot;*/

 ods graphics on;
 ods select lsmeans meanplot;  *print only stuff required for plot;

 proc glimmix data=test;
  class id a b;
  model y=a b a*b;
  random intercept / subject=id  type=un; * needed for individual SE estimates;
  lsmeans a*b / alpha=.32  plots=(meanplot(cl join sliceby=b));
 run;
 ods graphics off;
******************************************** 
 * kernel density estimate

  %Gpdf(figname= plot17_1.pdf, cex= 1)  

 proc sort data=pubh.brain_glucose;
   by diabetic;
 proc kde data=pubh.brain_glucose bwm=.6 out=g;
   by diabetic;
   var brain_glucose;

 proc gplot data=g;
  plot density * brain_glucose = diabetic;
  symbol1 value=none i=join line=2;
  symbol2 value=none i=join line=1;

  */
 ******************************************** 
  *  From SAS Lifetest help, example 1;

   symbol1 c=blue; symbol2 c=orange; symbol3 c=green;
   symbol4 c=red; symbol5 c=cyan; symbol6 c=black;
   title 'VA Lung Cancer Data';
   proc lifetest data=VALung plots=(s,ls,lls) outtest=Test maxtime=600;
      time SurvTime*Censor(1);
      id Therapy;
      strata Cell;
      test Age Prior DiagTime Kps Treatment;
      ods select Quartiles HomStats HomTests;

   run;

 ************************************************;
 * Simulation to check power for stratified logistic regression;

 data sim;
   do rep=1 to 500;
     do j=1 to 10; strata=1; trt="2Pbo"; y = (RANUNI(-1) < .08); output; end;
     do j=1 to 20; strata=1; trt="1Sit"; y = (RANUNI(-1) < .30); output; end;
     do j=1 to 20; strata=2; trt="2Pbo"; y = (RANUNI(-1) < .480); output; end;
     do j=1 to 40; strata=2; trt="1Sit"; y = (RANUNI(-1) < .800); output; end;
     end;
   keep rep trt strata y;

  
 proc logistic descending  data=sim;
   class trt strata;
   model y = strata trt / ctable pprob=.5;
   by rep;
   ods select OddsRatios;
   ods  output  OddsRatios=OR;
      
 data sig;
  set OR;
  if Effect="trt    1Sit vs 2Pbo";
  detect = (LowerCL > 1);

 proc freq data=sig;
   tables detect;

 

 Poisson rates from Lachin, formula (8.20)

  proc means n sum data=years1;
   var number_late_rej  loss_risk ;
   by cprop_6m_cat;
   output out=f sum= s_rej  s_loss_risk;

  data rates;
    set f;
    scale = 100.0;
    d = s_rej;
    m = s_loss_risk/scale;
    rate_rej = d/m;
    rate_rej_SE = sqrt(d)/m;

 **************  AUC code ***********************;
 data  hunger_AUC   missing;
  set VAS_hunger;
  array m[9] hunger1-hunger9;
  array t[9] time1-time9;
   time1=0; time2=0.25; time3=0.5; time4=0.75; time5=1; 
   time6=1.5; time7=2; time8=3; time9=4;
  AUC = 0;
  do j=1 to 8;
     next_trapezoid = 0.5 * (t[j+1] - t[j])*(m[j] + m[j+1]);
     if (next_trapezoid = .) then do;  
         output  missing;
         GOTO Duluth;  
         end;
     AUC= sum(AUC,next_trapezoid);
     end;
  output hunger_AUC;
  Duluth: 


 /*
 * Correct integer age (from SAS Dinosaur page 170);

 Age=intck('year',BirthDate,ActualDate)
     /* January 1 occurrences between the dates */
    -(put(BirthDate,mmddyy4.) gt put(ActualDate,mmddyy4.))
       /* Birth date is later than Actual date */
    +(put(BirthDate,mmddyy4.)||put(ActualDate,mmddyy4.)||
      put(ActualDate+1,mmddyy4.)='022902280301')
    /* Feb 29 birth date increments on Feb 28 if the following day is March 1 */;

 ***************************************
 Goptions reset=all vsize=5in ftext=simplex device=pdf lfactor=.5
   noborder gsfname=graphout gsfmode= replace;
   filename graphout "C:\Documents and Settings\Administrator\Desktop\xx.pdf";
   symbol1 interpol=std1mj value=none l=1 width=2;
   symbol2 interpol=std1mj value=none l=2 width=2;
 proc gplot data=mayod327.hamd;
  plot response*factor1 = factor2 ; *vaxis=40 to 75 by 5;

  %points(pdf_gif=gif ,figname=gm_muffin1.gif , cex=1 )
 symbol1 value=circle h=1 c=black;
 symbol2 value=dot h=1 c=green;
 symbol3 value=dot h=1 c=blue;
 symbol4 value=dot h=1 c=orange; 
 symbol5 value=dot h=1 c=red; 
 proc Gplot data=willis.AUC;
  plot q4_AUC * pleasant15a = trt/ vref=0;

  %Gpdf(figname= plot_nh1.pdf, height=4, cex= 1) 
 symbol1 value=dot;  symbol2 value=square;
 symbol3 value=triangle;  symbol4 value=circle;

*******************
 %out_plot(pdf_gif= ,figname= ,height=);
 symbol1 l=1 width=3;
 symbol2 l=3 width=3;


  symbol1 l=1 c=blue;
 symbol2 l=3 c=blue;
 symbol3 l=1 width=3 c=red;
 symbol4 l=3 width=3 c=red;
  
 proc lifetest data=pubh.stomach_cancer notable plots=(s) censoredsymbol="|";*(s,ls,lls)  maxtime;
   time days*death(0);
   strata group / test=(logrank tarone wilcoxon);;
*******************
 * jittered binary data plot;
 data plot;
   set mayod327.age_bmi;
   obese_jitter = obese + 0.1*(ranuni(-1) - .5);

 %Gpdf(figname=plot22_4.pdf, height=5, cex= 1)
   symbol1 interpol=sm50s  line=1 width=4;
   symbol2 interpol=sm50s  line=2 width=4;

 proc gplot data=plot;
*   plot obese_jitter * age / haxis=15 to 65 by 5;
   plot obese_jitter * age = female / haxis=15 to 65 by 5;

******************************************** 
  scatter plot matrix;
proc insight data=hsb2;
  scatter   list of Y var * list of X var;