Notes for Thursday, March 19, 1998

An Extended 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 cross-classification 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.

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

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

females. A 95% confidence interval for β_{2}, the 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).

402 > stroke <- read.table("stroke.dat", head=T)
402 > 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
402 > mymod <- glm(Status ~ Side + Gender + Age + LOSD +WAB1, family=binomial,
+ data=stroke)
402 > 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

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

     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
402 > 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
402 >

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 posthospital 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)

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.

402 > plot(WAB1, Status)

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

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

Figure 1: Residual plots

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.

402 > scatter.smooth(WAB1, Status,span=0.2)
402 > title("Span = 0.2")
402 > scatter.smooth(WAB1, Status,span=0.35)
402 > title("Span = 0.35")
402 > scatter.smooth(WAB1, Status,span=0.5)
402 > title("Span = 0.5")
402 > scatter.smooth(WAB1, Status,span=0.65)
402 > 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.

402 > 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,

402 > mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age + Age^2, family=bin>
402 > 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
402 > 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).

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

402 > WW <- rep(1, length(Side))
402 > WW[59] <- 0
402 > 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.

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

Now, what does the summary and anova look like?

402 > 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.

402 > formula(mymod)
Status ~ Side + Gender + LOSD + WAB1 + Age + Age^2
402 > 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?

402 > 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)
402 >

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.

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.

36-402: Advanced Data Analysis II Spring 1998 16 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.