Home / Shared Math Details / DOE / Logistic Estimation
Logistic Estimation¶
Quantum XL fits a binary output with one engine and a nominal output with another. Neither has a closed form solution, so both find their coefficients by iterating, and both then report a block of goodness of fit numbers built from the same fit. This page defines the log-likelihood the iteration maximises, the iteration itself, the test on a single coefficient, the two goodness of fit statistics, and the two per-pattern residuals the binary statistics are built from.
Where each number comes from. The two logistic engines are Quantum XL's own code. Every quantity on this page is formed by that code, including the degrees of freedom, and the sections describe the computation directly. Two p values are the exception in part: Quantum XL computes the statistic and chooses the degrees of freedom, then asks the regression library it uses for a tail probability of a named standard distribution. Those are the p value beside a coefficient, which comes from the standard normal distribution, and the two p values in the goodness of fit block, which come from the chi-square distribution.
The two engines are not the same, and the differences are not small. They sum the log-likelihood over different things, they build different information matrices, and they form the two goodness of fit statistics by different arithmetic. Only the binary engine produces per-pattern residuals, and only the binary engine looks for separation when a fit fails. Every section below says which engine it describes.
What the sheet writes is one block, not two. Both logistic reports carry a single block headed Goodness of Fit Test, with one row labelled Pearson and one labelled Deviance, under the columns Chi-Sq, DF and P. The two rows share the same degrees of freedom.
Notation¶
| Term | Description |
|---|---|
| \(X\) | the design matrix, one row per fitted row and one column per model term, the constant column included |
| \(\mathbf{x}_i^{\prime}\) | row \(i\) of the design matrix |
| \(p\) | the number of terms in the model, counting the constant |
| \(J\) | the number of distinct covariate patterns, that is, of distinct rows of the design matrix |
| \(K\) | the number of levels of a nominal output |
| \(\mathbf{b}\), \(b_j\) | the vector of fitted coefficients, and the coefficient of term \(j\) |
| \(SE(b_j)\) | the standard error of the coefficient of term \(j\) |
| \(\eta_i\), \(\eta_{ik}\) | the linear predictor \(\mathbf{x}_i^{\prime}\mathbf{b}\) at row \(i\); for a nominal fit, the one belonging to level \(k\) |
| \(m_i\) | the number of trials at covariate pattern \(i\) of a binary fit |
| \(y_i\) | the number of events at covariate pattern \(i\) of a binary fit |
| \(y_{ik}\) | \(1\) when run \(i\) of a nominal fit had outcome \(k\), otherwise \(0\) |
| \(w_i\) | the weight of run \(i\) of a nominal fit |
| \(\hat{\pi}_i\) | the fitted probability at pattern \(i\) of a binary fit |
| \(\hat{\pi}_{ik}\) | the fitted probability of level \(k\) at pattern \(i\) of a nominal fit |
| \(h_i\) | the leverage of pattern \(i\) of a binary fit |
| \(V\) | the diagonal weight matrix of a binary fit, entry \(m_i \hat{\pi}_i (1 - \hat{\pi}_i)\) for pattern \(i\) |
| \(\ln L\) | the log-likelihood the engine reports, which is the last entry of its log-likelihood history |
| \(\nabla \ln L(\mathbf{b})\) | the gradient of the log-likelihood with respect to the coefficients, at \(\mathbf{b}\) |
| \(I(\mathbf{b})\) | the observed information matrix at \(\mathbf{b}\) |
| \(u\) | the number of half steps already taken at one iteration |
| \(\varepsilon\) | the convergence criterion |
| \(Z_j\) | the Wald statistic of term \(j\) |
| \(p_j\) | the p value reported beside term \(j\) |
| \(\Phi\) | the standard normal cumulative distribution function |
| \(r_{P,i}\) | the Pearson residual at pattern \(i\) of a binary fit |
| \(d_i\) | the deviance residual at pattern \(i\) of a binary fit |
| \(\chi^2_P\), \(\chi^2_D\) | the Pearson and deviance goodness of fit statistics |
| \(O_{ik}\), \(E_{ik}\) | the observed and expected count of outcome \(k\) at pattern \(i\) of a nominal fit |
The log-likelihood¶
The log-likelihood is the quantity both logistic engines maximise. It is written \(\ln L\), it is always negative, and its value has no reading on its own: it exists to be compared, between nested models and inside the model fit statistics built from it.
The two engines do not sum it over the same things.
A binary fit sums over covariate patterns. Before it fits anything the binary engine collapses rows whose design rows are identical: their event counts and their trial counts are added onto the first such row, and the duplicates are then dropped along with any row whose trial count is zero. What is left is one row per covariate pattern, and every later quantity in the fit, this sum included, runs over those \(J\) patterns rather than over the original runs:
\(y_i\) is the number of events at pattern \(i\) and \(m_i\) the number of trials there. The binomial coefficient a full binomial log probability would also carry does not appear, because it does not depend on the coefficients and cannot change which coefficients maximise the sum. When every pattern carries a single trial the expression reduces to one term per run.
A term enters the sum only when the number inside its logarithm is strictly greater than zero. A fitted probability of exactly \(0\) or exactly \(1\) therefore contributes nothing at that end rather than contributing an infinite value.
A nominal fit sums over runs. The nominal engine does not collapse anything before fitting. It sums over the rows as they were supplied, each carrying its own weight \(w_i\):
\(y_{ik}\) is \(1\) when run \(i\)'s outcome is level \(k\) and \(0\) otherwise, and all \(K-1\) of a run's indicators are \(0\) when its outcome is the reference level. The inner sums run over the non-reference levels only, which is why the reference level has no coefficients and no term of its own.
Both engines keep every value the sum has taken, one per iteration plus the value before the first iteration, and report the last entry of that history as the fit's log-likelihood.
The model fit statistics built on this number are described where they are reported, and are not repeated here: McFadden's R squared, G, AIC and BIC on the binary results page and on the nominal results page, together with RMSE and the association measures, which the binary report carries and the nominal one does not.
The binary results page carries this same sum over covariate patterns, reducing to one term per run only when every pattern carries a single trial, which is what a design with one unit per run gives.
The grouping the binary engine performs before it fits is the same grouping several other reports use. See Covariate Pattern.
How the fit iterates¶
There is no closed form for the coefficients of a logistic model, so both engines find them by iterating. The method is Newton's method on the log-likelihood, which for this model is the same computation as iteratively reweighted least squares.
The coefficients start at zero unless a starting vector is supplied. One iteration forms the gradient of the log-likelihood and the observed information matrix at the current coefficients, inverts the information matrix, and steps:
\(u\) is the number of half steps already taken at this iteration and is \(0\) on the first attempt, so the first attempt at every iteration is a full Newton step.
The information matrix is not the same matrix in the two engines. For a binary fit it is \(X^{\prime} V X\), square of side \(p\), with \(V\) the diagonal matrix whose entry for pattern \(i\) is \(m_i \hat{\pi}_i (1 - \hat{\pi}_i)\). For a nominal fit it is square of side \((K-1)p\) and is accumulated one run at a time over a category-expanded design, because a nominal model carries a separate block of coefficients for each non-reference level.
When every row of a binary fit carries the same number of trials, the engine leaves that common count out of both the gradient and the information matrix. That does not change the step, because a factor common to both divides out of \(I(\mathbf{b})^{-1} \nabla \ln L(\mathbf{b})\).
Half stepping. After a step the log-likelihood is recomputed and compared with the last accepted value. If it is larger the step is kept and the iteration ends there. If it is not larger the coefficients are put back where they were, \(u\) is increased by one and the step is tried again at half its previous length, so the multipliers attempted are \(1\), then \(0.5\), \(0.25\), \(0.125\) and so on. The number of half steps allowed at one iteration is a parameter of the engine and defaults to \(7\), which makes \(0.5^7\), one part in \(128\) of the full Newton step, the shortest step ever attempted.
Reaching the half step limit ends the fit, and the fit is then reported as converged. The coefficients kept are the ones from before the rejected step. The line above the coefficient table is the ordinary converged message, the one naming the number of iterations, exactly as it would read for a fit that met the convergence test. Nothing on the sheet separates a fit that stopped at the half step limit from one that converged normally. Both engines behave this way. Each engine does store a distinct message in this branch, but the same branch also marks the fit converged, and a report writes a stored message only for a fit marked not converged, so that message is never written to a sheet.
The convergence test is on the coefficients, not on the log-likelihood. At the end of each iteration the largest relative change across the coefficients is compared with the convergence criterion:
Each difference is divided by that coefficient's own new value, so the test is relative rather than absolute and a large coefficient is allowed to move further in absolute terms than a small one. A coefficient whose new value is exactly zero contributes zero to the maximum rather than a division by zero. The criterion defaults to \(0.000001\). Both engines use this test.
The iteration ceiling defaults to 20, and the two engines do not report reaching it the same way.
A binary fit stops as soon as either the convergence test is met or the ceiling is reached. It is reported as converged whenever that test was met, the last iteration the ceiling allows included, and as not converged only when the ceiling was reached with the test still unmet.
A nominal fit is reported as not converged whenever its iteration count reaches the ceiling, even when the convergence test was met on that last allowed iteration. Its loop ends on either cause and the check that follows reads the iteration count alone, so it cannot tell the two apart. A nominal fit that settles on exactly its twentieth iteration therefore carries the failure warning while a binary fit in the same position does not.
In both engines a fit reported as not converged still returns the coefficients of its last iteration, and every statistic on the sheet is still computed from them, with the warning above the table saying so.
Only the binary engine looks for separation. When a binary fit fails to converge it forms one value per pattern from the observed proportion and the linear predictor,
and adds a sentence to the warning naming complete separation when every one of those values exceeds \(0.99\), or quasi-complete separation when more than half of them exceed \(0.95\). A nominal fit that does not converge carries the warning without that sentence.
The maximum number of iterations, the maximum number of half steps and the convergence criterion are all parameters of the engines. No Quantum XL dialog exposes them, so a logistic fit run from the product uses the three values named above.
The nominal engine's information matrix is described in full under the coefficient standard error on Regression Diagnostics, which is where its inverse is read for the standard errors.
The Wald test on a term¶
Each coefficient is tested on its own by the Wald test, which is the coefficient divided by its own standard error:
The reported p value is twice the upper tail of the standard normal distribution at the absolute value of that ratio:
Quantum XL forms the ratio itself and asks the regression library for the normal tail probability. The reference distribution is the normal one and not a \(t\) distribution, which is why the column is headed Z: the ratio's distribution is a large sample result here rather than an exact one, so there are no residual degrees of freedom to carry into it.
When the standard error is zero, neither number is reported. The engines guard the division rather than performing it, and both the Z cell and the P cell are written as unavailable. A standard error is reported as exactly zero when the matching diagonal entry of the inverted information matrix is not strictly positive.
Both engines do this, and they differ only in how many tests there are. A binary fit has one \(Z_j\) and one \(p_j\) per model term. A nominal fit has one of each per term and per non-reference level, so a model with \(p\) terms and \(K\) levels carries \((K-1)p\) of them, laid out as one coefficient table per non-reference level.
The Z and P columns are described where they sit on the sheet, on the binary results page and the nominal results page. The standard error the ratio divides by is defined on Regression Diagnostics.
The Pearson goodness of fit statistic¶
The Pearson statistic compares observed counts with the counts the fitted model expects, and it is referred to a chi-square distribution. Both engines compute the statistic and choose the degrees of freedom themselves, and ask the regression library for the upper tail probability of the chi-square distribution at that statistic on those degrees of freedom.
For a binary output it is the sum of the squared Pearson residuals over the \(J\) covariate patterns:
on \(J - p\) degrees of freedom, the number of patterns less the number of model terms counting the constant. A pattern whose residual is not a number is left out of the sum rather than making the whole sum not a number. That arises at a pattern the model reproduces exactly, where the fitted probability is exactly \(0\) or exactly \(1\) and the observed count agrees with it, so the residual is zero divided by zero. When the number of patterns equals the number of terms there are no degrees of freedom left and no p value is reported.
For a nominal output it is not built from residuals at all. The nominal engine first groups the runs into \(J\) distinct covariate patterns, which it does only at this point and not for the fit itself, and builds a table of \(J\) rows by \(K\) columns holding the observed count \(O_{ik}\) of each outcome at each pattern and the expected count \(E_{ik}\), which is the pattern's total weight times the fitted probability of that outcome. The statistic is the ordinary chi-square over that whole table, the reference level's column included:
on \((K-1)J - (K-1)p\) degrees of freedom, which is the count of free cells less the \((K-1)p\) coefficients the model estimates. A cell whose value is not a number, which is a cell whose observed and expected counts are both zero, is left out of the sum. When the degrees of freedom come out at zero or below, no p value is reported.
Both engines are affected by how much replication a design has, because \(J\) enters the degrees of freedom directly. A design in which no set of factor settings is repeated has \(J\) equal to its number of runs.
The deviance goodness of fit statistic¶
The deviance statistic asks the same question as the Pearson statistic by a different route, comparing the fitted model's log-likelihood with that of a model that reproduces every observed proportion exactly. It is referred to the same chi-square distribution on the same degrees of freedom as the Pearson statistic, and again the engines form the statistic and choose the degrees of freedom while the regression library supplies the upper tail probability.
For a binary output it is the sum of the squared deviance residuals over the \(J\) covariate patterns:
on \(J - p\) degrees of freedom, the same divisor the Pearson statistic uses. The sign each residual carries is lost in the squaring, so the two statistics share their degrees of freedom but not their arithmetic and they do not have to agree.
For a nominal output it is formed directly from the same table of observed and expected counts the nominal Pearson statistic uses, with no residual in between:
on \((K-1)J - (K-1)p\) degrees of freedom. A cell whose term is not a number, which is any cell holding an observed count of zero, is left out of the sum.
The Pearson residual¶
A logistic fit has no residual in the least squares sense, because the observation is a count of events and the prediction is a probability. The Pearson residual is one of the two standard replacements. It is the difference between the events observed at a covariate pattern and the events the model expects there, divided by the standard deviation of that count under the model:
The denominator is the binomial standard deviation at pattern \(i\), so the residual is in standard deviations of its own count and residuals from patterns of different sizes are comparable with each other. Squaring it and adding over the patterns gives the Pearson goodness of fit statistic.
This is a binary logistic quantity only. The nominal engine produces no Pearson residual, and no nominal output writes one.
The engine also forms a standardized version, which divides by the leverage factor in the same way the studentized residual of a least squares fit does:
Neither form is written to a sheet or drawn on a chart on its own. They are the pieces the reported per-pattern deletion measures are built from.
The deviance residual¶
The deviance residual is the other standard replacement for a least squares residual. It is the signed square root of one covariate pattern's contribution to the deviance, so squaring it and adding over the patterns gives the deviance goodness of fit statistic.
For a pattern where the events observed are neither none nor all of the trials it is
taking the negative sign when \(y_i - m_i \hat{\pi}_i\) is negative, so the residual is negative exactly when the model expected more events than were observed. If the quantity under the root is not strictly positive, which can only happen through rounding at a pattern the model fits almost exactly, the engine reports \(0.0000001\) in place of the root before applying the sign.
The two ends need their own expressions, because a logarithm in the form above would be taken of zero. Where a pattern had no events at all,
and where every trial was an event,
The sign is fixed by which end it is rather than being derived: a pattern with no events can only sit below its expectation, and a pattern of nothing but events can only sit above it.
This is a binary logistic quantity only. The nominal engine produces no deviance residual, and no nominal output writes one. Its deviance statistic is formed straight from observed and expected counts instead.
The per-pattern deletion measures a binary fit reports, and the columns they appear in, are described on Binomial Logistic Residuals.
See Also¶
References¶
- Hosmer, D. W., Lemeshow, S., and Sturdivant, R. X. (2013). Applied Logistic Regression, 3rd edition. Wiley.
- Agresti, A. (2013). Categorical Data Analysis, 3rd edition. Wiley.
- McCullagh, P., and Nelder, J. A. (1989). Generalized Linear Models, 2nd edition. Chapman and Hall.