Solution (source code)

= Solution

Use the <backfitting algorithm> to alternate conditional updates of the two additive functions. Initialize $\alpha=\overline Y$ and both centered function vectors at zero. With $S_e,S_d$ the <smoothing matrices> for the chosen experience and education spline penalties, one cycle is
$$
r_e=Y-\alpha\mathbf1-f_d^{\mathrm{old}},\qquad
\widetilde f_e=S_er_e,\qquad
f_e^{\mathrm{new}}=\widetilde f_e-\overline{\widetilde f_e}\mathbf1,
$$
followed by
$$
r_d=Y-\alpha\mathbf1-f_e^{\mathrm{new}},\qquad
\widetilde f_d=S_dr_d,\qquad
f_d^{\mathrm{new}}=\widetilde f_d-\overline{\widetilde f_d}\mathbf1.
$$
The centering transfers any constant component to the intercept; update $\alpha$ to the mean of $Y-f_e-f_d$ if necessary, which is $\overline Y$ 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
$$
\|Y-\alpha\mathbf1-f_e-f_d\|^2
+\lambda_e\int(f_e'')^2+\lambda_d\int(f_d'')^2,
$$
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).