Past exam of the mathematics course of the University of Cambridge / 2016 / iii / Paper 206 (Applied Statistics)

Use a common intercept and laboratory-specific slopes in the normal linear model. With $t_i$ denoting days and $I_{Bi},I_{Ci}$ the indicator variables for laboratories B and C, the design matrix has row $x_i^T=(1,t_i,t_i{I}_{Bi},t_i{I}_{Ci})$, andThe interaction terms change slopes, without introducing laboratory-specific intercepts. Thus the estimated daily growth rates arein millimetres per day, with fitted common intercept $0.6464$ mm. The residual degrees of freedom are $30-4=26$. The usual unbiased error variance estimate is the square of the residual standard error:The small difference from $5921.7/26\approx 227.758$ arises from rounding the standard error to two decimal places; the latter uses the more precise residual sum of squares supplied later. A maximum-likelihood estimate would instead divide that residual sum of squares by $30$.

Let $\hat{\beta }=(X^{T}X{)}^{-1}{X}^{T}Y$, $H=X(X^{T}X{)}^{-1}{X}^T$, ${\hat{y}}_i=x_i^T\hat{\beta }$, and $e_i=y_i-{\hat{y}}_i$. The horizontal axis of the residual-versus-fitted plot is ${\hat{y}}_i$, the fitted values, and the vertical axis is $e_i$, the regression residual, in millimetres. Under the normal linear model, $\mathrm{Var}(e_i)={\sigma }^2(1-h_{ii})$, where $h_{ii}$ is the regression leverage.

In the quantile-quantile plot, the vertical values are the ordered standardized regression residualsand the horizontal values are corresponding theoretical quantiles ${\Phi }^{-1}(p_k)$ of the standard normal distribution, with plotting positions such as $p_k=(k-1/2)/n$. The exact plotting-position convention has little practical effect here. Under Gaussian errors these points should approximately follow a straight line; they are not independent because fitting induces residual correlations.

The main visible concern is curvature in the conditional mean. The red smooth in the residual-versus-fitted plot descends from positive residuals at low fitted values, becomes negative in the middle, and rises again at high fitted values. This suggests the strictly linear time effects may be inadequate. There is no clear monotone widening of the residual scatter, so strong heteroscedasticity is not apparent. The quantile-quantile plot has modest tail deviations and a few labelled observations, but does not show a dramatic departure from normality. These plots cannot establish independence across days or laboratories; check residuals against day, laboratory, and sampling order as well. A large coefficient of determination does not remove the visible mean-model concern.

For a new independent tumour grown for ten days in B, use $x_0=(1,10,10,0{)}^T$. Its fitted mean isWrite $h_0=x_0^T(X^{T}X{)}^{-1}{x}_0$. The prediction interval in a normal linear model iswhere $s^2=\mathrm{RSS}/26$ and $t_{26,1-\alpha /2}$ is the indicated Student's t-distribution quantile. The design information needed to evaluate $h_0$ is not supplied, so a fully numerical interval cannot be recovered from the coefficient standard errors alone: the coefficient covariances are also needed.

To prove coverage, write $Y_0=x_0^T\beta +{\epsilon }_0$, with ${\epsilon }_0\sim N(0,{\sigma }^2)$ independent of the original sample. The error $Y_0-x_0^T\hat{\beta }$ has a normal distribution with mean zero and variance ${\sigma }^2(1+h_0)$. By the joint distribution of least-squares and variance estimators, $26s^2/{\sigma }^2\sim {\chi }_{26}^2$ and is independent of $\hat{\beta }$; it is also independent of ${\epsilon }_0$. Hence the pivotal quantityConsequently $\mathrm{Pr}(-t_{26,1-\alpha /2}\le T\le t_{26,1-\alpha /2})=1-\alpha$, which is exactly the claimed prediction interval coverage. This is repeated-sampling coverage conditional on the design, assuming the model; it need not hold conditional on the particular observed training responses. The extra $1$ represents new-observation noise and must not be omitted as in a mean confidence interval.

In the larger normal linear model, add laboratory intercept offsets ${\gamma }_B{I}_{Bi}+{\gamma }_C{I}_{Ci}$ to the mean in part (a). The nested-model F-test comparesThe smaller and larger models fit four and six mean coefficients respectively. Under $H_0$, independent Gaussian errors with common variance give independent quantities $({\mathrm{RSS}}_1-{\mathrm{RSS}}_2)/{\sigma }^2\sim {\chi }_2^2$ and ${\mathrm{RSS}}_2/{\sigma }^2\sim {\chi }_{24}^2$. Thus the test statistic isThe p-value is the upper-tail probability $\mathrm{Pr}(F_{2,24}\ge 1.9253)=0.1677$. At the conventional 5% significance level, do not reject the common-intercept restriction. Prefer the simpler model tumour1 on this comparison. This does not prove the intercepts coincide, nor address the curvature identified in part (b); the exact F-test calibration assumes the underlying error and mean assumptions. The one-decimal RSS values give $818.7$ when subtracted, whereas the table's difference uses unrounded values.

The fitted formula and coefficient table describe a random-slope linear mixed model:with random effects independent of the errors. There is a common fixed intercept, a population slope ${\beta }_1$, and a laboratory random slope deviation $b_j$. In tumour1, the three laboratory slopes are unrelated unknown fixed coefficients. Here they are modelled as draws from a common normal distribution; estimated laboratory deviations are shrunk toward zero. After integrating out $b_j$, observations in laboratory $j$ have covariance ${\tau }^2{t}_{ij}{t}_{kj}+{\sigma }^2{1}_{\{i=k\}}$, so they are marginally dependent even when their conditional errors are independent.

The displayed command contains a source inconsistency: size ~ 1 + (0 + days | lab) would impose zero population slope and could not produce the displayed fixed days coefficient. The output corresponds to size ~ days + (0 + days | lab). Its estimates are ${\hat{\beta }}_0=0.6464$, ${\hat{\tau }}^2=17.08$ and ${\hat{\sigma }}^2=227.76$. The estimated daily increase averaged over the laboratory population isbecause $\mathbb{E}{b}_j=0$.

A random slope is appropriate if these laboratories can reasonably be viewed as exchangeable representatives of a larger population and the goal includes population inference or predictions for other laboratories. If these particular three laboratories are the whole set of interest, or have systematic differences incompatible with exchangeability, laboratory-specific fixed slopes are preferable. Only three groups give weak information about between-laboratory variance; the mixed model does not itself resolve possible nonlinear growth. Its population mean is distinct from a realized finite-group average ${\beta }_1+(b_A+b_B+b_C)/3$.

Retain the random-slope linear mixed model of part (e), but replace the within-laboratory independent error covariance by continuous-time autoregressive residual correlation:Errors from different laboratories remain independent and are independent of the random slopes. The resulting marginal covariance within laboratory $j$ isThe exponential decay applies to the conditional errors, not the entire marginal correlation, which also contains the shared random slope.

Use nlme, since lme4::lmer does not fit this general residual correlation structure. With a data frame containing the observed variables, suitable code isHere corCAR1 estimates $\rho =e^{-\kappa }\in (0,1)$, the error correlation one day apart; correlation at separation $\Delta$ is ${\rho }^{|\Delta |}$. This works for noninteger or unequally spaced days. At equally spaced integer days, corAR1(form = ~ days | lab) gives the corresponding discrete-time construction when its parameter is positive. Distinct times within each error series are required: separate tumours with repeated day values should have separate series identifiers or an additional independent-error component, rather than being assigned perfectly correlated errors. The continuous-time autoregressive residual correlation documentation describes the grouping and continuous-time conventions: stat.ethz.ch/R-manual/R-patched/library/nlme/html/corCAR1.html. The fixed days term here follows the fitted output rather than the inconsistent command in part (e).

The Gaussian likelihood has log-likelihood, apart from no omitted parameter-dependent terms,For fixed ${\sigma }^2$, maximizing it is ordinary least squares. Differentiating with respect to $\beta$ gives the normal equation $X^{T}X\hat{\beta }=X^{T}Y$. Full column rank makes the quadratic strictly convex, soAt this value, differentiating with respect to $v={\sigma }^2$ gives $-n/(2v)+\mathrm{RSS}/(2v^2)=0$, soFor positive RSS this is the global likelihood maximum in $v>0$. RSS is positive almost surely under the model since $p<n$; an exactly zero observed RSS instead gives a boundary supremum as $v\downarrow 0$.

As a linear image of a multivariate normal vector, the coefficient estimate satisfiesIndeed, $\hat{\beta }-\beta =(X^{T}X{)}^{-1}{X}^T\epsilon$, and direct matrix multiplication gives its covariance matrix.

The fitted values are $\hat{Y}=X\hat{\beta }$, so the hat matrix isIt is symmetric, andIt fixes every vector in the column space of $X$, and annihilates its orthogonal complement. Thus $H$ is the orthogonal projection onto that column space, has matrix rank and trace $p$, and $I-H$ is the complementary orthogonal projection. ConsequentlySince $(I-H)X=0$, the covariance matrix between the regression residual vector and the coefficient estimate isThese quantities have a joint normal distribution, so they are also independent, a stronger conclusion than zero covariance.

Now $\hat{e}=(I-H)\epsilon$ has mean zero. Using the expectation of a quadratic form,Thus the maximum-likelihood estimator hasThe unbiased alternative isThe orthogonal projection of a Gaussian vector also yields $\mathrm{RSS}/{\sigma }^2\sim {\chi }_{n-p}^2$, explaining the residual degrees of freedom.

Let $s_i$ be sunny days and $l_i$ latitude. The R linear-model formula expands to a normal linear model with both main effects and their interaction term:The fitted conditional mean isThe residual standard error is $0.2836$, giving the usual variance estimate $s^2=0.08042896$, and $n=96+4=100$.

The estimated interaction term is tiny relative to its standard error: $t=-0.103$ and the two-sided p-value is $0.91811$. A sensible initial simplification is to remove the interaction and refit an additive linear model, lm(energy ~ sunny + latitude). Both main effects are significant in the supplied fit, and should initially be retained. Do not simply copy their old coefficient estimates into the new fit: removing a design column generally changes all estimates. Centering predictors makes intercepts and main-effect interpretations more useful, but does not change the fitted surface. The Box–Cox transformation in part (d) suggests a further response-scale improvement; after changing the response scale, reassess the interaction and diagnostics.

The Box–Cox transformation searches over transformations of the positive response,The plotted horizontal coordinate is $\lambda$ and the vertical coordinate is its maximized profile likelihood on the log scale. The likelihood includes the transformation Jacobian determinant; omitting it would compare response scales incorrectly. For the Gaussian regression it is, up to constants,The maximum occurs around $\lambda =2$. The horizontal 95% line is approximately ${\ell }_p(\hat{\lambda })-{\chi }_{1,0.95}^2/2$, and its intersections give an approximate interval roughly $1.2$ to $2.8$, read visually. In particular, $\lambda =1$ is below the line, and $\lambda =0$ is far below it.

Square the response: a convenient transformation is $Y_i^2$, equivalent up to affine rescaling to $T_2(Y_i)=(Y_i^2-1)/2$. Transform energy, not the sunny-days or latitude predictors. Refit and check the normal linear model on the transformed scale; the plot is evidence for a transformation, not proof that the transformed residual assumptions hold. Inference on transformed means cannot be inverted naively to obtain original-scale conditional means.

Check that the chosen mean structure adequately describes the transformed response. Examine the residual-versus-fitted plot and residuals versus each predictor: systematic curvature suggests omitted nonlinear effects or an interaction term. Plot residual spread against fitted values, for example using a scale-location plot, to assess homoscedasticity. Compare standardized residuals with a normal distribution using a quantile-quantile plot, especially for finite-sample Student's t-distribution and F-test inference.

Also check independence using the sampling design and residuals versus time, collection site, or other groups; spatially related photovoltaic systems may have correlated errors that are invisible in a residual-versus-fitted plot. Investigate large standardized regression residuals, high regression leverage, and influential observations using Cook's distance. Check that the design matrix has full rank and that severe multicollinearity is not making estimates unstable. Reassess whether the response transformation and retained predictors improve these diagnostics and prediction; a coefficient p-value alone cannot establish model adequacy. If observations are independent only conditional on site effects, use an appropriate dependence model rather than treating correlated systems as independent replicates.

Assume the known covariance matrix $\Sigma$ is positive definite. The Gaussian likelihood, up to constants, isIts gradient is $X^T{\Sigma }^{-1}(Y-X\beta )$. The negative Hessian matrix $X^T{\Sigma }^{-1}X$ is positive definite because $X$ has full column rank. Thus the unique maximum-likelihood estimator is the generalized least squares estimatorIt is unbiased, and multiplication of its linear coefficient matrix by $\Sigma$ givesIn fact ${\hat{\beta }}_{\mathrm{GLS}}\sim N_p(\beta ,(X^T{\Sigma }^{-1}X{)}^{-1})$. Equivalently, a whitening transformation premultiplies $Y$ and $X$ by ${\Sigma }^{-1/2}$, after which ordinary least squares applies. There is no extra ${\sigma }^2$ factor when $\Sigma$ is the complete known covariance; such a factor appears if only a shape matrix $\Omega$ is known and $\Sigma ={\sigma }^2\Omega$.

The initial Poisson regression contains altitude, temperature, and their interaction term, with a logarithmic link function. Its interaction Wald test has $p=0.521927$, suggesting removal of that term while initially preserving both main effects. The resulting additive model has a negligible temperature effect ($p=0.874$), suggesting a second simplification to altitude alone. The nested analysis of deviance for nested generalized linear models gives the same provisional decisions: deviance increases are $83.326-82.916=0.410$ and $83.351-83.326=0.025$, each for one removed coefficient, with approximate ${\chi }_1^2$ upper-tail probabilities $0.522$ and $0.874$. The Akaike information criterion also decreases at each step.

For independent island counts the final proposed model iswith fitted relationThis is an exponential mean relation on the count scale, not a linear relation between count and altitude. One additional metre multiplies the fitted expected count by $e^{0.01898}\approx 1.01916$. The null model's insignificance of individual main effects in the interaction fit does not require discarding every predictor: correlated design columns can make those conditional tests imprecise. These reductions concern relative mean complexity under the working Poisson regression; the large residual deviance in all three models warns that their variance or mean assumptions still need investigation.

For independent Poisson distributions, ignoring the factorial term only when taking likelihood differences,The saturated likelihood maximizes each mean separately at ${\tilde{\mu }}_i=y_i$, including its limiting value at zero. Therefore the Poisson deviance iswhere $0\log (0/{\hat{\mu }}_i)$ is defined to be zero and ${\hat{\mu }}_i=e^{3.24926+0.01898a_i}$. The expression is twice a saturated-model likelihood ratio difference, not twice a likelihood ratio itself.

With sufficiently large fitted counts and regularity, $D$ is approximately ${\chi }_{30-2}^2$ under a correctly specified Poisson regression. Here $D=83.351$ on 28 residual degrees of freedom, with upper-tail probability $2.0935\times 10^{-7}$. The Poisson model fits poorly on this diagnostic. Its deviance is nearly three times its degrees of freedom, consistent with overdispersion, although a wrong mean function, dependence, or influential islands could also cause the discrepancy. The smaller Akaike information criterion than its competitors does not make this an adequate absolute fit.

The Pearson dispersion estimator is $\hat{\varphi }=3.099079$, well above the Poisson value one. A Quasi-Poisson regression retains the logarithmic mean model but usesThis is a quasi-likelihood specification; it does not by itself identify a count probability distribution. For a common $\varphi$, the quasi-score equation is the Poisson score divided by $\varphi$, so its roots and fitted means are unchanged. ThusCoefficient standard errors are multiplied by $\sqrt{\hat{\varphi }}$. The altitude standard error becomes $0.00447\sqrt{3.099079}\approx 0.0078691$, giving the requested approximate normal confidence intervalIt excludes zero, so retain altitude on this approximate 5% test. Its test statistic is about $2.412$ after allowing for overdispersion, rather than the Poisson value $4.245$. Rescaling uncertainty cannot repair a systematically misspecified mean; the diagnostics in part (e) motivate reconsidering that mean as well.

Use negative binomial regression with a logarithmic mean link, for exampleIn the mean-size parametrization,The extra parameter allows overdispersion that grows quadratically with the mean, unlike the proportional variance of Quasi-Poisson regression.

The printed assertion needs qualification: a negative-binomial model with fixed size is a generalized linear model. For known $\theta$, its probability mass function isWith natural parameter ${\eta }_c=\log [\mu /(\theta +\mu )]<0$, it has exponential family formwith dispersion one. A logarithmic link is permitted even though it is not the canonical link function. This is the fixed-size negative binomial generalized linear model.

What prevents a single ordinary fixed-family GLM fit in glm.nb is joint estimation of the unknown size $\theta$: it changes the cumulant, carrier, and variance function, rather than merely multiplying a fixed variance function by a scalar dispersion. The routine alternates a GLM fit at fixed $\theta$ with size-parameter likelihood updates. Thus the intended distinction is between that enlarged fitting problem and a single fixed-family GLM, not a claim that negative-binomial regression can never be a GLM. The package documentation explicitly supports both conventions: stat.ethz.ch/R-manual/R-patched/library/MASS/html/glm.nb.html and stat.ethz.ch/R-manual/R-patched/library/MASS/html/negative.binomial.html.

The residual smooth rises from negative values, arches above zero, then falls at high fitted values. This suggests a nonlinear altitude effect on the log mean. Since the fitted log mean in the linear model increases with altitude, the horizontal fitted-value axis is also an increasing rescaling of altitude. A generalized additive model with Quasi-Poisson regression is therefore a reasonable candidate:The centering constraint identifies the intercept. Here bs="cr" uses a penalized cubic regression spline. The plotted count residuals are deviance residuals, with the Q-Q panel using leverage-standardized versions. The quantile-quantile plot also shows tail discrepancies, but normality of count residuals is not a distributional assumption of this model: its main requirements are an adequate mean, variance relation, and independence. Investigate the labelled islands and possible excess variation as well.

The estimated smooth in the second plot increases fairly rapidly at low altitude, flattens around altitude 20 to 25, and is nearly level at the upper end. The pointwise bands widen where there are fewer observations. A penalized spline can describe a rise followed by a plateau without forcing a global parabola. A quadratic log mean has slope ${\beta }_1+2{\beta }_{2}a$, so negative curvature eventually forces a decline and its curvature is constant everywhere. A cubic regression spline allows different curvature over different ranges, with a penalty controlling unnecessary wiggles and natural linear tails outside its boundary knots. The displayed smooth supports this greater flexibility, but it is not by itself a formal rejection of every quadratic model. Select complexity and compare predictive performance using cross-validation, then reassess residuals and the Pearson dispersion estimator.

For $0<p<1$, the Bernoulli distribution has massIn exponential dispersion family notation $\exp \{[y\theta -b(\theta )]/\varphi +c(y,\varphi )\}$, takeThe support is $\{0,1\}$; outside it the probability is zero. The natural parameter of an exponential family is the log odds, and the dispersion parameter is fixed rather than freely estimated. Differentiating the cumulant givesAt $p=0$ or $1$ the degenerate distributions are limiting boundary cases with natural parameter $-\infty$ or $+\infty$.

For independent patients, the Bernoulli logistic-regression model iswhere $d_i=1_{\{{\mathrm{frail}}_i=1\}}$ with nonfrail patients as the reference level in a regression factor. The fitted linear predictor isThe null model fits one intercept, so its 99 residual degrees of freedom imply 100 patients; the full model's $100-3=97$ confirms this. Unit weights correspond to individual binary responses, not grouped binomial counts.

The Akaike information criterion estimates expected relative out-of-sample discrepancy, using the maximized log-likelihood plus a complexity penalty:Here $k=3$, since the Bernoulli dispersion is fixed at one. The saturated model puts probability one on each observed binary outcome and has log likelihood zero, so the binomial deviance is $D=-2\ell (\hat{\beta })$. ConsequentlyThere is no estimated error-variance parameter to count. The printed name “Aikake” should be read as Akaike.

This is a sequential analysis of deviance for nested generalized linear models. Adding age to the intercept-only logistic regression reduces deviance by $6.5187$, with approximate ${\chi }_1^2$ p-value $0.01067$. Adding frailty to the age model reduces it by another $16.6474$, with p-value $4.501\times 10^{-5}$. Thus both variables have evidence of association with surgical success, under the model assumptions, and neither should be dropped solely on these results.

The second test assesses frailty adjusted for age, whereas the first assesses age before frailty is included. Sequential tests depend on term order. Evidence for age adjusted for frailty is additionally provided by its full-model Wald test, $p=0.01390$. Both fitted coefficients are negative: age and frailty are associated with lower fitted success probabilities. On the odds scale, one extra year at fixed frailty multiplies the odds by $e^{-0.08942}\approx 0.9145$, and frailty at fixed age multiplies them by $e^{-3.15336}\approx 0.04271$. These are odds ratios, not risk ratios, and an observational statistical association does not by itself establish a causal effect.

The proposed generalized additive model retains the binary response and logit link, but replaces the linear age effect by a smooth function:A penalized cubic regression spline represents $f$; frailty remains a parametric coefficient. The plotted age smooth is centered, so its vertical values represent an age contribution to log odds, not probabilities or the entire fitted predictor.

The plot gives no evidence that the nonlinear model improves the mean structure. The fitted curve is essentially a straight descending line, and the label s(age,1) reports effective degrees of freedom close to one, consistent with a linear effect. The widening bands near the ends reflect sparse data and increasing uncertainty, not detected curvature. A linear-age logistic regression is therefore adequate on this evidence and simpler to interpret. A formal comparative claim would need a suitable test or predictive/model-selection criterion, but the figure supplies no reason to retain extra nonlinear flexibility.

At a target $x$, approximate the unknown mean locally by $a+b(X_i-x)$ and solve the local linear regression weighted least-squares problemwhere $K$ is a nonnegative regression kernel and $h>0$ a smoothing bandwidth. A common factor $1/h$ in all weights makes no difference. Estimate $f(x)$ by the fitted intercept ${\hat{a}}_x$, and optionally $f^{\prime }(x)$ by ${\hat{b}}_x$. Repeat at each target rather than fitting a single global straight line.

For an explicit expression put $S_k={\sum }_i{w}_i(X_i-x{)}^k$ and $T_k={\sum }_i{w}_i(X_i-x{)}^k{Y}_i$. The normal equations yieldThe denominator is positive provided there are at least two distinct positively weighted predictor values. The weights sum to one and their first centered moment is zero, so local linear regression exactly reproduces affine functions and reduces boundary bias compared with a local constant fit. Smaller $h$ uses more local detail and increases variance; larger $h$ smooths more and can increase bias. Choose the bandwidth by cross-validation or another justified criterion. A nearest-neighbour span provides an adaptive bandwidth, as in the local fits used later.

Assume initially that the design points are distinct and ordered on $[a,b]$. The cubic smoothing spline minimizesover functions with an absolutely continuous first derivative and square-integrable second derivative. Dividing the residual term by $n$ simply changes the scaling assigned to $\lambda$. The minimum roughness property of the natural cubic spline interpolant reduces the problem to natural cubic splines with knots at the design points: among functions with any prescribed fitted values, that interpolant has the smallest integrated squared curvature. The minimizer is therefore a natural cubic smoothing spline, with natural second-derivative boundary conditions and linear tails when extended.

Choose a full basis $B_1,\dots ,B_n$ for this natural-spline space. Put $B_{ij}=B_j(X_i)$ and ${\Omega }_{jk}={\int }_a^b{B}_j^{\prime \prime }(t)B_k^{\prime \prime }(t)dt$. Writing $g(t)={\sum }_j{c}_j{B}_j(t)$ changes the objective toDifferentiation gives $(B^{T}B+\lambda \Omega )\hat{c}=B^{T}Y$. Although the curvature penalty has an affine nullspace, the full basis at distinct design points makes the combined matrix positive definite. HenceEquivalently, use the cardinal basis that interpolates the observations, making $B=I$ and $S_{\lambda }=(I+\lambda \Omega {)}^{-1}$. These are smoothing matrices, symmetric but generally not idempotent. Their effective degrees of freedom are $\mathrm{tr}(S_{\lambda })$. As $\lambda \downarrow 0$ the fit interpolates; as $\lambda \uparrow \infty$ it approaches ordinary least squares on the affine functions. With replicated predictors, use a design matrix for distinct knots and observation weights; the basis form remains valid when identifiable.

Useful choices include Leave-one-out cross-validation, generalized cross-validation, and estimation of smoothness by restricted maximum likelihood through a mixed-model representation. For fixed $\lambda$, the leave-one-out residual identity for a linear smoother givesMinimize these over a suitable range. Choose smoothing by predictive performance or a variance-component likelihood, not by arbitrarily forcing a visually appealing curve.

Both fits use the Gaussian generalized additive modelwhere $e_i$ is experience and $d_i$ education. Centering makes the additive decomposition identifiable. No experience-education interaction term is included.

The first fit uses local linear regression for each term. span=1 includes all observations in each nearest-neighbour neighbourhood, but distance weights and a moving target still vary: it does not make the fit exactly a global straight line. The second uses cubic smoothing splines, with separate tuning values spar=0.8 and spar=1.2. Increasing spar increases the curvature penalty and smooths more; spar is a monotone tuning parametrization, not the numerical $\lambda$ in part (b). The relation depends on the predictor scaling and design.

The plotted experience effect is increasing in both fits. The local fit gives a broadly smooth, nearly linear rise; the spline fit follows more of the steep low-experience rise and subsequent flattening, especially near the lower boundary. Education has a weaker increasing effect. The local estimate shows some curvature, while the spline fit at spar=1.2 is closer to a straight line. The clearest difference is the amount of smoothing and the experience curvature, rather than opposite effect directions. These are centered partial-effect plots with partial residuals, not scatterplots of wage against each predictor alone.

A smaller local span would allow the first fit to capture more experience curvature. A moderately smaller education spar could reveal nonlinear detail suppressed in the second fit; its small but significant nonparametric component in part (e) supports checking that possibility. Do not tune solely to follow every point, especially at sparse boundaries. Compare candidate spans or separate spline penalties using cross-validation or generalized cross-validation; increase smoothing if flexible fits introduce unsupported wiggles. See the documented monotone spar–penalty convention at stat.ethz.ch/R-manual/R-patched/library/stats/html/smooth.spline.html.

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 isfollowed byThe 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 ofsubject 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).

The null model's 99 residual degrees of freedom imply $n=100$. Each predictor has a one-degree-of-freedom linear contribution in the parametric table, plus the reported nonlinear contribution. Consequently the termwise effective degrees of freedom areInclude the intercept once, giving about $12.2$ total effective degrees of freedom. The more precise residual degrees of freedom in the output giveThere is no conflict: $8.8$ and $0.4$ are rounded, so they cannot recover the precise total. Individual unrounded termwise values are not provided.

The package's residual mean-square estimate of error variance ismatching the rounded $10328$ in the table. The nonlinear experience contribution is highly significant, while the smaller nonlinear education contribution has approximate $p=0.03923$.

This calculation uses the effective residual degrees of freedom reported by this fit. For a general non-idempotent linear smoother $S$, the noise contribution to expected RSS is ${\sigma }^2\mathrm{tr}[(I-S{)}^T(I-S)]$, and the expectation can additionally contain squared smoothing bias. Thus dividing by $n-\mathrm{tr}(S)$ is not a universal exact-unbiasedness theorem. The output's separation of linear and nonlinear term degrees of freedom is documented at cran.r-project.org/web/packages/gam/gam.pdf.

Write $t_{ij}$ for stirring rate of observation $i$ in furnace $j$. The independent random-intercept and random-slope model isErrors and random effects are independent. The separate R terms (1 | furnace) and (0 + stir | furnace) impose independent random intercepts and random slopes; (1 + stir | furnace) would instead estimate their covariance as well. The estimates areThe corresponding estimated standard deviations are $6.9813$, $0.1147$, and $6.0249$; they are not additional parameters.

For a furnace with predictor vector $t_j$, the marginal formulation, obtained by integrating the Gaussian random effects, isindependently across furnaces. In particular,In stacked notation $Y=X\beta +Zb+\epsilon$, $b\sim N(0,D)$, giving the Gaussian linear mixed model marginal law $N(X\beta ,V)$ with $V=ZD{Z}^T+{\sigma }^{2}I$.