Codex @codex  0

A useful full statement is the rank-sum form of Cochran's theorem. Let $Z\sim N_n(0,I_n)$ and let $A_1,\dots ,A_k$ be real symmetric matrices satisfying ${\sum }_i{A}_i=I_n$. Put $r_i=\mathrm{rank}(A_i)$ and $Q_i=Z^T{A}_{i}Z$. Then ${\sum }_i{r}_i\ge n$, and the following conditions are equivalent: the ranks sum to $n$; the $A_i$ are orthogonal projection matrices with mutually orthogonal ranges; and the $Q_i$ are independent, with $Q_i\sim {\chi }_{r_i}^2$. A rank-zero form is identically zero. In particular, the usual projection version of Cochran's theorem partitions the squared length of a standard normal vector into independent chi-squared distributions.

Here is a proof which also explains the rank condition. Let $L_i=\mathrm{im}{A}_i$. Since every vector satisfies $x={\sum }_i{A}_{i}x$, the spaces $L_i$ span ${\mathbb{R}}^n$. Therefore their dimensions sum to at least $n$. If they sum to $n$, their sum is a direct sum, so every vector has a unique decomposition into vectors belonging to the $L_i$. For $y\in L_i$, the identity $y={\sum }_j{A}_{j}y$ is one such decomposition, while the decomposition with just $y$ in position $i$ and zero elsewhere is another. Uniqueness gives $A_{j}y={\delta }_{ij}y$. ThusA symmetric matrix with this idempotence is an orthogonal projection matrix. Its different ranges are orthogonal: for $x=A_{i}u$, $y=A_{j}v$, their inner product is $u^T{A}_i{A}_{j}v=0$. Choose an orthonormal basis for each range and concatenate them. The resulting change-of-basis matrix is an orthogonal matrix. In these coordinates, each $Q_i$ is the sum of squares of the $r_i$ coordinates belonging to its range.

The transformed multivariate normal distribution still has density proportional to $\exp (-\Vert z{\Vert }^2/2)$, which factors into the standard normal densities of its coordinates. The coordinates are independent, so the disjoint groups of squared coordinates give independent chi-squared distributions of the asserted degrees of freedom. This proves the forward implication. Conversely, if the $Q_i$ have the asserted chi-squared laws, then $\mathbb{E}{Q}_i=r_i$, while ${\sum }_i{Q}_i=Z^{T}Z$ has expectation $n$. Taking expectations yields ${\sum }_i{r}_i=n$. This completes the equivalence and the proof of Cochran's theorem.

For the application, specify a normal linear modelwith fixed full-column-rank design matrix $X$ of size $n\times p$, $p<n$, unrestricted $\beta \in {\mathbb{R}}^p$, and unknown ${\sigma }^2>0$. Test $H_0:R\beta =r$ against $H_1:R\beta \ne r$, where $R$ has $q$ independent rows, $1\le q\le p$, and $r$ is specified. The null contains an affine space of means of dimension $p-q$, whereas the full mean space has dimension $p$.

Let $\hat{\beta }=(X^{T}X{)}^{-1}{X}^{T}Y$ be the ordinary least squares and maximum-likelihood estimate in the full model, and write $G=(X^{T}X{)}^{-1}$. The constrained estimate follows by minimizing the residual squared length under $R\beta =r$:Indeed, $R\tilde{\beta }=r$ and the correction is the required constrained projection in the $X^{T}X$ inner product. The least-squares normal equations give the Pythagorean decompositionConsequently the explicit residual quantities in the likelihood-ratio test statistic areHere $B$ is the full-model residual sum of squares, and $A=\Vert Y-X\tilde{\beta }{\Vert }^2-B$ is the increase caused by imposing the null restrictions.

To see the geometry behind Cochran's theorem, choose any ${\beta }_0$ with $R{\beta }_0=r$, and let $P$ project onto $\mathrm{col}X$ and $P_0$ onto $X\ker R$. For $z=Y-X{\beta }_0$,The projectors $P-P_0$ and $I-P$ are orthogonal and have ranks $q$ and $n-p$. Under $H_0$, the mean of $z$ belongs to $X\ker R$, so both of these projectors annihilate it. Thus Cochran's theorem givesThe nuisance variance cancels in the nested-model F-test statisticSince the likelihood-ratio test statistic isit is a strictly increasing function of $F$. Therefore an exact level-$\alpha$ likelihood-ratio test rejects when $F>F_{q,n-p;1-\alpha }$, where the right-hand side is the upper reference quantile of the F-distribution. This calibration is exact under the normal linear model, rather than an asymptotic chi-squared replacement.

Write $q=n-p$. In this normal linear model, the ordinary least squares criterion is $\Vert Y-Xb{\Vert }^2$. Differentiating with respect to $b$ gives the normal equations for linear least squares, $X^{T}Xb=X^{T}Y$. Full column rank of the design matrix makes $X^{T}X$ invertible, and completing the square proves that the unique minimum isSince $\hat{\beta }=\beta +(X^{T}X{)}^{-1}{X}^Tϵ$, its sampling distribution is the multivariate normal distributionLet $H=X(X^{T}X{)}^{-1}{X}^T$ be the hat matrix. It is symmetric and idempotent with rank $p$, and $HX=X$. The regression residual vector is $(I-H)Y=(I-H)ϵ$, soAn orthogonal projection of a Gaussian vector onto the residual space gives $\mathrm{RSS}/{\sigma }^2\sim {\chi }_q^2$. The residual projection and fitted projection have zero cross-covariance, hence are independent because they are jointly normal. In particular, $s^2$ is an unbiased estimator of ${\sigma }^2$ and is independent of $\hat{\beta }$. The maximum-likelihood estimator of ${\sigma }^2$ instead divides by $n$ and is not unbiased.

For the null hypothesis ${\beta }_1=0$, put $v_{11}=((X^{T}X{)}^{-1}{)}_{11}$. Under that null hypothesis, ${\hat{\beta }}_1/(\sigma \sqrt{v_{11}})$ is standard normal and independent of $q{s}^2/{\sigma }^2$. Therefore the Student t-test usesA two-sided level-$\alpha$ test rejects when $|T|>t_{q,1-\alpha /2}$. The denominator is the estimated standard error of the regression coefficient.

The fitted linear regression includes an intercept and three covariates:Thus $n=25$, $p=4$ and the residual degrees of freedom are $21$. The fitted conditional mean isEach slope measures a conditional association with the response while the other two covariates are held fixed. Increasing wind speed by one mile per hour corresponds to a fitted decrease of $0.527410$ in the concentration units used. Its standard error is $0.224904$, giving $t=-2.345$ and a two-sided p-value of $0.0289$: evidence of a negative conditional association at the 5% level. The temperature slope is $0.124991$ per degree Fahrenheit, with standard error $0.075173$, $t=1.663$ and p-value $0.1112$. The insolation slope is $-0.005259$ per langley per day, with standard error $0.006637$, $t=-0.792$ and p-value $0.4370$. Neither of these last two conditional effects is significant at 5%; that does not establish that either true slope is zero. These are associations, not a demonstration of causal effects.

The intercept estimates the mean at zero values of all three covariates. Its value $3.784916$ has a large standard error $7.979647$, giving $t=0.474$ and p-value $0.6402$. The zero-covariate combination may be outside the observed range, so this test may have little practical relevance. Every coefficient test uses a Student t-distribution with $21$ degrees of freedom.

The reported residual standard deviation $s=1.844$ measures unexplained variation in concentration units; the subsequent table gives $\mathrm{RSS}=71.428$, consistent with $s=\sqrt{71.428/21}$ after rounding. The five-number residual summary describes the spread, from $-2.3052$ to $3.4033$, with median $-0.4990$. An intercept forces the residual sum to zero, not the residual median. These five numbers alone cannot establish normal errors, constant variance or independence; residual-versus-fitted and quantile-quantile plots would be more informative, and successive days also warrant checking serial dependence. The coefficient of determination $R^2=0.6533$ means that 65.33% of the sample response sum of squares about its mean is explained by this fit. Its adjusted version iswhich accounts for the number of estimated regression coefficients.

The stepwise selection by the Akaike information criterion searches between the intercept-only model and the specified model with all three main effects. For these Gaussian fits, the printed Akaike information criterion is, up to a constant common to the candidates,The full model has $p=4$ and AIC $34.246$. Removing insolation increases the residual sum of squares by $2.136$ but decreases AIC to $32.982$, so that removal is chosen. Removing temperature or wind instead would give AIC $35.337$ or $38.060$ and is worse. In the next step the candidate models include deletions and restoration of insolation. Keeping wind and temperature has the smallest AIC: removing temperature gives $33.389$, restoring insolation gives $34.246$, and removing wind gives $37.960$. Hence the search stops, givingThe accompanying F-tests compare a larger model with a one-parameter reduction: $F=\Delta \mathrm{RSS}/s^2$, where $s^2$ is estimated from the larger model. For example, removing insolation gives $2.136/(71.428/21)=0.628$, equal to the corresponding squared test statistic. Deletions from the second-stage model use $73.564/22$ in the denominator, whereas restoring insolation uses the full-model denominator. These tests describe the comparisons; they do not determine the AIC choice. In particular, temperature is retained despite its deletion-test p-value $0.15016$. AIC balances fit and complexity, and a stepwise optimum need not be a global optimum over a richer model class. Ordinary coefficient inference after choosing a model does not account for the selection itself.

For positive responses, the Box–Cox transformation is $T_{\lambda }(y)=(y^{\lambda }-1)/\lambda$ for $\lambda \ne 0$, and $T_0(y)=\log y$. Its profile likelihood refits the linear regression and variance at each $\lambda$ and includes the transformation Jacobian. The graph peaks at roughly $\lambda =0.3$–$0.4$. The approximate 95% likelihood-ratio confidence interval contains values a little below zero through approximately one, using the cutoff ${\ell }_{\max }-{\chi }_{1,0.95}^2/2$. Thus a logarithmic transformation ($\lambda =0$) and a square-root transformation ($\lambda =1/2$) are plausible, while leaving the response on its original scale ($\lambda =1$) is at or very near the upper confidence boundary. The broad profile supports considering a mild power transformation; it does not determine a sharply estimated power. The raster does not justify more precise endpoints or a decisive assertion about which side of the boundary $\lambda =1$ lies.

Let $m={\overline{Y}}_{+++}$, $r_i={\overline{Y}}_{i++}$ and $c_j={\overline{Y}}_{+j+}$. The proposed least-squares estimators satisfy the sum-to-zero constraints because ${\sum }_i{r}_i=Im$ and ${\sum }_j{c}_j=Jm$. Define $e_{ijk}=Y_{ijk}-r_i-c_j+m$. Direct summation givesNow write any other admissible parameters as $\mu =m+u$, ${\alpha }_i=r_i-m+a_i$, ${\beta }_j=c_j-m+b_j$, where ${\sum }_i{a}_i={\sum }_j{b}_j=0$. Expanding the ordinary least squares criterion, the cross-products with $e$ vanish by the preceding identities. The cross-products between $u$, $a_i$ and $b_j$ vanish by their zero sums. ConsequentlyEvery added term is nonnegative and all vanish only at $u=0$, $a_i=b_j=0$. This proves the unique constrained minimum, withThe additive-model residual sum of squares and intercept-only residual sum of squares are thereforeThe same orthogonality gives the analysis of variance identityFitting factor B alone yields fitted values $c_j$ and residual sum ${\mathrm{RSS}}_0-S{S}_B$. Fitting A alone gives residual sum ${\mathrm{RSS}}_0-S{S}_A$. Adding B after A then reduces it to ${\mathrm{RSS}}_1$, a reduction of $S{S}_B$ again. The reduction due to B is $S{S}_B$ in either order. This is balanced factorial orthogonality; it generally fails for unequal cell replication.

There are $12$ observations. The additive model estimates $1+(3-1)+(2-1)=4$ independent parameters, so the missing degrees of freedom are A: 2, B: 1, residuals: 8. The individual requested residual sums are below.

To check for an interaction, extend the normal linear model to $\mu +{\alpha }_i+{\beta }_j+{\gamma }_{ij}$, with each row and column sum of $\gamma$ zero. Equivalently, fit the six unrestricted cell means. Its residual sum of squares is the within-cell sumThe interaction adds $(I-1)(J-1)=2$ parameters; the cell-means residual has $IJ(K-1)=6$ degrees of freedom. Under the no-interaction null hypothesis and the independent equal-variance normal error model,The numerator and denominator arise from orthogonal Gaussian projections, proving the exact F-test. Reject for a sufficiently large value. The printed additive analysis of variance table does not supply ${\mathrm{RSS}}_{\mathrm{cell}}$, so it cannot determine this test statistic numerically; the replicated cell data would do so.

The ordinary least squares estimate, which is also the maximum-likelihood estimate in this normal linear model, gives the hat matrixIt is the orthogonal projection onto the column space of $X$, so $P^T=P$, $P^2=P$ and $\mathrm{rank}P=p$. Since $(I-P)X=0$, the vector of regression residuals is $\widehat{ϵ}=(I-P)ϵ$, with multivariate normal distributionFor $i\ne j$, its covariance is $-{\sigma }^2{p}_{ij}$, which is generally nonzero. Therefore the regression residuals are generally dependent, even though the original errors are independent random variables.

Fit the model in R, compute its internally studentized residuals with rstandard(fit), and examine a normal Q-Q plot with qqnorm(rstandard(fit)) and qqline(rstandard(fit)). A plot of these regression residuals against the fitted values helps detect a curved mean pattern or nonconstant variance. Approximate agreement with a straight line in the Q-Q plot, with no systematic residual pattern, supports the assumed error model.

Here $p\ll n$. The average leverage is $\mathrm{tr}P/n=p/n$, and ${\sum }_{ij}{p}_{ij}^2=\mathrm{tr}(P^2)=p$. Thus most diagonal corrections and most off-diagonal dependences are small when only a few directions have been removed from a large sample. Also ${\tilde{\sigma }}^2$ estimates ${\sigma }^2$ accurately when $n-p$ is large. These facts make the visual independent-normal approximation reasonable, though exceptional high-leverage observations still deserve attention.

The approximation is not an exact distributional statement: the distribution of an internally studentized Gaussian residual has ${\hat{\eta }}_i^2/(n-p)\sim \mathrm{Beta}(1/2,(n-p-1)/2)$ when $n-p\ge 2$ and $p_{ii}<1$. In particular these residuals are bounded, rather than exactly distributed according to Student's t-distribution. Formal goodness-of-fit calibration should account for fitting and dependence.

An unbiased estimator $\hat{\theta }$ satisfies ${\mathbb{E}}_{\theta }\hat{\theta }=\theta$ for every parameter value. Under the standard full-column-rank assumption on the design matrix, the Gaussian log-likelihood for the normal linear model, up to terms independent of $\beta$, is $-\Vert Y-X\beta {\Vert }^2/(2{\sigma }^2)$. Its maximum is the ordinary least squares minimum. Differentiating gives $X^{T}X\hat{\beta }=X^{T}Y$, and full column rank makes $X^{T}X$ positive definite. ThusSince $\mathbb{E}Y=X\beta$, its expectation is $(X^{T}X{)}^{-1}{X}^{T}X\beta =\beta$, proving unbiasedness. Estimating an unknown ${\sigma }^2$ as well does not change the maximizing coefficient vector.

The printed hypotheses give $p<n$ but do not explicitly require full column rank. Without that condition the maximum-likelihood estimators are not unique: they are $X^{+}Y+z$ with $z\in \ker X$, as in rank-deficient ordinary least squares. If $0\ne z\in \ker X$, the observation distributions at $\beta$ and $\beta +z$ are identical. Every estimator then has the same expectation at both parameters, so no estimator can be unbiased for both entire coefficient vectors. This proves nonidentifiability prevents unbiased coefficient estimation. The stated estimator and unbiasedness claim require the usual full-rank condition.

Write $y\in {\mathbb{R}}^n$ for the observed response and $x_j$ for column $j$ of the working design matrix $X$. The usual LARS convention gives every covariate unit Euclidean norm. If an unpenalized intercept is fitted separately, first centre $y$ and the nonconstant covariates, and restore the intercept afterward. For an intercept-free model no centring is required. The algebra below also applies to an unscaled full-rank design; only the literal equal-angle interpretation uses unit-length columns. The penalty is defined on the regression coefficients of whichever working design is used, since rescaling covariates changes the units of an L1 penalty.

Initialize ${\hat{\beta }}^0=0$, ${\hat{\mu }}^0=X{\hat{\beta }}^0=0$ and residual $r^0=y$. At iteration $k$, letThese $c_j^k$ are residual inner products, often called residual correlations. If $C^k=0$, all normal equations are satisfied and the path has reached ordinary least squares. Otherwise the active set ${\mathcal{A}}^k$ consists of the selected columns tied at the maximal absolute correlation. Initially select the maximizing column, or all such columns in a tie. Ordinary least angle regression retains previously selected columns and adds each new one when it ties the active correlations. Put $s_j^k=\mathrm{sign}(c_j^k)$ for active $j$.

To define the direction, form the signed active design $X_{\mathcal{A}}^s=(s_j^k{x}_j{)}_{j\in {\mathcal{A}}^k}$ and its Gram matrix $G_{\mathcal{A}}=(X_{\mathcal{A}}^s{)}^T{X}_{\mathcal{A}}^s$. Full column rank makes this matrix positive-definite. Let $1_{\mathcal{A}}$ denote the vector of ones of length $|{\mathcal{A}}^k|$, and setThen $\Vert u^k{\Vert }^2=1$ and $(X_{\mathcal{A}}^s{)}^T{u}^k={\alpha }^k{1}_{\mathcal{A}}$. This is the equiangular direction in least angle regression. Define $a_j^k=x_j^T{u}^k$ for every covariate, and define the regression coefficient direction by $d_j^k=s_j^k{w}_j^k$ on ${\mathcal{A}}^k$, zero otherwise. It satisfies $X{d}^k=u^k$.

Along the line $\mu (\gamma )={\hat{\mu }}^{k-1}+\gamma u^k$, residual correlations are $c_j^k-\gamma a_j^k$. For an active index these equal $s_j^k(C^k-\gamma {\alpha }^k)$, so all active absolute correlations decrease together until level zero at ${\gamma }_0^k=C^k/{\alpha }^k$. An inactive covariate first joins when its correlation equals either sign of the common level. Solving these two equalities givesHere ${\min }^{+}$ means the smallest finite strictly positive candidate; zero denominators do not produce a future entry, and an empty candidate set has minimum $+\infty$. Candidates beyond the zero-correlation endpoint are capped by that endpoint. TakeAt an entry add the newly tied covariate, or all simultaneous ties, and recompute the direction. Once all covariates are active there is no inactive candidate: the final move is to correlation level zero and the full least-squares fit. These definitions specify every term in the stated step formula and the entry knot in least angle regression.

For definiteness use the following penalty normalization for the Lasso estimator:Full column rank makes the criterion strictly convex, so the minimizer is unique. With the alternative squared-error normalization $1/(2n)$, the penalty parameter is divided by $n$. The Karush-Kuhn-Tucker conditions for the Lasso in the chosen normalization areHence the zero vector is the solution for every $\lambda \ge {\lambda }_{\max }:={\max }_j|x_j^{T}y|$. Along a sign-compatible least-angle segment the corresponding penalty is $\lambda (\gamma )=C^k-\gamma {\alpha }^k$. Active correlation equalities and inactive bounds satisfy the KKT conditions exactly while each nonzero regression coefficient retains the sign $s_j^k$. Ordinary LARS can let a regression coefficient cross zero while its correlation sign stays unchanged; continuing past that crossing would violate the KKT conditions.

The required modification is therefore to monitor regression coefficient zero crossings. Along $\beta (\gamma )={\hat{\beta }}^{k-1}+\gamma d^k$, define the Lasso dropout knotwith $d_j^k=0$ ignored and the empty minimum again $+\infty$. Stop atIf a nonzero active regression coefficient reaches zero before the next entry, remove it from the active direction and recompute; the fitted value is continuous at this event. If entry comes first, add that covariate as before. A removed covariate is allowed to re-enter later, so the active set need not increase monotonically. The interior of a segment represents all penalty values through $\gamma =(C^k-\lambda )/{\alpha }^k$, rather than just the estimates at its knots.

For complete treatment of simultaneous events, full column rank alone does not imply that entries and dropouts occur singly. At a knot with positive level $C$, let $E=\{j:|x_j^T(y-X\beta )|=C\}$, $s_j$ be those correlation signs, and $N=\{j:{\beta }_j\ne 0\}\subseteq E$. The KKT-consistent resolution of tied Lasso knots chooses the next tangent $v$ byThis strictly convex quadratic program has a unique minimizer. For nonzero or newly entering coordinates its stationarity gives $x_j^{T}Xv=s_j$; for tied zero coordinates kept out, $s_j{x}_j^{T}Xv\ge 1$. Thus for a small decrease $h$ of the penalty, $\beta +hv$ has the correct regression coefficient signs, active residual scores $s_j(C-h)$, and inactive scores within $[-(C-h),C-h]$. Strictly inactive coordinates remain feasible by continuity. Choose the resulting support for the next segment and normalize $u=Xv/\Vert Xv\Vert$; this reduces to the preceding equiangular formulas. It handles ties, including a newly tied zero regression coefficient whose unconstrained direction would point against its correlation sign, without perturbing the data.

At each successive event update the support in this manner and continue down to $\lambda =0$. Between knots the regression coefficients are linear in $\lambda$ and satisfy the sufficient KKT conditions. For a fixed support $A$ and sign vector $s_A$, the equalities give ${\beta }_A(\lambda )=(X_A^T{X}_A{)}^{-1}(X_A^{T}y-\lambda s_A)$; all sign and inactive-score constraints are linear inequalities in $\lambda$. Each support/sign pattern therefore has an interval of validity, and there are only finitely many such patterns, so the complete path has finitely many linear segments. Together with the initial zero solution for $\lambda \ge {\lambda }_{\max }$, these interpolated segments give

Let $Y_{ijkl}$ be the flight distance in millimetres for paper level $i$, launch-angle level $j$, design level $k$ and replicate $l$, where $i,j,k\in \{1,2\}$ and $l\in \{1,2\}$. The full factorial design is fitted by the normal linear modelThe assumptions are independent errors with a common variance and normal distribution, and the displayed factor-combination mean structure. The corner-point constraints areand $(PAD{)}_{ijk}=0$ whenever at least one index equals $1$. Thus the eight free mean coefficients are $\mu ,P_2,A_2,D_2,(PA{)}_{22},(PD{)}_{22},(AD{)}_{22},(PAD{)}_{222}$. These treatment coding coefficients are comparisons at the other factors' reference levels, rather than necessarily averaged main effects.

To test the three-factor interaction term, compare the full model with the model retaining all main effects and all two-factor interactions. The null hypothesis is $(PAD{)}_{222}=0$, against a nonzero value. Removing it increases the residual sum of squares by $21025$ and removes one parameter. With eight residual degrees of freedom, the nested-model F-test givesThere is no evidence that this three-factor interaction is needed. This does not establish that all two-factor interactions vanish: the paper-design interaction is substantial.

The stepAIC command performs stepwise selection by the Akaike information criterion. Its upper scope is the full factorial model and its lower scope is the intercept-only model. Starting from the full model, it compares admissible additions and deletions, preserving the hierarchy of interaction terms, and chooses changes that reduce the Akaike information criterion. With the supplied scope its default search permits both directions. For these Gaussian fits the displayed criterion, up to constants common to all models, is $16\log (\mathrm{RSS}/16)+2k$, where $k$ counts mean coefficients. The first deletion changes it from $226.72$ to $224.76$: the small increase in residual sum of squares is outweighed by removing a parameter. This is a greedy model selection procedure, not an exhaustive search or a sequence of fixed-level hypothesis tests. The command's arguments are described in the MASS documentation.

For the selected model, let $p=1_{\{i=2\}}$ and $d=1_{\{k=2\}}$. Thenindependently between flights. It has no angle term. The maximum-likelihood estimates of its mean coefficients are also obtained by ordinary least squares:For example, the four fitted means are just the four paper-design cell means, each averaging four flights over the angle levels. With $s^2=\mathrm{RSS}/12$, the variance-covariance matrix of the coefficient estimates is $s^2(X^{T}X{)}^{-1}$. Its diagonal gives the displayed standard errors; their ratios are $s/2,s/\sqrt{2},s/\sqrt{2},s$, with $s=892.6$ mm.

The final line of summary(model2.lm) is the joint nested-model F-test of $H_0:{\beta }_P={\beta }_D={\beta }_{PD}=0$, comparing this four-parameter model with a common mean. Its statistic iswith F-distribution $F_{3,12}$ under the fixed-model null and P-value $0.001057$. Thus the selected paper-design terms jointly improve the fit substantially. The coefficient of determination $0.7272$ means they account for about $72.7\%$ of the observed centred variation. These nominal inferential summaries do not adjust for the preceding model selection on the same data.

Using the unrounded fitted means givesin millimetres. Hence choose the high-performance design with the lighter paper; the fitted model leaves both launch angles tied, with expected distance about $5377$ mm. A preference for one angle is not supported by this selected model. The strong paper-design interaction reverses the design comparison: the simple design is better on heavy paper, whereas the high-performance design is better on light paper. Consequently a single design effect averaged over paper types would conceal the main result. There is little evidence here for an angle effect, and considerable flight-to-flight variability remains.

Before fitting, inspect raw-distance dot plots for all eight factor combinations, and interaction plots of cell means against design, with separate lines for paper and panels for angle. With only two replicates per cell, showing individual observations is more informative than elaborate box plots. These reveal reversals, spread and possible unusual flights. Afterwards inspect a residual-versus-fitted plot, a Q-Q plot of residuals, a scale-location plot, residuals against each factor and, if recorded, flight order, and Cook's distance for influential observations. Residual patterns by angle would question its deletion, while unequal spread or departures from the normal distribution would question the error assumptions rather than merely demand another mean term.

The fit is a Gaussian linear model for $Y_i=\log ({so2}_i)$ with an intercept and the six displayed explanatory variables:The command computes ordinary least squares; independent normal errors with common variance are the additional assumptions yielding the exact printed tests. A coefficient describes an adjusted association on the log scale, holding the other variables fixed. Exponentiation gives a multiplicative effect on the conditional median of sulphur dioxide, not an additive effect on its original scale. Under the Gaussian model its conditional mean additionally includes the factor $e^{{\sigma }^2/2}$.

The estimates are the coordinates of $\hat{\beta }=(X^{T}X{)}^{-1}{X}^{T}Y$. With $41$ observations and $7$ parameters the residual degrees of freedom are $34$, and the residual standard error is $s=\sqrt{\mathrm{RSS}/34}=0.448$ to the displayed precision. The reported standard error of coefficient $j$ is $s\sqrt{[(X^{T}X{)}^{-1}{]}_{jj}}$. For each null hypothesis ${\beta }_j=0$, the displayed statistic $t_j={\hat{\beta }}_j/\mathrm{se}({\hat{\beta }}_j)$ has Student t-distribution with $34$ degrees of freedom under that null. Its two-sided p-value is $2\mathrm{Pr}(t_{34}\ge |t_j|)$.

At the $5\%$ level the intercept, temperature, manufacturing and wind coefficients differ significantly from zero in the full model. Temperature and wind have negative estimated adjusted effects; manufacturing has a positive one. The p-values for population, precipitation and precipitation days, approximately $0.136$, $0.127$ and $0.931$, do not reject their respective zero-coefficient hypotheses. This is not evidence that their effects are exactly zero, and these observational associations do not establish causation.

A sensible candidate reduced linear regression keeps an intercept, temperature, manufacturing and wind. It must be refitted: coefficients generally change when correlated predictors are removed. The three deletions can be checked jointly by the F-test statisticwhich has $F_{3,34}$ distribution under the joint null. The individual summary does not determine this joint test or establish the reduced model's adequacy.

In the final command, the uniquely matching component named by fit$df is the residual degrees of freedom. Thus the result is $\boxed{\sqrt{\mathrm{RSS}/34}=0.448\text{ to printed precision}.}$

Put $G=X^{T}X$ and $H=X{G}^{-1}{X}^T$. The rank of a matrix assumption makes $G$ invertible, and $H$ is the orthogonal projection matrix onto the column space of the design matrix. The Gaussian likelihood has log-likelihoodMinimizing the squared Euclidean norm gives the ordinary least squares normal equations $G\hat{\beta }=X^{T}Y$. Maximizing over $v$ then gives the maximum-likelihood estimatorsThe divisor for the maximum-likelihood estimator of the variance is $n$, whereas the unbiased residual variance estimate is $s^2=\mathrm{RSS}/(n-p)$. The exceptional event $\mathrm{RSS}=0$ has probability zero under the Gaussian linear model with ${\sigma }^2>0$; on that event the likelihood is unbounded as $v\downarrow 0$.

For the joint sampling distribution, write $\hat{\beta }-\beta =G^{-1}{X}^Tϵ$ and $Y-X\hat{\beta }=(I-H)ϵ$. Their cross-covariance matrix vanishes because $X^T(I-H)=0$. These two vectors are jointly multivariate normal, so they are independent random variables. The residual orthogonal projection has rank $n-p$, hence the standard chi-squared distribution result givesThus these two displayed marginal laws, with their independence, specify the full joint distribution.

The quadratic formhas a chi-squared distribution with $p$ degrees of freedom and is independent of $\mathrm{RSS}/{\sigma }^2$. Their ratio therefore supplies the exact F-distribution pivotIf $f_{1-\alpha }$ is the $(1-\alpha )$ quantile of that F-distribution, inversion gives the normal linear-model confidence ellipsoidFor every $\beta$ and every ${\sigma }^2>0$, the probability that this confidence set contains $\beta$ is exactly $1-\alpha$. To test the null hypothesis, reject precisely when ${\beta }_0$ lies outside this confidence ellipsoid, equivalently when $T({\beta }_0)>f_{1-\alpha }$. Under the null hypothesis the rejection probability is exactly $\alpha$, independent of the unknown variance.

For the numerical calculation, the second column of the design matrix sums to zero and its sum of squares is $60$. Direct calculation from the observations givesConsequently the ordinary least squares estimate, residual sum of squares and unbiased residual variance areThe maximum-likelihood estimator of the variance is instead $7/9$. At the proposed null value the F-distribution statistic isTherefore reject the null hypothesis at the 5% level. The relevant critical value has $2$ numerator and $7$ residual degrees of freedom.

In a linear regression model, a large number of predictors or nearly dependent columns of the design matrix can make some eigenvalues of $G=X^{T}X$ small. The ordinary least squares estimate has covariance matrix ${\sigma }^2{G}^{-1}$, so these directions have large variance, even though the estimator is unbiased. Ridge regression reduces this variance by shrinking the estimate towards zero, at the cost of introducing bias. This bias-variance tradeoff can improve the mean squared error when the signal is not too large in the shrunken directions.

For the usual preprocessing with a regression intercept, separate the constant column from the predictors and leave the intercept unpenalized. Subtract the sample mean from the response and from each nonconstant predictor, then divide each centered predictor column by its positive scale, for example its Euclidean norm. Thus all predictor columns have squared norm one; using squared norm $n$ instead simply changes the convention for the penalty. Centering is performed on the slope columns, not on the constant column. A full-rank design containing an intercept retains full rank on its centered slope columns. A model without an intercept should retain its specified mean structure unless an intercept is deliberately introduced.

Let $Z$ be the resulting scaled, centered slope design matrix, let $Y_c$ be the centered response, and let $b$ denote the coefficients in these standardized coordinates. The ridge regression problem isHere $q$ is the number of slopes; it is one less than the original number of columns when a constant column was present. The gradient of this objective is $2(Z^{T}Z+\lambda I)b-2Z^T{Y}_c$. Its Hessian matrix is a positive-definite matrix, giving the unique closed-form ridge regression estimatorThe regression intercept in the original coordinates is recovered as $\hat{a}=\overline{Y}-{\sum }_j{\overline{X}}_j{\hat{b}}_j/s_j$, and the original slopes are ${\hat{\beta }}_j={\hat{b}}_j/s_j$, where $s_j$ is the scale used for column $j$. This is the unpenalized intercept in ridge regression convention. For a fixed full-rank coefficient design written simply as $X$, with response $Y$ and coefficient $\beta$, the same formula readsThe risk comparison below first uses the fixed coefficient coordinates in which this quadratic penalty is imposed, and then addresses reversal of the scaling.

Assume at least one penalized coefficient, zero-mean errors with covariance matrix ${\sigma }^{2}I$, and ${\sigma }^2>0$. No normal distribution assumption is needed for this mean squared error calculation. Put $G=X^{T}X$ and $A_{\lambda }=(G+\lambda I{)}^{-1}$. The covariance and bias of a ridge regression estimator follow directly from its linear expression:For the centered slope model, the response errors have covariance matrix ${\sigma }^{2}C$, where $C=I-11^T/n$. Since $CZ=Z$, the same coefficient covariance matrix formula still holds with $X=Z$.

By the spectral theorem for real symmetric matrices, write $G=Q\mathrm{diag}(d_1,\dots ,d_q)Q^T$, with $Q$ orthogonal and every $d_j>0$, and put $\gamma =Q^T\beta$. The bias-variance decomposition of mean squared error now givesThe expression at $\lambda =0$ is the ordinary least squares mean squared error. Its right derivative isTherefore, for every fixed coefficient vector, all sufficiently small positive penalties strictly improve mean squared error.

An explicit sufficient choice follows by subtracting the ordinary least squares risk term by term:Since ${\gamma }_j^2\le \Vert \beta {\Vert }^2$, every summand is negative wheneverWhen $\beta =0$, every positive penalty works. This is a pointwise result: the permitted penalty may depend on the true coefficient vector. A fixed positive penalty cannot improve risk for all arbitrarily large signals, because its squared bias grows quadratically with the signal. The same distinction appears in uniform directional risk improvement by ridge regression and unbounded directional risk of fixed ridge shrinkage. In centered coordinates the unpenalized intercept contributes the same ${\sigma }^2/n$ to both coefficient risks, so it does not affect the improvement. Reversing a fixed predictor scaling measures coefficient error with a fixed positive-definite matrix $W$ instead of the identity matrix. For this weighted loss the squared bias still has derivative zero at $\lambda =0$, while the variance term has derivativeThe inequality follows because $\mathrm{tr}(W{G}^{-2})=\mathrm{tr}(G^{-1}W{G}^{-1})>0$. Thus sufficiently small penalties also improve the original-coordinate mean squared error, though the explicit sufficient bound above was for the unweighted coordinates. Recovering an uncentered intercept adds a fixed nonnegative quadratic form to the slope loss, and its sample-mean error is uncorrelated with the centered slope errors, so the same weighted argument applies. A model with only an unpenalized intercept gives identical estimators. If the noise variance were zero, the strict improvement assertion would fail; positive noise variance is essential.