Backfitting algorithm 2026-10-06
The backfitting algorithm estimates each additive mean function by smoothing its partial residual, obtained after subtracting the current fits of all other terms, and cycling through terms until convergence. Centering each smooth separates its constant part from the intercept. For fixed quadratic smoothness penalties it is block coordinate minimization of a penalized least-squares objective; identifiability and positive definiteness of the constrained problem ensure a unique converged fit. Non-Gaussian generalized additive models can use weighted backfitting inside iteratively reweighted least squares.
For the Poisson regression, the log-likelihood is
Since , the score function and its derivative are
Consequently a finite maximum-likelihood estimator satisfies
These nonlinear equations are normally solved by iteratively reweighted least squares. If has full column rank, the Hessian is negative definite and any finite solution is the unique maximizer. Full rank does not by itself rule out a boundary solution at infinite coefficients for every conceivable data set. The Poisson generalized linear model has dispersion one, so there is no separate scale parameter in this likelihood.
Use the backfitting algorithm to alternate conditional updates of the two additive functions. Initialize and both centered function vectors at zero. With the smoothing matrices for the chosen experience and education spline penalties, one cycle is
followed by
The centering transfers any constant component to the intercept; update to the mean of if necessary, which is under the stated constraints. Repeat until changes in the functions or objective are negligible. Each function is estimated by smoothing its partial residual, namely the response after subtracting all other currently fitted terms.
For fixed penalties these updates are block minimizations of
subject to the centering constraints. Each step cannot increase the objective; with identifiable additive terms and a positive-definite constrained quadratic, backfitting converges to its unique minimizer. Highly dependent predictors can make the decomposition poorly identified or slow convergence. The Gaussian identity-link model needs these least-squares updates directly; non-Gaussian generalized additive models use weighted backfitting inside iteratively reweighted least squares. The package implementation also separates each spline's unpenalized linear component from its nonlinear component, giving the two ANOVA tables displayed in part (e).
For EM for zero-inflated negative binomial regression, at iteration apply Bayes theorem to each zero in the zero-inflated negative binomial model. The E step computes the conditional expectation of its latent indicator:
Treating these values as fixed, the M step maximizes the expected complete-data log-likelihood. Its two parameter blocks separate:
The second block is a weighted negative binomial regression with fixed size . Writing and , its parameter-dependent objective is
Its score function and negative Hessian matrix are
The objective is concave, so a finite maximizer, when it exists, can be found by Newton method or iteratively reweighted least squares. The latter uses working response and weights . This is not an unweighted regression on only the positive counts: zeros retain fractional weight in the count component.
Start with and finite coefficients, alternate the two steps, and monitor the observed log-likelihood. Exact maximization makes the expectation-maximization algorithm nondecreasing in that likelihood, although the mixture may have multiple stationary points. All-zero data can lead to boundary fits and poor identification; multiple initializations and explicit checks for boundary parameters are useful.
The quasibinomial regression retains the logit link and mean function but uses the working variance function . Its quasi-score equation is
so multiplying by the common dispersion parameter does not change the coefficient estimates. Iteratively reweighted least squares therefore gives the same fitted means and coefficients as the binomial model. The usual Pearson dispersion estimator is
The estimated coefficient covariance matrix is , with . From the intercept standard errors, ; the rounded slope errors give approximately the same factor. The display uses approximate Student's t-tests with 118 residual degrees of freedom instead of fixed-dispersion normal tests.
For a genuinely individual binary response, the Bernoulli distribution identity forces . Thus is a working quasi-likelihood specification, not a different independent binary distribution with that mean. The mild estimated underdispersion does not establish an improved model, and the mean predictions are unchanged. In particular, a scalar rescaling does not model correlation among users, and usual likelihood-based Akaike information criterion comparisons are unavailable for a family without a specified probability likelihood. The output gives no convincing reason to prefer model3 to model2.
Write for the log-likelihood and for the score function. Newton method uses the Observed Fisher information ; Fisher scoring uses the Fisher information .
For independent responses in an exponential family, with fixed dispersion parameter and fixed weights , a canonical link function gives and
The last matrix depends on but not on the responses, so taking its expectation leaves it unchanged. Hence the two updates are identical, wherever this Fisher information matrix is invertible. Equivalently both solve the same iteratively reweighted least squares problem. This statement concerns the regression coefficients at a fixed dispersion parameter; a joint update of an unknown dispersion parameter need not have identical observed and expected full Fisher information matrices.
At a current value , put
and define the working response
Then the displayed formulas in the question can be written compactly as
and
Indeed, the th component of the latter is
which is precisely the stated score.
Substitution into the Fisher-scoring update gives
By part (b), this is the weighted least squares fit of the working response on . Both and are recomputed from the current fitted means after every update. Fisher scoring for a generalized linear model is therefore the iteratively reweighted least squares algorithm.
The estimates maximize the grouped-binomial likelihood, numerically obtained by iteratively reweighted least squares. At the fitted probabilities, the estimated covariance matrix is the inverse observed information, approximately
The Std. Error column is the square root of its diagonal. Each z value is Estimate/Std. Error, and the two-sided Wald test p-value is
for the null hypothesis that the corresponding coefficient is zero.
Thus is the fitted log odds for a nonmalignant patient at site A. The site coefficients compare B and C with A, while compares malignant with nonmalignant tumours after adjustment for site. In odds-ratio form,
At level , the B-versus-A contrast is not significant (), the C-versus-A contrast is not significant (), and malignancy has a significant negative association with survival (). The intercept test merely compares the baseline survival probability with .
The null deviance compares the intercept-only fit with the saturated model and has residual degrees of freedom. The residual deviance compares the four-parameter fitted model with the saturated model and has degrees of freedom. Finally, the Akaike information criterion is
Logistic regression assumes independent with
For the binomial exponential family, , so the mean is and the canonical parameter is its logit. It may be fitted by iteratively reweighted least squares.
Under LDA, Bayes' rule gives
which is a logistic model with an intercept.