Binary logistic regression

Predicting a binary outcome (participated/did not, recovered/did not, voted/did not) from one or more predictors, continuous or categorical: odds ratios with confidence intervals, the likelihood-ratio test for the model, pseudo-R², the Hosmer-Lemeshow test, a classification table with sensitivity and specificity, and diagnostics – multicollinearity, linearity of the logit and residuals.

Calculator

Loading the data

Drag an Excel file or an SPSS file (.sav) here, or click to choose
In the site's standard template: row = participant, column = variable, first row = headers. CSV is also accepted.
If the file includes a "legend" sheet, the value labels (0 = no, 1 = yes) will appear in the tables, and a "nominal" variable with three or more categories will automatically be marked as categorical (dummy coding).
Just want to try it out? Download a sample data file (320 residents: participation in community activity by length of residence, sense of belonging, age, gender and type of locality): English example
or

Paste columns directly from the sheet (including the header row) – select them in Excel, copy, and paste here:

A
B
C
D
E
F
You can also just paste – the data will be loaded and displayed automatically.

Variable selection

Guide

When to use logistic regression?

When the outcome is a category – participated or not, recovered or not, voted or not – and the goal is to predict it from predictors, continuous or categorical, and to know the contribution of each predictor given the others (Field, 2024, Section 20.2; Pallant, 2020, Chapter 14). It is an extension of the linear model to a binary outcome: the chi-square test (the Chi-square test of independence page) tests the association between two categorical variables; logistic regression predicts the outcome from several predictors together, including continuous ones. For an outcome with three or more categories – multinomial logistic regression, which is not on this page; for a continuous outcome – linear regression.

The model: from probability to log-odds

The linear model predicts a value; here the value is a probability, and a linear relationship between it and the predictor would give values outside 0-1 (Figure 1). The solution: predict not the probability but its log-odds – ln(p/(1 − p)) = b₀ + b₁X₁ + … + bₙXₙ (Field, 2024, Equation 20.13) – and get the probability back through p = 1/(1 + e−(b₀ + b₁X₁ + …)) (equation 20.5). The "odds" are the ratio of the probability of yes to the probability of no: a probability of .75 is odds of 3 (3:1). The parameters are estimated by maximum likelihood – the iterative search for the values that yield the highest probability for the observed data (Field, section 2.7 and tip 20.1) – so there is no "closed formula" here but a process that converges, and sometimes does not converge (the assumptions section).

Figure 1. Why not linear regression: the logistic curve stays between 0 and 1

Horizontal axis – a predictor (e.g. past success rate, Field, 2024, section 20.3.1); vertical axis – the probability of the outcome. The solid line – the logistic model p = 1/(1 + e−(−6 + .12x)); the dashed line – a straight line with the same slope in the middle

020406080100 0.25.5.751.0 The logistic model Straight line – leaves 0-1

Note. Computed from the formula, not simulated. A probability must stay between 0 and 1, and a linear relationship between a predictor and the probability violates this at the ends; the logistic model solves this by having the linear part predict the log-odds (logit) rather than the probability itself (Field, 2024, Sections 20.3.1-20.3.2, Equations 20.5 and 20.13) – and the result is an S-curve. The slope of the curve is not constant: the same B changes the probability a lot in the middle and little at the ends – which is why odds ratios are interpreted, not a "change in probability".

The odds ratio – Exp(B)

The coefficient B is the change in the log-odds for a one-unit increase in the predictor, with the others held constant – and "no one understands what a change in log-odds is. No one" (Field, 2024, Section 20.3.3). So eB, the odds ratio, is reported: the ratio of the odds after a one-unit increase in the predictor to the odds before it (equation 20.17). OR = 1 – no association; OR > 1 – the odds rise with the predictor; OR < 1 – they fall. Field's example (section 20.6.3, box 20.1): treatment versus no treatment, B = 1.229, OR = e1.229 = 3.42 – the odds of recovery are 3.42 times higher for those treated; the direct computation: the odds of recovery without treatment 24/32 = .75, with treatment 41/16 = 2.56, and the ratio 3.42. Pallant (2020, chapter 14): for a continuous predictor the OR is per unit (year, hour, scale point – and therefore depends on the unit); for a categorical predictor – the category versus the reference category; and an OR smaller than 1 is conveniently inverted (1/OR) – "for every hour less of sleep the odds increase by a factor of 1.56" – and then the confidence interval is inverted too. The confidence interval for the OR: if it does not include 1, the direction of the association is significant (Field, Section 20.6.3); if it includes 1, the direction in the population may be the opposite.

Model significance, pseudo-R² and Hosmer-Lemeshow

The log-likelihood and the deviance. Goodness of fit is measured by the log-likelihood (Field, 2024, Equation 20.19) – the counterpart of the residual sum of squares: how much information remains unexplained. It is more convenient to work with the deviance −2LL (deviance; Equation 20.20), because it has a chi-square distribution. The model test (equation 20.21): χ² = (−2LL)null − (−2LL)model, with df = the number of predictors – do the predictors together improve on a model that guesses the most frequent category for everyone. In Field: 154.084 − 144.158 = 9.926, p = .002. In SPSS this is Omnibus Tests of Model Coefficients (Pallant, 2020, chapter 14: "you want a significant value"). Wald for a single predictor (section 20.3.4): z = B/SE, which SPSS prints squared; but when B is large the standard error inflates and Wald shrinks (Menard, 2002, as cited in Field) – the likelihood-ratio test of adding the predictor (a model without and a model with) is more reliable. Pseudo-R². There is no true R²; Cox & Snell (equation 20.25) never reaches 1, and Nagelkerke (equation 20.26) corrects this – "pseudo measures" that are not a proportion of explained variance, but allow comparing models on the same outcome (Pallant). In Field: .084 and .113 for the same model. Hosmer-Lemeshow. A goodness-of-fit test that divides the participants into ten groups by predicted probability and compares observed with expected: here you actually want p above .05 – a low p means poor fit (Pallant, chapter 14: SPSS calls it "the most reliable test" of fit). Field mentions it only as an option (section 20.5.11). Gelman and Hill (2007, p. 100) give a rule of thumb for the deviance: a predictor that is pure noise reduces the deviance by 1 on average, and an informative predictor – by more than 1; k added predictors should reduce it by more than k.

Building the model: forced entry, hierarchy and parsimony

The page enters all the selected predictors in one block – forced entry (Enter), the default of SPSS and of Pallant (2020, chapter 14). Field (2024, section 20.3.8) prefers theory-driven hierarchical building: start with the main predictor, add blocks and check at each stage whether −2LL decreased significantly (the Block χ² in the SPSS output) – and if not, go back to the simpler model ("parsimony": do not enter a predictor without explanatory benefit). In his eel example, adding the duration of the problem changed χ² by .002 (p = .964) and the interaction by .061 (p = .805) – both were dropped. On the site this is done in two runs: a model with the basic predictors and a model with the addition, and the difference between the two χ² values is the test of the addition with df = the number of predictors added. Stepwise methods (Forward/Backward) – Pallant: "not recommended", because they are affected by random variation; Field: not for testing theory, and if used anyway – Backward with the likelihood-ratio criterion. An interaction in the model requires both of its main effects, even if they are not significant (Field).

How to read the output?

Table 1. First the model as a whole: χ² and p (Omnibus). Not significant – stop, the predictors do not improve the guess. Significant – continue to Nagelkerke's R² (comparison with other models) and to Hosmer-Lemeshow (p above .05 = no evidence of poor fit). Table 2. For each predictor: the sign of B gives the direction (Pallant: when everything is coded 0 = no, 1 = yes, a negative B = the predictor reduces the odds of 1); the Wald p – whether the unique contribution is significant; the OR and the confidence interval – the size. Pallant reads her example like this: "the odds of a person reporting a sleep problem are 7.27 times higher for those who have trouble staying asleep, with the others held constant"; and for hours of sleep OR = .64 – "for every additional hour of sleep the odds decrease by a factor of .64". For a categorical predictor with three or more categories, the row for the variable as a whole (Wald with df = the number of categories − 1) says whether the variable contributes, and the dummies say which category differs from the reference. Table 3. The classification table against the baseline (Block 0): if 60% of the sample is in one category, correct classification of 60% is no achievement. Pallant adds the medical measures – sensitivity, specificity, positive and negative predictive value – and the page computes them in the table note. Field (section 20.6.5) also reads the classification plot: in a good model the participants cluster at the ends, not around .5. Gelman and Hill (2007, pp. 99-100) put it as an error rate: the null model errs at p or 1 − p (whichever is lower), and "an error rate equal to that of the null model is terrible"; but the error rate does not distinguish between a prediction of .6 and one of .9 – a model most of whose participants are close to the mean can be very informative at the ends and still improve the overall classification only slightly. Table 4. The diagnostics – in the next section.

Assumptions and problems

Linearity of the logit and independence. The model assumes a linear relationship between each continuous predictor and the log-odds (Field, 2024, Section 20.4.1) – not the probability. The check of Hosmer et al. (2013, as cited in Field, Section 20.7.3): add to the model the interaction of each continuous predictor with its own log; a significant interaction = the assumption is violated. The page does this automatically for positive continuous predictors (Table 4). And the residuals are independent – the same person must not be measured twice.

Incomplete information from the predictors. If a certain combination of categories is empty (smokers who eat tomatoes, in Field's example, section 20.4.2), the model has nothing on which to base the probability there – and with continuous predictors too: "an 80-year-old, highly anxious, Buddhist, left-wing woman" is rare even in a sample of a thousand. The sign: huge standard errors. Field: careful researchers inspect crosstabulations of all the categorical predictors before the model (with the expected-frequencies condition from the chi-square page); "those who do neither – should expect trouble". Complete separation (section 20.4.3): when a predictor or combination of predictors predicts the outcome exactly – cats versus burglars by weight – the estimates run off to infinity, the standard errors explode, and the software may "quietly give up and hand you a completely wrong result". The page detects non-convergence and extreme coefficients and warns; the solution – a simpler model, merging categories, or more data. Sample size (Pallant, 2020, Chapter 14): a small sample with many predictors – convergence problems, especially with sparse categories; run descriptives for each predictor and merge or remove sparse categories.

Multicollinearity. SPSS does not print collinearity diagnostics in logistic regression ("which creates the illusion that it doesn't matter" – Field, section 20.7.4); run a linear regression on the same variables and read only Tolerance and VIF: Tolerance below .1 or VIF above 10 – a problem (Menard, 2002, as cited in Field; Pallant). The page computes the VIF directly (Table 4). Outliers and influential cases (Field, Section 20.6.7; Pallant): the standardized residuals (ZResid) – no more than 5% above 2, no more than 1% above 2.5, and above 3 is an outlier; Pallant prints the cases above 2 (Casewise List) and suggests examining them – perhaps a group the model does not fit (shift workers). Field: an outlier is not a reason to delete – but to examine why; delete only with a documented reason. Influence measures (Cook, DFBeta, leverage) are not on the page.

What the page does not compute. Hierarchical blocks (run twice and subtract the χ² values); stepwise methods; interactions between predictors (Field, Section 20.7); bootstrap for the coefficients (section 20.5.12); influence measures; the SPSS R and Hosmer-Lemeshow's R²L (section 20.3.7); multinomial logistic regression; Bayesian / penalized-likelihood estimation for separation (Gelman and Hill, in the next section).

Three additions from Gelman and Hill: divide-by-four, binned residuals and separation

Gelman and Hill (2007, chapter 5) teach logistic regression from a different direction – not SPSS output but understanding the model – and three of their tools made it onto the page. The "divide-by-four" rule (p. 82). The logistic curve is steepest at its center, where p = .5, and there its slope is exactly B/4; therefore B/4 is the upper bound on the change in probability corresponding to a one-unit increase in the predictor – a reasonable approximation when the probabilities in the data are close to .5, and only a bound when they are close to 0 or 1. In their example, B = 0.33 for income (on a 1-5 scale) gives 0.33/4 = .08: one more income category raises the probability of voting for Bush by no more than 8 percentage points – and the derivative at the mean of the data (pp. 81-82: "evaluation at and near the mean of the data") comes out at .13, close. Table 2 prints B/4 for each predictor, and the table note – the slope at the mean of the data, B·p̂(1 − p̂). Gelman and Hill admit that "the concept of odds can be hard to understand, and odds ratios are even more obscure", and therefore prefer to interpret on the original probability scale (p. 83); on the site the OR is reported (Field; Pallant; APA) and B/4 is added as a translation into probability. Another rule of theirs (p. 83): a coefficient more than two standard errors away from zero is significant – this is the Wald test of Table 2 in different clothing.

Binned residuals (pp. 97-98, figures 5.13-5.14). The residual in logistic regression is y − p̂, and since y is 0 or 1 it is discrete: for p̂ = .7 it is −.7 or +.3; a plot of raw residuals shows two bands and teaches nothing. The solution: divide the participants into equal-sized groups by predicted probability, and for each group plot the mean residual against the mean probability. If the model is correct, about 95% of the points should fall within the bounds ±2√(p̄(1 − p̄)/n) (n – the number of participants in the group), and there is no pattern. In their Bangladesh wells example (N = 3,020, 40 groups of 75) one point out of 40 fell outside – "not a surprise". The same plot by a continuous predictor (figure 5.14) reveals what the overall plot hides: by arsenic level, a bottom group was found in which only 32% switched wells versus the 49% the model predicted, and a rising-and-falling pattern along the axis – "a signal to the experienced modeler" to take the logarithm of the predictor (p. 98), or to add a quadratic term. The page draws the binned-residuals plot by predicted probability and by each continuous predictor (a selector in the chart), in √N groups; Table 4 counts the groups outside the bounds. The complement to the Hosmer-Lemeshow test: the test gives one p, the plot shows where the fit fails.

Separation and nonidentifiability (p. 104). Two reasons why a logistic regression is "not identified" – that the parameters cannot be estimated from the data: collinearity, as in linear regression; and separation, which stems from the discreteness of the outcome. If a predictor is perfectly aligned with the outcome – y = 1 for everyone above a threshold and y = 0 for everyone below it – the best estimate of the coefficient is infinity (figure 5.16: a logistic line with infinite slope); and a linear combination of predictors that is aligned with the outcome also sends at least one coefficient to infinity. This is Field's complete separation (section 20.4.3), and its sign is the same – a huge coefficient with a huge standard error. Gelman and Hill's solution differs from Field's: instead of simplifying the model, Bayesian or penalized-likelihood estimation (bayesglm in the arm package in R, or brlr) that "provides a little information about all the coefficients, even those not identified by the data alone", and returns a finite estimate. The page checks for each predictor whether it is perfectly aligned with the outcome, or whether at one of its values all participants are in the same category (an empty cell), and warns (Table 4); it does not estimate Bayesianly – the R section demonstrates.

And how do you do it in SPSS?

Analyze › Regression › Binary Logistic… (Pallant, 2020, Chapter 14; Field, 2024, Section 20.5.5): the outcome to Dependent, the predictors to Covariates, Method: Enter. Under the Categorical… button – the categorical predictors, Contrast: Indicator, Reference Category: First and click Change (Pallant; Field also demonstrates Last – depending on your coding). Under the Options… button: Classification plots, Hosmer-Lemeshow goodness-of-fit, Casewise listing of residuals, CI for exp(B) and Iteration history (for the −2LL of the null model). Pallant's syntax, with the history added:

LOGISTIC REGRESSION VARIABLES outcome
  /METHOD=ENTER sex age fallsleep staysleep hrswknight
  /CONTRAST(sex)=Indicator(1)
  /CONTRAST(fallsleep)=Indicator(1)
  /CONTRAST(staysleep)=Indicator(1)
  /CLASSPLOT
  /CASEWISE OUTLIER(2)
  /PRINT=GOODFIT CI(95) ITER(1)
  /CRITERIA=PIN(.05) POUT(.10) ITERATE(20) CUT(.5).

Coding the outcome: SPSS internally codes the low value as 0 and the high one as 1 (Dependent Variable Encoding); if the category to be predicted is in fact the low one, create a 1/0 variable – the page's syntax does this with COMPUTE when needed. Reading the output: Block 0 – the classification table of the null model (the baseline) and Iteration History with the initial −2LL; Block 1 – Omnibus Tests of Model Coefficients (the Model row = the χ² of Table 1; the Block row = the addition of the last block), Model Summary (−2LL, Cox & Snell, Nagelkerke), Hosmer and Lemeshow Test, Classification Table = Table 3, Variables in the Equation = Table 2 (B, S.E., Wald, df, Sig., Exp(B), 95% CI), Casewise List – the participants with ZResid above 2 = Table 4. VIF: Analyze › Regression › Linear… with the same variables and Statistics › Collinearity diagnostics – read only the Tolerance and VIF columns (Field, Section 20.7.4). Binned residuals: SPSS does not draw them; /SAVE=PRED RESID saves the predicted probability and the residual, RANK /NTILES bins and MEANS computes means for each group – the bounds and the plot are on the page. The "Download SPSS syntax" button in the output generates all the steps – loading, dictionary, coding the outcome if needed, LOGISTIC REGRESSION with /SAVE, REGRESSION for VIF, binning the residuals, and the check of linearity of the logit.

And how do you do it in R?

glm(…, family = binomial) in base R; a categorical predictor as a factor receives dummy coding with the first level as the reference – like Indicator(first). exp(coef) and exp(confint.default) give the OR and the confidence interval on the log scale (as in SPSS; confint without .default gives a profile-likelihood interval, slightly different); hoslem.test from the ResourceSelection package – the Hosmer-Lemeshow test; binnedplot and bayesglm from Gelman and Hill's arm package – the binned-residuals plot and Bayesian estimation for separation. The script was run as is on the sample data file (in Rscript, with png() around the plot). The site rule for R: the file is saved as UTF-8, every comment line starts with an English word, and Hebrew names require a UTF-8 environment (on Windows: R 4.2 or later, or Sys.setlocale('LC_ALL', 'Hebrew')).

# Binary logistic regression – Professor Stat. This file is saved as UTF-8; every comment line starts with an English word
library(readxl)             # reading Excel
library(ResourceSelection)  # hoslem.test
library(arm)                # binnedplot, bayesglm (Gelman & Hill)

# Loading: the standard template – row = participant, column = variable; the sheet "data"
d <- read_excel("C:/Rdata/DATA_logistic_example_eng.xlsx", sheet = "data")
d$"settlement" <- factor(d$"settlement")        # categorical predictor – dummy coding, the first category (1 = city) is the reference

# The model: participation (0/1) from length of residence, sense of belonging, age, gender and type of locality; listwise deletion as in SPSS
m <- glm(participation ~ tenure + belonging + age + gender + settlement, data = d, family = binomial)
summary(m)$coefficients          # B, standard error, z (Wald in SPSS = z²), p
exp(coef(m))                     # odds ratios Exp(B): residence 1.11, belonging 2.81, kibbutz vs city 2.22
exp(confint.default(m))          # 95% confidence interval for the OR (on the log scale, as in SPSS)
nobs(m)                          # N = 294

# The model versus a null model: likelihood-ratio χ², and pseudo-R²
mf <- model.frame(m); m0 <- glm(participation ~ 1, data = mf, family = binomial)   # null model on the same 294 participants
chi2 <- m0$deviance - m$deviance; df <- m0$df.residual - m$df.residual
c(chi2 = chi2, df = df, p = pchisq(chi2, df, lower.tail = FALSE))       # 59.13, 6, p < .001
n <- nobs(m); cs <- 1 - exp(chi2 / -n); nag <- cs / (1 - exp(-m0$deviance / n))
c(cox_snell = cs, nagelkerke = nag)                                    # .182, .243

# Classification table at the .5 cutoff
table(observed = m$y, predicted = as.integer(fitted(m) >= 0.5))     # 99 47 / 50 98 – 67.0% correct

# Hosmer-Lemeshow test (10 groups) – hoslem.test bins by quantiles of p, so it may differ slightly from the page and from SPSS when there are repeated values
hoslem.test(m$y, fitted(m), g = 10)                                    # 15.53, df 8, p = .050 (the page: 15.05, p = .058)

# Multicollinearity check: VIF from a linear regression on the same predictors (Field, 2024, section 20.7.4)
lm_vif <- lm(participation ~ tenure + belonging + age + gender + settlement, data = d)
1 / (1 - summary(lm(tenure ~ belonging + age + gender + settlement, data = model.frame(lm_vif)))$r.squared)   # VIF of residence = 1.01

# Binned residuals (Gelman & Hill, 2007, pp. 97-98) – binnedplot from the arm package; nclass = the number of groups (the page: √N = 17)
binnedplot(fitted(m), residuals(m, type = "response"), nclass = 17, xlab = "Estimated Pr(participation)")
# The values themselves: mean p̂, mean residual and 2·sd/√n in each group (arm computes the bound from the standard deviation of the residuals in the group, not from p̄(1 − p̄) as in the book – the page follows the book)
br <- arm:::binned.resids(fitted(m), residuals(m, type = "response"), nclass = 17)$binned
round(br[, c("xbar", "ybar", "n", "2se")], 4)                          # group 1: p̂ .135, residual −.135, n 17 (the page: identical)

# The "divide-by-four" rule (Gelman & Hill, 2007, p. 82): B/4 – the upper bound on the change in probability per predictor unit
round(coef(m)[-1] / 4, 3)                                              # residence .026, belonging .258, kibbutz vs city .199

# Complete separation (Gelman & Hill, 2007, p. 104): glm warns "fitted probabilities numerically 0 or 1 occurred"; bayesglm gives a finite estimate
x <- c(-3, -2.5, -1, 0, 0.5, 1, 1.5, 2.5, 3, 3.5, 4, 5, 6); y <- as.integer(x > 2)
coef(summary(glm(y ~ x, family = binomial)))        # B = 44.3, SE = 55,296 – do not interpret
coef(summary(bayesglm(y ~ x, family = binomial)))   # B = 1.86, SE = .88

Reporting the results

Field (2024, section 20.8): like any linear model – a table with B, the standard error and p, the odds ratio with its confidence interval, the constant (so that the reader can reconstruct the equation), and overall model statistics (R², χ²) in the table note; also the nonsignificant predictors, "because that can be as valuable as knowing which predictors are significant". According to APA (2020, section 6.36), R² and p without a leading zero; the OR can exceed 1 and is therefore written with a zero. Pallant (2020, chapter 14) demonstrates the paragraph: "direct logistic regression was performed to assess the impact of… on the likelihood that… The full model was statistically significant, χ²(5, N = 241) = 76.02, p < .001, indicating that it was able to distinguish between… The model correctly classified 75.1% of cases. As shown in Table 1, only three of the variables made a unique statistically significant contribution… The strongest predictor was…, with an odds ratio of 7.27" – followed by a B / SE / Wald / df / p / OR / 95% CI table. The page's report sentences are built on Pallant's paragraph, and Table 2 is her table.

Common mistakes

Interpreting B as in linear regression. B is a change in log-odds; report the OR. And the OR is not a "change in probability" either: 3 times the odds is not 3 times the probability – Gelman and Hill's "divide-by-four" rule gives the bound on the change in probability.

Getting confused about the direction of the coding. An OR of .64 and one of 1.56 are the same finding; SPSS predicts the higher value, and the page – the category you selected. Say in words whose odds of what, and invert the confidence interval together with the reciprocal.

Entering a categorical variable with three categories as a number. "City = 1, moshav = 2, kibbutz = 3" as a continuous predictor says that a kibbutz is "three times a city". Mark it "categorical" and get dummies.

Reading pseudo-R² as explained variance. Cox & Snell and Nagelkerke compare models; they do not say "the model explains 24%".

Being impressed by a classification percentage without a baseline. 70% correct when the baseline is 68% – nothing. And in the opposite direction: a significant model can classify poorly, and is not invalid because of that.

Ignoring huge standard errors. An SE of 40 with an OR of a million is not a finding – it is complete separation or an empty cell (Field, Section 20.4).

Letting SPSS choose the predictors. Stepwise – not for testing hypotheses (Pallant; Field). Build hierarchically, according to theory.

Saying "causes". An odds ratio in correlational research is "associated with"; causality requires a design (the Variables and research designs page).

References

American Psychological Association. (2020). Publication manual of the American Psychological Association (7th ed.). https://doi.org/10.1037/0000165-000

Field, A. (2024). Discovering statistics using IBM SPSS statistics (6th ed.). SAGE.

Gelman, A., & Hill, J. (2007). Data analysis using regression and multilevel/hierarchical models. Cambridge University Press.

Pallant, J. (2020). SPSS survival manual: A step by step guide to data analysis using IBM SPSS (7th ed.). Routledge.

Hosmer et al. (2013), Menard (2002), Cox and Snell (1989), Nagelkerke (1991), Agresti (2013) and Lo et al. (2004) – as cited in Field (2024); Tabachnick and Fidell (2013), Hosmer and Lemeshow (2000), Peat (2001) and Wright (1995) – as cited in Pallant (2020); Landwehr et al. (1984), Albert and Anderson (1984), Firth (1993) and Heinze and Schemper (2003) – as cited in Gelman and Hill (2007, p. 105). By the site rule, the reference list includes only sources read in the original.