Worked example 1
Our first worked example uses raw binary data where each binary observation has unique values of one or more explanatory variables for each individual case. In this situation we must feed in each observation with a 1 or a 0 along with its values for the explanatory variables. We use data from Crawleyon costs of reproduction in a perennial plant. The first column (State) indicates survival through the winter (1) or death between growing seasons (0). Explanatory variables are seed production this summer (Flowers) and the size of its rootstock (Root).
Examine relationships between variables
Plot out data to get a visual assessment of the relationships between variables.
As you can see the raw binary data (black circles) are not very easy to interpret in this form as they are all either 0 or 1. However, one can observe that none of the plants that produced the largest number of flowers survived (there is none in the top right hand corner of the left hand plot). On the contrary not many of the plants with the biggest roots died (bottom right of the right hand plot).
We can get a much better idea of what is going on if we fit simple univariate logistic regressions to each graph, and the look at the fitted relationships (shown as red points in the flip). The probability of survival is negatively related to flower production and positively related to root size.
Fit full model with interaction
Below we have fitted the full model which contains both flower production and root size and their interaction.
As expected we obtain a positive relationship with root size and a negative relationship with flower production. But in addition we have a significant interaction between these two factors. Since this complicates interpretation, we check whether we can drop the interaction (albeit its significance level suggests it cannot be dropped).
Compare models with and without interaction
We use anova to provide an analysis of deviance table for the two models. Note that we have used both chi square and Mallow's Cp as test statistics.
Clearly we cannot drop the interaction term - the difference between the two models is highly significant using chi square, and the model without the interaction has a much larger value for Cp.
Examination of the model coefficients shows that a given level of flowering has a bigger impact on reducing survival on plants with small roots than on larger plants.
Note that we cannot use the residual deviance to tell us anything about goodness of fit of the model because we have entered the data in raw binary format.
Worked example 2
Our second example uses a result from a cross-sectional survey on the prevalence of dystocia in cats. We previously looked at this work in relation to the confidence intervals attached to the odds and risk ratios in relation to breed of cat.
|Prevalence of dystocia in cats in relation to breed
|Breed||No. positive||No negative||% positive
Fit logistic regression model
We have highlighted the key parameters in red. Note that we cannot use the residual deviance to tell us anything about goodness of fit of the model because the residual deviance is (effectively) equal to zero, as are its degrees of freedom. This is because we are fitting the saturated model. The natural log of the odds ratio is 2.7528 with a standard error of 1.1804, exactly as we obtained previously. We can detransform the odds ratio as exp(2.7528) to give 15.686.
Obtain likelihood confidence interval
The P-value of 0.0197 for breedManx suggests we can accept this odds ratio differs significantly from 1. But this should not be accepted at face value. Instead work out the confidence interval for the odds ratio. You cannot use a normal approximation confidence interval because this is invalid when some frequencies are very low. You need to estimate a likelihood confidence interval:
The (detransformed) interval for the odds ratio is 0.755 - 129.85 - in other words it overlaps 1, so you cannot be confident that the effect is real. You may note that if you just plugged in your result to a logistic regression package (such as here) you get no warning about the problem of very small frequencies. Just because a regression package will do a logistic regression on your 2 × 2 table (or any other set of data), does not mean that the result will necessarily be reliable.
Worked example 3
|Probability of developing|
We use data from Agresti
on the proportion of people developing AIDS symptoms (given as number positive and number negative) in relation to taking the antiretroviral drug azt (azt = yes or no) and race (col = black or white). Hence we have two binary explanatory variables and a binary response variable.
Fit full model
We first fitted the full model with the interaction between race and drug use. The interaction proved non-significant (in fact all terms were non-significant), so we dropped the interaction term to give the following 'simple main effects' model.
This suggests that taking AZT significantly reduces the chance of developing AIDS symptoms (P = 0.0099). The effect of race (col) (P = 0.848) is not significant and can possibly be dropped from the model. Note that the residual deviance (at 1.3835) is similar to the degrees of freedom (1), suggesting there is no problem of overdispersion.
Test effect of dropping race from model
We drop race (col) from the model using the function 'update' in base R, and test whether the change in deviance is significant using 'anova' to compare the two models.
The change in deviance is not significant (P = 0.8473) so race is dropped from the model.
Summarize final model and obtain likelihood CI
The likelihood CI to the log(odds ratio) for developing AIDS symptoms for taking versus not taking AZT is obtained using the function 'confint' (in the package MASS but now usually included in base R).
The (detransformed) odds ratio for developing AIDS symptoms for taking versus not taking AZT is therefore e-0.7218 = 0.486 (95% CI: 0.278 to 0.833). The residual deviance is 1.4206 on 2 degrees of freedom which is not significant suggesting a good fit to the model.
One could proceed to estimate probabilities for each of the four predictor combinations (see p87 in Laura Thompson (2009) for this example) but this seems unjustified if one is dropping race from the model.
We have used P-values in the model simplification process, but better approach would be to use AIC. We have done this below using the R 'step' function on the full model (including the interaction) to demonstrate that in this case it gives the same result (key values in red).
Check goodness of fit
We have already been given the deviance residual goodness of fit statistic above as the residual deviance (1.421) with 2 degrees of freedom, but we will print out this again along with the Pearson goodness of fit statistic, and (for good measure) two measures of pseudo-R-square, the Cox & Snell R square and Nagelkerke's R square (thanks to Sören Vogel.)
The deviance goodness of fit statistic is 1.4206 with 2 df, so P = 0.491.
The Pearson goodness of fit statistic is 1.414 with 2 df, so P = 0.493.
Both these are compatible with a model which explains most of the variation. The Cox & Snell R square (at 0.823) and Nagelkerke's R square (at 0.940) also indicate a very acceptable model.