Advanced Data Analysis, II

36-402  December 19, 1997  Spring, 1996

Logistic Regression—Again

The Challenger Data

You are probably all familiar with the Challenger Space Shuttle accident of January 20, 1986. I will leave out all the interesting substantive details about the accident and focus instead on the data about O-rings that was available at the time of the accident.

Well, I do need a little background, so here goes. The space shuttle uses two booster rockets and each rocket is made of several pieces whose joints are sealed by rubber O-rings. The O-rings are supposed to stop the hot gases that are inside the booster from escaping. Each booster contains 3 primary O-rings. There was data about the O-rings on 23 previous flights of space shuttles. Here is that data.

Temp Fail          Temp Fail
53   2             70   1
57   1             70   1
58   1             72   0
63   1             73   0
66   0             75   0
67   0             75   2
67   0             76   0
67   0             76   0
68   0             78   0
69   0             79   0
70   0             81   0
70   0

On the day of the fated flight, the forecast launch temperature was 31^{∘}F. The O-rings were known to be sensitive to temperature. The previous coldest launch temperature was 53^{∘}. There was discussion about the cold temperature before launch, and the engineers did look at the data, but only the data for the failures. They produced a plot that looked something like this.

> chal <- read.table("chal.dat", head=T)
> attach(chal, pos=1)
> ls()
[1] ".Last.value" "Fail"        "Temp"
> n <- rep(6, length(Fail))
> plot(Temp[Fail >0], Fail[Fail>0])

Based on the erroneous assumption that the flights with no failures did not add any information, they noted that there did not appear to be any relationship between failures and temperature.

Had they produced a plot of the full data, which follows, they might have behaved differently.

> plot(Temp, Fail)

In addition before each flight the O-rings were given a pressure test. The pressure at which that test was conducted increased over time, and it is possible that the higher pressure test could have caused more failures. That data, and a lot more discussion, is available in the paper Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure, by S.R Dalal, E.B. Fowlkes, and B. Hoadley, in V. 84, pp.945–957, of the Journal of the American Statistical Association (JASA). You should be able to read most of that article. It is well written and easy to follow.

How to look at these data

Obviously, we are most interested how probability of failure relates to temperature, and other quantities known before a flight. Given our data, it is not immediately obvious what response variable we should be looking at. For example, should we look at the number of failures of O-rings in a flight, or just the fact that an O-ring failed? These might lead to different results (but fortunately for us they can be looked at using the same analyses).

We have in the past looked at models for count data that were of the form

log m = some combination of factor effects

Although these loglinear models did not require discrete (categorical) predictors, the only examples we saw were of that form. The current data can also be modelled using the standard loglinear model, but an equivalent and more natural approach is to set up the model as a logistic regression.

In this case, allow p(t) to be the probability of failure given the temperature t. If one wanted to be really careful, one might write this as

p(F |t) = probability of failure given temperature was t

but the shorthand notation will probably suffice.

A first approach is to consider that the failures are a binomial process with the probability of failure dependent on t, but with the events themselves independent. Thus using the first form of the data (where we explicitly model the number of O-rings that failed) we are assuming that

• The 6 O-rings on any flight are fail independently (but with a common probability), and

• The 23 flights are independent.

The second assumption is more tenable than the first, but even the first seems quite reasonable. In order to fret about the first assumption, we will also analyze a version of the data that just records whether there was a failure or not.

The Logistic model

Clearly, we can not get away with modelling

p(t) = some arbitrary function of t

unless we are very careful about the arbitrary function. One would quickly find oneself in the position where the models produced negative probabilities, or probabilities greater than 1.

One might think of the model (based on our loglinear experience)

log p(t) = some arbitrary function of t

but it suffers from the same problems.

Lets go back to first principles and think about the data we are trying to model.

We believe (at least I am suggesting) that the data can be modelled as follows. Let n be the number of trials that are done for each run of the experiment (in the Challenger case this is either 6 if each O-ring is considered or 1 if just overall success/failure is considered) and let r be the number of failures in those n trials. Let R be the corresponding random variable. Then we believe that

P (R = r|t) = \binom{n}{r}p(t)^{r}(1 − p(t))^{n−r}
= \binom{n}{r}(\frac{p(t)}{1 − p(t)})^{r}(1 − p(t))^{n}
= Fn of r × Fn of p and r × Fn of only p

Likelihood theory (and common sense) suggests that the only thing of any importance here is the part that depends on both r and p. In other words, if you want to maximize the likelihood only the middle bit counts.

If we go one step further, and consider the log-likelihood

log \binom{n}{r} + r × log (\frac{p(t)}{1 − p(t)}) + n × log (1 − p(t))

then we see that modelling might be easy and natural (and the asymptotic statistics might work) if we form models that are linear in

log (\frac{p(t)}{1 − p(t)})

That is exactly what we will do.

The function

log (\frac{p}{1 − p})

is called the logistic transformation. If

x = log (\frac{p}{1 − p})

then the inverse function is

p = \frac{exp(x)}{1 + exp(x)}

We will (obviously, from the name) model the logit of the probabilities as a function of the predictors.

Back to the Challenger Data

How do we want to model the challenger data? We (with the current data) only have temperature as a predictor, so we are stuck with models of the form,

logit(p(t)) = some function oft

The obvious first step is to plot the data (how can we do that?), then fit some simple models, then think about transformations. In this example we are particularly interested in extrapolations (dangerous though that is!) to temperatures that are less than we have seen.

Plotting the data

This is harder than one might imagine. One (mostly) can not just plot the logits, because much of the data is at the extremes where p = 0 or (1 − p) = 0. Often the best we can do is the plots that we did earlier, perhaps adding some sort of smooth line to the data. (Next example, we will try that). We will be able to do better once we have a model and can consider residuals—but be careful which residuals we think about.

The model

The S-Plus logistic regression command needs to know both the number of successes and the number of failures, so we need to start out by creating such a variable,

> FN <- cbind(Fail, n-Fail)

Now, just use the glm command, with the family=binomial option, to fit the model, viz.

> mymod <- glm( FN ~ Temp, family=binomial)
> mymod <- glm( cbind(Fail, n-Fail) ~ Temp, family=binomial)

Either of the above forms will work, but the first is easier in the long run. Now look at the results of the fit,

> summary(mymod)
Call: glm(formula = FN ~ Temp, family = binomial)
Deviance Residuals:
        Min         1Q     Median          3Q      Max
 -0.9522773 -0.7829968 -0.5411832 -0.04378972 2.651492
Coefficients:
                 Value Std. Error   t value
(Intercept)  5.0848620  3.0477855  1.668379
       Temp -0.1155992  0.0469317 -2.463136
(Dispersion Parameter for Binomial family taken to be 1 )
    Null Deviance: 24.23036 on 22 degrees of freedom
Residual Deviance: 18.08633 on 21 degrees of freedom
Number of Fisher Scoring Iterations: 4
Correlation of Coefficients:
     (Intercept)
Temp -0.9932182

> anova(mymod)
Analysis of Deviance Table
Binomial model
Response: FN
Terms added sequentially (first to last)
     Df Deviance Resid. Df Resid. Dev
NULL                    22   24.23036
Temp  1 6.144035        21   18.08633

The model fits well (although, just like R^{2} is regression, one should be cautious about this, and the temperature coefficient is certainly significant.

At this point we need to

• Get predictions for the low-temperature extrapolations, and

• Evaluate the fit of the model

Extrapolating Lets work on that, in that order (which is really back-to-front, but useful for teaching purposes).

The fitted values (on the probability scale) are obtained via,

> fitted(mymod)
         1         2         3          4          5          6          7
 0.2607849 0.1817866 0.1652202 0.09994062 0.07278457 0.06535816 0.06535816
          8          9         10         11         12         13         14
 0.06535816 0.05864156 0.05257638 0.04710711 0.04710711 0.04710711 0.04710711
         15         16         17         18         19         20        21
 0.03775043 0.03376847 0.02698616 0.02698616 0.02411117 0.02411117 0.0192299
         22         23
 0.01716665 0.01367159

The predicted values (on the logit scale) are,

> predict(mymod)
         1        2         3         4         5         6         7
 -1.041893 -1.50429 -1.619889 -2.197885 -2.544682 -2.660281 -2.660281
         8        9        10        11        12        13        14
 -2.660281 -2.77588 -2.891479 -3.007079 -3.007079 -3.007079 -3.007079
        15        16        17        18        19        20        21
 -3.238277 -3.353876 -3.585074 -3.585074 -3.700673 -3.700673 -3.931872
        22        23
 -4.047471 -4.278669
> fitted(mymod) - unlogit(predict(mymod))
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

(unlogit is the following simple function),

> unlogit
function(x)
exp(x)/(1 + exp(x))

For communication purposes, it makes sense to work on the probability scale (and hence use fitted), but for all other purposes the predicted values are what we want. In particular, the predict function has an option that will give the Standard Errors of the predicted values (on the logit scale), viz.,

> predict(mymod,se=T)
$fit:
         1        2         3         4         5         6         7
 -1.041893 -1.50429 -1.619889 -2.197885 -2.544682 -2.660281 -2.660281
         8        9        10        11        12        13        14
 -2.660281 -2.77588 -2.891479 -3.007079 -3.007079 -3.007079 -3.007079
        15        16        17        18        19        20        21
 -3.238277 -3.353876 -3.585074 -3.585074 -3.700673 -3.700673 -3.931872
        22        23
 -4.047471 -4.278669
$se.fit:
         1         2        3         4         5         6         7
 0.7682786 0.5943294 0.556385 0.4298909 0.4298809 0.4441498 0.4441498
         8         9        10        11        12        13        14
 0.4441498 0.4647341 0.4908397 0.5216385 0.5216385 0.5216385 0.5216385
        15        16        17        18        19        20        21
 0.5942932 0.6348851 0.7222049 0.7222049 0.7682357 0.7682357 0.8637766
        22       23
 0.9129226 1.013299
$residual.scale:
[1] 1.189907
$df:
[1] 21

So, a 95% confidence interval for the first predicted value on the logit scale would be

−1.041893 ± 2 × 0.7682786

Transferred to the probability scale, this is an interval (0.070, 0.621), with a center at 0.260. But what about predictions nearer the take-off temperature?

 predict(mymod, data.frame(Temp=c(30,31,32,33)), se=T)
$fit:
        1        2        3       4
 1.616887 1.501288 1.385689 1.27009
$se.fit:
        1        2        3       4
 1.972255 1.917739 1.863302 1.80895

(That was pretty sneaky, but the predict function requires a data frame predictor values).

The bottom line is that at 31^{∘}F our model extrapolates to p = 0.83, with a 95% confidence interval of (0.089, 0.996), which is VERY wide. Even with this wide interval, the point estimate is something to worry about.

Evaluating the Model

> plot(mymod)

Notice that the residual plot has at least two recognizable components, but it couldn’t help but look like this. There is perhaps some suggestion that a function of temperature might be needed, like perhaps t^{2}, but don’t be fooled. The places where there are big residuals are for temperatures where there were both no failures, and failures. I tried a number of additional functions of temperature, and nothing much worked.

The Zero-One version of the model

Very quickly, lets look at the model when we just look at a failure or not.

> Bn <- cbind(sign(Fail), 1 - sign(Fail))
> Bn
   [,1] [,2]
 1    1    0
 2    1    0
 3    1    0
 4    1    0
 5    0    1
 6    0    1
 7    0    1
 8    0    1
 9    0    1
10    0    1
11    0    1
12    0    1
13    1    0
14    1    0
15    0    1
16    0    1
17    0    1
18    1    0
19    0    1
20    0    1
21    0    1
22    0    1
23    0    1

then

> mymod <- glm(Bn ~ Temp, family=binomial)
> summary(mymod)
Call: glm(formula = Bn ~ Temp, family = binomial)
Deviance Residuals:
       Min       1Q     Median        3Q      Max
 -1.061112 -0.76128 -0.3782822 0.4523828 2.217506

Coefficients:
                  Value Std. Error   t value
(Intercept)  15.0422911  7.3366528  2.050293
       Temp  -0.2321537  0.1076141 -2.157279
(Dispersion Parameter for Binomial family taken to be 1 )
    Null Deviance: 28.26715 on 22 degrees of freedom
Residual Deviance: 20.31519 on 21 degrees of freedom
Number of Fisher Scoring Iterations: 4
Correlation of Coefficients:
     (Intercept)
Temp -0.99716
> anova(mymod)
Analysis of Deviance Table
Binomial model
Response: Bn
Terms added sequentially (first to last)
     Df Deviance Resid. Df Resid. Dev
NULL                    22   28.26715
Temp  1  7.95196        21   20.31519

So far things look very similar, including the residual plot (believe me).

Another Example, Stroke

Background

The data are from a prospective study whose purpose was to characterize the spontaneous recovery of language functioning in patients during the first 3 to 4 months following a stroke. Details of the design and results of this study may be found in Holland et al. (1985) and Holland, Greenhouse, Fromm, and Swindell (1989). Briefly, patients admitted to a hospital with the diagnosis of stroke and who met the study inclusion criteria were seen daily (6 days/week) for 15 minutes beginning 48 to 72 hours poststroke and continuing until hospital discharge. The daily visits consisted of semistructured conversational interactions and systemic scoring of the patient’s communication behavior. At hospital discharge, patients were formally tested using the Western Aphasia Battery (WAB). If the WAB indicated aphasic dysfunction (AQ < 93.8), the patient was seen and assessed at 1 month postdischarge. If this second WAB demonstrated continued difficulty the patient was assessed again at 2 months postdischarge. The data set to be discussed here consists of 61 patients who had unilateral (right- or left-hemisphere) stroke, ischemic or hemorrhagic etiology, and impaired language functioning at time of hospital discharge. We will use these data to investigate the correlates of language recovery during the two month period post-hospital discharge. Definitions and the codes for the variables of interest are presented in Table 2.

The multiple logistic regression model will allow us to investigate simultaneously the effect of a number of different variables on the likelihood of spontaneous language recovery. A feature of the logistic regression model is that it allows us to use explanatory variables that are either quantitative or discrete. Using the variables defined and coded in Table 2 for the stroke recovery study, the specification of the deterministic component of the multiple linear logistic regression model on the logit scale is

ln(\frac{π}{1−π}) = β_{0} + β_{1} · (SIDE) + β_{2} · (GENDER) + β_{3} · (AGE) + β_{4} · (LOSD) + β_{5} · (W AB1). (1)

Note that if we were to dichotomize length of stay, age and WAB1 then a multiple crossclassification table of these data would consist of 2^{6} = 64 cells or subgroups and that with 61 patients such a table would have numerous sparse or empty cells which would be a problem for a chi-square analysis. An advantage of a statistical model such as the one given in equation (1), however, is that the model provides a structure that smooths the fit of the probability of a response to all the cells in such a table.

Side  Gender  Age  LOSD  WAB1  Status
  0     1      74   17   92.5    1
  1     1      75   30   87.8    0
  0     1      72   17   65.9    0
  0     1      60   21   86.4    1
  0     0      69   11   90.0    1
  0     0      62   13   6.2     0

Table 1: Example of what the data file ’stroke.dat’ looks like

Interpretation of Regression Coefficients

The output below presents part of the S-Plus output listing the estimates of the regressions coefficients and their standard errors. We see that the coefficients for gender, length-of-stay, and the WAB score at discharge are all significantly related to the probability of spontaneous language recovery at the 0.05 level. The coefficient for length-of-stay is negative indicating that the longer patients are in the hospital the smaller the chance of spontaneous language recovery relative to patients who have shorter lengths-of-stay even after controlling for the effect of the WAB1 score. The coefficient for WAB1 is positive as we would expect. Males, after controlling for the other variables in the model, are exp(2.88) = 17.8 times as likely to recover spontaneous language function as females. A 95% confidence interval for β_{2}, the

Variable Description  Values/Codes
STATUS Status at last  0 = Not Recovered
observation  1 = Recovered
SIDE  Side of lesion  0 = Left

1 = Right

GENDER  0 = Female

1 = Male

AGE  Age at stroke  In years

onset

LOSD  Length of hospital In days

stay

WAB1  WAB at  Range: 0 - 100

hospital discharge

Table 2: Variables and Codes from the Stroke Recovery Study

adjusted log odds ratio of recovery for males relative to females, however, is large and ranges from 0.47 to 5.29. We also note that side of stroke and subject’s age do not have a statistically significant association with spontaneous language recovery after controlling for the effects of the other variables. Details about these tests can be found in Hosmer and Lemeshow (1989, Chapter 5).

> stroke <- read.table("stroke.dat", head=T)
> stroke[1:10,]
  Side Gender Age LOSD WAB1 Status
1    0      0  79   27 18.2      0
2    0      0  82   17 86.0      0
3    1      0  72   42 46.3      0
4    0      0  59   23 35.6      0
5    0      0  60   28  0.0      0
6    1      1  75   12 81.9      1
7    0      0  75   53 56.0      0
8    0      0  85   46 17.3      0
9    0      0  79   58 74.9      0
10    0      1  82   34 18.2      0

> mymod <- glm(Status ~ Side + Gender + Age + LOSD +WAB1, family=binomial,
+ data=stroke)
> summary(mymod)
Call: glm(formula = Status ~ Side + Gender + Age + LOSD + WAB1, family = binomial,
        data = stroke)
Deviance Residuals:
       Min         1Q       Median            3Q      Max
 -1.317524 -0.2999725 -0.003778562 -0.0002212266 2.862154
Coefficients:
                    Value Std. Error    t value
(Intercept) -11.072759221 6.28633143 -1.7614024
       Side  -1.369026127 1.10318267 -1.2409786
     Gender   2.876569185 1.25325078  2.2952862
        Age  -0.006938071 0.03062231 -0.2265691
       LOSD  -0.196137069 0.09595028 -2.0441532
       WAB1   0.164467235 0.07172681  2.2929674
(Dispersion Parameter for Binomial family taken to be 1 )
    Null Deviance: 65.71913 on 60 degrees of freedom
Residual Deviance: 28.9878 on 55 degrees of freedom
Number of Fisher Scoring Iterations: 7
Correlation of Coefficients:
       (Intercept)       Side     Gender        Age       LOSD
  Side  0.2553975
Gender -0.5007016  -0.2793548
   Age -0.3373404  -0.1265137  0.2006563
  LOSD  0.1199385   0.2046374 -0.5309336 -0.1894954
  WAB1 -0.9292116  -0.3055970  0.5278826  0.0616042 -0.2969773
> anova(mymod)
Analysis of Deviance Table
Binomial model
Response: Status
Terms added sequentially (first to last)
       Df Deviance Resid. Df Resid. Dev

  NULL                    60   65.71913
  Side  1  0.10011        59   65.61902
Gender  1  3.36010        58   62.25892
   Age  1  1.45449        57   60.80444
  LOSD  1 12.25425        56   48.55018
  WAB1  1 19.56238        55   28.98780
>

Analysis of the Fit of the Model

Figure 1 presents a residual analysis for this logistic regression.

One relatively easy method for the analysis of the fit of the logistic regression model is based on calculating the predicted probability of spontaneous language recovery based on the model in (1) compared to the observed probability of language recovery Specifically, substitute the estimated regression coefficients given in above into the multiple logistic regression equation in (1). For each patient in the study calculate his or her individual predicted probability of language recovery by calculating the linear combination of each patient’s individual set of explanatory variables with the estimated regression coefficients and use the inverse of the logit transformation to solve for the probability of language recovery, π̂. For example, consider a patient whose stroke was on the left side (SIDE=0), was a 61 year old male (GENDER=1, AGE=61), stayed in the hospital 19 days (LOSD=19), and whose WAB at discharge was 81.6 (WAB1=81.6). The estimated probability of spontaneous language recovery based on the model for such a patient is

π̂ = \frac{1}{1 + e^{−(−11.1−1.37·(0)+2.88·(1)−0.0069·(61)−0.196·(19)+0.165·(81.6))}} = 0.77.

This value can be interpreted to mean that 77 patients out of 100 with the above characteristics can be expected to recover normal language functioning spontaneously within the first 2 months post-hospital discharge. To study the goodness of fit of the model, we obtain the frequency distribution of the values of π̂, and determine, say, the quintiles of the distribution. The observed number of patients who had normal language functioning are tallied within each quintile, and the expected number of patients who had normal language functioning (predicted by the model) is determined by summing the π̂’s for all the patients in the quintile. Table 3 presents the expected and observed number of patients who had normal language functioning for each quintile of risk of language recovery. Comparisons of the expected number predicted by the model to the observed number suggests that this model fits the data reasonably well.

A test for whether there is a significant difference between the observed and expected numbers in Table reftab:fit can be constructed based on a Pearson chi-square test. This chi-square test would have degrees-of-freedom equal to the number of groups (e.g., 5 for the number of quintiles) minus 2 (Hosmer and Lemeshow, 1989, chapter 5). In our example, this chi-square statistics is equal to 3.60 with 3 degrees-of-freedom and has a non-significant p-value=0.69, indicating that the model adequately fits the data. A small p-value would indicate a significant lack of fit of the model.

Expected Observed
Quintile of risk Number Number
Highest  5  10.7  11
4  3.1  2
3  0.2  1
2  0.0  0
Lowest  1  0.0  0
Total  14.0  14

Table 3: Analysis of the fit of the logit model in Equation (1) Figure 1: Residual plots

More Logistic Regression

How to plot these data

Look at the stroke data. This is explicitly 0–1 data, and it is very hard to think about plotting it. A simple plot of, say, Status against WAB1, looks like the following, and it not all that informative.

> plot(WAB1, Status)

Even if I jitter the points a little, it is still not very easy to see any trends or patterns.

> plot(WAB1, Status+0.1*runif(length(Status)))

A somewhat more useful approach is to plot the data and add some sort of smooth curve that will allow you to at least look at some of the trends in the data. The particular approach I took with these data was to superimpose a (black-box) smooth curve generated by loess. The various pictures result from different smoothing parameters.

> scatter.smooth(WAB1, Status,span=0.2)
> title("Span = 0.2")
> scatter.smooth(WAB1, Status,span=0.35)

> title("Span = 0.35")
> scatter.smooth(WAB1, Status,span=0.5)
> title("Span = 0.5")
> scatter.smooth(WAB1, Status,span=0.65)
> title("Span = 0.65")

This is still not perfect, but you certainly get the sense that the probability goes up, probably more or less monotonely with WAB1.

Compare this to the similar picture(s) for Age.

> for (i in c(0.2, 0.35, .5, 0.65)){
+ scatter.smooth(Age, Status, span=i)
+ title(paste("Span = ", i))
+ }

If anything, there is a bit of evidence for some type of quadratic effect of Age, but the probable story is that Age is not that important, on its own. We should always remember that variables that do not seem to be important (from these marginal plots), may be quite important when paired with other variables (that is viewed conditionally).

With that hint about the Age variable, I re-ran the logistic regression model for these data and obtained,

> mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age + Age^2, family=bin>
> summary(mymod)
Call: glm(formula = Status ~ Side + Gender + LOSD + WAB1 + Age + Age^2, family =
        binomial)
Deviance Residuals:
       Min         1Q        Median            3Q      Max
 -1.317304 -0.2637303 -0.0002506098 -6.014239e-07 2.770847
Coefficients:
                   Value   Std. Error   t value
(Intercept) -42.13652444 17.240544854 -2.444037
       Side  -1.70612965  1.320176507 -1.292350
     Gender   4.03384776  1.847166535  2.183803

       LOSD  -0.28388626  0.130491084 -2.175522
       WAB1   0.24891576  0.104667163  2.378165
        Age   0.92453559  0.421765743  2.192059
   I(Age^2)  -0.00791642  0.003578441 -2.212254
(Dispersion Parameter for Binomial family taken to be 1 )
    Null Deviance: 65.71913 on 60 degrees of freedom
Residual Deviance: 21.7621 on 54 degrees of freedom
> anova(mymod)
Analysis of Deviance Table
Binomial model
Response: Status
Terms added sequentially (first to last)
         Df Deviance Resid. Df Resid. Dev
    NULL                    60   65.71913
    Side  1  0.10011        59   65.61902
  Gender  1  3.36010        58   62.25892
    LOSD  1 13.46451        57   48.79441
    WAB1  1 19.75552        56   29.03889
     Age  1  0.05110        55   28.98780
I(Age^2)  1  7.22570        54   21.76210

This is somewhat surprising (because I didn’t think there was an Age effect). Could we have seen the same trend from a residual analysis of these data?

Let us go back to the original model and look at residuals.

There is at least one large (Pearson) residual, and it is the same point as the large (deviance) residual. What about a plot of (deviance) residual against Age (when Age is in the model).

> mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age, family=binomial)
> scatter.smooth(Age, (residuals(mymod)))
> identify(Age, (residuals(mymod)))
[1] 59
> scatter.smooth(Age[-59], (residuals(mymod))[-59])

We see that observation 59 is the one with the large residual. If I plot the data (the lower plot) without that observation, the quadratic trend seems to go away, at least a bit.

Observation 59 is now suspect, so let us remove that observation, and see what the resulting model looks like.

> WW <- rep(1, length(Side))
> WW[59] <- 0
> mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age,
+ weight=WW, family=binomial)
Warning messages:
  linear convergence not obtained in  10  iterations. in: glm.fitter(x = X, y =
        Y, w = w, start = start, offset = offset, family = family, ....

Oh-oh, what does that mean? (Basically, the optimization algorithm is converging slowly, which is an indication that there are some fitted probabilities that are very close to 0 or 1). We can try a few more iterations, via.

> mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age,
+ weight=WW, maxit=30,family=binomial)

Now, what does the summary and anova look like?

> summary(mymod)

Call: glm(formula = Status ~ Side + Gender + LOSD + WAB1 + Age, family = binomial,
        weights = WW, maxit = 30)
Deviance Residuals:
      Min           1Q        Median            3Q      Max
 -1.21685 -0.006138687 -2.107342e-08 -2.107342e-08 2.120409
Coefficients:
                    Value  Std. Error     t value
(Intercept) -46.034665461 23.60416082 -1.95027757
       Side  -7.479606389  4.79366143 -1.56031178
     Gender   7.889726328  4.28891089  1.83956406
       LOSD  -0.677382165  0.43418335 -1.56012928
       WAB1   0.662724777  0.34770876  1.90597667
        Age  -0.001289562  0.04096666 -0.03147832
(Dispersion Parameter for Binomial family taken to be 1 )
    Null Deviance: 62.71879 on 60 degrees of freedom
Residual Deviance: 13.32727 on 54 degrees of freedom
Number of Fisher Scoring Iterations: 10
Correlation of Coefficients:
       (Intercept)       Side     Gender       LOSD       WAB1
  Side  0.8469550
Gender -0.9187997  -0.8860396
  LOSD  0.7944690   0.9174434 -0.8555061
  WAB1 -0.9831416  -0.9061715  0.9317554 -0.8758313
   Age -0.0387868   0.0615410  0.0243350 -0.0022418 -0.0596445
Warning messages:
  1 rows with zero weights not counted in: summary(mymod)

If I now start adding the quadratic term in Age, I get myself into all sorts of trouble.

> formula(mymod)
Status ~ Side + Gender + LOSD + WAB1 + Age + Age^2
> anova(mymod)
Analysis of Deviance Table
Binomial model
Response: Status

Terms added sequentially (first to last)
         Df Deviance Resid. Df Resid. Dev
    NULL                    60   62.71879
    Side  2  0.62858        58   62.09021
  Gender  1  2.34693        57   59.74328
    LOSD  1 13.67051        56   46.07277
    WAB1  1 32.74451        55   13.32826
     Age  1  0.00099        54   13.32727
I(Age^2)  1 13.32727        53    0.00000
Warning messages:
1: linear convergence not obtained in  10  iterations. in: glm.fit(x, y, w,
        offset = offset, family = family, maxit = control$maxit,  ....
2: linear convergence not obtained in  10  iterations. in: glm.fit(x, y, w,
        offset = offset, family = family, maxit = control$maxit,  ....

Whoa. What is going on here? We suddenly have a model that fits perfectly?

The summary looks even weirder?

> summary(mymod)
Call: glm(formula = Status ~ Side + Gender + LOSD + WAB1 + Age + Age^2, family =
        binomial, weights = WW, maxit = 50)
Deviance Residuals:
         Min            1Q        Median            3Q          Max
 -3.9065e-05 -2.107342e-08 -2.107342e-08 -2.107342e-08 3.281089e-05
Coefficients:
                    Value  Std. Error       t value
(Intercept) -2858.6085408 2227320.001 -0.0012834297
       Side  -164.0618144  165491.082 -0.0009913635
     Gender   273.5752941  237195.052  0.0011533769
       LOSD   -16.2464405   13633.473 -0.0011916582
       WAB1    24.1122280   17153.825  0.0014056473
        Age    39.8782980   35880.527  0.0011114190
   I(Age^2)    -0.3549147     302.566 -0.0011730157
(Dispersion Parameter for Binomial family taken to be 1 )
    Null Deviance: 62.71879 on 60 degrees of freedom
Residual Deviance: 0 on 53 degrees of freedom
Number of Fisher Scoring Iterations: 29
Correlation of Coefficients:

         (Intercept)       Side     Gender       LOSD       WAB1        Age
    Side  0.3661039
  Gender -0.8397332  -0.5132934
    LOSD  0.7676113   0.5867586 -0.7405370
    WAB1 -0.9678428  -0.5023038  0.8384746 -0.7822973
     Age -0.9537356  -0.2802197  0.7986261 -0.8085521  0.8643352
I(Age^2)  0.9585630   0.3372172 -0.8129406  0.8428222 -0.8863173 -0.9963403
Warning messages:
1: fitted values close to 0 or 1 in: family$deviance(mu, y, w, residuals = T)
2: 1 rows with zero weights not counted in: summary(mymod)
>

On reflection, it is clear what is going on. Look at the first plot of Status against WAB1. Here we see, that absent the one observation in the middle (which is 59), we have an almost perfect discimination between the two levels of Status. That is the clue to what is “broken”. With the 7 parameters (anything like WAB1 squared would have done), we can find a function that perfectly discriminates between the two levels of status, and then maximum likelihood modelling breaks down. There is a lot more to be said here, I just don’t know what it is.

End of Stroke example

Influence of a point

There are good measures of influence for glm models, but they don’t seem to be implemented (i.e., I could not find them) in the standard packages. Rather than just give up, I quickly hacked together a “Poor Man’s” influence function. The basic idea of the following code is to fit the glm model with each data point in turn removed. The summary is then a plot of each of the coefficients, with as many points on the plot as there are data points. This plot is surprisingly effective and quickly shows us where the pathology in the data lies.

# Heavily annotated, for class
# Generate a vector of all ones, of length the number of observations.
W <-  rep(1,length(Side))
# Generate current weights.
WW <- W
# Fit the base model
mymod <- glm(Status ~ Gender + LOSD + Side + Age + WAB1, weight=WW,
   family=binomial)
# Save the coefficients, just in case
base.coef <- coefficients(mymod)
# Set up a matrix to hold the coefficients for each model with one
# point deleted.
infs <- matrix(0, nrow=length(Side), ncol=length(coefficients(mymod)))
# Loop through the data, setting one weight equal to 0.
# Save the coefficients in the infs matrix
for (i in 1:length(Side)){
        WW<-W; WW[i] <- 0;
        infs[i,] <- coefficients(update(mymod, . ~ .))
}
## Can't think of an easy, smart summary, so just plot the
## coefficients.
#Set up a 2 x 3 plotting frame
par(mfrow=c(2,3))
#Do the 6 plots.  After each plot, allow the user to
# identify the points which look interesting.
for(i in 1:6){
        plot(1:length(Side), infs[,i])
        title(names(coefficients(mymod))[i])

        identify(1:length(Side), infs[,i])
}

After running that code, and interacting with the various plots, you end up with the following plot.

There are other things we could have plotted (like following the deviance as the models changed), but this plot alone tells us much of the story.

More Examples

Here I will pull together a number of examples, drawn directly from the Book by Venables and Ripley.

First, the budworm data. The following data is the result of an experiment on the toxicity of the tobacco budworm to doses of particular drug, that the moths are beginning to show resistance to. Batches of 20 moths of each sex were exposed for 3 days to the chemical and the number of dead (or disabled) moths in each group were recorded. The data is,

Dose Response Gender
1       1       M
2       4       M
4       9       M
8       13      M
16      18      M
32      20      M
1       0       F
2       2       F
4       6       F
8       10      F
16      12      F
32      16      F

The doses were in micro-grams, and the dose level when up by powers of two. It is a judgement call as to whether you will try to model with dose as the response or log_{2}(dose. We will try both.

> bud <- read.table("budworm.dat", head=T)
> bud
   Dose Response Gender
 1    1        1      M
 2    2        4      M
 3    4        9      M
 4    8       13      M
 5   16       18      M
 6   32       20      M
 7    1        0      F
 8    2        2      F
 9    4        6      F
10    8       10      F
11   16       12      F
12   32       16      F
> attach(bud,pos=1)
> N <- 20

First, lets plot the data, and see if anything obvious shows up.

> plot( Dose, Response/N)
> points( Dose[Gender=="F"], (Response/N)[Gender=="F"], pch="F")
> points( Dose[Gender=="M"], (Response/N)[Gender=="M"], pch="M")

The same plot on the log(Dose) scale, is

> plot( Dose, Response/N, log="x")
> points( Dose[Gender=="F"], (Response/N)[Gender=="F"], pch="F")
> points( Dose[Gender=="M"], (Response/N)[Gender=="M"], pch="M")

Given that we have more than 0-1 data, it makes a lot of sense to at least try to plot the (observed) logits,

> par(mfrow=c(2,1))
> plot( Dose, logit(Response/N), log="x")
> title("Logit vs Log(Dose)")
> plot( Dose, logit(Response/N))
> title("Logit vs Dose")

This plot makes it quite clear, that the log scale has a much better chance of working. So, proceed to fit the model. We surely want a model with different slopes for the two

sexes

> ldose <- log(Dose)
> SF <- cbind(Response, N - Response)
> mymod <- glm(SF ~ Gender * ldose, family=binomial)
> summary(mymod)
 ..........
Coefficients:
                  Value Std. Error    t value
 (Intercept) -2.9935414  0.5525295 -5.4178852
      Gender  0.1749865  0.7781556  0.2248733
       ldose  1.3071341  0.2410133  5.4234939
Gender:ldose  0.5091459  0.3894474  1.3073547

    Null Deviance: 124.8756 on 11 degrees of freedom
Residual Deviance: 4.993727 on 8 degrees of freedom
> anova(mymod,test="Chisq")
Analysis of Deviance Table

Binomial model
Response: SF
Terms added sequentially (first to last)
             Df Deviance Resid. Df Resid. Dev   Pr(Chi)
        NULL                    11   124.8756
      Gender  1   6.0770        10   118.7986 0.0136955
       ldose  1 112.0415         9     6.7571 0.0000000
Gender:ldose  1   1.7633         8     4.9937 0.1842088
>

Now, what is going on here? Gender seems to have no effect (based on the t-value), but the deviance is quite large? Even more, our plots seem to imply what looks like a gender effect. The best way to find out is to distrust the parameter stuff (it can get easily confused), and try fitting a model without the interaction. This makes sense anyway, because there is little/no evidence of an interaction.

> anova(mymod <- glm(SF ~ Gender + ldose, family=binomial))
Terms added sequentially (first to last)
       Df Deviance Resid. Df Resid. Dev
  NULL                    11   124.8756
Gender  1   6.0770        10   118.7986
 ldose  1 112.0415         9     6.7571
> summary(mymod)
Coefficients:
                Value Std. Error   t value
(Intercept) -3.473154  0.4682939 -7.416612
     Gender  1.100743  0.3557226  3.094385
      ldose  1.535336  0.1890119  8.122959
    Null Deviance: 124.8756 on 11 degrees of freedom
Residual Deviance: 6.757064 on 9 degrees of freedom

That is a little more like it, and seems to be a VERY good fitting model. We should at least look at some diagnostics, and some other plots, but unless they show something very odd, we are done.

> plot(mymod)

That looks very good, we should finish by plotting the data and fitted curves on the same graph, and that is left as an exercise for you (or look it up in the book). The book also goes on to talk about LD50, and how to estimate it