*****************************************************;
* Examples of logistic regression (and stepwise logistic regression):
*
* Suppose that we have Bernoulli trials (that is, observations that
*   yield Yes/No results, or values Y=1 or Y=0). The Bernoulli trials
*   can depend on a covariate or set of covariates, such as age,
*   health, income, etc. 
*
* For definiteness, suppose that Y=1 means Yes (or success) and Y=0
*   means No (or failure). These are the other possible values of Y.
*   Using data about the outcomes of Bernoulli trials with values of
*   their covariates X, we want we want to estimate the probability
*   of success as a function of X:
*
*      p = P(Y=1|X) = p(X) = a function of a+bX
*
* Assuming a linear relationship
*
*      p(X) = P(Y=1|X) = a + b X
*
*   would not work very well, since  0 < a + bX < 1  would imply
*   awkward conditions on a, b, and X. The function
*
*      p(X) = P(Y=1|X) = exp(a + bX)
*
*   is always positive, but can be greater than one. A better choice
*   is
*
*      Prob(Y=1|X) = p(X) = exp(a+bX)/(1 + exp(a+bX))	  (*)
*
*   which maps the line `a+bX' into values between 0 and 1.
*
* The relation (*) is equivalent to p(X)/(1-p(X)) = exp(a+bX), or
*
*      log (p(X)/(1-p(X))) = a + bX   (**)
*
*   Since
*
*      log (p(X)/(1-p(X))) = log (P(Y=1|X)/P(Y=0|X))
*
*   the function log(p/(1-p)) is also called the ``log odds ratio''.
*
* Note that (*) and (**) imply that X affects P(Y=1|X) monotonically.
*   That is, if X increases, then P(Y=1|X) increases if b>0 and
*   P(Y=1|X) decreases if b<0.
*
* While (*) or (**) resembles a classical linear regression
*
*      Y = a + bX + Error,   Error is N(0,sigma^2)  (***)
*
*   it is somewhat different: For Bernoulli (Yes/No) data , the 
*   left-hand side of (***) always has Y=0 or Y=1, and there is no
*   reasonable way to model the error distribution. The relation (**)
*   is more analogous to writing a normal regression as
*
*      E(Y|X) = a + bX	   (****)
*
*   A second difference is that, in (**), we need a LINK FUNCTION
*   that maps the range of the parameter p (0<p<1) to the range of
*   the linear function  a+bX  in (**) and (****), namely the entire
*   real line.
*
* The basic idea of the regression (*) or (**), and indeed on the
*   classical linear regression (***), is that we can use a set of
*   pairs of outcomes and covariates (Y,X) to estimate a linear
*   function of X that estimates the parameters of the distribution
*   of Y. In this way, we can estimate the distribution of Y as
*   a function of the covariates X for new observations of Y.
*
* The function on the left-hand side of (**) is also called the LOGIT
*   function (logit(p)=log(p/(1-p))). Thus we can write (**) as
*
*      logit[p(X)] = a + bX
*
*   The name logit() comes from the fact that its inverse function
*
*      Logistic(y) = exp(y)/(1-exp(y))
*
*   is the LOGISTIC FUNCTION of y. Thus (*) can also be written
*
*      p(X) = Prob(Y=1|X) = Logistic(a+bX).
*
*   For this reason, the model (**) is called the LOGISTIC REGRESSION
*   of the Bernoulli variable Y on the covariates X.
*
* The most commonly used link functions for Bernoulli outcomes, along
*   with their inverse functions, are
*
*     logit(p) = log(p/(1-p))       Logistic(y) = exp(y)/(1+exp(y))
*
*     probit(p) = Phi^{-1}(p)       Phi(y) = Pr[N(0,1)<=y]
*
*     cloglog(p) = log(-log(1-p))   Gompertz(y) = 1-exp(-exp(y))
*
* where N(0,1) means a standard normal random variable. Thus probit(p)
*   is the standard normal quantile for p: That is, probit(p)=y  if
*   and only if  Prob[N(0,1)<=y] = p. The term normit(p) is also used
*   for probit(p). The term LOGISTIC REGRESSION is used even with these
*   alternative link functions.
*
* SAS allows you to enter data for a logistic regression in either of
*   two different ways. The standard way is to enter (for example)
*
*       1  X1
*       1  X2
*       1  X3
*       0  X4
*       0  X5
*       0  X6
*
*   where 1 or 0 in the first column denotes success or failure and
*   the second column are the values of a covariate.
*
* If there are many ties in covariate values, for example if X is
*   always a small integer, then you can also enter the data as
*
*       3   5  X1
*       7  11  X2
*       4  22  X3
*
*   This syntax says (for example) that there were 5 events observed
*   with covariate X1, of which 3 were successes (Y=1) and 2 were
*   failures (Y=0). Similarly, there were 11 events observed with X=X2
*   of which 7 were successes (Y=1), as so forth.
*
* With data in either format, SAS's `proc logistic' uses maximum
*   likelihood methods to estimate the parameters a and b in (*) or
*   (**), where b may be vector valued. The output is similar to that
*   for AFT or linear regressions.
*
* S. Sawyer - Washington University in St.Louis - November 2, 2005
*****************************************************;


title 'EXAMPLES OF LOGISTIC REGRESSION -- YOUR NAME';
options ls=75 ps=60 pageno=1 nocenter;

*****************************************************;
* For each of 10 different values of a covariate, read the results of
*   a number of Bernoulli trials for that covariate. The table is a
*   summary of 133 Bernoulli trials for 10 different covariate values.
* For the logistic regression model, does the covariate have a
*   significant effect on the probability of success?
* The numbers below suggest that increasing the value of X increases
*   the numbers of successes per trial, but this may not be significant.
*****************************************************;

data log1;
     input rr nn xx;
datalines;
    3    10     1
    4     9     2
    5    11     3
    4    15     4
    7    17     5
   11    20     6
    5     9     7
    8    13     8
    7    12     9
   12    17    10
  run;


proc logistic;
     title1 'THE FIRST EXAMPLE: SUCCESSES/FAILURES WITH ONE COVARIATE';
     model rr/nn = xx;
     run;


*****************************************************;
* The following has the results of 32 Bernoulli trials with 12
*   success and 20 failures, along with the values of 4 covariates.
*   The model is
*
*   logit[Prob(Y=1|X)] = a + b_1 X_1 + b_2 X_2 + b_3 X_3 + b_4 X_4
*
*****************************************************;

title1 'SECOND EXAMPLE: OBSERVATIONS WITH FOUR COVARIATES';

data log2;
     input obsno yy x1 x2 x3 x4;
     * For plotting and charting;
     if yy=1 then ysymbol='Y';
     else if yy=0 then ysymbol='N';
     else ysymbol='?';
datalines;
   1     1      48      94      38      93
   2     1      53      79      58      56
   3     1      42     123      51      66
   4     1      70      96      41      55
   5     1      57     102      51      69
   6     1      69     112      36      53
   7     1      41      79      45      67
   8     1      30      86      40      67
   9     1      64     117      40      60
  10     1      54     114      61      88
  11     1      49      79      36      58
  12     1      46      95      45      76
  13     0      49      57      39      65
  14     0      28      74      40      67
  15     0      44      49      29      50
  16     0      43      70      62     112
  17     0      52      97      60      86
  18     0      57      71      42     104
  19     0      39      83      52      86
  20     0      63     100      41      90
  21     0      65     101      49      67
  22     0      48      88      51      86
  23     0      48      74      42      68
  24     0      40      62      54      73
  25     0      50      73      60      89
  26     0      45      61      49      64
  27     0      46      83      33      99
  28     0      47      77      18      73
  29     0      50      88      35      95
  30     0      56      91      40      69
  31     0      38      72      57      97
  32     0      60     105      42      88
     run;



*****************************************************;
* Let's visualize the data by charting the covariate values for
*   the two groups `successes' and `failures'.
*****************************************************;

proc chart;
     title2 'Successes/Failures for different covariate values';
     vbar x1-x4 / subgroup = ysymbol;
     run;


*****************************************************;
* By default, if the dependent value is entered as Y-1 or Y=0, SAS
*   models the probability of the LOWER value of Y. That is, SAS
*   SAS regresses Prob(Y=0|X) on the covariates and not Prob(Y=1|X),
*   as if SAS considers 0 to be a success and 1 to be a failure. 
*
* To reverse this and have SAS model Prob(Y=1) in terms of the
*   covariates (with Y=0 as the baseline) and not vice versa, use
*   the option `descending'.
*
* Since Prob(Y=1) = 1-Prob(Y=0) and logit(p) = -logit(1-p), replacing
*   p by 1-p just changes the sign of logit(p), so that modeling 
*   Prob(Y=0) instead of Prob(Y=1) just changes the signs of the
*   estimated coefficients. However, it sometimes can be difficult to
*   tell whether a particular covariate makes things better or worse
*   (that is, makes Y=1 more or less likely), so that it can be
*   important to keep track of the signs in the output.
* 
* The reason for this odd convention is that SAS also allows Y to be an
*   ordinal variable, for example Y=0,1,2,3,4 for different levels of
*   success. In this case SAS uses the model (for example with two
*   covariates)
*
*     logit[Prob(Y <= y)] = A(y) + b_1 X_1 + b_2 X_2
*
*   so that P(Y<=y) has the same slopes with respect to the covariates
*   but different intercepts. This makes it more natural for SAS to
*   model Prob(Y=0) than Prob(Y=1) when Y is a binary variable.
*
* In the following analysis, only X2 and X4 have significant effects
*   on Prob(Y=1|X). The effect of X2 is positive (that is, larger values
*   of X2 increases Prob(Y=1)) and the effect of X4 is negative. This
*   is quantified by the `odds-ratios' field in the output, which
*   means `odds-ratio per unit of increased covariate'.
*
*****************************************************;


proc logistic descending;
     title2 'LOGISTIC REGRESSION FOR FULL MODEL';
     model yy = x1 x2 x3 x4;
     run;


*****************************************************;
* Let's see if SAS can find the best covariate subset on its own.
*   Whether SAS starts with no covariates and adds them out one
*   at a time, or else starts with all covariates and throws them
*   out one at a time, it ends up with the two covariates x2 x4.
*****************************************************;

proc logistic descending;
     title2 'STEPWISE LOGISTIC REGRESSION';
     model yy = x1 x2 x3 x4 / selection=stepwise;
     run;

proc logistic descending;
     title2 'BACKWARDS STEPWISE LOGISTIC REGRESSION';
     model yy = x1 x2 x3 x4 / selection=backwards;
     run;