/* $Id: simple6.sas,v 1.5 2009/06/23 19:04:52 kkleinma Exp $ */


/*
These examples use the HELP data sets, available from http://www.math.smith.edu/sasr/datasets.php.  
The code provided in these code files reads the data sets in the various formats provided there.  
For smooth running, we suggest downloading all of the files available.  For Windows based systems, 
the files should be stored in c:\data.  The code in the book was run on a windows system in this 
configuration.  For other operating systems, you will need to modify the code below to reflect 
different directory structures and references.  
*/

proc power;
   twosamplemeans groupmeans=(0 0.5) stddev=1 power=0.9 ntotal=.;
run;

proc power;
   twosamplemeans groupmeans=(0 0.5) stddev=1 power=. ntotal=200;
run;

proc power;
   twosamplefreq test=pchi ntotal=. groupproportions=(.1 .2) power=0.9;
run;

proc power;
   twosamplefreq test=pchi ntotal=200 groupproportions=(.1 .2) power=.;
run;

data simpower1;
   effect = 0.35;  /* effect size */
   corr = 0.4;     /* desired correlation */   
   covar = (corr)/(1 - corr);  /* implied covariance given variance = 1*/   
   numsim = 1000;  /* number of datasets to simulate */   
   numfams = 100;  /* number of families in each dataset */   
   numkids = 3;    /* each family */   
   do simnum = 1 to numsim;   /* make a new dataset for each simnum */      
      do famid = 1 to numfams; /* make numfams families in each dataset */         
         inducecorr = normal(42)* sqrt(covar);
                        /* this is the mechanism to achieve the desired	 
                           correlation between kids within family */         
         do kidnum = 1 to numkids;   /* generate each kid */            
            exposed = ((kidnum eq 1) or (famid le numfams/2)) ; 
                        /* assign kid to be exposed */            
            x = (exposed * effect) + (inducecorr + normal(0))/sqrt(1 + covar);            
            output;         
         end;      
      end;   
   end;
run;  

*** This generates too much output: use ODS even in this simplified code.;
ods select none;

options nonotes;  *** Don't clutter the log with reports on each model; 
ods output solutionf=simres;
proc mixed data=simpower1 order=data;
by simnum;
   class exposed famid;   
   model x = exposed / solution;   
   random int / subject=famid;
run;
ods select all;
options notes;

data powerout;
set simres;   
   where exposed eq 1;
   reject=(probt lt 0.05);
run;

proc freq data=powerout;
   tables reject / binomial (level='1');
run;

data sim;   
   sigbsq=4; beta0=-2; beta1=1.5; beta2=0.5; beta3=-1; n=1500;
   do i = 1 to n;
      x1 = (i lt (n+1)/2);
      randint = normal(0) * sqrt(sigbsq);
      do x2 = 1 to 3 by 1;
          x3 = uniform(0);
          linpred = beta0 + beta1*x1 + beta2*x2 + beta3*x3 + randint;
          expit = exp(linpred)/(1 + exp(linpred));
          y = (uniform(0) lt expit);
          output;
      end;
   end;
run;

proc nlmixed data=sim qpoints=50;
   parms b0=1 b1=1 b2=1 b3=1;
   eta = b0 + b1*x1 + b2*x2 + b3*x3 + bi1;
   mu = exp(eta)/(1 + exp(eta));
   model y ~ binary(mu);
   random bi1 ~ normal(0, g11) subject=i;
   predict mu out=predmean;
run;

proc glimmix data=sim order=data;
   nloptions maxiter=100 technique=dbldog;
   model y = x1 x2 x3 / solution dist=bin;
   random int / subject=i;
run;

data test;
   p1=.15; p2=.25;
   corr=0.4;
   p1p2=corr*sqrt(p1*(1-p1)*p2*(1-p2)) + p1*p2;
   do i = 1 to 10000;
      cat=rand('TABLE', 1-p1-p2+p1p2, p1-p1p2, p2-p1p2);
      y1=0;
      y2=0;
      if cat=2 then y1=1;
      else if cat=3 then y2=1;
      else if cat=4 then do;
         y1=1;
         y2=1;
      end;
   output;
   end;
run;


filename census1 url "http://www.math.smith.edu/sasr/datasets/co25_d00.dat";

data pcts cents;
infile census1;
retain cntyid;
input @1 endind $3. @;  /* the trailing '@' means to hold onto this line */
   if endind ne 'END' then do;   
   input @7 neglat $1. @;  /* if this line does not say 'END', then
                              check to see if the 7th character is '-' */   
   if neglat eq '-' then do;   /* if so, it has a boundary point */     
      input @7 x y;

/* NOTE: the following statement idavertently omitted from the text */
 
	  const = 1; 

/* See comment above */

      output pcts;             /* write out to boundary dataset */   
      end;   
   else if neglat ne '-' then do;  /* if not, it must be the centroid */      
      input @9 cntyid 2. x y ;  	  
      output cents;                /* write it to the centroid dataset */   
      end;
   end;
run;

filename census2 url "http://www.math.smith.edu/sasr/datasets/co25_d00a.dat"; 
data cntynames;
infile census2 DSD;   
format cntyname $17. ;   
input cntyid 2. cntyname $;
run;

proc sort data=cntynames; by cntyid; run;
proc sort data=cents; by cntyid; run;

data nameloc;   
   length function style color $ 8 position $ 1 text $ 20;   
   retain xsys ysys "2" hsys "3" when "a";
   merge cntynames cents;
   by cntyid;   
   function="label"; 
   style="swiss"; 
   text=cntyname; 
   color="black";   
   size=3; position="5";    
   output;
run;

ods pdf file="map_plot.pdf";

pattern1 value=empty;
proc gmap map=pcts data=pcts;
   choro const / nolegend coutline=black annotate=nameloc;
   id cntyid;
run; quit;

ods pdf close;

filename myurl url "http://www.math.smith.edu/sasr/datasets/helpmiss.csv"    lrecl=704;

proc import replace datafile=myurl out=help dbms=dlm;
   delimiter=',';
   getnames=yes;
run;

ods select misspattern;
proc mi data=help nimpute=0;
   var homeless female i1 sexrisk indtot mcs pcs;
run;

ods select all;

/* Note: error in the code below in text: should read data=help  */

proc mi data=help nimpute=20 out=helpmi20 noprint;
   mcmc chain=multiple;
   var homeless female i1 sexrisk indtot mcs pcs;
run;

ods select none;
ods output parameterestimates=helpmipe covb=helpmicovb;
proc logistic data=helpmi20 descending;by _imputation_;
   model homeless=female i1 sexrisk indtot / covb;
run;
ods select all;

proc mianalyze parms = helpmipe covb=helpmicovb;
modeleffects intercept female i1 sexrisk indtot;
run;

proc import datafile='c:/book/help.csv' out=help dbms=dlm;
   delimiter=',';
   getnames=yes;
run;

proc genmod data=help;
   class substance;
   model i1 = female substance age / dist=poisson;
*   bayes;
run;

*** Code chunk number 1 ***;
libname k "c:/book";

data help;
set k.help;
run;

*** Code chunk number 2 ***;
proc corr data=help alpha nomiss;
   var f1a -- f1t;
run;

*** Code chunk number 3 ***;
*** Code chunk number 4 ***;
data helpcc;
set help;
   n = n(of f1a--f1t);
   if n=20;
run;

*** Code chunk number 5 ***;
proc factor data=helpcc nfactors=3 method=ml rotate=varimax;
   var f1a--f1t;
run;

*** Code chunk number 6 ***;
*** Code chunk number 7 ***;
proc discrim data=help out=ldaout;
   class homeless;
   var age cesd mcs pcs;
run;

*** Code chunk number 10 ***;
axis1 label=("Prob(homeless eq 1)");
proc univariate data=ldaout;
   class homeless;
   var _1;
   histogram _1 / nmidpoints=20 haxis=axis1;
run;

*** Code chunk number 11 ***;
proc varclus data=help outtree=treedisp centroid;
   var mcs pcs cesd i1 sexrisk;
run;

proc tree data=treedisp nclusters=5;
   height _varexp_; 
run;