Detailed Outline 287 D. Modeling strategy: All Es, no Cs q qq Model : Logit PðXÞ ¼ a þ ~ biEi þ ~ ~ d*ii0 EiEi0 i¼1 i¼1 i0¼1 i6¼i0 Step 1: Define initial model (above) Step 2: Assess interaction involving Es. Option A*: Overall chunk test for EEs, followed by backward elimination of EEs Option B*: Skip chunk test for EEs; start with backward elimination of EEs Skip previous Step 3 Step 4: Test for nonsignificant Es if not components of significant EEs E. How causal diagrams can influence choice of initial model? III. Screening Variables (pages 263–270) A. Problem Focus: Model contains one E, and a large number of Cs and E Â Cs, but computer program does not run or fitted model unreliable (“large” p) B. Screening: Exclude some Cj one-at-a-time; fit reduced model C. Method 0: Consider predictors one-at-a-time; screen-out those Xi not significantly associated with the outcome (D) D. Questions and Brief Answers about Method 0: 1. Any criticism? Yes: does not consider confounding or interaction involving Cs 2. Depends on types of Xs? Yes: use if only Es and no Cs. 3. How large k compared to n? No good answer. 4. Other ways than Method 0? Yes: evaluate confounding and/or interaction for Cs. 5. Collinearity and/or screening? Consider collinearity prior to and following screening. E. Assessing Confounding and Interaction when Screening C variables. Confounding: Compare Logit P(X) ¼ a þ bE with Logit P*(X) ¼ a* þ b*E þ g*C Does OdRDE ¼ eb^ 6¼ OdRDEjC ¼ eb^* ? Interaction: Test H0: d ¼ 0 for the model Logit P(X) ¼ a þ bE þ gC þ dEC F. How to proceed if several Es and several Cs: It depends! G. How to proceed if several Es and no Cs: Use method 0.
288 8. Additional Modeling Strategy Issues IV. Collinearity (pages 270–275) A. The Problem: If predictors (Xs) are “strongly” related, then b^j unreliable, Va^rb^j high, or model may not run. B. Diagnosing Collinearity. Use condition indices (CNIs) and variance decomposition proportions (VDPs). Collinearity detected if: largest CNI is large (>30?) and at least 2 VDPs are large (!0.5?) C. Collinearity for Logistic Regression Requires computer macro (program not available in popular computer packages). CNIs derived from inverse of Information Matrix (I21) D. Difficulties How large is large for CNIs and VDPs? Guidelines provided are “soft.” How to proceed? We recommend sequential procedure: fix one collinearity problem at a time. How to fix problem? Usual approach: drop one of the collinear variables; or, define new variable. E. Example using MRSA data V. Influential Observations (pages 275–279) A. The Problem: Does removal of subject from the data result in “significant” change in b^j or OdR? B. Measures: Delta-betas (Dbs) and Cook’s distance- type measures (Cs). C. Computer packages: provide plots of Dbs and Cs for each subject. D. What to do with influential observations: Not easy to decide whether or not to drop subject from the data. Conservative approach: drop subjects only if their data is incorrect and cannot be corrected. VI. Multiple Testing (pages 280–282) A. The Problem: should you adjust a when performing several tests? B. Bonferroni approach: Use a ¼ a0/T, where a0 ¼ family-wise error rate (FWER) and T ¼ number of tests. C. Criticisms of Bonferroni approach: low power; based on unrealistic “universal H0”; other. D. Model building problem: number of tests (T) not known in advance; therefore, no foolproof approach. VII. Summary (pages 283–285)
Practice Practice Exercises 289 Exercises 1. Consider the following logistic regression model in which all predictors are (0,1) variables: logit PðXÞ ¼ a þ b1E1 þ b2E2 þ b3E3 þ g1C1 þ g2C2 þ g3C3 þ g12C1C2 þ d11E1C1 þ d21E2C1 þ d12E1C2 þ d22E2C2 þ d33E3C3 þ de13E1E3 þ de23E2E3 þ d112E1C1C2 þ d212E2C1C2 For the above model, determine which of the following statements are True or False. (i.e., Circle T or F) T F a. The above model is hierarchically well-formu- lated. T F b. Suppose the chunk test for the null hypothesis H0: d112 ¼ d212 ¼ 0 is found to be significant and backward elimination involving these two three- factor product terms results in only E2C1C2 remaining in the model. Then based on the “hier- archy principle,” the final model must contain the variables C1C2, C1, C2, E2C1C2, E2C1, E2C2, E2E3, and E2. T F c. Suppose the chunk test for the null hypothesis H0: d112 ¼ d212 ¼ 0 is found to be significant and, as in the previous question, backward elimina- tion involving these two three-factor product terms results in only E2C1C2 remaining in the model. Then, based on the hierarchy principle and the hierarchical backward elimination approach, the only variables that remain as can- didates for being dropped from the model at this point are E1, E3, E1E3, E2E3, E1C1, E3C3, and C3. T F d. Suppose that after the interaction assessment stage, the only terms remaining in the model are E2C1C2, E2C1, E2C2, E3C3, C1C2, C1, C2, C3, E1, E2, and E3. Then, at this point, the odds ratio formula for comparing a person for whom E1 ¼ E2 ¼ E3 ¼ 1 to a person for whom E1 ¼ E2 ¼ E3 ¼ 0 is given by the expression OR ¼ exp[b1 þ b2 þ b3 þ d21C1 þ d22C2 þ d33C3 þ d212C1C2] where the coefficients in the formula are estimated from the reduced model obtained after interaction assessment. T F e. Suppose that neither E1C1C2 nor E2C1C2 remains in the model after interaction assess- ment of these two three-factor products (but prior to interaction assessment of two-factor products). Suppose further that separate Wald (and corresponding likelihood ratio) tests for
290 8. Additional Modeling Strategy Issues H0: d11 ¼ 0, H0: d21 ¼ 0, and H0: d33 ¼ 0 are non- significant in the reduced model (without E1C1C2 and E2C1C2). Then, as the next step in interaction assessment, the model should be further reduced to logit PðXÞ ¼ a þ b1E1 þ b2E2 þ b3E3 þ g1C1 þ g2C2 þ g3C3 þ g12C1C2 þ d12E1C2 þ d22E2C2 þ de13E1E3 þ de23E2E3: prior to the assessment of confounding. 2. Suppose that after interaction assessment involving both EVi and Ei Ej terms in the initial model stated in question 1, the following reduced model is obtained: logit PðXÞ ¼ a þ b1E1 þ b2E2 þ b3E3 þ g1C1 þ g2C2 þ g3C3 þ g12C1C2 þ d11E1C1 þ d22E2C2 Suppose further that the assessment of confounding will only consider changes in the odds ratio that com- pares a person for whom E1 ¼ E2 ¼ E3 ¼ 1 to a per- son for whom E1 ¼ E2 ¼ E3 ¼ 0. Based on the recommended guidelines for the assess- ment of confounding described in this chapter: a. What is the formula for the estimated odds ratio in the gold standard model that should be used to assess confounding? 2 b. Assuming that you will need to consider tables of odds ratios that consider different subsets of potential con- founders, describe what a table of odds ratios would look like for the gold standard model using the rectangle shown below for the (outside) borders of the table. In your answer, make sure to state the formulae for the odds ratios that will go into the different boxes in the table. Hint. You will need to draw horizontal and ver- tical lines to subdivide the rectangle and label the different row and column categories of the table, recognizing that the odds ratios being represented reflect the interaction effects that are present in the gold standard model.
Practice Exercises 291 c. Considering only those variables that are candidates for being assessed as nonconfounders, i. How many subsets of these variables need to be considered to address confounding? ii. List the variables that are contained in each of the above subsets. d. Suppose that the following results are obtained when comparing tables of ORs for different subsets of V variables in the previously stated (question 2) reduced model obtained after interaction assessment: Model Variables dropped Table of ORs Within 10% # from model of Gold Standard Table? 1 C1 No 2 C2 Yes 3 C3 No 4a C1C2 Yes 4b C1 and C3 Yes 5 C1 and C1C2 No 6 C2 and C1C2 No 7 C3 and C1C2 Yes 8 C1 and C2 Yes 9 C1 and C3 Yes 10 C2 and C3 Yes 11 C1, C2 and C3 Yes 12 C1, C2 and C1C2 Yes 13 C1, C3 and C1C2 Yes 14 C2, C3 and C1C2 Yes 15 C1, C2, C3, and C1C2 No 16 None Yes (GS model) Based on the above results, what models are eligible to be considered as final models after confounding assess- ment? e. Based on your answer to part d, how would you address precision? 3. In addition to the variables E1, E2, E3, C1, C2, C3 consid- ered in questions 1 and 2, there were 25 other variables recorded on study subjects that were identified from the literature review and conceptualization of the study as potential control variables. These variables were screened out by the investigators as not being neces- sary to include in the multivariable modeling analyses that were carried out. Assume that screening was carried out by putting all 25 variables in a logistic regression model together with the variables E1, E2, E3, C1, C2, and C3 and then using a back- ward elimination to remove nonsignificant variables.
292 8. Additional Modeling Strategy Issues State or describe at least three issues/problems that would not have been addressed by the above screening approach. 4. Consider the following logistic regression model in which all predictor variables are (0, 1) variables: logit PðXÞ ¼ a þ bE þ g1C1 þ g2C2 þ g3C3 þ g12C1C2 þ d1EC1 þ d2EC2 þ d3EC3 þ d12EC1C2 Determine which of the following statements are True or False (Circle T or F). T F a. Suppose that using the collinearity macro described in the text results in the three highest condition indices (CNIs) being 77, 56, and 47. Then, using a sequential approach for diagnos- ing collinearity, there is at least 1 and possibly 3 collinearity problems associated with fitting this model. T F b. Suppose that using the collinearity macro described in the text, the highest condition index (CNI) is found to be 77 and the only VDPs that are determined to be high for this CNI are associated with the variables E and C3. Then it is reasonable to drop C3 from the model, and recompute collinearity diagnostics to fur- ther reassess collinearity for the reduced model. T F c. Suppose that using the collinearity macro described in the text, the three highest condition indices (CNIs) are 77, 56, and 47. For the highest CNI of 77, suppose there are exactly four VDPs other than the intercept that are larger than .5. Suppose also that the variables associated with these four VDPs are C1, C3, EC3, and EC1C2. Then this collinearity problem can be addressed by dropping either EC3 or EC1C2 from the model, but not C1, C2 or C3, after which a reduced model is fit to see if there are additional collinearity problems. T F d. Suppose that for the same situation as described in part c above, the collinearity problem is addressed by dropping EC3 from the model, after which a reduced model without EC3 is fit. Then the collinearity diagnostics obtained for the reduced model will indicate another collin- earity problem that involves the variable EC1C2. 5. Suppose for the data analyzed using the model stated in question 4, it was desired to evaluate whether or not any subjects were influential observations. a. Assuming that the study objective was to assess the effect of E controlling for C1, C2, and C3, how can
Test Test 293 you criticize the use of a Cooks-distance-type mea- sure to identify influential subjects? b. Suppose you identified five influential subjects in the dataset. How should you decide whether or not to drop any of these subjects from the data? c. How can you criticize using the model stated in question 4 to identify influential subjects? 6. In attempting to determine a best model, starting with the model stated in question 4, suppose you wanted to use a Bonferroni-type adjustment to correct for multi- ple testing of interaction terms. What problem will you have in trying to determine the significance level for individual tests in order to achieve a family-wise error rate (FWER) not higher than 0.05? THE SCENARIO. A firm mattress is often believed to be beneficial for low back pain, although evidence supporting this recommendation is lacking. A randomized, double- blind, controlled, multicenter trial was carried out to deter- mine the effect of mattress firmness on the clinical course of patients with chronic, nonspecific low back pain. A series of 313 adults with chronic, nonspecific low back pain, without referred pain, who complained of backache while lying in bed and upon rising were randomly assigned to four types of mattresses: 3 ¼ firm, 2 ¼ medium firm, 1 ¼ medium soft, 0 ¼ soft. Clinical assessments were carried out at baseline and at 90 days. The endpoint (i.e., the health outcome variable) was improvement in pain while in bed and upon rising from bed between 0 and 90 days. The data set included the following variables: D ¼ improvement score for pain in bed from baseline to 90 days (0 ¼ no improvement or improvement in bed only, 1 ¼ improvement upon rising only or improvement both in bed and upon rising) F ¼ firmness of mattress (3 ¼ firm, 2 ¼ medium firm, 1 ¼ medium soft, 0 ¼ soft) BASE ¼ type of base of bed (1 ¼ firm, 0 ¼ not firm) POST ¼ posture while sleeping at baseline (1 ¼ supine or fetal, 0 ¼ other) PF ¼ subjective perception of firmness of mattress (3 ¼ firm, 2 ¼ medium firm, 1 ¼ medium soft, 0 ¼ soft) OCC ¼ type of occupation (1 ¼ sedentary, 0 ¼ not sedentary) AGE ¼ age in years of subject GEN ¼ gender (0 ¼ male, 1 ¼ female)
294 8. Additional Modeling Strategy Issues Questions about how to analyze these data now follow: 1. In addition to the variables listed above, there were 12 other variables also listed in the dataset and identified from the literature review and conceptualization of the study as potential control variables. These variables were screened out by the investigators as not being necessary to include in the multivariable modeling ana- lyses that were carried out. a. Assume that screening was carried out one variable at a time using tests of significance for the relation- ship between each potential control variable and the outcome variable, so that those potential con- trol variables not found significantly associated with the outcome variable were not included in any modeling analysis. How can you criticize this approach to screening? b. Why was some kind of screening likely necessary for this analysis? Suppose that the logistic regression treated the variables mattress type (F) and perceived mattress type (PF) as ordi- nal variables. Suppose also that the variable BASE is also considered to be an exposure variable (even though it was not involved in the randomization) in addition to the vari- able mattress type (F). Suppose further that this model allows for two-way inter- actions (i.e., products of two variables) between mattress type (F) and each of the other independent variables (POST, BASE, PF, OCC, AGE, and GEN) and two-way interactions between BASE and each of the control vari- ables PF, POST, OCC, AGE, and GEN. 2. State the logit formula for the logistic regression model just described. Make sure to consider both F and BASE as exposure variables. 3. For each of the following product terms, state whether the product term is an EE variable, an EV variable, or neither: F 3 POST F 3 BASE POST 3 BASE PF 3 POST BASE 3 PF (To answer this part, simply write EE, EV, or neither next to the variable regardless of whether the product term given is contained in the revised model described in question 2.)
Test 295 4. Suppose that in carrying out interaction assessment for the model of question 2, a chunk test for all two-way product terms is significant. Suppose also that upon further interaction testing: a chunk test for two-way product terms involving F with POST, OCC, PF, AGE, and GEN has a P-value of 0.25 and a chunk test for two-way product terms involving BASE with POST, OCC, PF, AGE, and GEN has a P-value of 0.70. Which of the following choices are “reasonable” as the next step in the assessment of interaction? Circle as many choices as you consider reasonable. a. All product terms except F 3 BASE should be dropped from the model. b. The two-way product terms involving F with POST, OCC, PF, AGE, and GEN should be dropped from the model and backward elimination should be car- ried out involving the two-way product terms involving BASE with POST, OCC, PF, AGE, and GEN. c. The two-way product terms involving BASE with POST, OCC, PF, AGE, and GEN should be dropped from the model and backward elimination should be carried out involving the two-way product terms involving F with POST, OCC, PF, AGE, and GEN. d. Carry out a backward elimination of all two-way product terms. Note: In answering the above question, the word “reason- able” should be interpreted in the context of the hierarchi- cal backward elimination strategy described in this text. Also, recall that F and BASE are exposure variables. 5. Describe how you would test for the significance of the two-way product terms involving F with POST, OCC, PF, AGE, and GEN using the model in question 2. In answering this question, make sure to state the null hypothesis in terms of model parameters, describe the formula for the test statistic, and give the distribution and degrees of freedom of the test statistic under the null hypothesis. 6. Suppose that at the end of the interaction assessment stage, it was determined that the variables F 3 BASE, F 3 POST, and BASE 3 OCC need to remain in the model as significant interaction effects. Based on the hierarchical backward elimination strategy described in Chap. 7, what V variables are eligible to be dropped
296 8. Additional Modeling Strategy Issues from the model as possible nonconfounders? Briefly explain your answer. 7. Based on the interaction assessment results described in question 6, is it appropriate to test the significance for the main effect of POST and/or OCC? Explain briefly. 8. a. State the logit formula for the reduced model obtained from the interaction results described in question 6. b. Based on your answer to 8 a, give a formula for the odds ratio that compares the odds for improvement of pain both in bed and upon rising to no improve- ment for a subject getting a firm mattress (F ¼ 3) and a firm base (BASE ¼ 1) to the corresponding odds for a subject getting a medium firm (F ¼ 2) mattress and an infirm base (BASE ¼ 0), controlling for POST, PF, OCC, AGE, and GEN. 9. Assume that the odds ratio formula obtained in ques- tion 8 represents the gold standard odds ratio for describing the relationship of mattress type (F) and mattress base (BASE) to pain improvement controlling for POST, OCC, PF, AGE, and GEN, i.e., your only interest is the OR comparing (F ¼ 3, BASE ¼ 1) with (F ¼ 2, BASE ¼ 0). One way to assess confounding among the variables eligible to be dropped as noncon- founders is to compare tables of odds ratios for each subset of possible confounders to the gold standard odds ratio. a. How many subsets of possible confounders (other than the set of possible confounders in the gold standard odds ratio) need to be considered? b. Describe what a table of odds ratios would look like for any of the subsets of possible confounders, i.e., draw such a table and specify what quantities go into the cells of the table. c. How would you use the tables of odds ratios described above to decide about confounding? In your answer, describe any difficulties involved. d. Suppose you decided that PF and GEN could be dropped from the model as nonconfounders, i.e., your reduced model now contains F, BASE, POST, OCC, AGE, F 3 BASE, F 3 POST, and BASE 3 OCC: Describe how you would determine whether precision was gained when dropping PF and GEN from the model. 10. Suppose that after all of the above analyses described in the previous questions, you realized that you
Test 297 neglected to check for the possibility of collinearity in any of the models you considered so far. a. Assuming that all the models you previously con- sidered ran using the SAS’s LOGISTIC procedure (i.e., they all produced output without any error or warning messages), why should you still be concerned about the possibility of collinearity? Suppose, further, that you use SAS’s LOGISTIC procedure to fit a model containing the independent variables F, BASE, POST, OCC, PF, AGE, GEN, all two-way product terms involving F with the control variables (POST, OCC, PF, AGE, GEN), and all two-way product terms involving BASE with these same control variables. b. You now run the collinearity macro with the above logistic model, and you find that there are three condition indices with values 97, 75, and 62, with all other condition indices less than 25. How would you proceed “sequentially” to use this information to assess collinearity? c. Suppose the condition index of 97 has “high VDP values” on the variables F, BASE, PF, and F 3 BASE. Does this result cause you difficulty in accepting your previous assessment of interaction in which you found that the variables F 3 BASE, F 3 POST, and BASE 3 OCC needed to remain in the model as significant interaction effects? Explain. d. Based on the collinearity results in part 10 b and c, which of the following choices is an appropriate next step? (Circle the “best” choice) i. Drop the product term F 3 BASE from the polytomous logistic regression model and redo the hierarchical backward elimination strategy without further consideration of collinearity. ii. Determine whether the next highest condition index of 75 corresponds to high VDP loadings two or more predictors. iii. Drop the product term F 3 BASE from the logistic regression model, and apply collinearity diagnostics to the reduced model to determine if there is an additional collinearity problem. iv. Ignore the collinearity diagnostics results and use the model obtained from the hierarchical backward elimination strategy previously used. 11. a. Assuming that mattress type (F) and type of base (BASE) are the only two exposures of interest, with PF, GEN, POST, OCC, and AGE considered as control variables, briefly outline how you would
298 8. Additional Modeling Strategy Issues assess whether or not any subjects in the dataset are influential observations. Make sure to indicate whether you would prefer to use DeltaBeta mea- sures or Cook’s distance-type measures, both types of measures, or other measures. b. Why would it be questionable to automatically drop from your dataset any subjects that you find to be influential observations? 12. Suppose in your analysis strategy for determining a “best” model, you want to reduce the number of sta- tistical tests that you perform. What approach(es) can you use? Why are you not able to adequately carry out a Bonferroni-type of adjustment procedure? Answers to 1. a. True Practice b. False: E2E3 not needed. Exercises c. False: E1C2 also a candidate d. True e. False: Incorrect use of backward elimination. 2. a. OdR ¼ exp½b^1 þ b^2 þ b^3 þ d^11C1 þ d^22C2 b. C2 ¼ 1 C2 ¼ 0 C1 ¼ 1 OdR11 OdR10 C1 ¼ 0 OdR01 OdR00 OdRab ¼ exp½b^1 þ b^2 þ b^3 þ d^11ðC1Þ þ ^d22ðC2Þ c. i. 4 ii. {C3, C1C2} {C3} {C1C2} {Neither C3 nor C1C2} d. Models 16, 7, 4a e. Obtain tables of confidence intervals for the odds ratios for each of the three models stated in part d. Choose as the best model either: i. the model with the narrowest confidence interval ii. the gold standard model (16) if all four models have approximately the same width. iii. the more parsimonious model 3. i. Confounding for individual Cs not addressed by statistical testing. ii. Interaction of individual Cs with each Ei not addressed iii. Interaction of Ei with Ej not addressed iv. Screening of Cs were not distinguished from screening of Es.
Answers to Practice Exercises 299 4. a. False: Can’t tell until you check VDPs. Possible that all VDPs are not high (i.e., much less than 0.5) b. False: Model won’t be HWF if C3 is dropped. c. True d. False: May not be any remaining collinearity prob- lem once EC3 is dropped. 5. a. A Cook’s distance-type measure combines the infor- mation from all estimated regression coefficients in one’s model, whereas it would be preferable to con- sider either the Db or Dexp[b] for the E variable alone, since the E variable is the primary variable of interest. b. You should not automatically drop a subject from the dataset just because you have identified it as influential. A conservative approach is to drop only those subjects whose data are clearly in error and cannot be corrected. c. The model of question 4 may not be the best model, so that different conclusions might result about which subjects are influential if a different (“best”) model were used instead. 6. The number of tests to be performed cannot be deter- mined in advance of the modeling process, i.e., it is not clear what T will be for the individual significance level of 0.05/T.
9 Assessing Goodness of Fit for Logistic Regression n Contents Introduction 302 Abbreviated Outline 302 Objectives 303 342 Presentation 304 Detailed Outline 329 Practice Exercises 334 Test 338 Answers to Practice Exercises D.G. Kleinbaum and M. Klein, Logistic Regression, Statistics for Biology and Health, 301 DOI 10.1007/978-1-4419-1742-3_9, # Springer ScienceþBusiness Media, LLC 2010
302 9. Assessing Goodness of Fit for Logistic Regression Introduction Regression diagnostics are techniques for the detection and assessment of potential problems resulting from a Abbreviated fitted regression model that might either support, compro- Outline mise, or negate the assumptions made about the regres- sion model and/or the conclusions drawn from the analysis of one’s data. In this chapter, we focus on one important issue for evalu- ating binary logistic regression results, namely, goodness of fit (GOF) measurement. Although examination of data for potential problems, such as GOF, has always been con- sidered a requirement of the analysis, the availability of computer software to efficiently perform the complex cal- culations required has contributed greatly to fine-tuning the diagnostic procedures and the conclusions drawn from them. The outline below gives the user a preview of the material to be covered by the presentation. A detailed outline for review purposes follows the presentation. I. Overview (pages 304–305) II. Saturated vs. Fully Parameterized Models (pages 305–312) III. The Deviance Statistic (pages 312–317) IV. The HL Statistic (pages 318–320) V. Examples of the HL Statistic (pages 320–325) VI. Summary (page 326) VII. Appendix: Derivation of SS Deviance Formula (pages 327–328)
Objectives Objectives 303 Upon completing this chapter, the learner should be able to: 1. Explain briefly what is meant by goodness of fit. 2. Define “perfect prediction.” 3. Distinguish between “events–trials” format and “subject-specific” format for a dataset. 4. Define and illustrate the covariate patterns for a specific logistic model. 5. State or recognize the distinction between a fully parameterized and a saturated binary logistic model. 6. Given a specific binary logistic model, state or recognize the deviance formula for the model. 7. Explain briefly why the deviance statistic is not used for assessing goodness of fit when fitting a binary logistic regression model. 8. Given a printout of the results of a binary logistic regression: a. State or recognize the Hosmer–Lemeshow statistic b. Carry out a test of hypothesis for goodness of fit using the Hosmer–Lemeshow statistic
304 9. Assessing Goodness of Fit for Logistic Regression Presentation I. Overview This presentation describes methods for asses- sing the extent to which a logistic model esti- Focus Does estimated mated from a dataset predicts the observed logistic model predict outcomes in the dataset. The classical term for this topic is goodness of fit (GOF). observed outcomes in data? GOF is an issue that considers how well a given model, considered by itself, fits the data, rather Considers a given model than whether or not the model is more appro- priate than another model. Does not consider comparing models In most epidemiologic analyses, the primary goal is to assess an exposure–disease relation- Primary analysis goal: ship, so that we are usually more interested in Assess E–D relationship to derive deriving the “best” model for the relationship “best” model (which typically involves a strategy requiring the comparison of various models) than in GOF goal: using a GOF procedure. Nevertheless, once Determine how well final (“best”) we have obtained a final (i.e., best”) model, model fits the data we would also like this model to fit the data well, thus justifying a GOF procedure. Assume: Y is binary (0,1) Units of analysis: individual subjects GOF: Summary measure that com- Assuming that the outcome (Y) is binary, say pares coded as 0 or 1, and the unit of analysis is an individual subject, GOF typically requires a Yi to Y^i; where summary measure over all subjects that com- Yi ¼ observed response for pares the observed outcome (Yi) for subject i to the predicted outcome (Y^i) for this subject subject i obtained from the fitted model, i.e., Y^i ¼ P^ðXiÞ. Y^i ¼ P^ðXiÞ ¼ predicted response If when determined collectively over all sub- for subject i jects, the GOF measure is “small” or “nonsig- i ¼ 1, 2, . . . , n nificant,” we say the model has good fit, although, technically, we mean there is not Good fit: GOF measure “small” or “n.s.” sufficient evidence to conclude a bad fit. Otherwise, we say that the model has evidence Not sufficient evidence to conclude of lack of fit. a bad fit Lack of fit: Otherwise
Presentation: II. Saturated vs. Fully Parameterized Models 305 Widely used GOF measure: devi- A widely used GOF measure for many mathe- ance matical models is called the deviance. However, as we describe later, for a binary logistic But: deviance problematic for regression model, the use of the deviance for binary logistic regression assessing GOF is problematic. Popular alternative: A popular alternative is the Hosmer–Lemeshow Hosmer–Lemeshow (HL) (HL) statistic. Both the deviance and the HL statistic will be defined and illustrated in this statistic chapter. II. Saturated vs. Fully As stated briefly in the previous overview sec- Parameterized Models tion, a measure of goodness of fit (GOF) pro- vides an overall comparison between observed GOF: Overall comparison between (Yi) and predicted values ðY^iÞ of the outcome observed and predicted out- variable. comes We say there is perfect fit if Yi À Y^i ¼ 0 for all i. Perfect fit: Yi À Y^i ¼ 0 for all i. Although the mathematical characteristics of any logistic model require that the predicted Rarely happens value for any subject must lie between or Typically 0 < Y^i < 1 for most i including 0 or 1, it rarely happens that pre- dicted values are either 0 or 1 for all subjects in a given dataset. In fact, the predicted values for most, if not all, subjects will lie above 0 and below 1. Perfect fit: Not practical goal Thus, achieving “perfect fit” is typically not a Conceptual ideal practical goal, but rather is a conceptual ideal Saturated model when fitting a model to one’s data. Neverthe- (a reference point) less, since we typically want to use this “ideal” model as a reference point for assessing the fit EXAMPLE of any specific model of interest, it is conve- nient to identify such a model as a saturated model. SBP n=2 FOOT A trivial example of a saturated regression Subject # SBP model is obtained if we have a dataset contain- 1 115 9 170 •2 170 11 ing only n ¼ 2 subjects, as shown on the left. 115 • Here, the outcome variable, SBP, is continu- FOOT ous, as is the (nonsense) predictor variable foot length (FOOT). A “perfect” straight line fits the 9 11 data.
306 9. Assessing Goodness of Fit for Logistic Regression EXAMPLE (continued) The linear regression model here involves only two parameters, b0 and b1, whose estimates SB^ P ¼ b^0 þ b^1ðFOOTÞ; yield predicted values equal to the two where b^0 ¼ À132:5 and b^1 ¼ 27:5 observed values of 115 and 170. so SB^ P ¼ À132:5 þ 27:5ð9Þ ¼ 115 and SB^ P ¼ À132:5 þ 27:5ð11Þ ¼ 170 Thus, in this example, a saturated model is obtained when the number of model para- Linear model example: meters (k þ 1 ¼ 2) is equal to the number of k þ 1 (¼ n) ¼ 2, subjects in the dataset. (Note: k ¼ # of vari- ables in the model, and the “1” refers to the where k þ 1 ¼ # of parameters intercept parameter.) (including intercept) and n ¼ # of subjects More generally, the saturated model for a given dataset is defined as any model that con- Saturated Model (general): tains as many parameters as the number of kþ 1¼ n “observations” in the dataset, i.e., the sample size. where k þ 1 ¼ # of parameters (including intercept) n ¼ sample size EXAMPLE Observed Cohort Data To illustrate GOF assessment when using binary logistic regression, consider the follow- V=1 V=0 ing observed data from a cohort study on E=1 E=0 E=1 E=0 40 subjects. The outcome variable is called D, D = 1 3 7 10 there is one binary exposure variable (E), and D=1 6 4 10 there is one binary covariate (V). D=0 4 6 10 D = 0 7 3 10 These data indicate that V is an effect modifier of the E, D relationship, since the odds ratios of OR V=1 = 2.250 OR V=0 = 0.184 2.250 and 0.184 are very different and are on opposite sides of the null value of 1. >1 <1 Three models that may be fit to these data are Very different shown at the left. In model 1, E is the only ⇓ predictor. In model 2, both E and V are predic- tors. Model 3 includes the product term E Â V V is effect modifier of ORE,D in addition to both E and V main effect terms. Model 1: logit P(X) ¼ a þ bE Model 2: logit P(X) ¼ a þ bE þ gV Model 3: logit P(X) ¼ a þ bE þ gV þ dEV n ¼ 40 but k þ 1 ¼ 2, 3, or 4 Since the total sample size is 40, whereas each i.e., k þ 1 < n in all 3 models of these models contains 2, 3, and 4 para- i.e., no model is saturated meters, respectively, none of these three mod- els are saturated because k þ 1 < n for each (for predicting individual outcome) model.
Presentation: II. Saturated vs. Fully Parameterized Models 307 EXAMPLE (continued) Note, however, that concluding that “none of the 8 models are saturated” is based on the following assumption: the unit of analysis is the subject. >>>>>>>>>>>>>><>>>>>>>>>>>>>>>SKDMueaoUytbdanaejlisietlnsco2uet)fmsa(lpini)staitoleyndDs:ibsyisEsuthbejVescutb(jee.gct., This assumption is equivalent to listing the >>>>>>>>:>>>>>>>>>>>>>>>>>>>>>Ithf esunbM43j12...e09ocdtsela2reisunn1100... iotts 11 dataset using 40 lines of data, one for each 1... 1... subject. For Model 2, therefore, each dataline 00 would contain for a given subject, the values of 00 the outcome (D) and the predictor variables (E and V) in the model. of analysis, saturated When each subject is the unit of analysis, we can therefore claim that Model 2 is not (k þ 1 ¼ 3 < n ¼ 40) saturated because the number of parameters (3) in the model is less than total units in the 8 dataset (n ¼ 40). >>>>>>>>>>>>>>>>>>>>>>>>><>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>AMDGlartoUEetodvanrp(uneelsi3241aitlunapntotb2tteiesfj)(vser–geanctLn)rt)asiiassawsltldyues7463isgdmtfihsobpsrni1111tyasmi0000gomGaanertg:Eo0110rcuoopuvVa0101p(reia.gte., However, there is another way to view the dataset we have been considering: the unit of >>>>>>>>>>>>>>>>>>>:>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>“G((n(ErdMgGovgr¼oue)oofnapiodu#nltr-(epspsd^noln–oagifg2ptvft)c¼rmu)iroioderabdouvadlgssadtai=eerelcfnrdioltpvga”irpraotmperetneidraopodirtncvaia:ctsittidtthoi¼eeonerdnr)n4n:tsh¼ank þ 1 analysis is a group of subjects, all of whom have (i.e., perfectly predicts p^g) the same covariate pattern within a group. This assumption is equivalent to listing the dataset using only four lines of data, one for each covariate pattern. Each dataline would contain for a given group of subjects, the num- ber of cases (D ¼ 1) in each group (dg), the number of subjects in each group (ng), and the values of each predictor variable being modeled, as shown at the left for our dataset. This type of data layout is called an events–trials format, where there are dg events and ng trials. Using events–trials format, we can argue that the number of observations (n) consists of the total number of datalines (4 in our example), rather than the number of subjects (40). Here, the goal of model prediction no longer is to predict an individual’s (0 or 1) outcome, but rather to predict the observed proportion of persons in a group (the unit of analysis) that has the outcome, i.e., p^g ¼ dg=ng. Thus, using events–trials format, we can declare a model to be “group-saturated” if for each covar- iate pattern listed in the dataset, the model per- fectly predicts the observed proportion p^g.
308 9. Assessing Goodness of Fit for Logistic Regression EXAMPLE (continued) Let us now focus on Model 3. This model is the Model 3: logit P(X) ¼ a þ bE þ gV largest possible model that can be defined con- þ dEV taining the two “basic” variables, E and V. Largest possible model containing Since E and V are both (0,1) variables, E2 ¼ E binary E and binary V and V2 ¼ V, so we cannot add to the model any Note: E2 ¼ E and V2 ¼ V higher order polynomials in E or V or any product terms other than E Â V. Fully parameterized model: Model 3 is an example of a fully parameterized Contains maximum # of model, which contains the maximum number covariates defined from the of covariates that can be defined from the main-effect covariates, main-effect covariates in the model. No. of parameters (k þ 1) ¼ G Equivalently, the number of parameters in such covariate patterns, where a model must equal the number of covariate G ¼ # of covariate patterns patterns (G) that can be defined from the covari- ates in the model. Covariate patterns (i.e., subgroups): In general, for a given model with covariates Distinct specifications of X X ¼ (X1, . . . , Xk), the covariate patterns are EXAMPLE defined by the distinct values of X. Model 3: logit P(X) = α + βE + γV + δEV 4 covariate patterns k+ 1=4=G For Model 3, which contains four parameters, there are four distinct covariate patterns, i.e., X1: E = 1, V = 1 fully subgroups, that can be defined from the covari- X2: E = 0, V = 1 parameterized ates in the model. These are shown at the left. X3: E = 1, V = 0 X4: E = 0, V = 0 Model 1, which contains only binary E, is also fully parameterized, providing E is the only Model 1: logit P(X) = α + βE k+ 1=2=G basic predictor of interest. No other variables defined from E, e.g., E2 ¼ E, can be added to 2 covariate patterns fully model 1. Furthermore, Model 1 contains two parameterized parameters, which correspond to the two X1: E = 1 covariate patterns derived from E. X2: E = 0 Model 2: logit P(X) = a + bE + g V 4 covariate patterns k + 1 =3 ≠G=4 (same as Model 3) not fully parameterized However, Model 2 is not fully parameterized, since it contains three parameters and four covariate patterns. Assessing GOF: Thus, we see that a fully parameterized model has a nice property (i.e., k þ 1 ¼ G): it is the Saturated ðk þ 1 ¼ nÞ largest possible model we can fit using the vari- vs: ables we want to allow into the model. Such a model might alternatively be used to assess fully parameterized ðk þ 1 ¼ GÞ? GOF rather than using the saturated model as the (gold standard) referent point.
Presentation: II. Saturated vs. Fully Parameterized Models 309 Model A: X ¼ (X1, X2, X3) fully Note, however, even if a model A is fully para- parameterized but meterized, there may be a larger model B con- taining covariates not originally considered in Model B: X ¼ (X1, X2, X3, X4, X5, X6) model A that provides “better fit” than model A. “better fit” than Model A EXAMPLE e.g., 8 For instance, although Model 1 is the largest model that can be defined when only binary E ▸tLeRstÀÀÀÀ>>>>>><Mode¼l 1a:þlobgEit PðXÞ is considered, it is not the largest model that :>>>>>>Mode¼l 3a:vþlso:bgEitþPgðXV Þþ can be defined when binary V is also consid- dEV ered. Nevertheless, we can choose between Model 1 and Model 3 by performing a standard LR ¼ À2 ln L^Model 1 À ðÀ2 ln L^Model 3Þ likelihood ratio (LR) test that compares the ¼ 55:051 À 51:355 ¼ 3:696 two models; the (nonsignificant) LR results (P > 0.10) are shown at the left. $ w2ð2 dfÞ under H0: g ¼ d ¼ 0 ðP > 0:10Þ n:s: Fitted Model 3: Focusing now on Model 3, we show the fitted logit P^ðXÞ ¼ a^ þ b^E þ ^gV þ ^dEV; model at the left. Note that these same esti- where ^a ¼ 0:8473; b^ ¼ À1:6946; mated parameters would be obtained whether we input the data by individual subjects (40 ^g ¼ À1:2528; ^d ¼ 2:5055 datalines) or by using an events–trials format (4 datalines). However, using events–trials format, we lose the ability to identify which subjects become cases. Covariate Obs. Pred. The predicted risks obtained from the fitted pattern risk risk model for each covariate pattern are also shown at the left, together with their cor- X1: E ¼ 1, p^1 ¼ 0:6 P^ðX1Þ ¼ 0:6 responding observed risks. Notice that these V¼1 p^2 ¼ 0:4 predicted risks are equal to their corresponding p^3 ¼ 0:3 P^ðX2Þ ¼ 0:4 observed proportions computed from the X2: E ¼ 0, p^4 ¼ 0:7 observed (stratified) data. V¼1 P^ðX3Þ ¼ 0:3 X3: E ¼ 1, P^ðX4Þ ¼ 0:7 V¼0 X4: E ¼ 0, V¼0 E = 1, V = 1: Pˆ (X1) = 0.6 ≠ 0 or 1 Nevertheless, none of these predicted risks are E = 0, V = 1: Pˆ (X2) = 0.4 ≠ 0 or 1 either 0 or 1, so the fitted model does not per- fectly predict each subject’s observed outcome, Model 3: – E = 1: some D = 1, which is either 0 or 1. This is not surprising, some D = 0 since some exposed subjects develop the dis- No perfect fit E = 0: some D = 1, ease and some do not. some D = 0
310 9. Assessing Goodness of Fit for Logistic Regression Model 3. Fully parameterized model However, we will show below that because + Model 3 is a fully parameterized model, it per- fectly predicts the number of cases actually # of expected cases observed for each covariate pattern. That is, ¼ # of observed cases for the expected number of cases for each pattern, each covariate pattern based on the fitted model, equals the observed number of cases for each pattern. Covariate pattern Xg: ng subjects More specifically, suppose ng subjects have Xg ¼ values of X in group g covariate pattern Xg, and dg denotes the observed number of cases in group g. Then, dg ¼ observed cases in group g (binomial) since dg has the binomial distribution, the expected number of cases in group g is P^ðXgÞ: predicted risk in group g d^g ¼ ngP^ðXgÞ, where P^ðXgÞ is the predicted risk d^g ¼ ngP^ðXgÞ: expected cases in for any subject in that group. group g EXAMPLE X: EV Exp. Cases Obs. Cases Thus, for Model 3, we expect d^1 ¼ 10ð0:6Þ ¼ 6 cases among subjects with E ¼ 1 and V ¼ 1, X1:11 d^1 ¼ 10ð0:6Þ ¼ 6 d1 ¼ 6 d^2 ¼ 10ð0:4Þ ¼ 4 cases among subjects with X2:01 d^2 ¼ 10ð0:4Þ ¼ 4 d2 ¼ 4 E ¼ 0 and V ¼ 1, and so on for the other X3:10 d^3 ¼ 10ð0:3Þ ¼ 3 d3 ¼ 3 two covariate patterns. The corresponding X4:00 d^4 ¼ 10ð0:7Þ ¼ 7 d4 ¼ 7 observed number of subjects are also shown at the left. Notice that the corresponding observed Model 3. Perfect and expected cases are equal for Model 3. Yi À Y^i ¼6 0 for Prediction? So, even though Model 3 does not provide “per- all i (not fect prediction” in terms of individual out- Individuals: No comes, it does provide “perfect prediction” in saturated) terms of group outcomes. dg À d^g ¼ 0 for Groups: Yes In other words, although Yi À Y^i ¼6 0 for all sub- all g (fully (Patterns) jects, dg À d^g ¼ 0 for all covariate patterns. parameterized) Model 3 is “group-saturated”: Another way of saying this is that Model 3 is perfectly group outcomes “group-saturated” in the sense that Model 3 perfectly predicts the group outcomes corres- Two GOF approaches: ponding to the distinct covariate patterns. Compare fitted model to: 1. Saturated model: Provides Thus, we see that an alternative gold standard perfect individual prediction model for assessing GOF is a fully parameter- 2. Fully parameterized model: ized (group-saturated) model containing the Provides perfect group covariates of interest rather than a (subject- prediction (based on specific) saturated model that can rarely if covariate patterns) ever be achieved using these covariates.
Presentation: II. Saturated vs. Fully Parameterized Models 311 Classical GOF approach: The traditional (i.e., “classical”) GOF approach Saturated model gives perfect considers the saturated model as the ideal for fit for individual “perfect fit.” This makes sense when the units subjects of analysis are individual subjects, since their actual observed outcomes are 0 or 1, rather Why? than some value in between. However, as we will explain further (later below), use of the Yi ¼ 0 or 1 only possible saturated model to assess GOF for logistic outcomes for regression is problematic. subject i However, problematic for logistic regression Saturated model: k þ 1 ¼ n Recall that we originally defined the saturated model as that model for which the number of EXAMPLE parameters (k þ 1) equals the sample size (n). Previous example (n ¼ 40, 4 covariate For our example involving four covariate pat- patterns): terns for 40 subjects, the subject-specific (SS) Model 4 (SS saturated model) saturated model is shown at the left. This model does not have an intercept term, but logit PðXÞ ¼ o1Z1 þ o2Z2 þ o3Z3 does contain 40 parameters, as defined by the ( þ Á Á Á þ o40Z40 oi. The Zi are dummy variables that distinguish 1 if subject i; i ¼ 1; 2; . . . ; 40 the 40 subjects. Zi ¼ 0 otherwise Y40 The likelihood function for this (SS) saturated LSS ¼ PðXiÞYi ð1 À PðXiÞÞ1ÀYi model is shown at the left. In this formula, Yi denotes the observed value (either 0 or 1) for i¼1 the ith individual in the dataset. where Xi denotes the values of X for Note that for subject i, Zi ¼ 1 and Zk ¼ 0 for subject i k 6¼ i, so P(Xi) can be written in terms of the regression coefficient oi that involves only that Subject i : Zi ¼ 1; other Zs ¼ 0 one subject. + logit PðXiÞ ¼ oi and PðXiÞ ¼ 1=½1 þ expðÀoiÞ Saturated model Furthermore, since the saturated model per- + fectly fits the data, it follows that the maximum likelihood (ML) estimate Y^i; which equals P^ðXiÞ Y^i def P^ðXiÞ ¼ Yi; i ¼ 1; 2; . . . ; n by definition, must be equal to the observed Yi for each subject i.
312 9. Assessing Goodness of Fit for Logistic Regression Y40 It follows that the formula for the ML value of L^SS max ¼ YiYi ð1 À YiÞ1ÀYi the SS saturated model (L^SS max) involves sub- i¼1 stituting Yi for P(Xi) in the above formula for For any binary logistic model: the likelihood, as shown at the left. Yi ¼ 1 : YiYi ð1 À YiÞ1ÀYi ¼ 11ð1 À 1Þ1À1 ¼ 1 From simple algebra, it also follows that the Yi ¼ 0 : YiYi ð1 À YiÞ1ÀYi ¼ 00ð1 À 0Þ1À0 ¼ 1 expression YiYi ð1 À YiÞ1ÀYi will always be equal to one when Yi is a (0,1) variable. L^SS max 1 always Consequently, the maximized likelihood will Important implications for GOF always equal 1. This result has important implications when attempting to assess GOF using a saturated model as the gold standard for comparison with one’s current model. III. The Deviance Statistic As mentioned at the beginning of this chapter, a widely used measure of GOF is the deviance. Deviance: The general formula for the deviance (for any regression model) is shown at the left. In this Dev(βˆ ) = −2ln(Lˆ c / Lˆ max) formula, b^ denotes the collection of estimated regression coefficients in the current model b^ ¼ ðb^0; b^1; b^2; . . . ; b^pÞ being evaluated, L^c denotes the maximized likelihood for the current model, and L^max L^c ¼ ML for current model denotes the maximized likelihood for the saturated model. L^max ¼ ML for saturated model (Note: If subjects are the unit of Thus, the deviance contrasts the likelihood of the current model with the likelihood of the analysis, model that perfectly predicts the observed out- L^max L^SS max) comes. The closer are these two likelihoods, the better the fit (and the smaller the deviance). Lˆ c closer to Lˆ max ⇓ better fit (smaller deviance) perfect fit Lˆ c= Lˆ max ⇒ –2 ln( Lˆ c / Lˆ max) In particular, if L^c ¼ L^max, then the deviance is = –2 ln(1) = 0 0, its minimum value. In contrast, if L^c is much smaller than L^max, then the ratio L^c=L^max is a Lˆ c << Lˆ max ⇒ Lˆ c / Lˆ max small fraction small fraction, so that the logarithm of the ratio is a large negative number and À2 times poor fit ⇒ ln(Lˆ c / Lˆ max) this large negative number will be a large posi- large,negative tive number. ⇒ –2 ln(Lˆ c / Lˆ max) large,positive
Presentation: III. The Deviance Statistic 313 Properties of deviance These properties of the deviance, i.e., its values similar to range from zero to larger and larger positive numbers, correspond to the properties of a chi- properties of w2 statistic square statistic as used in a likelihood ratio test. Common to test for GOF In fact, when using the deviance to test for by comparing deviance GOF, it is common, though not strictly legiti- with w2nÀkÀ1 value: mate, to compare the deviance to chi-square values with n À k À 1 degrees of freedom when ðquestionably legitimateÞ the current model contains k þ 1 parameters. GOF H0: model fits The GOF null hypothesis is that the “model HA: model does not fit fits,” and the alternative hypothesis is that the “model does not fit.” Thus, if the deviance is Deviance “significant” ) poor fit “significantly” large, the model is considered to have poor fit. Devðb^Þ ¼ À2 lnðL^c=L^maxÞ Note that the deviance statistic is, by definition, a ¼ À2 ln L^c À ðÀ2 ln L^maxÞ likelihood ratio (LR) statistic for comparing one’s current model to the saturated model. Thus, the Recall use of the chi-square distribution to test for the LR ¼ À2 ln L^R À ðÀ2 ln L^FÞ $ w2 significance of the deviance appears justified R ¼ reduced model because the LR statistic has an approximate F ¼ full model chi-square distribution under H0 when compar- ing full vs. reduced (nonsaturated) models that are fitted using ML estimation. DevRðb^Þ À DevFðb^Þ In particular, it can be shown from simple ¼ ½À2 lnðL^R=L^maxÞ À ½À2 lnðL^F=L^maxÞ algebra that the LR test for comparing two hier- = [–2 ln Lˆ R – (–2 ln Lˆ max)] – [–2 ln Lˆ F – (– 2 ln ˆLmax)] archical nonsaturated regression models is equivalent to the difference in deviances ¼ À2 ln L^R À ðÀ2 ln L^FÞ LR between the two models. This follows, as shown at the left, because the maximized like- lihood for the saturated model drops out of the difference in deviance scores. Nevertheless for the logistic model: Nevertheless, as we shall soon explain, when w2 approximation of using the deviance statistic to assess GOF for a single (nonsaturated) logistic model, the deviance statistic is chi-square approximation for the LR test is questionable (see below) questionable.
314 9. Assessing Goodness of Fit for Logistic Regression Alternative (Events–Trials) Devi- First, we present an alternative formula for the ance formula: for the logistic deviance in a logistic model that considers the co- model: variate patterns defined by one’s current model. G covariate patterns We assume that this model contains G covariate Xg ¼ ðXg1; Xg2; . . . ; XgpÞ; g ¼ 1; 2; . . . ; G patterns Xg ¼ (Xg1, Xg2, . . . , Xgp), with ng subjects d^g ¼ ngP^ðXgÞ ¼ expected cases having pattern g. As defined earlier, d^g ¼ ngP^ðXgÞ denotes the expected cases, where P^ðXgÞ is dg ¼ observed cases the predicted risk for X ¼ Xg, and dg denotes the observed number of cases in subgroup g. DevET ðb^ Þ The alternative deviance formula is shown here at the left. This formula corresponds to the data- ¼ À2 ln L^\"c;ET À ðÀ2 ln L^max;ETÞ set listed in events–trials format, where there ! are G datalines, dg and ng denote the number G of events and number of trials, respectively, on dg ln dg ¼ À2 ~ d^g !# g¼1 þ ðng À dgÞ ln ng À dg ; the gth dataline. ng À d^g where À2ln L^c,ET and À2ln L^max,ET are defined using events–trials (ET) format EXAMPLE Model 3: logit P(X) ¼ a þ bE þ gV Recall that for the data set on n ¼ 40 subjects þ dEV described above, Model 3 has G ¼ 4 covariate patterns. For each pattern, the corresponding Exp. Obs. values for ng, d^g, and dg are shown at the left. X: EV ng Cases Cases X1 : 11 10 d^1 ¼ 6 d1 ¼ 6 X2 : 01 10 d^2 ¼ 4 d2 ¼ 4 X3 : 10 10 d^3 ¼ 3 d3 ¼ 3 X4 : 00 10 d^4 ¼ 7 d4 ¼ 7 DevET ðb^ Þ for Model 3: ! Substituting the values in the table into the alternative (events–trials) deviance formula, ¼ À2 6 ln 6 þ 4 ln 4 we find that the resulting deviance equals zero. 6 4 ! À 2 4 ln 4 þ 6 ln 6 How can we explain this result, since we know 4 6! that the fully parameterized Model 3 is not À 2 3 ln 3 þ 7 ln 7 saturated in terms of perfectly predicting the 3 7! 0 or 1 outcome for each of the 40 subjects in the À 2 7 ln 7 þ 3 ln 3 dataset? The answer is that since the ET devi- 73 ance formula corresponds to an events–trials format, the units of analysis being considered ¼0 by the formula are the four groups of covariate patterns rather than the 40 subjects. Deviance formula DevETðb^Þ uses events–trials format + Units of analysis are groups ðnot subjectsÞ
EXAMPLE Presentation: III. The Deviance Statistic 315 “group-saturated”: Consequently, Model 3 is “group-saturated” in that the number of parameters in this model is # of parameters = # of groups equal to the total number of groups being con- sidered. That is, since the units of analysis are sample size = 4 the four groups, the sample size corresponding to the alternative deviance formula is 4, rather than 40. Model 3: logit P(X) ¼ a þ bE þ gV So, then, how can we compare Model 3 to the þ dEV saturated model we previously defined (Model 4) that perfectly predicts each subject’s 0 or 1 vs. outcome? Model 4 (perfectly predicts 0 or 1 outcome): logit PðXÞ ¼ o1Z1 þ o2Z2 þ o3Z3 þ Á Á Á þ o40Z40 Subject-specific deviance formula: n This requires a second alternative formula for Yi the deviance, as shown at the left. Here, DevSS ðb^ Þ ¼ À2 ~ Yi ln Y^i the summation (i) covers all subjects, not all groups, and Yi and Y^i denote the observed and i¼1 ! predicted response for the ith subject rather than the gth covariate pattern/group. þ ð1 À YiÞ ln 1 À Yi 1 À Y^i This alternative (subject-specific) formula is not identical to the events–trials formula Yi ¼ observed (0, 1) response for given above unless G ¼ n. Moreover, since G << n for Model 3, the SS deviance will be subject i nonzero for this model. Y^i ¼ predicted probability for the subject i DevSSðb^Þ ¼6 DevETðb^Þ unless G ¼ n DevSSðb^Þ > 0 since G << n for Model 3 Equivalent formula for DevSS(b^): We now show how to compute the subject- (from calculus and algebra) specific (SS) deviance formula for Model 3. However, we first provide an equivalent SS DevSSðb^Þ \" P^ðXiÞ ! deviance formula that can be derived from n À P^ðXiÞ both algebra and some calculus (see this P^ðXiÞ ln 1 chapter’s appendix for a proof ). ¼ À2 ~ # i¼1 þ lnð1 À P^ðXiÞÞ
316 9. Assessing Goodness of Fit for Logistic Regression EXAMPLE Obs. Pred. To compute this formula for Model 3, we need risk risk to provide values of P^ðXiÞ for each of the 40 Covariate subjects in the data set. Nevertheless, this cal- pattern p^1 ¼ 0:6 P^ðX1Þ ¼ 0:6 culation can be simplified since there are only four distinct values of P^ðXiÞ over all 40 sub- X1: E ¼ 1, P^2 ¼ 0:4 P^ðX2Þ ¼ 0:4 jects. These correspond to the four covariate V¼1 patterns of Model 3. p^3 ¼ 0:3 P^ðX3Þ ¼ 0:3 X2: E ¼ 0, The calculation now shown at the left, where V¼1 p^4 ¼ 0:7 P^ðX4Þ ¼ 0:7 we have substituted each of the four distinct values of P^ðXiÞ 10 times in the above formula. X3: E ¼ 1, V¼0 We have thus seen that the events–trial and subject-specific deviance values obtained for X4: E ¼ 0, Model 3 are numerically quite different. V¼0 The reason why DevETðb^Þ is zero is because the DevSS ðb^ Þ ET formula assumes that the (group-) saturated ¼ À 2ð10Þ½0:6 lnð0:6=0:4Þ þ lnð0:4Þ model is the fully parameterized Model 3 and þ 0:4 lnð0:4=0:6Þ þ lnð0:6Þ the current model being considered is also þ 0:3 lnð0:3=0:7Þ þ lnð0:7Þ Model 3. So the values of their corresponding þ 0:7 lnð0:7=0:3Þ þ lnð0:3Þ log likelihood statistics (À 2 ln L^) are equal and ¼ 51:3552 their difference is zero. DevETðb^Þ ¼ 0:0 ¼6 DevSSðb^Þ ¼ 51:3552 In contrast, the reason why DevSSðb^Þ is differ- ent from zero is because the SS formula DevETðb^Þ ¼ 0:0 because assumes that the saturated model is the “clas- À 2 ln L^ET saturated ¼ À2 ln L^Model 3 sical” (SS) saturated model that perfectly pre- dicts each subject’s 0 or 1 outcome. As ¼ À2 ln L^C mentioned earlier, À 2 ln L^ is always zero for the SS saturated model. Thus, DevSSðb^Þ simpli- so fies to À 2 ln L^C for the current model (i.e., Model 3), whose value is 51.3552. DevETðb^Þ ¼ À 2 ln L^C À ðÀ2 ln L^ET saturatedÞ Mathematically, the formula for À 2 ln L^C differs with the (ET or SS) format being used to specify ¼ À 2 ln L^Model 3 the data. However, these formulae differ by a À ðÀ2 ln L^Model 3Þ constant K, as we show on the left and illustrate for Model 3. The formula for K, however, does ¼ 0:0 not involve b. Thus, the ML estimate b^ will DevSSðb^Þ ¼6 0:0 because be the same for either format. Consequently, À 2 ln L^SS saturated ¼ 0:0 some computer packages (e.g., SAS) present the same value (i.e., À 2 ln L^C;SS) regardless of so the data layout used. DevSSðb^Þ ¼ À 2 ln L^C À ðÀ2 ln L^SS saturatedÞ ¼ À 2 ln L^Model 3 À 0 ¼ 51:3552 Note: À 2 ln L^C;ET ¼6 À2 ln L^C;SS –2 ln Lˆ C,ET = –2 ln Lˆ C,SS −2K , where K = ln G ng! K does not involve β, so ∑ βˆ is same for SS or ET g = 1 dg!(ng – dg)! Model 3: –2 ln Lˆ C,ET = 10.8168 À 2 ln L^C;SS ¼ 51:3552; K ¼ 20:2692
Presentation: III. The Deviance Statistic 317 Deviance not always appropriate We are now ready to discuss why the use of the for logistic regression GOF deviance formula is not always appropriate for assessing GOF for a logistic regression model, Alternative approach: and we will describe an alternative approach, Hosmer–Lemeshow using the Hosmer–Lemeshow statistic, which is typically used instead of the deviance. When G<<n, we can assume DevETðb^Þ is approximately w2nÀkÀ1 When the number of covariate patterns (G) is considerably smaller than the number of obser- under H0: good fit vations (n), the ET deviance formula can be assumed to have an approximate chi-square dis- EXAMPLE tribution with n À k À 1 degrees of freedom. Previous data: n ¼ 40 For the data we have illustrated above involv- G ¼ 2; 4; 4 for Models 1; 2; 3; ing n ¼ 40 subjects, G ¼ 2 for Model 1, and respectively G ¼ 4 for Models 2 and 3. So a chi-square test + for GOF is appropriate using the ET deviance formula. w2 test for GOF is OK However, when G % n, we cannot However, when G is almost as large as n, in assume particular, when G equals n, then the deviance cannot be assumed to have a chi-square distri- DevETðb^Þ is approximately w2nÀpÀ1 bution (Collett, 1991). This follows from large- under H0: good fit sample statistical theory, where the primary problem is that in this situation, the number (Statistical theory: ng small % 1 as of subjects, ng, for each covariate pattern n ! 1) remains small, e.g., close to 1, as the sample size increases. Xi continuous, e.g., Xi ¼ AGE ) G % n Note that if at least one of the variables in the Many situations where predictors model, e.g., AGE, is continuous, then G will are continuous tend to be close to n whenever the age range in the sample is reasonably wide. Since logistic + regression models typically allow continuous Cannot use deviance to test for GOF variables, there are many situations in which the chi-square distribution cannot be assumed when using the deviance to test for GOF. G ¼ n: When G ¼ n, each covariate pattern involves Each covariate pattern: only one subject, and the SS deviance formula 1 subject, is equivalent to the ET deviance formula, which nevertheless cannot be assumed to ng 1 for all g, g ¼ 1, . . . , n have a chi-square distribution under H0. DevSSðb^Þ ¼ DevETðb^Þ Moreover, the SS deviance formula, shown but not w2 under H0: again here contains only the predicted values P^ðXiÞ for each subject. Thus, this formula tells GOF adequate nothing about the agreement between observed \"! (0,1) outcomes and their corresponding pre- n P^ðXiÞ dicted probabilities. DevSS ðb^ Þ ¼ À2 P^ðXi Þ ln À P^ðXi Þ ~ 1 i¼1 # þ ln 1 À P^ðXiÞ Provides P^ðXiÞ but not observed Yi
318 9. Assessing Goodness of Fit for Logistic Regression IV. The Hosmer– Lemeshow (HL) Statistic Alternative to questionable use To avoid questionable use of the deviance to pro- of deviance vide a significance test for assessing GOF, the Hosmer–Lemeshow (HL) statistic has been deve- Available in most computer loped and is available in most computer packages. packages HL widely used regardless of The HL statistic is widely used regardless of whether G << n or G % n: whether or not the number of covariate patterns (G) is close to the number of observations. Nev- Requires G > 3 ertheless, this statistic requires that the model Rarely significant when G < 6 considers at least three covariate patterns, rarely Works best when G % n (e.g., results in significance when G is less than 6, and works best when G is close to n (the latter occurs some Xs are continuous) when some of the predictors are continuous). HL 0 for fully parameterized Moreover, the HL statistic has the property model that it will always be zero for a fully parame- terized model. In other words, the “saturated” “Saturated” model is fully model for the HL statistic is essentially a fully parameterized in ET format parameterized (group-saturated) model coded in events–trials format. Steps for computing HL statistic: The steps involved in computing the HL statis- 1. Compute P^ðXiÞ for all n subjects tic are summarized at the left. Each step will 2. Order P^ðXiÞ from largest to then be described and illustrated below, fol- lowing which we will show examples obtained smallest values from different models with differing numbers 3. Divide ordered values into Q of covariate patterns. percentile groupings (usually Q ¼ 10 ) deciles of risk) 4. Form table of observed and expected counts 5. Calculate HL statistic from table 6. Compare computed HL to w2 with Q À 2 df Step 1: Compute P^ðXiÞ; i ¼ 1; 2; . . . ; n At the first step, we compute the predicted risks n ¼ 200 ) 200 values for P^ðXiÞ P^ðXiÞ for all subjects in the dataset. If there are, although some values are identi- say, n ¼ 200 subjects, there will be 200 pre- dicted risks, although some predicted risks cal if Xi Xj for subjects i and j. will be identical for those subjects with the same covariate pattern. Step 2: Order values of P^ðXiÞ: e.g., n ¼ 200 At the second step, we order the predicted risks from largest to smallest (or smallest to largest). Order # P^ðXiÞ Again, there may be several ties when doing 1 0.934o(largest) this if some subjects have the same covariate pattern. 2 0.901 tie 3... 0.9... 01 199 0.123 200 0.045 (smallest)
Presentation: IV. The Hosmer–Lemeshow (HL) Statistic 319 Step 3: Form Q percentile group- At the third step, we divide the ordered pre- ings. dicted risks into Q percentile groupings. The typical grouping procedure involves Q ¼ 10 Typically, Q ¼ 10, i.e., deciles of deciles. Thus, if the sample size is 200, each risk e.g., n ¼ 200 decile will contain approximately 20 subjects. Henceforth, we will assume that Q ¼ 10. Decile No. of subjects 1 %20 ..2. %2...0 9 %20 10 %20 Total 200 Ties ) # of subjects ¼6 exactly Note, however, because some subjects may 20 (¼ n/Q) in all deciles have identical predicted risks (i.e., ties), the number of subjects per decile may vary some- Must keep subjects with identical what to keep subjects with identical predicted values of P^ðXiÞ in the same decile risks in the same decile. Step 4: Form table of observed At the fourth step, we form (typically using and expected cases and a convenient computer program) the table, noncases shown at the left, that contains observed and expected cases and noncases within each dec- Deciles Obs. Exp. Obs. non Exp. non ile. In this table, the values Ocq, Ecq, Oncq, and of risk cases cases cases cases Encq, q ¼ 1, 2, . . . , 10 are defined as follows: Oc1 Ec1 Onc1 Enc1 Ocq ¼ # of observed cases in the qth decile 1 Oc2 Ec2 Onc2 Enc2 2 Ecq ¼ # of expected cases in the qth decile Oc3 Ec3 Onc3 Enc3 3 Oncq ¼ # of observed noncases in the qth decile 10 Oc10 Ec10 Onc10 Enc10 Encq ¼ # of expected noncases in the qth decile Observed cases and noncases: The observed cases (Ocq) and noncases (Oncq) in each decile are obtained by simply counting Ocq counts # of cases (Yi ¼ 1) in qth decile the numbers of subjects in that decile who are Oncq counts # of noncases (Yi ¼ 1) in qth decile cases (i.e., Yi ¼ 1) and noncases (i.e., Yi ¼ 0), respectively. Note that once we count Ocq, we Note: Oncq ¼ nq À Ocq can obtain Oncq by subtraction from nq, the total number of subjects in the qth decile. Expected cases and noncases: The expected cases (Ecq) in each decile are obtained by summing the predicted risks nq P^ðXiÞ for all subjects in that decile. The expected number of noncases (Encq) are Ecq ¼ ~ P^ðXiqÞ and Encq ¼ nq À Ecq, obtained by subtraction from nq. i¼1 where Xiq ¼ covariate values for ith subj in qth decile
320 9. Assessing Goodness of Fit for Logistic Regression EXAMPLE e.g., q ¼ 3, n3 ¼ 4, P^ðXiÞ: 0.30, 0.35, For example, if the third decile contains four 0.40, 0.45 subjects with predicted risks of 0.30, 0.35, 0.40, and 0.45, then the expected number of cases Ec3 ¼ n3 P^ðXi3Þ ¼ 0:30 þ 0:35 þ 0:40 (Ec3) would be their sum 0.30 þ 0.35 þ 0.40 þ 0.45 ¼ 1.50 (regardless of whether or not a ~ subject is an observed case). The expected noncases in the same decile (Enc3) would be i¼1 4 À 1.50 ¼ 2.50. þ 0:45 ¼ 1:50 and Enc3 ¼ n3 À Ec3 ¼ 4 À 1.50 ¼ 2.50 Step 5: In step 5, the HL statistic is calculated using the formula at the left. This formula involves HL ¼ Q ðOcq À Ecq Þ2 þ Q ðOncq À EncqÞ2 summing Q values of the general form (Oq À Ecq Encq Eq)2/Eq for cases and another Q values for non- ~ ~ cases. When Q ¼ 10, the HL statistic therefore involves 20 values in the summation. q¼1 q¼1 Q ¼ 10 ) 20 values in summation Step 6: In step 6, the HL statistic is tested for signifi- HL ~ approx c2Q–2 under H0: Good fit cance by comparing the computed HL value to a percentage point of w2 with Q À 2 degrees of Q = 10 ⇒ df = Q – 2 = 8 i.e., not enough evidence to freedom. When Q ¼ 10, therefore, the HL sta- tistic is approximately w2 with 8 df. indicate lack of fit V. Examples of the HL We now illustrate the use of the HL statistic Statistic with the Evans County data (n ¼ 609). This dataset has been considered in previous chap- EXAMPLE ters, and is described in detail in the Computer Appendix. SASs Logistic procedure was used Evans County Data (n ¼ 609) for the computations. (see previous chapters and Computer Appendix) In our first illustration, we fit the two models shown at the left. The outcome variable is CHD Model EC1 (no interaction): status (1 ¼ case, 0 ¼ noncase), and there are logit P(X) ¼ a þ bCAT þ g1AGEG three basic (i.e., main effect) binary predictors, CAT (1 ¼ high, 0 ¼ low), AGEG (1 ¼ age 55, þ g2ECG 0 ¼ age > 55), and ECG (1 ¼ abnormal, 0 ¼ normal). Recall that the Evan County dataset Model EC2 (fully parameterized): is described in the Computer Appendix. logit PðXÞ ¼ a þ bCAT þ g1AGEG þ g2ECG þ g1AGEG Â ECG þ d1CAT Â AGE þ d2CAT Â ECG þ d3CAT Â AGE Â ECG
Presentation: V. Examples of the HL Statistic 321 EXAMPLE (continued) The data layout used to fit both models in events–trials format is now shown at the left. Datalines in events trials format Model EC1 is a no-interaction model involving (n 5 609) only the main effects CAT, AGEG, and ECG as Group predictors. Model EC2 is a fully parameterized model since there are eight model parameters (g) dg ng CAT AGEG ECG as well as eight covariate patterns, i.e., p þ 1 ¼ 8 ¼ G. 1 17 274 0 0 0 2 15 122 0 1 0 3 7 59 0 0 1 4 5 32 0 1 1 5 1 81 0 0 6 9 39 1 1 0 7 3 17 1 0 1 8 14 58 1 1 1 proc logistic data ¼ evans2; Here, we provide the computer code using model cases/total ¼ cat ageg ecg/ scale ¼ none aggregate ¼ (cat SASs PROC LOGISTIC used to fit Model EC1 and provide HL, deviance, and À 2 ln L^C ageg ecg) lackfit; statistics, as well as predicted risk values (i.e., output out ¼ pred p ¼ phat predprob ¼ (individual); run; “phat” in the code at the left) for each covariate proc print data ¼ pred; run; pattern. EXAMPLE Edited output for Model EC1 is shown here. Edited Output (Model EC1): The table of observed and expected cases (Variables – CAT, AGE, ECG) and noncases has divided the data into Q ¼ 6 percentile groups rather than 10 deciles. The 1 274 17 18.66 257 255.34 number of covariate patterns is G ¼ 8, so the number of percentile groups allowable is less 2 59 7 5.60 52 53.40 than 10. Also, from the phat values provided (below left), two pairs {0.11913, 0.11984} and 3 122 15 14.53 107 107.47 {0.16264, 0.16355} are essentially identical and should not be separated into different groups. 4 57 9 8.94 48 48.06 Since Q ¼ 6, the df for the HL test is Q À 2 ¼ 4. 5 39 9 7.85 30 31.15 The HL test statistic (0.9474) is not significant. Thus, there is not enough evidence to indicate 6 58 14 15.41 44 42.59 that Model EC1 has lack of fit. No evidence 0.9474 4 0.9177 The output here also gives the Deviance that Model EC1 (0.9544) and Pearson (0.9793) chi-square sta- has lack of fit 0.9544 4 0.9166 tistics as well as the log likelihood statistic 0.9793 4 0.9129 (418.181) for Model EC1. – 2 Log L 418.181 1 0.06810 0.06204 2 0.11913 0.12295 3 0.09497 0.11864 4 0.16264 0.15625 5 0.11984 0.12500 6 0.20128 0.23077 7 0.16355 0.17647 8 0.26574 0.24138
322 9. Assessing Goodness of Fit for Logistic Regression EXAMPLE (continued) The Pearson statistic is another GOF statistic that is similar to the Deviance in that it is not Since G ¼ 8 << n ¼ 609, recommended when the number of covariate Pearson statistic and Dev statistic patterns is close to the sample size. However, approx w2 under H0 since the number of covariate patterns (G ¼ 8) for Model EC1 is much less than the sample size (n ¼ 609), both statistics can be assumed to be approximately chi square under H0. Notice also that the Pearson and Deviance values are very close to the HL value (0.9474). À2 ln L^C;SS ¼ 418:181 The À2 log L value (418.181) in the output is ¼ DevSSðb^Þ for Model EC1 the SS statistic À 2 ln L^C;SS, where C is Model EC1. This statistic is equivalent to DevSSðb^Þ for Model EC1, since L^max;SS is always one. DevET(βˆ ) = 0.9544 418,181 The Deviance in the output (0.9544) is computed using the DevETðb^Þ formula based on the ET data = –2 ln Lˆ EC1,SS layout. This formula is also equivalent to the difference between SS log-likelihood statistics –(–2 ln Lˆ EC2,SS) for Model EC1 (418.181) and the (fully parame- terized) Model EC2 (417.226, in output below). 417,226 Also, from the table of probabilities (above left), Table of Probabilities ðphat vs: pÞ the observed and predicted probabilities are dif- + ferent, so Model EC1 does not provide perfect group (i.e., covariate pattern) prediction. No perfect group prediction e.g., group 3: phat ¼ 0.09497 p ¼ 0.11864 EXAMPLE We now show edited ouput for the fully para- meterized Model EC2. Edited Output (Model EC2): (Variables – CAT, AGE, ECG, AGE × ECG, CAT × AGE, As with Model EC1, Model EC2 has eight covar- AGE × ECG, and CAT × AGE × ECG) iate patterns, and only Q ¼ 6 percentile groups are obtained in the table of observed and 1 274 17 17.00 257 257.00 expected cases and noncases. However, since Model EC2 is fully parameterized (k þ 1 ¼ 8 ¼ 2 59 7 7.00 52 52.00 G), corresponding observed and expected cases and noncases are identical throughout the table. 3 122 15 15.00 107 107.00 Consequently, the HL test statistic is zero, as 4 57 9 9.00 48 48.00 are both the Deviance and Pearson statistics. 5 39 9 9.00 30 30.00 The log likelihood statistic of 417.226 is equiv- alent to the SS deviance (i.e., DevSSðb^Þ) for 6 58 14 14.00 44 44.00 Model EC2. Since this deviance value is differ- ent from 0, we know that Model EC2 is not the 0.0000 4 1.0000 (SS) saturated model that perfectly predicts the 0 or 1 outcome for each of the 609 subjects 0.0000 4 in the dataset. 0.0000 4 – 2 Log L 417.226 1 0.06205 0.06205 2 0.12295 0.12295 3 0.11864 0.11864 4 0.15625 0.15625 5 .012500 0.12500 6 0.23079 0.23079 7 0.17647 0.17647 8 0.24138 0.24138
Presentation: V. Examples of the HL Statistic 323 Pred ¼ obs for each g Yet, as the table of probabilities indicates, ðphat ¼ pÞ Model EC2 perfectly predicts the observed probabilities obtained for each covariate pat- + tern. Equivalently, then, this model perfectly d^g ¼ dg cases predicts the observed number of cases (dg) corresponding to each covariate pattern. EXAMPLE When none of the predictors are continuous, as Fully parameterized model with these data, these results support the use of (w/o conts. Xs): a fully parameterized model defined from all Above results support its use as covariate patterns as the gold standard model for assessing GOF. In particular, the HL statistic gold standard for assessing GOF reflects this framework, since the value of HL HL statistic always 0 when will always be zero whenever there is perfect group, rather than subject-specific, prediction. perfect group prediction We now provide a second illustration of GOF Second illustration assessment using the Evans County data (see (Evans County data): Computer Appendix) with models that involve Model EC3 (no interaction): continuous variables. In particular, we con- logit PðXÞ ¼ a þ bCAT þ g1AGE sider two previously considered models (see Chap. 7) shown on the left that involve the þ g2ECG þ g3SMK predictors CAT, AGE, ECG, SMK, CHL, and þ g4CHL þ g5HPT HPT. Here, AGE and CHL are continuous, whereas CAT, ECG, SMK, and HPT are binary Model EC4: variables. logit PðXÞ ¼ a þ bCAT þ g1AGE þ g2ECG þ g3SMK þ g4CHL þ g5HPT þ d1CAT Â CHL þ d2CAT Â HPT EXAMPLE Since both models EC3 and EC4 contain con- tinuous variables AGE and CHL, there are not AGE and CHL continuous likely to be many of the 609 subjects with iden- + tically the same values for these two variables. Consequently, the number of covariate pat- Few subjects with identical values for terns (G) for each model should be close to both AGE and CHL the sample size of 609. + G % nð¼ 609Þ SASs Proc Logistic automatically In fact, SASs Logistic procedure automatically outputs “Number of Unique outputs this number (identified in the output Profiles” (G) as the “number of unique profiles”), which turns out to be 599 (see output, below left), Here, G ¼ 599 i.e., G ¼ 599.
324 9. Assessing Goodness of Fit for Logistic Regression EXAMPLE Since G % n in both models, we therefore can- not assume that the Deviance or Pearson sta- G ¼ 599 % n (¼ 609): tistics are approximately chi square. However, Deviance nor Pearson not we can use the Hosmer–Lemeshow (HL) statis- tic to carry out a test for GOF. approximately $ w2 HL statistic approximately $ w2 On the left, we now show edited output for Model EC3, which provides the HL information (can use HL for GOF test) as well as Deviance, Pearson, and À2 log L statistics. Edited Output (Model EC3): (Variables – CAT, AGE, ECG, SMK, CHL, and HPT) From the output for Model EC3, we see that the predicted risks have been divided into (Q ¼ 10) 1 61 0 1.72 61 59.28 deciles, with about 61 subjects in each decile. Also, the observed and expected cases are 2 61 2 2.65 59 58.35 somewhat different within each decile, and the observed and expected noncases are some- 3 61 5 3.46 56 57.54 what different. 4 61 6 4.22 55 56.78 The HL statistic of 5.1028 has Q À 2 ¼ 8 degrees of freedom and is nonsignificant (P ¼ 0.7465). 5 61 7 5.21 54 55.79 Thus, there is no evidence from this test that the no-interaction Model EC3 has lack of fit. 6 61 6 6.10 55 54.90 Both the Deviance (400.3938) and Pearson 7 61 5 7.50 56 53.50 (589.6446) statistics are very different from each other as well as from the HL statistic 8 61 9 9.19 52 51.81 (5.1028). This is not surprising since the Devi- ance and Pearson statistics are not appropriate 9 61 12 11.86 49 49.14 for GOF testing here. 10 60 19 19.08 41 40.92 Note that since the log likelihood for the (SS) saturated model is always zero, the log likeli- No evidence 5.1028 8 0.7465 hood (À2 Log L) value of 400.394 is identical to Model EC3 the Deviance value (i.e., DevSSðb^ÞÞ for Model has lack of fit 400.3938 592 1.0000 EC3. We will use this value of À2 Log L in an 589.6446 592 0.5196 LR statistic (below) that compares the no-inter- action Model E3 to the interaction Model EC4. –2 Log L 400.394 À2 Log L ¼ DevSSðb^Þ ¼ À2 lnðL^SS;EC3=L^SS;maxÞ since Log L^SS;max 0
Presentation: V. Examples of the HL Statistic 325 EXAMPLE Edited output for the interaction model EC4 is now shown here at the left. This model con- Edited Output (Model EC4): tains CAT, AGE, ECG, SMK, CHL, and HPT, (Variables – CAT, AGE, ECG, SMK, CHL, HPT and the product terms CAT  CHL and CAT  HPT in addition to the six main effects. CAT × CHL, and CAT × HPT) From this output, as with Model EC3, we see 1 61 2 0.94 59 60.06 that the predicted risks have been divided into (Q ¼ 10) deciles, with about 61 subjects 2 61 1 1.96 60 59.04 in each decile. Also, the observed and expected cases are somewhat different within each dec- 3 61 2 2.68 59 58.32 ile, and the observed and expected noncases are somewhat different. 4 61 5 3.37 56 57.63 The HL statistic of 7.9914 has Q À 2 ¼ 8 5 61 4 4.07 57 56.93 degrees of freedom and is nonsignificant (P ¼ 0.4343). Thus, there is not sufficient evi- 6 61 2 4.78 59 56.22 dence from this test to conclude that interac- tion Model EC4 has lack of fit (as concluded for 7 61 4 5.77 57 55.23 Model EC3). 8 61 12 7.66 49 53.34 As with Model EC3, both the Deviance (347.2295) and Pearson (799.0580) statistics 9 61 10 11.05 51 49.95 for Model EC4 are very different from each other as well as from the HL statistic (7.9914). 10 60 29 28.71 31 31.29 Although we can conclude from the HL test No evidence 7.9914 8 0.4343 that both Models EC3 and EC4 do not indicate Model EC4 lack of fit, we can decide between these two has lack of fit 347.2295 590 1.0000 models by performing an LR test that com- 799.0580 590 <.0001 pares corresponding log likelihood statistics – 2 Log L for the two models. 347.230 The null hypothesis is that the coefficients of Deviance statistic (347.2295) the two product terms in Model EC4 are both very different from zero. As previously seen in Chap. 7, the test statistic is approximately chi square with Pearson statistic (799.0580) 2 degrees of freedom under the null hypothesis. Both Models EC3 and EC4 The resulting LR value is 53.164, which is do not have lack of fit: highly significant (P < 0.0001). Use LR test to compare Models EC3 Consequently, the interaction Model EC4 is vs. EC4 preferred to the no-interaction Model EC3. H0: d1 ¼ d2 ¼ 0 in Model EC4 LR ¼ À2 ln L^EC3 À ðÀ2 ln L^EC4Þ $ w22 df under H0 Model EC3 Model EC4 À2 Log L 400.394 347.230 LR ¼ 400.394 À 347.230 ¼ 53.165 (P < 0.001): Model EC4 preferred over Model EC3
326 9. Assessing Goodness of Fit for Logistic Regression VI. SUMMARY 3 Chapter 9: Assessing Goodness This presentation is now complete. We have of Fit for Logistic described how to assess the extent to which a Regression binary logistic model of interest predicts the observed outcomes in one’s dataset. Saturated model: We have identified two alternative models, a saturated model and a fully parameterized Contains as many parameters model, that can be used as possible gold stan- (p þ 1) as the number of dard referent points for evaluating the fit of a subjects (n) in the dataset given model. Provides perfect prediction of We have also distinguished between two alter- the observed (0, 1) outcomes on native data layouts that can be used – subject each subject specific (SS) vs. events–trials (ET) formats. Fully parameterized model: A widely used GOF measure for many mathe- matical models is called the deviance. How- Contains the maximum ever, the deviance is not recommended for a number of covariates that can binary logistic regression model in which the be defined from the basic number of covariate patterns (G) is close to predictors (X) being considered the number of subjects (n). for the model In the latter situation, a popular alternative is the Hosmer–Lemeshow (HL) statistic, which Provides perfect prediction of is computed from a table of observed and the observed proportion of cases expected cases and noncases categorized by within subgroups defined by percentile subgroups, e.g., deciles of pre- distinct covariate patterns of X dicted probabilities. Subject-specific (SS) format: Datalines listed by subjects Used for GOF measure of model fit for (0, 1) outcomes Events–trials (ET) format: Datalines listed by subgroups based on (G) covariate patterns Used for GOF measure of model fit for subgroup proportions Deviance: Likelihood ratio (LR) statistic for comparing one’s current model to the saturated model Not recommended when G % n Hosmer–Lemeshow (HL) statistic: GOF statistic appropriate when G%n Computed using O and E cases and noncases in percentile subgroups
Presentation: VII. Appendix: Derivation of the Subject-pecific (SS) 327 We suggest that you review the material cov- ered here by reading the detailed outline that follows. Then do the practice exercises and test. Chapter 10: Assessing Discrimi- In the next chapter (Chap. 10), we describe natory Performance of a Binary methods for assessing the discriminatory per- Logistic Model formance is of a binary logistic model using misclassification tables and ROC curves. VII. Appendix: Derivation \"! of the Subject- n P^ðXkÞ Specific (SS) Deviance DevSS ðb^ Þ ¼ À2 P^ðXi Þ ln À P^ðXk Þ Formula ~ 1 k¼1 # þ ln 1 À P^ðXkÞ Proof. We first write the Deviance formula in a convenient form as follows: DevSSðb^Þ ¼ À2 lnL^ML^CAX ¼ À2 ln L^C since ln L^MAX 0 À2 ln; L^Cdefi¼nition À 2 n h ln P^ðXkÞ: Yk ~ k¼1 i þ ð1 À YkÞ lnð1 À P^ðXkÞ ! n P^ ðXk Þ algebra ¼ À2 ~ Yk ln 1 À P^ðXk Þ k¼1 ! þ lnð1 À P^ðXkÞ : We now write the log of the logistic likelihood function in a convenient form and take its derivative: definition Y P^ðXk ÞYk ð1 À P^ðXk ÞÞ1ÀYk LðbÞ ¼ so k ln LðbÞ ¼ ~ Yk ln P^ðXkÞ þ ð1 À YkÞ lnð1 À P^ðXkÞÞ: k
328 9. Assessing Goodness of Fit for Logistic Regression Taking derivatives, we obtain () @ ln LðbÞ ¼ ~ Yk À 1 À Yk P^ðXkÞ @bj P^ðXk À P^ðXk k Þ 1 Þ Â ð1 À P^ðXkÞÞXjk; which can be rewritten as @ ln LðbÞ ¼ ~ n ð1 À P^ðXk ÞÞ À ð1 À Yk ÞP^ðXk o @bj Yk Þ Xjk k and further simplied to @ ln LðbÞ ¼ ~ ðYk À P^ ðXk ÞÞXjk : @bj k We can then write ~ bj @ ln LðbÞ ¼ ~ À ~ bj Xjk @ bj Yk P^ðXkÞ j k j ! P^ ðXk Þ ¼ ~ Yk À P^ðXkÞ ln 1 À P^ðXkÞ k ¼ ~ Yk À P^ðXkÞ logit P^ðXkÞ: k Since @ ~ln@bðLjYðb^kÞÀ¼P^0ðXfokrÞÞtlhoegiMt PL^ðeXsktÞim¼a0t.e b^ , we can write k ~ Yk logitðP^ðXkÞÞ ¼ It then follows that k ~ P^ðXkÞ logitðP^ðXkÞkÞ. k We then replace ~ Yk logitðP^ðXkÞÞ by ~ P^ðXkÞ logitðP^ðXkÞÞ in k above simplified the k formula for the deviance to obtain n ÂP^ðXkÞ logitðP^ðXkÞÞ DevSS ðb^ Þ ¼ À2 à ~ ÞÞ : þ k¼1 À P^ðXk lnð1
Detailed Detailed Outline 329 Outline I. Overview (pages 304–305) A. Focus: Goodness of fit (GOF) – assessing the extent to which a logistic model estimated from a dataset predicts the observed outcomes in the dataset. B. Considers how well a given model, considered by itself, fits the data. C. Provides a summary measure over all subjects that compares the observed outcome (Yi) for subject i to the predicted outcome ðY^iÞ for this subject obtained from the fitted model. D. Widely used measure is the deviance; however, for binary logistic regression, use of deviance is problematic. Alternative measure: Hosmer–Lemeshow (HL) statistic. II. Saturated vs. Fully Parameterized Models (pages 305–312) A. Saturated model i. Provides perfect prediction of the (0, 1) outcome for each subject in the dataset ii. Contains as many parameters as the number of “subjects” in the dataset iii. Uses data layout in subjects-specific (SS) format iv. Classical model used as gold standard for assessing GOF B. Fully parameterized model i. Contains the maximum number of covariates that can be defined from the basic predictors (X) being considered for the model. ii. The number of parameters (k þ 1) equals the number (G) of distinct covariate patterns (or subgroups) that can be defined from the basic predictors. iii. Uses data layout in events–trials (ET) format. iv. Provides perfect prediction of the observed proportion of cases within subgroups defined by distinct covariate patterns of X. v. An alternative gold standard model for determining GOF. C. Example: n ¼ 40, 2 basic predictors: E (0, 1), V(0, 1) i. Fully parameterized model (G ¼ 4 covariate patterns, k ¼ 3 variables): logit PðXÞ ¼ a þ bE þ gV þ dEV
330 9. Assessing Goodness of Fit for Logistic Regression ii. Saturated model (n ¼ 40 parameters) logit PðXÞ ¼ o1Z1 þ o2Z2 þ o3Z3 ( þ Á Á Á þ o40Z40 # 1 if subject i; i ¼ 1; 2 . . . ; 40 Zi ¼ 0 otherwise : III. The Deviance Statistic (pages 312–317) A. Formula: Devðb^Þ ¼ À2 lnðL^c=L^maxÞ, where b^ ¼ ðb^0; b^1; b^2; . . . ; b^kÞ L^c ¼ ML for current model L^max ¼ ML for saturated model B. Contrasts the likelihood of the current model with the likelihood of the model that perfectly predicts the observed outcomes C. The closer L^c and L^max are to one another, the better the fit (and the smaller the deviance D. Common to test for GOF by comparing deviance with wn2ÀpÀ1 value, but questionably legitimate E. There are two alternative formulae for the deviance: i. Dev ETðb^Þ ! G ln dg þ ðng À dgÞ ln ng Àdg d^g ng Àd^g ¼ À2 ~ dg g¼1 uses events–trials format, where d^g ¼ ngP^ðXgÞ ¼ # of expected cases, dg ¼ # of observed cases, G ¼ # of covariate patterns ii. Dev SSðnb^hÞ i Yi 1ÀYi ¼ À2 ~ Yi ln Y^i þ ð1 À YiÞ ln 1ÀY^i i¼1 uses subject-specific format, where Yi ¼ observed (0, 1) response for subject i and Y^i ¼ predicted probability for subject i ¼ P^ðXiÞ iii. DevSSðb^Þ 6¼ DevETðb^Þ unless G ¼ n iv. Fully parameterized model: DevETðb^Þ ¼ 0 always but DevSSðb^Þ is never 0
Detailed Outline 331 F. Using DevETðb^Þ to test for GOF: i. When G << n, can assume DevETðb^Þ is approximately w2nÀpÀ1 under H0: good fit ii. However, when G % n, cannot assume DevETðb^Þ is approximately w2nÀpÀ1 under H0: good fit iii. Xi continuous, e.g., Xi ¼ AGE ) G % n, so cannot test for GOF G. Why DevSSðb^Þ cannot be used to test for GOF: i. Alternative formula for DevSSðb^Þ : DevSSðb^Þ ¼ À2 n h Þ ln1ÀP^ðP^XðXiÞi þ À P^ðXi i P^ðXi ln 1 Þ ~ Þ i¼1 ii. The above formula contains only the predicted values P^ðXiÞ for each subject; tells nothing about the agreement between observed (0, 1) outcomes and their corresponding predicted probabilities IV. The HL Statistic (pages 318–320) A. Used to provide a significance test for assessing GOF: i. Avoids questionable use of the deviance when G % n ii. Available in most computer procedures for logistic regression iii. Requires that the model considers at least three covariate patterns, rarely results in significance when G is less than 6, and works best when G is close to n, e.g., with continuous predictors B. Steps for computation: 1. Compute P^ðXiÞ for all n subjects 2. Order P^ðXiÞ from largest to smallest values 3. Divide ordered values into Q percentile groupings (usually Q ¼ 10, i.e., deciles) 4. Form table of observed and expected counts 5. Calculate HL statistic from table 6. Compare computed HL to w2 with Q À 2 df
332 9. Assessing Goodness of Fit for Logistic Regression C. Table of observed and expect counts (Step 4) Deciles Obs. Exp. Obs. non Exp. non of risk cases cases cases cases Oc1 Ec1 Onc1 Enc1 1 Oc2 Ec2 Onc2 Enc2 2 Oc3 Ec3 Onc3 Enc3 3 10 Oc10 Ec10 Onc10 Enc10 D. HL Statistic formula (Step 5): HL ¼ Q ðOcq À EcqÞ2 þ Q ðOncq À EncqÞ2 Ecq Encq ~ ~ q¼1 q¼1 V. Examples of the HL Statistic (pages 320–325) A. Two examples, each using the Evans County data (n ¼ 609). B. Example 1 uses two models involving three binary predictors with data layout in events–trials format (G ¼ 8). i. The models Model EC1 (no interaction): logit PðXÞ ¼ a þ bCAT þ g1AGEG þ g2ECG Model EC2: logit PðXÞ ¼ a þ bCAT þ g1AGEG þ g2ECG þ g3AGEG Â ECG þ d1CAT Â AGE þ d2CAT Â ECG þ d3CAT Â AGE Â ECG ii. Model EC2 is fully parameterized, which, as expected, perfectly predicts the observed number of cases (dg) corresponding to each covariate pattern. iii. The HL test statistic for Model EC2 is zero. C. Example 2 uses two models that involve continuous variables. i. The models: Model EC3 (no interaction): logit PðXÞ ¼ a þ bCAT þ g1AGE þ g2ECG þ g3SMK þ g4CHL þ g5HPT Model EC4: logit PðXÞ ¼ a þ bCAT þ g1AGE þ g2ECG þ g3SMK þ g4CHL þ g5HPT þ d1CAT Â CHL þ d2CAT Â HPT. ii. The number of covariate patterns (G) for each model is 599, which is quite close to the sample size (n) of 609.
Detailed Outline 333 iii. The HL statistics computed for both models are not significant, i.e., there is no evidence from this test that either model has lack of fit. iv. Can decide between Models EC3 and EC4 by performing an LR test that corresponding log likelihood statistics or deviances for the two models. VI. Summary (page 326) VII. Appendix: Derivation of SS Deviance Formula (pages 327–328)
334 9. Assessing Goodness of Fit for Logistic Regression Practice The following questions and computer information con- Exercises sider the Evans Country dataset on 609 white males that has been previously discussed and illustrated in earlier chapters of this text. Recall that the outcome variable is CHD status (1 ¼ case, 0 ¼ noncase), the exposure variable of interest is CAT status (1 ¼ high CAT, 0 ¼ low CAT). In this example, we consider only two categorical control variables AGEG (1 ¼ age > 55, 0 ¼ age 55) and ECG (1 ¼ abnormal, 0 ¼ normal). The dataset involving the above variables is given as follows: Cases Total CAT AGE ECG 17 274 0 0 0 15 122 0 1 0 7 59 0 0 1 5 32 0 1 1 1 8100 9 39 1 1 0 3 17 1 0 1 14 58 1 1 1 The SAS output provided below was obtained for the fol- lowing logistic model: Logit PðXÞ ¼ a þ b1CAT þ g1AGE þ g2ECG Deviance and Pearson Goodness-of-Fit Statistics Criterion Value DF Value/DF Pr > ChiSq Deviance 0.9544 4 0.2386 0.9166 Pearson 0.9793 4 0.2448 0.9129 Number of unique profiles: 8 Model Fit Statistics Criterion Intercept Only Intercept and Covariates À2 Log L 438.558 418.181 Analysis of Maximum Likelihood Estimates Parameter DF Estimate Std Wald Intercept 1 À2.6163 Error Chi-Sq Pr > ChiSq cat 1 0.6223 0.2123 151.8266 <.0001 age 1 0.6157 0.3193 ecg 1 0.3620 0.2838 3.7978 0.0513 0.2904 4.7050 0.0301 1.5539 0.2126
Practice Exercises 335 Partition for the Hosmer and Lemeshow Test Event Nonevent Group Total Observed Expected Observed Expected 1 274 17 18.66 257 255.34 2 59 7 5.60 52 53.40 3 122 15 14.53 107 107.47 4 57 9 8.94 48 48.06 5 39 9 7.85 30 31.15 6 58 14 15.41 44 42.59 Hosmer and Lemeshow Goodness-of-Fit Test Chi-Square DF Pr > ChiSq 0.9474 4 0.9177 Questions about the above output begin on the follow- ing page. 1. Is data listing described above in events trials (ET) format or in subject-specific (SS) format? Explain briefly. 2. How many covariate patterns are there for the model being fitted? Describe them. 3. Is the model being fitted a fully parameterized model? Explain briefly. 4. Is the model being fitted a saturated model? Explain briefly. 5. a. Is the deviance value of 0.9544 shown in the above output calculated using the deviance formula Devðb^Þ ¼ À2 lnðL^c=L^maxÞ; where L^c ¼ ML for current model and L^max ¼ ML for saturated model? Explain briefly. b. State the logit form of two logistic models that can be used to calculate the deviance value of 0.9544. Hint: One of these models is the model being fitted. c. How can the deviance value of 0.9544 be calculated using the difference between two log likelihood values obtained from the two models stated in part b? What are the values of these two log likelihood functions? d. What is actually being tested using this deviance statistic? Explain briefly. e. How can you justify that this deviance statistic is approximately chi-square under the null hypothesis that the fitted model has adequate fit to the data?
336 9. Assessing Goodness of Fit for Logistic Regression 6. a. What can you conclude from the Hosmer– Lemeshow statistic provided in the above output about whether the model has lack of fit to the data? Explain briefly. b. Why does the output shown under “Partition for the Hosmer and Lemeshow Test” involve only 6 groups rather than 10 groups, and why is the degrees of freedom for the test equal to 4? Explain briefly. c. What two models are actually being compared by the Hosmer–Lemeshow statistic of 0.9474? Explain briefly. d. How can you choose between the two models described in part c? e. Does either of the two models described in part c perfectly fit the data? Explain briefly. Additional questions using the same Evans County data described at the beginning of these exercises consider SAS output provided below for the following (interac- tion) logistic model: Logit PðXÞ ¼ a þ b1CAT þ g1AGE þ g2ECG þ g3AGE Â ECG þ d1CAT Â AGE þ d2CAT Â ECG þ d3CAT Â AGE Â ECG Deviance and Pearson Goodness-of-Fit Statistics Criterion Value DF Value/DF Pr > ChiSq Deviance 0.0000 0 · · Pearson 0.0000 0 · · Number of unique profiles: 8 Model Fit Statistics Intercept and Criterion Intercept Only Covariates À2 Log L 438.558 417.226 Analysis of Maximum Likelihood Estimates Std Wald Parameter DF Estimate Error Chi-Sq Pr > ChiSq Intercept 1 À2.7158 0.2504 117.6116 <.0001 cat 1 0.7699 1.0980 0.4917 0.4832 age 1 0.7510 0.3725 4.0660 0.0438 ecg 1 0.7105 0.4741 2.2455 0.1340 catage 1 À0.00901 1.1942 0.0001 0.9940 catecg 1 À0.3050 1.3313 0.0525 0.8188 ageecg 1 À0.4321 0.7334 0.3471 0.5557 cae 1 0.0855 1.5245 0.0031 0.9553
Practice Exercises 337 Partition for the Hosmer and Lemeshow Test Group Total Event Nonevent 1 274 Observed Expected Observed Expected 2 59 3 122 17 17.00 257 257.00 4 57 7 7.00 52 52.00 5 39 6 58 15 15.00 107 107.00 9 9.00 48 48.00 9 9.00 30 30.00 44 44.00 14 14.00 Hosmer and Lemeshow Goodness-of-Fit Test Chi-Square DF Pr > ChiSq 0.0000 4 1.0000 7. Is the model being fitted a fully parameterized model? Explain briefly. 8. Is the model being fitted a saturated model? Explain briefly. 9. a. Is the deviance value of 0.0000 shown in the above output calculated using the deviance formula Devðb^Þ ¼ À2 lnðL^c=L^maxÞ; b. where L^c ¼ ML for current model and L^max ¼ ML c. for saturated model? Explain briefly. 10. a. How can the deviance value of 0.0000 be calculated b. using the difference between two log likelihood c. functions? What is actually being tested using this deviance statistic? Explain briefly. What can you conclude from the Hosmer– Lemeshow statistic provided in the above output about whether the interaction model has lack of fit to the data? Explain briefly. What two models are actually being compared by the Hosmer–Lemeshow statistic of 0.0000? Explain briefly. Does the interaction model perfectly fit the data? Explain briefly.
Search
Read the Text Version
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
- 108
- 109
- 110
- 111
- 112
- 113
- 114
- 115
- 116
- 117
- 118
- 119
- 120
- 121
- 122
- 123
- 124
- 125
- 126
- 127
- 128
- 129
- 130
- 131
- 132
- 133
- 134
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 147
- 148
- 149
- 150
- 151
- 152
- 153
- 154
- 155
- 156
- 157
- 158
- 159
- 160
- 161
- 162
- 163
- 164
- 165
- 166
- 167
- 168
- 169
- 170
- 171
- 172
- 173
- 174
- 175
- 176
- 177
- 178
- 179
- 180
- 181
- 182
- 183
- 184
- 185
- 186
- 187
- 188
- 189
- 190
- 191
- 192
- 193
- 194
- 195
- 196
- 197
- 198
- 199
- 200
- 201
- 202
- 203
- 204
- 205
- 206
- 207
- 208
- 209
- 210
- 211
- 212
- 213
- 214
- 215
- 216
- 217
- 218
- 219
- 220
- 221
- 222
- 223
- 224
- 225
- 226
- 227
- 228
- 229
- 230
- 231
- 232
- 233
- 234
- 235
- 236
- 237
- 238
- 239
- 240
- 241
- 242
- 243
- 244
- 245
- 246
- 247
- 248
- 249
- 250
- 251
- 252
- 253
- 254
- 255
- 256
- 257
- 258
- 259
- 260
- 261
- 262
- 263
- 264
- 265
- 266
- 267
- 268
- 269
- 270
- 271
- 272
- 273
- 274
- 275
- 276
- 277
- 278
- 279
- 280
- 281
- 282
- 283
- 284
- 285
- 286
- 287
- 288
- 289
- 290
- 291
- 292
- 293
- 294
- 295
- 296
- 297
- 298
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 310
- 311
- 312
- 313
- 314
- 315
- 316
- 317
- 318
- 319
- 320
- 321
- 322
- 323
- 324
- 325
- 326
- 327
- 328
- 329
- 330
- 331
- 332
- 333
- 334
- 335
- 336
- 337
- 338
- 339
- 340
- 341
- 342
- 343
- 344
- 345
- 346
- 347
- 348
- 349
- 350
- 351
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 363
- 364
- 365
- 366
- 367
- 368
- 369
- 370
- 371
- 372
- 373
- 374
- 375
- 376
- 377
- 378
- 379
- 380
- 381
- 382
- 383
- 384
- 385
- 386
- 387
- 388
- 389
- 390
- 391
- 392
- 393
- 394
- 395
- 396
- 397
- 398
- 399
- 400
- 401
- 402
- 403
- 404
- 405
- 406
- 407
- 408
- 409
- 410
- 411
- 412
- 413
- 414
- 415
- 416
- 417
- 418
- 419
- 420
- 421
- 422
- 423
- 424
- 425
- 426
- 427
- 428
- 429
- 430
- 431
- 432
- 433
- 434
- 435
- 436
- 437
- 438
- 439
- 440
- 441
- 442
- 443
- 444
- 445
- 446
- 447
- 448
- 449
- 450
- 451
- 452
- 453
- 454
- 455
- 456
- 457
- 458
- 459
- 460
- 461
- 462
- 463
- 464
- 465
- 466
- 467
- 468
- 469
- 470
- 471
- 472
- 473
- 474
- 475
- 476
- 477
- 478
- 479
- 480
- 481
- 482
- 483
- 484
- 485
- 486
- 487
- 488
- 489
- 490
- 491
- 492
- 493
- 494
- 495
- 496
- 497
- 498
- 499
- 500
- 501
- 502
- 503
- 504
- 505
- 506
- 507
- 508
- 509
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 519
- 520
- 521
- 522
- 523
- 524
- 525
- 526
- 527
- 528
- 529
- 530
- 531
- 532
- 533
- 534
- 535
- 536
- 537
- 538
- 539
- 540
- 541
- 542
- 543
- 544
- 545
- 546
- 547
- 548
- 549
- 550
- 551
- 552
- 553
- 554
- 555
- 556
- 557
- 558
- 559
- 560
- 561
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 571
- 572
- 573
- 574
- 575
- 576
- 577
- 578
- 579
- 580
- 581
- 582
- 583
- 584
- 585
- 586
- 587
- 588
- 589
- 590
- 591
- 592
- 593
- 594
- 595
- 596
- 597
- 598
- 599
- 600
- 601
- 602
- 603
- 604
- 605
- 606
- 607
- 608
- 609
- 610
- 611
- 612
- 613
- 614
- 615
- 616
- 617
- 618
- 619
- 620
- 621
- 622
- 623
- 624
- 625
- 626
- 627
- 628
- 629
- 630
- 631
- 632
- 633
- 634
- 635
- 636
- 637
- 638
- 639
- 640
- 641
- 642
- 643
- 644
- 645
- 646
- 647
- 648
- 649
- 650
- 651
- 652
- 653
- 654
- 655
- 656
- 657
- 658
- 659
- 660
- 661
- 662
- 663
- 664
- 665
- 666
- 667
- 668
- 669
- 670
- 671
- 672
- 673
- 674
- 675
- 676
- 677
- 678
- 679
- 680
- 681
- 682
- 683
- 684
- 685
- 686
- 687
- 688
- 689
- 690
- 691
- 692
- 693
- 694
- 695
- 696
- 697
- 698
- 699
- 700
- 701
- 702
- 703
- 704
- 705
- 706
- 707
- 708
- 709
- 1 - 50
- 51 - 100
- 101 - 150
- 151 - 200
- 201 - 250
- 251 - 300
- 301 - 350
- 351 - 400
- 401 - 450
- 451 - 500
- 501 - 550
- 551 - 600
- 601 - 650
- 651 - 700
- 701 - 709
Pages:
