While the goodness-of-fit statistics can tell us how well a particular model fits the data, they don't tell us much about why a model may fit poorly. Do the model assumptions hold? To assess the lack of fit we need to look at regression diagnostics.

Let us begin with a bit of a review of Linear Regression diagnostics. The standard linear regression model is given by:

\(y_i \sim N(\mu_i,\sigma^2)\)

\(\mu_i =x_i^T \beta\)

The two crucial features of this model are

1. the assumed mean structure, \(\mu_i =x_i^T \beta\) , and
2. the assumed constant variance \(\sigma^2\) (homoscedasticity).

The most common diagnostic tool is the residuals, the difference between the estimated and observed values of the dependent variable. There are other useful regression diagnostics, e.g. measures of leverage and influence, but for now, our focus will be on the estimated residuals.

The most common way to check these assumptions is to fit the model and then plot the residuals versus the fitted values \(\hat{y}_i=x_i^T \hat{\beta}\) .

• If the model assumptions are correct, the residuals should fall within an area representing a horizontal band, like this,

• If the residuals have been standardized in some fashion (i.e., scaled by an estimate of \(\sigma\)), then we would expect most of them to have values within \(\pm2\) or \(\pm3\); residuals outside of this range are potential outliers.
• If the plot reveals some kind of curvature—for example, like this,

  it suggests a failure of the mean model; the true relationship between \(\mu_i\), and the covariates might not be linear.

• If the variability in the residuals is not constant as we move from left to right—for example, if the plot looks like this,

  or like this,

  then the variance \(V(Y_i)\) is not constant but changes as the mean \(\mu_i\) changes.

The logistic regression model says that the mean of \(Y_i\) is

\(\mu_i=n_i \pi_i\)

where

\(\log\left(\dfrac{\pi_i}{1-\pi_i}\right)=x_i^T \beta\)

and that the variance of \(Y_i\) is

\(V(Y_i)=n_i \pi_i(1-\pi_i)\).

After fitting the model, we can calculate the Pearson residuals:

\(r_i=\dfrac{y_i-\hat{\mu}_i}{\sqrt{\hat{V}(Y_i)}}=\dfrac{y_i-n_i\hat{\pi}_i}{\sqrt{n_i\hat{\pi}_i(1-\hat{\pi}_i)}}\)

If the \(n_i\)s are relatively large, these act like standardized residuals in linear regression. To see what's happening, we can plot them against the linear predictors:

\(\hat{\eta}_i=\log\left(\dfrac{\pi_i}{1-\pi_i}\right)=x_i^T \hat{\beta}_i\)

which are the estimated log-odds of success, for cases \(i = 1, \ldots , N\).

• If the fitted logistic regression model was true, we would expect to see a horizontal band with most of the residuals falling within \(\pm 3\):

• If the \(n_i\)s are small, then the plot might not have such a nice pattern, even if the model is true. In the extreme case of ungrouped data (with all \(n_i\)s equal to 1), this plot becomes uninformative. From now on, we will suppose that the \(n_i\)s are not too small so that the plot is at least somewhat meaningful.

• If outliers are present—that is, if a few residuals or even one residual is substantially larger than \(\pm 3\), then \(X^2\) and \(G^2\) may be much larger than the degrees of freedom. In that situation, the lack of fit can be attributed to outliers, and the large residuals will be easy to find in the plot.

Curvature

• If the plot of Pearson residuals versus the linear predictors reveals curvature—for example, like this,

then it suggests that the mean model

\(\log\left(\dfrac{\pi_i}{1-\pi_i}\right)=x_i^T \beta\) (2)

has been misspecified in some fashion. That is, it could be that one or more important covariates do not influence the log-odds of success in a linear fashion. For example, it could be that a covariate \(x\) ought to be replaced by a transformation such as \(\sqrt{x}\), \(\log x\), \(x^2\), etc.

To see whether individual covariates ought to enter in the logit model in a non-linear fashion, we could plot the empirical logits

\(\log\left(\dfrac{y_i+1/2}{n_i-y_i+1/2}\right)\) (3)

versus each covariate in the model.

• If one of the plots looks substantially non-linear, we might want to transform that particular covariate.

• If many of them are nonlinear, however, it may suggest that the link function has been mis-specified—i.e., that the left-hand side of (2) should not involve a logit transformation but some other function.

Changing the link function will change the interpretation of the coefficients entirely; the \(\beta\)s will no longer be log-odds ratios. But, depending on what the link function is, they might still have a nice interpretation. For example, in a model with a log link, \(\log\pi_i=x_i^T \beta\), an exponentiated coefficient \(\exp(\beta_j)\) becomes a relative risk.

Test for correctness of the link function

Hinkley (1985) suggests a nice, easy test to see whether the link function is plausible:

\(\hat{\eta}_i=x_i^T \hat{\beta}\)

• Fit the model and save the estimated linear predictors
• Add \(\eta^2_i\) to the model as a new predictor and see if it's significant.

A significant result indicates that the link function is misspecified. A nice feature of this test is that it applies even to ungrouped data (\(n_i\)s equal to one), for which residual plots are uninformative.

Non-constant Variance

Suppose that the residual plot shows non-constant variance as we move from left to right:

Another way to detect non-constancy of variance is to plot the absolute values of the residuals versus the linear predictors and look for a non-zero slope:

Non-constant variance in the Pearson residuals means that the assumed form of the variance function,

\(V(y_i)=n_i \pi_i (1-\pi_i)\)

is wrong and cannot be corrected by simply introducing a scale factor for overdispersion. Overdispersion and changes to the variance function will be discussed later.

Regression diagnostics can also tell us how influential each observation is to the fit of the logistic regression model. We can evaluate the numerical values of these statistics and/or consider their graphical representation (like residual plots in linear regression).

Some measures of influence:

• Pearson and Deviance Residuals identify observations not well explained by the model.
• Hat Matrix Diagonal detects extremely large points in the design space. These are often labeled as "leverage" or "\(h_i\)" and are related to standardized residuals. A general rule says that if \(h_i > 2p/n\) or \(> 3p/n\), the points is influential. Here "p" is the number of parameters in the model and "n" is the number of observations. For the donner data example with the main effects model, \(2(3)/45=0.13\). There are two points possibly influential.
• The regression diagnostics work is very similar to linear regression.

      „ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ†
H 0.2 ˆ                                                                                                                            ˆ
a     ‚                                                                                                                            ‚
t     ‚                                                                                *                                           ‚
    H ‚                                                      *                                                                     ‚
D     ‚                              *     *                                                                                       ‚
i     ‚                                                                                                                            ‚
a 0.1 ˆ                    *                                                         *                                             ˆ
g     ‚                                      *       * *   *     *             * *       *     *                 *                 ‚
o     ‚                                * *                     *   *   *   *                 *                                     ‚
n     ‚                  *   *     *               *                 *   *                 *       * * *   * * *                   ‚
a     ‚                        * *             * *       *                   *     *             *       *                         ‚
l     ‚                                                                                                                            ‚
  0.0 ˆ                                                                                                                            ˆ
      ŠƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒŒ
                       0         5        10        15        20        25        30        35        40        45                  
                                                                                                                                    
                                                           Case Number   INDEX                                                      

• DFBETA assesses the effect of an individual observation on the estimated parameter of the fitted model. It’s calculated for each observation for each parameter estimate. It's a difference in the parameter estimated due to deleting the corresponding observation. For example, for the main effects model there are 3 such plots (one for intercept, one for age, and one for sex). Here is the one for AGE, showing that observation #35 is potentially influential.

          „ƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ†
      0.5 ˆ                                                                                                                        ˆ
a         ‚                                                                                    *                                   ‚
g         ‚                                                                            *                                           ‚
e DFBETA1 ‚                  * *                                                                                                   ‚
          ‚                                                                                                                        ‚
D         ‚                *                   *     *                 *       * *           *   *   * * *       *                 ‚
f     0.0 ˆ                      * *             * *       *                         *             *       *                       ˆ
B         ‚                          *   * *           * *       * *     * * *     *       *                 * *   *               ‚
e         ‚                            *     *               * *     *                   *                                         ‚
t         ‚                                                                                                                        ‚
a         ‚                                                                                                                        ‚
          ‚                                                                                                                        ‚
     -0.5 ˆ                                                                                                                        ˆ
          ŠƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒŒ
                         0         5        10        15        20        25        30        35        40        45                
                                                                                                                                    
                                                             Case Number   INDEX                                                    

• C and CBAR provide scalar measures of the influence of the individual observations on the parameter estimates, similar to the Cook distance in linear regression. The same rules in linear regression hold here too. For example, observation #35 again seems influential (see the rest of donner.lst).
• DIFDEV and DIFCHISQ detect observations that heavily add to the lack of fit, i.e. there is a large disagreement between the fitted and observed values. DIFDEV measures the change in the deviance (\(G^2\)) as an individual observation is deleted. DIFCHISQ measures the change in \(X^2\).

In SAS, under PROC LOGISTIC, INFLUENCE option and IPLOTS option will output these diagnostics. In R see, influence() function.

As pointed out before, the standardized residual values are considered to be indicative of lack of fit if their absolute value exceeds 3 (some sources are more conservative and take 2 as the threshold value).

The Index plot, in which the residuals are plotted against the indices of the observations, can help us determine outliers, and/or systematic patterns in variation and thus assess the lack of fit. In our example, all Pearson and Deviance residuals fall within \(-2\) and \(+2\); thus, there does not seem to be any outliers or particular heavy influential points; see donner.html for the IPLOTS.

For our example, here is the IPLOT of the deviance residuals vs. the case number.

         „ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ†
D      2 ˆ                     *                                                               *                                   ˆ
e        ‚                                                                     *                                                   ‚
v        ‚                                       *         *               *         *                         *                   ‚
i RESDEV ‚                   *                                       *       *   *               *                                 ‚
a        ‚                                             * *         *               *       *                       *               ‚
n        ‚                                                   *                                                                     ‚
c      0 ˆ                                                                                                                         ˆ
e        ‚                               * *                     *       *                                                         ‚
         ‚                           *                                                                       *                     ‚
R        ‚                 *     * *   *     *     * *         *       *                 *   *     * * * * *     *                 ‚
e        ‚                                                                             *                                           ‚
s        ‚                                     *                                                                                   ‚
i     -2 ˆ                                                                                                                         ˆ
d        ŠƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒŒ
u                        0         5        10        15        20        25        30        35        40        45                
a                                                                                                                                   
                                                             Case Number   INDEX                                                    

Overdispersion is an important concept in the analysis of discrete data. Many times data admit more variability than expected under the assumed distribution. The extra variability not predicted by the generalized linear model random component reflects overdispersion. Overdispersion occurs because the mean and variance components of a GLM are related and depend on the same parameter that is being predicted through the predictor set.

Overdispersion is not an issue in ordinary linear regression. In a linear regression model

\(y_i \sim N(x_i^T \beta, \sigma^2)\)

the variance \(σ^2\) is estimated independently of the mean function \(x_i^T \beta\). With discrete response variables, however, the possibility for overdispersion exists because the commonly used distributions specify particular relationships between the variance and the mean; we will see the same holds for Poisson.

For the binomial response, if \(Y_i\sim Bin(n_i, \pi_i)\), the mean is \(\mu_i=n_i\pi_i\), and the variance is \(\mu_i(n_i - \mu_i) / n_i\).

• Overdispersion means that the variance of the response \(Y_i\) is greater than what's assumed by the model.
• Underdispersion is also theoretically possible but rare in practice. More often than not, if the model's variance doesn't match what's observed in the response, it's because the latter is greater.

In the context of logistic regression, overdispersion occurs when the discrepancies between the observed responses \(y_i\) and their predicted values \(\hat{\mu}_i=n_i \hat{\pi}_i\) are larger than what the binomial model would predict. Overdispersion arises when the \(n_i\) Bernoulli trials that are summarized in a line of the dataset are

• not identically distributed (i.e., the success probabilities vary from one trial to the next), or
• not independent (i.e., the outcome of one trial influences the outcomes of other trials).

In practice, it is impossible to distinguish non-identically distributed trials from non-independence; the two phenomena are intertwined.

The problem of overdispersion may also be confounded with the problem of omitted covariates. If we have included all the available covariates related to \(Y_i\) in our model and it still does not fit, it could be because our regression function \(x_i^T \beta\) is incomplete. Or it could be due to overdispersion. Unless we collect more data, we cannot do anything about omitted covariates. But we can adjust for overdispersion.

The most popular method for adjusting for overdispersion comes from the theory of quasi-likelihood. Quasilikelihood has come to play a very important role in modern statistics. It is the foundation of many methods that are thought to be "robust" (e.g. Generalized Estimating Equations (GEE) for longitudinal data) because they do not require the specification of a full parametric model.

In the quasilikelihood approach, we must first specify the "mean function" which determines how \(\mu_i=E(Y_i)\) is related to the covariates. In the context of logistic regression, the mean function is

\(\mu_i=n_i \exp(x_i^T \beta)\),

which implies

\(\log\left(\dfrac{\pi_i}{1-\pi_i}\right)=x_i^T \beta\)

Then we must specify the "variance function," which determines the relationship between the variance of the response variable and its mean. For a binomial model, the variance function is \(\mu_i(n_i-\mu_i)/n_i\). But to account for overdispersion, we will include another factor \(\sigma^2\) called the "scale parameter," so that

\(V(Y_i)=\sigma^2 \mu_i (n_i-\mu_i)/n_i\).

• If \(\sigma^2 \ne 1\) then the model is not binomial; \(\sigma^2 > 1\) corresponds to "overdispersion", and \(\sigma^2 < 1\) corresponds to "underdispersion."

• If \(\sigma^2\) were known, we could obtain a consistent, asymptotically normal and efficient estimate for \(\beta\) by a quasi-scoring procedure, sometimes called "estimating equations." For the variance function shown above, the quasi-scoring procedure reduces to the Fisher scoring algorithm that we mentioned as a way to iteratively find ML estimates.

Note that no matter what \(\sigma^2\) is assumed to be, we get the same estimate for \(\beta\). Therefore, this method for overdispersion does not change the estimate for \(\beta\) at all. However, the estimated covariance for \(\hat{\beta}\) changes from

\(\hat{V}(\hat{\beta})=(x^T W x)^{-1}\)

to

\(\hat{V}(\hat{\beta})=\sigma^2 (x^T W x)^{-1}\)

That is, the estimated standard errors must be multiplied by the factor \(\sigma=\sqrt{\sigma^2}\).

 How do we estimate \(\sigma^2\)?

One recommendation is to use

\(\hat{\sigma}^2=X^2/(N-p)\)

where \(X^2\) is the usual Pearson goodness-of-fit statistic, \(N\) is the number of sample cases (number of rows in the dataset we are modeling), and \(p\) is the number of parameters in the current model (suffering from overdispersion).

• If the model holds, then \(X^2/(N - p)\) is a consistent estimate for \(\sigma^2\) in the asymptotic sequence \(N\rightarrow\infty\) for fixed \(n_i\)s.

• The deviance-based estimate \(G^2/(N - p)\) does not have this consistency property and should not be used.

This is a reasonable way to estimate \(\sigma^2\) if the mean model \(\mu_i=g(x_i^T \beta)\) holds. But if important covariates are omitted, then \(X^2\) tends to grow and the estimate for \(\sigma^2\) can be too large. For this reason, we will estimate \(\sigma^2\) under a maximal model, a model that includes all of the covariates we wish to consider.

The best way to estimate \(\sigma^2\) is to identify a rich model for \(\mu_i\) and designate it to be the most complicated one that we are willing to consider. For example, if we have a large pool of potential covariates, we may take the maximal model to be the model that has every covariate included as a main effect. Or, if we have a smaller number of potential covariates, we decide to include all main effects along with two-way and perhaps even three-way interactions. But we must omit at least a few higher-order interactions, otherwise, we will end up with a model that is saturated.

In an overdispersed model, we must also adjust our test statistics. The statistics \(X^2\) and \(G^2\) are adjusted by dividing them by \(\sigma^2\). That is, tests of nested models are carried out by comparing differences in the scaled Pearson statistic, \(\Delta X^2/\sigma^2\), or the scaled deviance, \(G^2/\sigma^2\) to a chi-square distribution with degrees of freedom equal to the difference in the numbers of parameters for the two models.

If the data are overdispersed — that is, if

\(V(Y_i) \approx \sigma^2 n_i \pi_i (1-\pi_i)\)

for a scale factor \(\sigma^2 > 1\), then the residual plot may still resemble a horizontal band, but many of the residuals will tend to fall outside the \(\pm 3\) limits. In this case, the denominator of the Pearson residual will tend to understate the true variance of the \(Y_i\), making the residuals larger. If the plot looks like a horizontal band but \(X^2\) and \(G^2\) indicate lack of fit, an adjustment for overdispersion might be warranted. A warning about this, however: If the residuals tend to be too large, it doesn't necessarily mean that overdispersion is the cause. Large residuals may also be caused by omitted covariates. If some important covariates are omitted from \(x_i\), then the true \(\mu_i\)s will depart from what your model predicts, causing the numerator of the Pearson residual to be larger than usual. That is, apparent overdispersion could also be an indication that your mean model needs additional covariates. If these additional covariates are not available in the dataset, however, then there's not much we can do about it; we may need to attribute it to overdispersion.

Note, there is no overdispersion for ungrouped data. McCullagh and Nelder (1989) point out that overdispersion is not possible if \(n_i=1\). If \(y_i\) only takes values 0 and 1, then it must be distributed as Bernoulli(\(\pi\)), and its variance must be \(\pi_i(1-\pi_i)\). There is no other distribution with support {0,1}. Therefore, with ungrouped data, we should always assume scale=1 and not try to estimate a scale parameter and adjust for overdispersion.

The usual way to correct for overdispersion in a logit model is to assume that:

where \(\sigma^2\) is a scale parameter. Under this modification, the Fisher-scoring procedure for estimating \(\beta\) does not change, but its estimated covariance matrix becomes \(\sigma^2(x^TWx)^{-1}\)—that is, the usual standard errors are multiplied by the square root of \(\sigma^2\). This will make the confidence intervals wider.

Let's get back to our example and refit the model, making an adjustment for overdispersion.

model y/n= logconc / scale=pearson; 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -15.834      4.923  -3.216 0.001298 ** 
logconc        5.578      1.680   3.319 0.000902 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 4.08)

    Null deviance: 83.631  on 8  degrees of freedom
Residual deviance: 29.346  on 7  degrees of freedom
AIC: 62.886

There are other corrections that we could make. If we were constructing an analysis-of-deviance table, we would want to divide \(G^2\) and \(X^2\) by \(\hat{\sigma}^2\) and use these scaled versions for comparing nested models. Moreover, in reporting residuals, it would be appropriate to modify the Pearson residuals to

\( r_i^\ast=\dfrac{y_i-n_i\hat{\pi}_i}{\sqrt{\hat{\sigma}^2n_i\hat{\pi}_i(1-\hat{\pi}_i)}}\);

that is, we should divide the Pearson residuals (or the deviance residuals, for that matter) by \(\sqrt{\hat{\sigma}^2}\).

A Receiver Operating Characteristic Curve (ROC) is a standard technique for summarizing classifier performance over a range of trade-offs between true positive (TP) and false positive (FP) error rates (Sweets, 1988). ROC curve is a plot of sensitivity (the ability of the model to predict an event correctly) versus 1-specificity for the possible cut-off classification probability values \(\pi_0\).

For logistic regression we can create a \(2\times 2\) classification table of predicted values from your model for the response if \(\hat{y}=0\) or 1 versus the true value of \(y = 0\) or 1. The prediction if \(\hat{y}=1\) depends on some cut-off probability, \(\pi_0\). For example, \(\hat{y}=1\) if \(\hat{\pi}_i>\pi_0\) and \(\hat{y}=0\) if \(\hat{\pi}_i \leq \pi_0\). The most common value for \(\pi_0 = 0.5\). Then \(sensitivity=P(\hat{y}=1|y=1)\) and \(specificity=P(\hat{y}=0|y=0)\).

The ROC curve is more informative than the classification table since it summarizes the predictive power for all possible \(\pi_0\).

The position of the ROC on the graph reflects the accuracy of the diagnostic test. It covers all possible thresholds (cut-off points). The ROC of random guessing lies on the diagonal line. The ROC of a perfect diagnostic technique is a point at the upper left corner of the graph, where the TP proportion is 1.0 and the FP proportion is 0.

The Area Under the Curve (AUC), also referred to as index of accuracy (A), or concordance index, \(c\), in SAS, and it is an accepted traditional performance metric for a ROC curve. The higher the area under the curve the better prediction power the model has. \(c = 0.8 \) can be interpreted to mean that a randomly selected individual from the positive group has a test value larger than that for a randomly chosen individual from the negative group 80 percent of the time.

options nocenter nodate nonumber linesize=72;
data assay; 
input logconc y n; 
cards; 
2.68  10  31 
2.76  17  30 
2.82  12  31 
2.90   7  27 
3.02  23  26 
3.04  22  30 
3.13  29  31 
3.20  29  30 
3.21  23  30 
; 
run; 
                                                                                 
proc logistic data=assay; 
  model y/n= logconc / scale=pearson outroc=roc1; 
  output out=out1 xbeta=xb reschi=reschi; 
  run; 
                                                                                 
axis1 label=('Linear predictor');
axis2 label=('Pearson Residual');
proc gplot data=out1; 
  title 'Residual plot'; 
  plot reschi * xb / haxis=axis1 vaxis=axis2; 
run; 
symbol1 i=join v=none c=blue;
proc gplot data=roc1;
  title 'ROC plot';
  plot  _sensit_*_1mspec_=1 / vaxis=0 to 1 by .1 cframe=ligr ;
run;

                                  #### ROC curve
                                  #### sensitivity vs 1-specificity
                                  lp = result$linear.predictors
                                  p = exp(lp)/(1+exp(lp))
                                  cbind(yes,no,p)
                                  p0 = 0
                                  sens = 1
                                  spec = 0
                                  total = 100
                                  for (i in (1:total)/total)
                                  { 
                                  yy = sum(r*(p>=i))
                                  yn = sum(r*(p=i))
                                  nn = sum(n*(p<i))
                                

In the past two lessons, we've explored ways to fit and evaluate logistic regression models for a binary response. These results generalize what we saw earlier for two and three-way tables to explain associations between two variables while controlling for additional variables and allowing for both categorical and continuous types. One shortcoming, however, would be in accommodating response variables with more than two levels (something more general than "success" or "failure"). For this, we can utilize a multinomial random component; the link and systematic components need not change. The challenge, as we'll see, is in defining and interpreting odds when more than one complementary category is involved.