Past exam of the mathematics course of the University of Cambridge / 2005 / iii / Paper 41 (Applied Statistics)

scan reads the numerical prices, including the three missing-value markers, into a vector. The data must be read row by row for the subsequent factor construction to match them. The sample size available for analysis is therefore $37$, rather than $40$. The numerical summary describes the observed prices: the sample median is £43.54 and the sample mean is £46.94. The upper tail, particularly the boots, raises the sample mean above the sample median. These summaries mix very different products and are not adjusted country comparisons.

The two character scans create labels. gl(8,5,length=40,labels=item) makes a regression factor with each of eight item labels repeated five times. gl(5,1,length=40,labels=country) cycles through the five country labels once for each item. For a regression factor and numerical response, the two plot commands produce grouped box plots. The country box plots suggest lower US prices, while the item box plots show particularly expensive boots and inexpensive cardigans. They are descriptive displays of this small basket, with different item compositions where prices are missing; they cannot separate country effects from item effects.

The first lm fits an additive two-way analysis of variance:The constraints are R's usual treatment coding. The normal linear model used for the printed tests assumes independent errors with common variance and approximately a normal distribution, and no item-by-country interaction in the mean. Missing prices are omitted, leaving $37$ responses and $1+4+7=12$ independent mean coefficients. Thus the residual degrees of freedom are $37-12=25$. This deletion is defensible for the conditional price model when missingness does not depend on the unobserved error after conditioning on the two factors; the output itself does not establish such a missing-data mechanism.

The anova tables report sequential sums of squares, rather than order-independent adjusted effects. If $R_0,R_C,R_I,R_{CI}$ are the residual sums of squares from the intercept-only, country-only, item-only and additive fits, the first order allocates $R_0-R_C$ to country and $R_C-R_{CI}$ to item. The second order allocates $R_0-R_I$ to item and $R_I-R_{CI}$ to country. Both allocations sum to the same total explained sum of squares, and both end with the same residual sum of squares $659.4$. The three missing combinations destroy balanced factorial orthogonality, so the allocation changes with order. To test country after allowing for items, use the country row in the second table; to test item after allowing for countries, use the item row in the first table. The first-entered rows do not supply those adjusted tests.

Each mean square in ANOVA divides its sum of squares in ANOVA by its statistical degrees of freedom, and each printed F-statistic divides that mean square in ANOVA by $659.4/25=26.38$. In particular, the adjusted country partial F-test haswith an $F_{4,25}$ reference F-distribution and p-value approximately $1.86\times 10^{-6}$. The adjusted item partial F-test gives $F=91.58$ on $7,25$ statistical degrees of freedom. Both factors contribute strongly under the additive normal linear model. The larger item sum of squares in ANOVA measures the much greater differences between these products, but is not evidence that every pair of items differs significantly.

The intercept $50.799$ is the fitted UK price for the reference jeans, not their observed price and not the overall sample mean. Each country regression coefficient is its adjusted price difference from the UK, assumed common across items: Sweden is £5.224 lower, France £10.708 lower, Germany £7.676 lower and the US £21.328 lower. Each item regression coefficient is its adjusted difference from the reference jeans, assumed common across countries. Thus the fitted UK boots price is $50.799+52.854=103.653$, and the fitted US boots price is $103.653-21.328=82.325$.

The coefficient standard errors come from $s^2(X^{T}X{)}^{-1}$, where $s^2=659.4/25$ and $X$ is the design matrix. Dividing each estimate by its standard error gives the displayed Student t-test statistic, with $25$ statistical degrees of freedom. The stars summarize the corresponding two-sided p-values. The individual tests compare each indicated country or item with its reference level, conditional on all other additive terms; they are neither tests of every pairwise comparison nor a simultaneous confidence interval procedure. In particular, Sweden's p-value $0.0626$ provides weaker evidence than the other country contrasts, and does not establish equality with the UK. The intercept test has little substantive relevance, and interpreting many individual stars also raises multiple hypothesis testing issues.

The residual standard error $s=5.136$ measures unexplained prices in pounds. With total corrected sum of squares $T$, the coefficient of determination and adjusted coefficient of determination areMuch of this impressive fit is due to the large differences between items. The overall F-test compares the additive fit with a constant mean: $F=((T-659.4)/11)/(659.4/25)=62.12$ on $11,25$ statistical degrees of freedom. It tests that all eleven non-intercept coefficients vanish, not that the additive model is adequate. Regression diagnostics should examine regression residuals, changing price variability across items and possible interaction. Percentage differences may be more plausible than common pound differences, motivating a natural logarithm of price. No random sample of countries or products, and no independent replication of a price within a cell, is supplied; broad claims about all designer goods or the causes of country differences require more than these conditional tests.

Finally, Item*Country expands to both main effects and their full interaction. A cell-means parameterization has one separate fitted mean for every observed cell. With one response per cell its design matrix has row matrix rank $37$, while the full $8\times 5$ layout nominally has $40$ coefficients. Three coefficient combinations are aliased because their cells are missing. This is saturation of an unreplicated two-factor regression: every observed price is fitted exactly and the residual degrees of freedom are zero. The final fit has zero residuals, $R^2=1$ and no residual-based estimate of the error variance or ordinary interaction significance test. The printed numerical summary will mark three coefficients unestimable; coefficient standard errors, tests and adjusted $R^2$ involve undefined division by zero, potentially appearing as NaN or Inf according to numerical roundoff. Perfect interpolation here is a consequence of a saturated statistical model, not evidence for a superior predictive model.

Write $a$ for the scalar intercept denoted $\mu$ in the question, to distinguish it from the individual means ${\mu }_i$. The log-likelihood of the independent Poisson observations isDifferentiating gives the score equations for an interior maximum-likelihood estimator:In particular, when $m={\sum }_i{y}_i>0$,The left side is a weighted covariate mean whose derivative is its weighted variance, so it is strictly increasing when the covariates are not all equal. An interior solution is unique if the observed count-weighted mean is strictly between the smallest and largest covariates. At an extreme it may require a limiting slope, and if all counts are zero the maximum occurs as $a\to -\infty$, rather than at finite parameters. Constant covariates make the slope and intercept unidentifiable.

Put $S_j={\sum }_i{x}_i^j{\mu }_i$. The negative Hessian matrix is deterministic at given parameters and is also the Fisher information matrix:Under the usual interior, increasing-information regularity for asymptotic normality of a maximum likelihood estimator, replace the unknown means by their fitted values and read the lower-right entry:Equivalently, Poisson slope information after eliminating an intercept is ${\sum }_i{\mu }_i(x_i-{\bar{x}}_{\mu }{)}^2$, with ${\bar{x}}_{\mu }=S_1/S_0$. This explains both the positive denominator for varying covariates and the loss of slope information when their spread is small.

The Poisson deviance compares the fitted Poisson regression with a saturated statistical model:using $0\log 0=0$. If the usual deviance goodness-of-fit test approximation is adequate, compare $27.2$ with a chi-squared distribution on $29$ statistical degrees of freedom. Its upper-tail p-value is about $0.561$, and $D/29=0.938$. There is no detected lack of fit or evidence of substantial overdispersion from this deviance. This neither proves the Poisson regression correct nor establishes that the slope is nonzero. Small expected counts, dependence or inadequate covariate specification can undermine the reference approximation, so deviance residuals and other regression diagnostics still matter.

The Wald test of $H_0:\beta =0$ useswith an approximate standard normal distribution under the null. Its two-sided p-value is $0.403$, and the approximate $95\%$ Wald confidence interval isThere is no convincing evidence of a covariate effect, and the slope is very imprecisely estimated. The model can fit acceptably while its slope is statistically indistinguishable from zero. In the Poisson regression, a unit covariate increase multiplies the mean by $e^{\beta }$; exponentiating the two interval limits shows how poorly this multiplier is determined on the stated scale. One might compare with an intercept-only model using a likelihood-ratio test, but its numerical deviance difference is not supplied. Failure to reject does not establish no effect.

Expanding the ordinary least squares objective givesFor a nonzero vector $v$, $v^T{X}^{T}Xv=\Vert Xv{\Vert }^2>0$, since the design matrix has matrix rank $p$. The Hessian matrix is therefore positive-definite, and $Q$ is strictly convex. Solving its normal equations gives the unique least-squares estimator:Indeed, expansion about this stationary point and $X^T(Y-X\hat{\beta })=0$ give $Q(\beta )=Q(\hat{\beta })+(\beta -\hat{\beta }{)}^T{X}^{T}X(\beta -\hat{\beta })$, which also proves the global minimum directly.

Define the hat matrix $H=X(X^{T}X{)}^{-1}{X}^T$. It is a symmetric matrix and an idempotent linear map, since $H^T=H$ and $H^2=H$, and its image is the column space of $X$. Thus $H$ is the orthogonal projection onto the fitted mean space. The fitted values are $HY$ and the regression residual vector is $(I-H)Y$. Because $I-H$ is also a symmetric matrix and an idempotent linear map,These minimization and projection identities do not need errors with a normal distribution; the normal linear model becomes relevant to the distributional calculations.

Put $G=I-H$. Since $GX=0$, the regression residual vector is $\hat{\epsilon }=G\epsilon$. A linear transformation of a multivariate normal distribution is again multivariate normal, soThis is a singular multivariate normal distribution supported on $\ker X^T$, with dimension $n-p$. In particular,Consequently regression residuals are neither independent nor identically distributed even under the stated model. Their linear constraints include $X^T\hat{\epsilon }=0$, and, if the design matrix contains an intercept, their sum is zero.

Choose an orthonormal basis of the residual space. The coordinates of $G\epsilon$ in that basis are $n-p$ independent $N(0,{\sigma }^2)$ variables, provingFor $n>p$, $s^2$ is an unbiased estimator of ${\sigma }^2$. Orthogonality of $H$ and $G$ gives zero cross-covariance, so the Gaussian fitted and residual vectors are independent. This explains why an appropriate residual-versus-fitted diagnostic has no systematic mean relationship under the model. If $n=p$, no residual information remains for these diagnostics or the variance estimate.

For $h_{ii}<1$, use standardized regression residuals $r_i={\hat{\epsilon }}_i/[s\sqrt{1-h_{ii}}]$ to account for the different residual variances. The diagonal $h_{ii}$ is the regression leverage; $h_{ii}=1$ makes that residual identically zero and cannot be diagnosed by this standardization. With known $\sigma$, the corresponding individually standardized residual has an exact standard normal distribution. With estimated $s$, the internally studentized residual is not exactly normal or a Student t-distribution, and the residuals remain correlated; a normal reference plot is a diagnostic approximation.

Plot the regression residuals or standardized regression residuals against fitted values and each covariate. A persistent curve suggests a misspecified mean, while a fan-shaped spread suggests heteroscedasticity; a scale-location plot helps separate these. A quantile-quantile plot against normal quantiles checks tails, skewness and possible regression outliers. Plot residuals against observation order or time to look for unexplained serial patterns, recognizing that projection already induces the displayed residual covariances. Inspect regression leverage together with Cook's distance for influential observations. Such patterns motivate revising the model; absence of a visible pattern does not prove all its assumptions.

The response errors have standard deviations proportional to $|{\mu }_i|$. Thus, when fitted means approximate the true means, increasing absolute fitted values tends to accompany increasing residual spread. A residual-versus-fitted plot shows a fan for positive means, or spread increasing away from zero for means of either sign. Standardization by one common residual standard deviation, even with the usual regression leverage correction, does not remove this systematic heteroscedasticity.

Under independent errors the exact covariance matrix after fitting iswhere $k>0$ is the common proportionality constant. Hence the simple fan description is a diagnostic tendency, not an assertion that each fitted residual has exactly variance $k{\mu }_i^2$.

For positive responses, the intended variance-stabilizing transformation is the natural logarithm. The delta method givesUse $\log Y_i$ for positive data with modest relative error, then refit and check the transformed mean and residual structure. In particular, if the original mean was $X_i^T\beta$, the transformed mean is approximately $\log (X_i^T\beta )-k/2$, not automatically a linear function of the original covariates.

There is a domain qualification in the literal printed model: a nondegenerate normal distribution has negative support. For positive ${\mu }_i$, the probability of $Y_i\le 0$ is $\Phi (-1/\sqrt{k})$, so a real $\log Y_i$ is not defined throughout the stated sampling model. When $k$ is small this is a negligible-tail approximation, not an exact transformation to normal errors. If signed values must be retained and all ${\mu }_i\ne 0$, write $Y_i={\mu }_i(1+\sqrt{k}{Z}_i)$ with standard normal $Z_i$. ThenThe second term has a common distribution with finite variance, so log absolute value stabilization of a normal scale family is exact, but loses the sign of the mean and does not produce normal errors. Zeros and a zero mean are excluded from this transform. Alternatively retain the response scale and fit the stated variance model with weighted least squares weights proportional to $1/{\hat{\mu }}_i^2$ where the means stay away from zero. No arbitrary addition of a constant before taking logarithms exactly solves the given variance relation.

read.table(...,header=T) reads four grouped rows with column names. The displayed outcome counts give the observed risks $0.1519,0.2145,0.2611,0.5111$ in the indicated factor order. rbind in the first test pools over predisposition: the two rows represent no cannabis use and some use, with event/non-event counts $(341,1775)$ and $(82,238)$. Thus the pooled risks are $341/2116=0.1612$ and $82/320=0.2563$. This test concerns the marginal statistical association between use and outcome, not an association adjusted for predisposition.

Under independence, expected counts in the two-by-two contingency table are $E_{ij}=O_{i+}{O}_{+j}/N$. The ordinary Pearson chi-squared statistic is $17.5183$. The printed $16.8618$ is reproduced by the default two-by-two Yates continuity correction:All absolute differences here exceed one half. Comparison with a chi-squared distribution on one statistical degree of freedom gives $p=4.02\times 10^{-5}$. There is strong evidence of marginal association. The simplified output omits the usual continuity-correction wording, but its numerical value shows which statistic was used.

attach makes the data columns available by name, and tot=with+without obtains group sizes $1936,275,180,45$. The grouped binomial regression fits $Y_g\sim \mathrm{Bin}(n_g,{\pi }_g)$, with response $Y_g/n_g$ and prior weights $n_g$. These weights tell R the binomial denominator; treating the four proportions as four equally precise individual observations would be wrong. The default logit link isLet $A$ indicate some use and $B$ indicate predisposition. The additive logistic regression specifies $\eta =a+bA+cB$. Its log-likelihood, including the binomial constants, isDifferentiation gives score equations $X^T(Y-n\pi )=0$ and Fisher information matrix $X^{T}WX$, where $W_{gg}=n_g{\pi }_g(1-{\pi }_g)$. Fisher scoring solves these equations, and the inverse fitted Fisher information matrix gives the reported approximate coefficient variances. The number of scoring iterations records numerical convergence, not degrees of freedom or statistical evidence.

The intercept corresponds to no use and no predisposition and gives baseline fitted risk ${\mathrm{logit}}^{-1}(-1.73881)=0.1495$. The additive coefficients imply common adjusted odds ratios:They multiply the odds, not the risks. The Wald tests divide the estimates by their standard errors and use a standard normal distribution; both predictors have small two-sided p-values. Approximate $95\%$ confidence intervals for these odds ratios are $(1.296,2.266)$ and $(1.685,3.110)$, obtained by exponentiating the log-scale intervals. The intercept test against zero asks whether the reference risk equals one half, and is not a test of either exposure association. The binomial dispersion is assumed to be one because $\mathrm{Var}(Y_g)=n_g{\pi }_g(1-{\pi }_g)$; the output does not demonstrate independence or rule out extra-binomial variation.

There are four grouped probabilities and three fitted coefficients, leaving one residual degree of freedom. The fitted risks are approximately $(0.1495,0.2314,0.2869,0.4080)$. The binomial deviance compares this fit with the four unrestricted observed proportions:The deviance goodness-of-fit test against ${\chi }_1^2$ has p-value $0.0796$. The common-odds-ratio fit is not rejected at $5\%$, although the doubly exposed group's observed risk exceeds its fitted risk appreciably and suggests possible interaction.

The second formula adds $\delta AB$. Its four coefficients fit all four probabilities, making a saturated statistical model with zero grouped residual degrees of freedom. The near-zero deviance is roundoff from exact interpolation; it supplies no remaining grouped lack-of-fit test. Unlike Q1's unreplicated Gaussian fit, however, these grouped responses represent many independent Bernoulli trials, and their binomial variances are specified by their sizes and probabilities. Coefficient standard errors remain available from binomial information without estimating dispersion from the zero residual degrees of freedom.

Indeed, let ${\eta }_{ab}=\log [Y_{ab}/(n_{ab}-Y_{ab})]$ be the fitted cell logit link value. The saturated coefficients areThese calculations yield the printed estimates. The delta method gives $\mathrm{Var}({\hat{\eta }}_{ab})\approx 1/Y_{ab}+1/(n_{ab}-Y_{ab})$, and independent groups make the variance of $\hat{\delta }$ the sum of all eight reciprocal counts. Its square root is $0.37857$, reproducing the printed standard error.

The saturated logistic regression describes the conditional odds ratios as follows:Thus the main coefficients now concern effects at the other factor's reference level. Logistic interaction as a ratio of odds ratios gives $e^{\delta }=1.939$, with approximate $95\%$ confidence interval $(0.923,4.073)$. The interaction Wald test has $p=0.0802$, while the nested likelihood-ratio test uses deviance difference $3.0733$ on one statistical degree of freedom and has $p=0.0796$. The data suggest a larger odds association in predisposed subjects, but do not establish this interaction at the conventional $5\%$ level. This conclusion is specific to the odds scale; absolute-risk differences can vary even with no logit-scale interaction.

The Akaike information criterion is $-2\hat{\ell }+2k$. The additive fit has $k=3$ and AIC $31.766$; the interaction fit has $k=4$ and AIC $30.692$. The extra coefficient improves twice the log-likelihood by $3.0733$ but incurs penalty $2$, so AIC falls by only $1.0733$. It weakly favors the larger fit rather than decisively establishing an interaction, and AIC is not a p-value.

The marginal odds ratio $1.793$ and the fitted conditional odds ratios are different quantities. Predisposition is associated with both use frequency and outcome, so confounding is a substantive possibility; noncollapsibility of the odds ratio also means marginal and conditional odds ratios need not agree. This is an observational association analysis with only one adjustment factor. It does not by itself establish a causal effect, exclude other confounding, or validate the assumed independent binomial sampling.

Let $\Delta C$ and $\Delta E$ be the new-minus-standard mean cost and survival increments. The incremental cost-effectiveness ratio isUsing $365$ days per year, this is £14,600 per additional life-year. It is not cost per quality-adjusted life-year because no quality-of-life adjustment is supplied. Both point-estimated increments are positive, so the new intervention exchanges higher cost for longer survival; this is not dominance.

Assume the two trial arms are independent samples, their cost and survival means have a valid joint normal approximation, and the supplied standard errors and correlation coefficient estimate the joint uncertainty in their mean differences. This could follow from a central limit theorem for patient-level cost/survival pairs, with finite second moments and suitably handled survival follow-up. The supplied covariance matrix for $(\widehat{\Delta C},\widehat{\Delta E})$ isThe off-diagonal entry has units pounds times days, and must not be discarded. Independence of trial arms does not mean cost and survival increments are independent: patient-level cost and survival are paired within each arm.

For a proposed ratio $r$, the null assertion is $\Delta C-r\Delta E=0$. Its estimated variance is $90000-4500r+225r^2$. Fieller's theorem inverts the associated two-sided normal tests. With $q={\Phi }^{-1}(0.975)\approx 1.96$, retain $r$ whenExpanding givesThe leading coefficient is positive, since $45>q(15)$; the effect increment is separated from zero at this level. With $q^2=3.841459$, the polynomial is approximatelyIt is nonpositive between its two real roots. The resulting approximate $95\%$ Fieller confidence set is therefore a bounded confidence interval:On the $365$-day life-year scale the endpoints are approximately £9,133 and £36,376. Exact coverage would need an exact pivot with its corresponding critical value; using estimated covariance and a central limit theorem makes this an asymptotic interval. In other data, a poorly determined denominator can give unbounded or disconnected Fieller confidence sets, which should not be replaced by an artificially bounded interval.

An alternative is a paired bootstrap: independently resample patients within each randomized arm, keeping each patient's cost and survival together; recompute both mean differences and their ratio in each replicate. The $2.5\%$ and $97.5\%$ empirical ratio quantiles give a percentile bootstrap confidence interval. The raw paired data are needed to calculate it, so the numerical endpoints cannot be obtained from this summary table alone. Near-zero bootstrap effect increments or sign changes warn that ordinary ratio-percentile intervals may be unreliable; one should retain the joint cost/effect uncertainty or invert a suitable bootstrap test instead.

A simpler first-order alternative is the delta method. For $r=\Delta C/\Delta E$, the gradient is $(1/\Delta E,-\Delta C/\Delta E^2)$, givingThe resulting approximate Wald confidence interval is $40\pm 1.96\sqrt{133.333}=(17.37,62.63)$ pounds per day. Its symmetric first-order form misses the denominator's nonlinear effect, explaining its marked difference from the Fieller confidence set; it is not the preferred ratio interval here.

The incremental cost-effectiveness ratio must be read with the signs of both increments. A negative ratio can mean either greater effectiveness at lower cost, which is favorable, or lower effectiveness at greater cost, which is unfavorable. A positive ratio can likewise lie in either the more-costly/more-effective or cheaper/less-effective quadrant. Thus ordering ratio values without locating the cost-effectiveness plane quadrant can reverse a decision.

The incremental cost-effectiveness ratio is undefined at a zero effect increment and extremely unstable near it. Its sampling distribution can be highly skewed and its mean may fail to exist; a usual symmetric confidence interval or an interval spanning positive and negative ratios can conceal that instability. Cost-effect correlation must be retained, as the Fieller confidence set calculation illustrates.

Interpretation also requires a specified comparator, effect unit, costing perspective, follow-up horizon and willingness-to-pay threshold. Different comparator choices give different incremental ratios, and lower ratios alone do not settle decisions across differently scaled outcomes. Longer survival need not mean a proportional quality-of-life gain; extrapolation, censored outcomes and omitted costs can affect the increments. Point estimates conceal both sampling uncertainty and uncertainty in these assumptions.

For a positive effect increment, comparison with a threshold $\lambda$ is equivalent to $\Delta C/\Delta E<\lambda$. For negative effect increments that inequality reverses when multiplied by $\Delta E$, so a rule that blindly compares all ratios with $\lambda$ is wrong. Use the signs and joint cost/effect uncertainty, or the incremental net monetary benefit $\lambda \Delta E-\Delta C$, rather than the ratio alone. No particular threshold is given in the paper, so the point estimate does not establish whether reimbursement should occur.

One presentation is the cost-effectiveness plane, with $\Delta E$ horizontally and $\Delta C$ vertically. Plot the point estimate $(45,1800)$ with a joint confidence region or paired bootstrap cloud. Its position shows the tradeoff quadrant, while the cloud shows uncertainty, dependence and possible dominance. A willingness-to-pay threshold $\lambda$ is the line $\Delta C=\lambda \Delta E$ through the origin. Points below that line have positive incremental net monetary benefit, including points in quadrants where a ratio rule would fail. The point-estimated switch occurs at $\lambda =40$ pounds per day.

A second presentation is a cost-effectiveness acceptability curve, giving an uncertainty-based measure of cost-effectiveness for a range of thresholds. For paired bootstrap increments $(\Delta E^{\ast (b)},\Delta C^{\ast (b)})$, plotThis counts positive net benefits, not merely ratios below the threshold, so it handles all quadrants. Using Bayesian inference the curve is the posterior probability that incremental net monetary benefit is positive; a frequentist bootstrap curve is a resampling summary, not literally a probability that a fixed parameter is positive.

Here the estimated incremental net monetary benefit and its standard error areThus a normal approximation uncertainty curve isIt passes through one half at £40 per day. Alternatively display the incremental net monetary benefit curve itself with approximate pointwise $95\%$ bands $45\lambda -1800\pm 1.96s_{\mathrm{INMB}}(\lambda )$. The plane preserves the joint uncertainty; the acceptability curve shows how the threshold changes the cost-effectiveness assessment. Neither presentation chooses society's willingness-to-pay threshold.