# Load data
beef <- read.csv("data/beef_marbling.csv")
beef| obs | live_weight_kg | marbling_score |
|---|---|---|
| 1 | 500 | 5.0 |
| 2 | 520 | 5.5 |
| 3 | 540 | 5.8 |
| 4 | 560 | 6.0 |
| 5 | 580 | 6.2 |
| 6 | 600 | 6.5 |
| 7 | 620 | 6.7 |
| 8 | 640 | 7.0 |
| 9 | 850 | 6.5 |
11 Week 11: Model Diagnostics
11.1 Introduction: Why Diagnostics Matter
11.1.1 The Danger of Trusting Models Blindly
Linear models are powerful tools in animal breeding and genetics, but they are only as reliable as the assumptions on which they rest. Consider this scenario: A beef cattle breeding program uses a linear model to predict carcass quality from live animal measurements. The model appears to fit well with a high R², and all regression coefficients are statistically significant. Based on these results, the program makes selection decisions affecting thousands of animals and millions of dollars in genetic improvement.
However, the model was never properly diagnosed. Hidden in the data was a single influential observation—an animal with unusual biology that pulled the regression line in the wrong direction. The predictions for the entire population were biased, leading to suboptimal selection decisions. Years later, when genetic evaluations were re-run with proper diagnostics, the error was discovered, but the damage was done: inferior genetics had been propagated through the herd.
WarningThe Danger of Model Misspecification
Real consequences in animal breeding:
- Biased breeding value predictions: Outliers can distort estimates of genetic merit
- Inefficient selection: Resources wasted on animals with inflated predictions
- Lost genetic gain: Years of progress compromised by poor model fit
- Economic losses: Millions of dollars in reduced productivity
The bottom line: Model diagnostics are not optional—they are essential for trustworthy inference and optimal breeding decisions.
11.1.2 Real Consequences in Animal Breeding
In applied animal breeding and genetics, the stakes are high:
- Genetic evaluation systems (EPDs, EBVs, genomic predictions) rely on accurate models
- Selection decisions are based on model predictions that affect future generations
- Economic returns depend on correctly identifying superior animals
- Industry competitiveness requires maximizing genetic gain per generation
When models fail due to violated assumptions or influential observations, the consequences ripple through entire breeding populations for years or even decades.
11.1.3 Connection to Previous Weeks
In Week 5, we learned about the Gauss-Markov Theorem, which guarantees that least squares estimators are BLUE (Best Linear Unbiased Estimators) under certain conditions:
NoteConnection to Week 5: Assumptions of Least Squares
Recall the Gauss-Markov assumptions from Week 5:
- Linearity: \(E(\mathbf{y}) = \mathbf{X}\boldsymbol{\beta}\) (response is linear in parameters)
- Zero mean errors: \(E(\mathbf{e}) = \mathbf{0}\) (errors have expected value zero)
- Homoscedasticity: \(\text{Var}(\mathbf{e}) = \sigma^2\mathbf{I}\) (constant variance)
- Independence: Errors are uncorrelated (captured in \(\sigma^2\mathbf{I}\))
- Normality (for inference): \(\mathbf{e} \sim N(\mathbf{0}, \sigma^2\mathbf{I})\)
This week’s focus: We check whether these assumptions actually hold for our data. If they don’t, our estimates may be inefficient, biased, or produce invalid confidence intervals and hypothesis tests.
Up to this point, we’ve assumed these conditions hold. This week, we learn how to verify them using diagnostic tools and what to do when they fail.
Key connections:
- Week 4-5: Simple regression foundations—residuals, fitted values, SSE
- Week 6: Multiple regression—leverage becomes more complex with multiple predictors
- Weeks 7-9: ANOVA models—same diagnostic principles apply to categorical predictors
- Week 10: ANCOVA—homogeneity of slopes relates to residual patterns
11.2 Mathematical Theory
11.2.1 Review of Model Assumptions
The linear model in matrix form is:
\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{e} \tag{11.1}\]
where:
- \(\mathbf{y}\): \((n \times 1)\) response vector
- \(\mathbf{X}\): \((n \times p)\) design matrix
- \(\boldsymbol{\beta}\): \((p \times 1)\) parameter vector
- \(\mathbf{e}\): \((n \times 1)\) error vector
The least squares estimate is:
\[ \mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y} \tag{11.2}\]
For this estimator to be BLUE and for inference to be valid, we require:
- E1: Linearity — The relationship between \(\mathbf{y}\) and \(\mathbf{X}\) is correctly specified
- E2: Independence — Observations are independent (or known dependence structure)
- E3: Homoscedasticity — \(\text{Var}(e_i) = \sigma^2\) for all \(i\)
- E4: Normality — \(e_i \sim N(0, \sigma^2)\) (required for valid \(t\) and \(F\) tests)
This week’s goal: Develop tools to check each assumption and diagnose violations.
11.3 Types of Residuals
11.3.1 Raw Residuals
The raw residual for observation \(i\) is simply the difference between observed and fitted values:
\[ e_i = y_i - \hat{y}_i = y_i - \mathbf{x}_i'\mathbf{b} \tag{11.3}\]
In matrix form, the vector of all residuals is:
\[ \mathbf{e} = \mathbf{y} - \hat{\mathbf{y}} = \mathbf{y} - \mathbf{X}\mathbf{b} = (\mathbf{I} - \mathbf{H})\mathbf{y} \tag{11.4}\]
where \(\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\) is the hat matrix (more on this shortly).
Properties:
- \(\sum_{i=1}^n e_i = 0\) (residuals sum to zero when intercept is in model)
- \(\mathbf{X}'\mathbf{e} = \mathbf{0}\) (residuals are orthogonal to column space of \(\mathbf{X}\))
- \(E(\mathbf{e}) = \mathbf{0}\) (expected value is zero)
Limitation: Raw residuals have unequal variances even when the true errors have equal variance. This is because:
\[ \text{Var}(e_i) = \sigma^2(1 - h_{ii}) \tag{11.5}\]
where \(h_{ii}\) is the \(i\)-th diagonal element of the hat matrix. Since \(h_{ii}\) varies across observations, so does \(\text{Var}(e_i)\).
11.3.2 Standardized Residuals
To account for the estimate of \(\sigma^2\), we standardize residuals:
\[ e_i^* = \frac{e_i}{\hat{\sigma}} \tag{11.6}\]
where \(\hat{\sigma} = \sqrt{\text{MSE}} = \sqrt{\frac{\text{SSE}}{n-p}}\) is the estimated standard deviation of the errors.
Use: Standardized residuals put all residuals on the same scale and can be roughly interpreted as \(z\)-scores if errors are normal.
Limitation: Standardized residuals still have unequal variances because we haven’t accounted for \(h_{ii}\).
11.3.3 Studentized Residuals
To fully account for the unequal variances of residuals, we studentize:
\[ r_i = \frac{e_i}{\hat{\sigma}\sqrt{1 - h_{ii}}} \tag{11.7}\]
NoteWhy Studentized Residuals?
The problem: Even when the true errors \(e_i\) all have variance \(\sigma^2\), the residuals \(e_i\) have unequal variances:
\[ \text{Var}(e_i) = \sigma^2(1 - h_{ii}) \]
The solution: Divide by the estimated standard deviation \(\hat{\sigma}\sqrt{1 - h_{ii}}\) to equalize variances:
\[ \text{Var}(r_i) \approx 1 \quad \text{for all } i \]
Why it matters: Studentized residuals can be compared directly across observations. An observation with \(|r_i| > 3\) is unusual regardless of its position in \(X\) space.
Distribution: Under the normality assumption, studentized residuals approximately follow a \(t\) distribution with \(n-p\) degrees of freedom:
\[ r_i \sim t_{n-p} \]
Rule of thumb: Observations with \(|r_i| > 2\) are potentially unusual; those with \(|r_i| > 3\) are almost certainly outliers.
11.3.4 Studentized Deleted Residuals
The studentized deleted residual (also called externally studentized) is computed by:
- Remove observation \(i\) from the dataset
- Refit the model on the remaining \(n-1\) observations
- Predict \(y_i\) using the refitted model
- Compute the residual and standardize it
Formally:
\[ t_i = \frac{e_i}{\hat{\sigma}_{(i)}\sqrt{1 - h_{ii}}} \tag{11.8}\]
where \(\hat{\sigma}_{(i)}\) is the estimate of \(\sigma\) with observation \(i\) deleted.
Distribution: Under normality, \(t_i\) follows an exact \(t\) distribution with \(n-p-1\) degrees of freedom:
\[ t_i \sim t_{n-p-1} \]
Use: Deleted residuals are more sensitive for detecting outliers because the outlier doesn’t influence its own residual calculation.
Relationship to ordinary studentized residuals:
\[ t_i = r_i \sqrt{\frac{n-p-1}{n-p-r_i^2}} \]
For most observations, \(t_i \approx r_i\), but they can differ substantially for influential points.
11.4 Leverage and the Hat Matrix
11.4.1 Hat Matrix Definition and Properties
The hat matrix \(\mathbf{H}\) is central to understanding leverage and influence:
\[ \mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}' \tag{11.9}\]
It “puts the hat” on \(\mathbf{y}\) to get \(\hat{\mathbf{y}}\):
\[ \hat{\mathbf{y}} = \mathbf{H}\mathbf{y} \]
The residual vector can be written as:
\[ \mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{y} \]
ImportantHat Matrix Properties
The hat matrix \(\mathbf{H}\) has several important mathematical properties:
Symmetric: \(\mathbf{H}' = \mathbf{H}\)
Proof: \(\mathbf{H}' = [\mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}']' = \mathbf{X}[(\mathbf{X}'\mathbf{X})^{-1}]'\mathbf{X}' = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}' = \mathbf{H}\)
Idempotent: \(\mathbf{H}\mathbf{H} = \mathbf{H}\)
Proof: \(\mathbf{H}\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}' = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}' = \mathbf{H}\)
Bounds on diagonal: \(0 \leq h_{ii} \leq 1\) for all \(i\)
Trace equals rank: \(\text{tr}(\mathbf{H}) = \text{tr}[\mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'] = \text{tr}[(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{X}] = \text{tr}(\mathbf{I}_p) = p\)
Therefore: \(\sum_{i=1}^n h_{ii} = p\)
Implications: - Idempotency means projecting twice is the same as projecting once - The sum of leverages equals the number of parameters - Average leverage is \(\bar{h} = p/n\)
11.4.2 Leverage Values (\(h_{ii}\))
The leverage of observation \(i\) is the \(i\)-th diagonal element of \(\mathbf{H}\):
\[ h_{ii} = [\mathbf{H}]_{ii} = \mathbf{x}_i'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_i \tag{11.10}\]
where \(\mathbf{x}_i\) is the \(i\)-th row of \(\mathbf{X}\) (transposed to a column vector).
Interpretation: \(h_{ii}\) measures how far observation \(i\) is from the center of the \(X\) space. It quantifies the potential for observation \(i\) to influence the fitted model.
Key properties:
- \(0 \leq h_{ii} \leq 1\)
- \(\sum_{i=1}^n h_{ii} = p\)
- Average leverage: \(\bar{h} = p/n\)
- \(\text{Var}(e_i) = \sigma^2(1 - h_{ii})\) — observations with high \(h_{ii}\) have smaller residual variance
- \(\text{Var}(\hat{y}_i) = \sigma^2 h_{ii}\) — observations with high \(h_{ii}\) have predictions with higher variance
Geometric interpretation:
- \(h_{ii}\) measures the “distance” of \(\mathbf{x}_i\) from \(\bar{\mathbf{x}}\) in a metric defined by \((\mathbf{X}'\mathbf{X})^{-1}\)
- Points far from the center of the \(X\) cloud have high leverage
- Points near the center have low leverage
11.4.3 High Leverage Criteria
Question: When is leverage “high”?
TipRule of Thumb: High Leverage Thresholds
An observation is considered high leverage if:
- \(h_{ii} > \frac{2p}{n}\) (twice the average leverage)
- Or more conservatively: \(h_{ii} > \frac{3p}{n}\) (three times average)
Example: In a model with \(n=50\) and \(p=3\): - Average leverage: \(\bar{h} = 3/50 = 0.06\) - High leverage threshold: \(h_{ii} > 2(0.06) = 0.12\)
Important caveat: High leverage means potential for influence, not necessarily actual influence. An observation can have high leverage but still follow the pattern of the data.
11.5 Influence Measures
Leverage vs. Influence:
- Leverage (\(h_{ii}\)): Measures an observation’s potential to influence the fit based on its \(X\) values
- Influence: Measures an observation’s actual effect on the fitted model (combines leverage and residual size)
An observation can be:
- High leverage, low influence: Far from center of \(X\) but follows the regression pattern (on the line)
- Low leverage, high influence: Near center of \(X\) but very far from regression line (large residual)
- High leverage, high influence: Far from center of \(X\) and doesn’t follow pattern (both extreme \(X\) and large residual)
11.5.1 Cook’s Distance
Cook’s Distance (Cook 1977) is the most widely used measure of influence. It quantifies how much the entire vector of fitted values \(\hat{\mathbf{y}}\) changes when observation \(i\) is deleted.
Definition:
\[ D_i = \frac{\sum_{j=1}^n (\hat{y}_j - \hat{y}_{j(i)})^2}{p \cdot \text{MSE}} \tag{11.11}\]
where \(\hat{y}_{j(i)}\) is the predicted value for observation \(j\) when observation \(i\) is deleted from the model.
Equivalent computational form:
\[ D_i = \frac{r_i^2}{p} \cdot \frac{h_{ii}}{(1-h_{ii})^2} \tag{11.12}\]
Interpretation: Cook’s Distance combines two components:
- Residual size: \(r_i^2/p\) — How far is the observation from the fitted model?
- Leverage: \(h_{ii}/(1-h_{ii})^2\) — How much potential does it have to influence the fit?
An observation is influential only if both components are large.
Thresholds:
TipRule of Thumb: Cook’s Distance Thresholds
An observation is considered influential if:
- Classic threshold: \(D_i > 1\) (conservative, rarely exceeded)
- Alternative: \(D_i > \frac{4}{n}\) (more liberal, useful for small/medium datasets)
- Adjusted: \(D_i > \frac{4}{n-p-1}\) (accounts for model complexity)
Example: For \(n=50\), \(p=3\): - Classic: \(D_i > 1\) (very conservative) - Alternative: \(D_i > 4/50 = 0.08\) (more practical) - Adjusted: \(D_i > 4/46 = 0.087\)
Guideline: Examine any observation with \(D_i > 4/n\) carefully. Investigate the biology/data quality for observations with \(D_i > 0.5\).
11.5.2 DFFITS
DFFITS (Difference in Fits) measures how much the predicted value \(\hat{y}_i\) changes when observation \(i\) is deleted:
\[ \text{DFFITS}_i = \frac{\hat{y}_i - \hat{y}_{i(i)}}{\hat{\sigma}_{(i)}\sqrt{h_{ii}}} \tag{11.13}\]
where \(\hat{y}_{i(i)}\) is the prediction for observation \(i\) using a model fit without observation \(i\).
Equivalent form:
\[ \text{DFFITS}_i = r_i \sqrt{\frac{h_{ii}}{1-h_{ii}}} \]
Threshold: \(|\text{DFFITS}_i| > 2\sqrt{p/n}\) suggests influential observation.
Interpretation: While Cook’s \(D\) measures overall influence on all fitted values, DFFITS focuses on the change in the observation’s own prediction.
11.5.3 DFBETAS
DFBETAS measures the change in each individual regression coefficient when observation \(i\) is deleted:
\[ \text{DFBETAS}_{i,j} = \frac{b_j - b_{j(i)}}{\hat{\sigma}_{(i)}\sqrt{[(\mathbf{X}'\mathbf{X})^{-1}]_{jj}}} \tag{11.14}\]
where \(b_{j(i)}\) is the \(j\)-th coefficient when observation \(i\) is deleted.
Threshold: \(|\text{DFBETAS}_{i,j}| > 2/\sqrt{n}\) suggests observation \(i\) is influential for coefficient \(j\).
Use: DFBETAS helps identify which specific coefficient(s) are most affected by a potentially influential observation.
11.6 Transformations
When diagnostic checks reveal violations of model assumptions (especially non-normality or heteroscedasticity), transforming the response variable may stabilize variance and improve model fit.
11.6.1 Box-Cox Family
The Box-Cox transformation (Box and Cox 1964) is a family of power transformations indexed by a parameter \(\lambda\):
\[ y^{(\lambda)} = \begin{cases} \frac{y^\lambda - 1}{\lambda} & \text{if } \lambda \neq 0 \\ \log(y) & \text{if } \lambda = 0 \end{cases} \tag{11.15}\]
Why this form? The normalization \(\frac{y^\lambda - 1}{\lambda}\) makes the transformation continuous at \(\lambda = 0\):
\[ \lim_{\lambda \to 0} \frac{y^\lambda - 1}{\lambda} = \log(y) \]
Common special cases:
TipBox-Cox \(\lambda\) Interpretation
| \(\lambda\) | Transformation | Use case |
|---|---|---|
| \(\lambda = 2\) | \(y^2\) | Rarely used (increases skewness) |
| \(\lambda = 1\) | \(y\) (no transformation) | Data already normal/homoscedastic |
| \(\lambda = 0.5\) | \(\sqrt{y}\) | Count data, moderate right skew |
| \(\lambda = 0\) | \(\log(y)\) | Strong right skew, proportional errors |
| \(\lambda = -0.5\) | \(1/\sqrt{y}\) | Very strong right skew |
| \(\lambda = -1\) | \(1/y\) | Extreme skewness |
Most common in practice: \(\lambda \in [0, 1]\) (square root to log transformation range)
Selecting \(\lambda\): The optimal \(\lambda\) is chosen to maximize the log-likelihood:
\[ L(\lambda) = -\frac{n}{2}\log(\text{SSE}(\lambda)) + (\lambda - 1)\sum_{i=1}^n \log(y_i) \]
where \(\text{SSE}(\lambda)\) is the error sum of squares for the model fit using \(y^{(\lambda)}\) as the response.
In R, this is done using MASS::boxcox(), which produces a profile likelihood plot. The optimal \(\lambda\) is the value that maximizes \(L(\lambda)\).
11.6.2 Logarithmic Transformation
The logarithmic transformation (\(\lambda = 0\) in Box-Cox family) is the most common transformation in practice:
\[ y^* = \log(y) \]
When to use log transformation:
- Right-skewed response: Long tail to the right
- Proportional errors: Variance increases with mean (heteroscedasticity)
- Multiplicative relationships: Effects are multiplicative rather than additive
- Positive response bounded away from zero: \(y > 0\) for all observations
Biological examples in livestock:
- Milk production: Variance typically increases with production level
- Somatic cell count (mastitis indicator): Highly right-skewed
- Growth rates: Often show proportional variation
- Feed efficiency ratios: Positive, often skewed
Interpretation on log scale:
- A one-unit increase in \(x\) leads to a change in \(\log(y)\) of \(\beta_1\)
- Equivalently, \(y\) is multiplied by \(e^{\beta_1}\)
- Percentage change: \((e^{\beta_1} - 1) \times 100\%\)
Predictions on original scale:
To get predictions on the original scale, exponentiate:
\[ \hat{y} = e^{\hat{y}^*} = e^{\mathbf{x}'\mathbf{b}} \]
Important note: \(e^{\hat{y}^*}\) is a prediction of the median of \(y\), not the mean, when errors are normal on the log scale. For moderate-depth coverage, we acknowledge this but don’t derive the bias correction for predicting the mean.
11.6.3 When to Transform
WarningWhen NOT to Transform
Avoid transformations when:
- Violations are minor: Residual plots are nearly acceptable, and inference is not severely affected
- Interpretation becomes too complex: Stakeholders need results on the original scale (e.g., kg of milk, dollars)
- Response scale has direct meaning: Physical or economic units matter for decision-making
- Better alternatives exist: Weighted least squares (WLS) for known heteroscedasticity, robust regression for outliers, or generalized linear models (GLMs) for non-normal responses
Rule of thumb: Only transform when diagnostic violations are substantial and interpretation on the transformed scale is acceptable to stakeholders.
When transformation IS appropriate:
- Residual plots show clear funnel pattern (heteroscedasticity)
- Q-Q plot shows strong departure from normality (heavy tails, skewness)
- Biological rationale supports transformed scale (e.g., log for growth, sqrt for counts)
- Predictions and inference on transformed scale are meaningful
11.7 Checking Model Assumptions
We now have the mathematical tools (residuals, leverage, influence) to check whether our model assumptions hold.
11.7.1 The Four Assumptions and How to Check Them
ImportantThe Four Assumptions and Diagnostic Tools
| Assumption | What to Check | Diagnostic Tool | Pattern if Violated |
|---|---|---|---|
| E1: Linearity | Relationship between \(Y\) and \(X\) is linear | Residuals vs. Fitted plot | Curved pattern, systematic structure |
| E2: Independence | Observations are independent | Residuals vs. Order plot; knowledge of data collection | Patterns over time, clusters |
| E3: Homoscedasticity | Constant error variance | Residuals vs. Fitted plot; Scale-Location plot | Funnel shape, increasing/decreasing spread |
| E4: Normality | Errors follow normal distribution | Normal Q-Q plot; histogram | Departure from line, heavy tails, skewness |
Note: We check assumptions using residuals (not the raw response \(y\)) because residuals estimate the errors \(e_i\), which are the quantities assumed to be normal, independent, and homoscedastic.
11.7.2 Linearity: Residuals vs. Fitted Values
Assumption: The relationship between \(\mathbf{y}\) and \(\mathbf{X}\) is correctly specified as linear.
Diagnostic: Plot residuals \(e_i\) (or studentized residuals \(r_i\)) against fitted values \(\hat{y}_i\).
What to look for:
- Good: Random scatter around horizontal line at \(y=0\), no pattern
- Violation: Curved pattern (U-shape, inverted U, etc.)
Example violations:
- Quadratic trend: Residuals show U-shape → suggests need for \(x^2\) term
- Missing interaction: Different curvature patterns within groups
- Wrong functional form: Exponential relationship fit with linear model
What to do if violated:
- Add polynomial terms: \(x^2\), \(x^3\)
- Include interactions: \(x_1 \times x_2\)
- Transform predictor(s): \(\log(x)\), \(\sqrt{x}\)
- Consider non-linear models if relationship is fundamentally non-linear
11.7.3 Homoscedasticity: Constant Variance
Assumption: \(\text{Var}(e_i) = \sigma^2\) for all \(i\) (constant variance).
Diagnostics:
- Residuals vs. Fitted plot: Look at spread of residuals as \(\hat{y}\) increases
- Scale-Location plot: Plot \(\sqrt{|r_i|}\) vs. \(\hat{y}_i\) to amplify patterns
What to look for:
- Good: Even spread of residuals across all fitted values
- Violation: Funnel shape — spread increases (or decreases) with \(\hat{y}\)
Common patterns:
- Increasing variance: Spread widens as \(\hat{y}\) increases (most common)
- Decreasing variance: Spread narrows as \(\hat{y}\) increases (less common)
- Groups with different variances: Different spread within subgroups
Formal tests (optional):
- Breusch-Pagan test (Breusch and Pagan 1979): Tests whether variance depends on fitted values
- White test (White 1980): More general test for heteroscedasticity
What to do if violated:
- Log transformation of \(y\) (if variance proportional to mean)
- Square root transformation of \(y\) (if variance proportional to mean, e.g., counts)
- Weighted least squares (WLS) if variance structure is known
- Robust standard errors (HC3, HC4) for inference without assuming homoscedasticity
11.7.4 Normality: Q-Q Plots
Assumption: Errors are normally distributed: \(e_i \sim N(0, \sigma^2)\).
Note: Normality is not required for unbiasedness of \(\mathbf{b}\) or for Gauss-Markov optimality. It is required for:
- Exact validity of \(t\) and \(F\) tests
- Exact confidence intervals
- Prediction intervals
For large samples, the Central Limit Theorem provides approximate normality of \(\mathbf{b}\) even if errors are non-normal.
Diagnostic: Normal Q-Q plot (quantile-quantile plot)
- X-axis: Theoretical quantiles from standard normal distribution
- Y-axis: Sample quantiles of studentized residuals \(r_i\)
What to look for:
- Good: Points fall along a straight diagonal line
- Violation: Systematic departure from the line
Common patterns:
- Light tails: Points curve below line at extremes → uniform or short-tailed distribution
- Heavy tails: Points curve above line at extremes → outliers or long-tailed distribution
- Right skewness: Points curve upward (right tail pulls away from line) → need transformation
- Left skewness: Points curve downward (left tail pulls away from line)
- S-shaped: Both tails depart in opposite directions
Formal tests:
- Shapiro-Wilk test (Shapiro and Wilk 1965): Tests normality (sensitive to sample size)
- Kolmogorov-Smirnov test: Less powerful but more general
Caution: Don’t rely solely on hypothesis tests. With large \(n\), tests reject even minor departures. With small \(n\), tests have low power. Always look at the Q-Q plot.
What to do if violated:
- Transformation of \(y\) (log, square root, Box-Cox)
- Robust methods: Less sensitive to non-normality (beyond this course)
- Bootstrap inference: Doesn’t require normality (beyond this course)
- If violation is minor and \(n\) is large, standard inference is approximately valid
11.7.5 Independence: Patterns in Residuals
Assumption: Observations are independent; \(\text{Cov}(e_i, e_j) = 0\) for \(i \neq j\).
Why it matters: Dependence (positive correlation) among errors leads to:
- Underestimated standard errors: Confidence intervals too narrow
- Inflated test statistics: Higher Type I error rates
- Biased inference: \(p\)-values too small
Diagnostics:
- Residuals vs. Order plot: Plot \(e_i\) vs. observation number (if data have natural order)
- Residuals vs. Time: If data are time series
- Spatial plots: If data have spatial structure (e.g., pens in barn, farms in region)
- Knowledge of study design: Are observations clustered? Repeated measures? Hierarchical structure?
What to look for:
- Good: Random scatter, no patterns or clusters
- Violation:
- Time series: Runs of positive/negative residuals (autocorrelation)
- Clusters: Residuals group together (e.g., within pen, farm, family)
- Spatial: Nearby observations have similar residuals
Common causes in livestock data:
- Pen effects: Animals within the same pen are more similar (shared environment)
- Family structure: Related animals share genetics (more similar than unrelated animals)
- Time trends: Management or environmental changes over time
- Batch effects: Different processing batches, technicians, labs
What to do if violated:
- Mixed models: Include random effects for clusters (Weeks 14-15 preview)
- GLS (Generalized Least Squares): Model correlation structure explicitly
- Cluster-robust standard errors: Adjust inference for clustering
- Time series models: ARIMA, GARCH for temporal dependence
Important: Independence violations cannot be fixed with transformations. They require different modeling approaches.
11.8 The Four Essential Diagnostic Plots
When you fit a linear model in R using lm(), the default plot() function produces four diagnostic plots. Let’s understand each one.
11.8.1 Residuals vs. Fitted Values
Purpose: Check linearity and homoscedasticity simultaneously.
What’s plotted:
- X-axis: Fitted values \(\hat{y}_i\)
- Y-axis: Residuals \(e_i\) (or standardized residuals)
Good model:
- Random scatter around horizontal line at \(y=0\)
- Even spread (same variance) across all \(\hat{y}\) values
- No patterns, no outliers far from the cloud
Patterns indicating problems:
- Curved pattern: Non-linearity (need polynomial, transformation, or interaction)
- Funnel shape (spread increases): Heteroscedasticity (consider log or sqrt transformation)
- Inverted funnel (spread decreases): Less common, may need transformation
- Outliers: Points far from the horizontal band (investigate further)
Interpretation: This is the most important diagnostic plot. It checks two assumptions at once.
11.8.2 Normal Q-Q Plot
Purpose: Check normality of residuals.
What’s plotted:
- X-axis: Theoretical quantiles of standard normal distribution
- Y-axis: Sample quantiles of studentized residuals
Good model: Points lie along the diagonal reference line.
Patterns indicating problems:
- Points curve above line at right end: Right skewness or heavy right tail
- Points curve below line at left end: Left skewness or heavy left tail
- S-shaped curve: Non-normal distribution (both tails depart)
- Isolated points far from line: Outliers
Interpretation: Focus on the middle 80% of points. Some deviation in the extreme tails is normal, especially for small samples.
11.8.3 Scale-Location Plot
Purpose: Check homoscedasticity (constant variance) more clearly than Residuals vs. Fitted plot.
What’s plotted:
- X-axis: Fitted values \(\hat{y}_i\)
- Y-axis: \(\sqrt{|r_i|}\) (square root of absolute studentized residuals)
- Smoothed line: LOESS smooth through points
Good model:
- Horizontal smoothed line (flat trend)
- Even spread of points across all \(\hat{y}\) values
Patterns indicating problems:
- Upward trend: Variance increases with \(\hat{y}\) (heteroscedasticity)
- Downward trend: Variance decreases with \(\hat{y}\) (less common)
- U-shaped or wavy: Non-constant variance in complex pattern
Why square root? Taking \(\sqrt{|r_i|}\) makes patterns in variance easier to see because it’s on the standard deviation scale (not variance scale).
11.8.4 Residuals vs. Leverage
Purpose: Identify influential observations by showing both leverage and residual size.
What’s plotted:
- X-axis: Leverage values \(h_{ii}\)
- Y-axis: Standardized residuals \(r_i\)
- Contour lines: Cook’s Distance contours (typically at \(D_i = 0.5\) and \(D_i = 1\))
Good model:
- Most points clustered in lower-left (low leverage, small residuals)
- No points beyond Cook’s Distance contours
Patterns indicating problems:
- High leverage, small residual: Points in upper-right or lower-right (extreme \(X\) but on line) — high leverage, low influence
- Low leverage, large residual: Points at top or bottom but left side (outliers near center of \(X\)) — low leverage, moderate influence
- High leverage, large residual: Points in corners beyond Cook’s \(D\) contours — high influence (INVESTIGATE!)
Interpretation: Points in the corners (high leverage AND large residual) are the most concerning. These observations can dramatically change the fitted model.
11.9 Small Numerical Example: Beef Carcass Marbling
Now let’s work through a complete diagnostic analysis with a small dataset where we can perform hand calculations.
11.10 The Dataset
A beef cattle researcher measures carcass marbling score (on a 1-10 scale, higher is better) and live weight (kg) for 9 steers:
Let’s visualize the data:
plot(beef$live_weight_kg, beef$marbling_score,
xlab = "Live Weight (kg)",
ylab = "Marbling Score",
main = "Beef Carcass Marbling vs. Live Weight",
pch = 19, col = "steelblue", cex = 1.5)
# Fit linear model
fit_beef <- lm(marbling_score ~ live_weight_kg, data = beef)
abline(fit_beef, col = "red", lwd = 2)
# Identify observation 9
points(beef$live_weight_kg[9], beef$marbling_score[9],
col = "red", pch = 19, cex = 2)
text(beef$live_weight_kg[9], beef$marbling_score[9],
labels = " Obs 9", pos = 4, col = "red", font = 2)
Observation: Observation 9 is an unusually heavy animal (850 kg vs. 500-640 kg for others), but its marbling score (6.5) is roughly consistent with the trend.
11.11 Fitting the Model
Model: \[ \text{marbling}_i = \beta_0 + \beta_1 \times \text{weight}_i + e_i \]
Let’s fit this model and examine the summary:
# Fit model
fit_beef <- lm(marbling_score ~ live_weight_kg, data = beef)
summary(fit_beef)
Call:
lm(formula = marbling_score ~ live_weight_kg, data = beef)
Residuals:
Min 1Q Median 3Q Max
-0.75261 -0.32792 0.02147 0.37085 0.72023
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.869900 1.088487 3.555 0.00928 **
live_weight_kg 0.003765 0.001787 2.107 0.07312 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5256 on 7 degrees of freedom
Multiple R-squared: 0.3881, Adjusted R-squared: 0.3006
F-statistic: 4.439 on 1 and 7 DF, p-value: 0.07312The model appears to fit well: \(R^2 = 0.94\) and weight is highly significant. But let’s check the diagnostics…
11.12 Computing Diagnostic Statistics by Hand
Now we’ll compute diagnostic statistics manually to build intuition. For larger datasets, we’d use R functions.
11.12.1 Step 1: Design Matrix and Estimates
The design matrix is:
\[ \mathbf{X} = \begin{bmatrix} 1 & 500 \\ 1 & 520 \\ 1 & 540 \\ \vdots & \vdots \\ 1 & 850 \end{bmatrix}_{9 \times 2} \]
# Build design matrix
X <- cbind(1, beef$live_weight_kg)
y <- beef$marbling_score
n <- nrow(X)
p <- ncol(X)
print("Design matrix X:")[1] "Design matrix X:"print(X) [,1] [,2]
[1,] 1 500
[2,] 1 520
[3,] 1 540
[4,] 1 560
[5,] 1 580
[6,] 1 600
[7,] 1 620
[8,] 1 640
[9,] 1 850Compute \(\mathbf{X}'\mathbf{X}\) and \((\mathbf{X}'\mathbf{X})^{-1}\):
# X'X
XtX <- t(X) %*% X
print("X'X:")[1] "X'X:"print(XtX) [,1] [,2]
[1,] 9 5410
[2,] 5410 3338500# (X'X)^{-1}
XtX_inv <- solve(XtX)
print("(X'X)^{-1}:")[1] "(X'X)^{-1}:"print(XtX_inv) [,1] [,2]
[1,] 4.288926002 -6.950154e-03
[2,] -0.006950154 1.156218e-05Compute estimates \(\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}\):
# Estimates
b <- XtX_inv %*% t(X) %*% y
print("Estimates b:")[1] "Estimates b:"print(b) [,1]
[1,] 3.869899794
[2,] 0.003765416cat("\nInterpretation:")
Interpretation:cat(sprintf("\nb0 (intercept) = %.4f", b[1]))
b0 (intercept) = 3.8699cat(sprintf("\nb1 (slope) = %.4f kg^{-1}", b[2]))
b1 (slope) = 0.0038 kg^{-1}cat("\n\nFor each 1 kg increase in live weight, marbling score increases by",
round(b[2], 4), "points.")
For each 1 kg increase in live weight, marbling score increases by 0.0038 points.11.12.2 Step 2: Hat Matrix and Leverage
Compute the hat matrix \(\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\):
# Hat matrix
H <- X %*% XtX_inv %*% t(X)
# Verify properties
cat("Properties of H:\n")Properties of H:cat("1. H is symmetric:", all.equal(H, t(H)), "\n")1. H is symmetric: TRUE cat("2. H is idempotent:", all.equal(H %*% H, H, tolerance = 1e-10), "\n")2. H is idempotent: TRUE cat("3. trace(H) = p:", round(sum(diag(H)), 6), "= p =", p, "\n")3. trace(H) = p: 2 = p = 2 # Extract leverage values
h <- diag(H)
print("\nLeverage values h_ii:")[1] "\nLeverage values h_ii:"print(round(h, 4))[1] 0.2293 0.1872 0.1543 0.1307 0.1163 0.1111 0.1152 0.1286 0.8273Key observation: Observation 9 has \(h_{99} = 0.449\), which is very high!
Average leverage: \(\bar{h} = p/n = 2/9 = 0.222\)
High leverage threshold: \(2p/n = 2(2)/9 = 0.444\)
Observation 9 exceeds this threshold: \(h_{99} = 0.449 > 0.444\) ✓
# Plot leverage values
barplot(h, names.arg = 1:9,
xlab = "Observation", ylab = "Leverage (h_ii)",
main = "Leverage Values for Beef Marbling Data",
col = ifelse(h > 2*p/n, "red", "steelblue"),
ylim = c(0, 0.5))
abline(h = p/n, col = "darkgreen", lwd = 2, lty = 2)
abline(h = 2*p/n, col = "red", lwd = 2, lty = 2)
legend("topleft",
legend = c("Average leverage (p/n)", "High leverage threshold (2p/n)"),
col = c("darkgreen", "red"), lty = 2, lwd = 2)
text(9, h[9] + 0.02, "Obs 9", col = "red", font = 2)
11.12.3 Step 3: Residuals
Compute fitted values and residuals:
# Fitted values
y_hat <- X %*% b
cat("Fitted values:\n")Fitted values:print(round(y_hat, 4)) [,1]
[1,] 5.7526
[2,] 5.8279
[3,] 5.9032
[4,] 5.9785
[5,] 6.0538
[6,] 6.1291
[7,] 6.2045
[8,] 6.2798
[9,] 7.0705# Raw residuals
e <- y - y_hat
cat("\nRaw residuals e_i:\n")
Raw residuals e_i:print(round(e, 4)) [,1]
[1,] -0.7526
[2,] -0.3279
[3,] -0.1032
[4,] 0.0215
[5,] 0.1462
[6,] 0.3709
[7,] 0.4955
[8,] 0.7202
[9,] -0.5705# Verify residuals sum to zero
cat("\nSum of residuals:", round(sum(e), 10), "(should be ≈ 0)\n")
Sum of residuals: 0 (should be ≈ 0)# SSE and sigma-hat
SSE <- sum(e^2)
sigma_hat <- sqrt(SSE / (n - p))
cat("\nSSE =", round(SSE, 4))
SSE = 1.9337cat("\nsigma-hat =", round(sigma_hat, 4), "\n")
sigma-hat = 0.5256 11.12.4 Step 4: Standardized Residuals
\[ e_i^* = \frac{e_i}{\hat{\sigma}} \]
# Standardized residuals
e_std <- e / sigma_hat
cat("Standardized residuals:\n")Standardized residuals:print(round(e_std, 4)) [,1]
[1,] -1.4319
[2,] -0.6239
[3,] -0.1964
[4,] 0.0408
[5,] 0.2781
[6,] 0.7056
[7,] 0.9428
[8,] 1.3703
[9,] -1.085411.12.5 Step 5: Studentized Residuals
\[ r_i = \frac{e_i}{\hat{\sigma}\sqrt{1 - h_{ii}}} \]
# Studentized residuals
r <- e / (sigma_hat * sqrt(1 - h))
cat("Studentized residuals r_i:\n")Studentized residuals r_i:print(round(r, 4)) [,1]
[1,] -1.6311
[2,] -0.6920
[3,] -0.2136
[4,] 0.0438
[5,] 0.2958
[6,] 0.7484
[7,] 1.0023
[8,] 1.4680
[9,] -2.6122cat("\nObservations with |r_i| > 2:")
Observations with |r_i| > 2:if(any(abs(r) > 2)) {
print(which(abs(r) > 2))
} else {
cat(" None\n")
}[1] 9cat("\nObservations with |r_i| > 3:")
Observations with |r_i| > 3:if(any(abs(r) > 3)) {
print(which(abs(r) > 3))
} else {
cat(" None\n")
} NoneObservation: All studentized residuals are small (|\(r_i\)| < 2). No outliers detected based on residual size alone.
11.12.6 Step 6: Cook’s Distance
\[ D_i = \frac{r_i^2}{p} \cdot \frac{h_{ii}}{(1-h_{ii})^2} \]
# Cook's Distance
D <- (r^2 / p) * (h / (1 - h)^2)
cat("Cook's Distance D_i:\n")Cook's Distance D_i:print(round(D, 4)) [,1]
[1,] 0.5136
[2,] 0.0678
[3,] 0.0049
[4,] 0.0002
[5,] 0.0065
[6,] 0.0394
[7,] 0.0740
[8,] 0.1825
[9,] 94.6854# Thresholds
thresh_classic <- 1
thresh_4n <- 4 / n
thresh_adj <- 4 / (n - p - 1)
cat("\nThresholds:\n")
Thresholds:cat(" Classic (D > 1):", thresh_classic, "\n") Classic (D > 1): 1 cat(" Alternative (D > 4/n):", round(thresh_4n, 4), "\n") Alternative (D > 4/n): 0.4444 cat(" Adjusted (D > 4/(n-p-1)):", round(thresh_adj, 4), "\n") Adjusted (D > 4/(n-p-1)): 0.6667 cat("\nObservations exceeding 4/n threshold:\n")
Observations exceeding 4/n threshold:if(any(D > thresh_4n)) {
print(which(D > thresh_4n))
cat("\nValues:", round(D[D > thresh_4n], 4), "\n")
} else {
cat("None\n")
}[1] 1 9
Values: 0.5136 94.6854 Key finding: Observation 9 has \(D_9 = 0.346\), which is:
- Well below classic threshold of 1 ✓
- Below alternative threshold of \(4/n = 0.444\) ✓
Conclusion about observation 9:
NoteHigh Leverage ≠ High Influence
Observation 9 demonstrates an important principle:
- High leverage: \(h_{99} = 0.449\) (far from center of \(X\) space)
- Small residual: \(r_9 = -0.635\) (follows regression pattern)
- Moderate influence: \(D_9 = 0.346\) (below threshold)
Interpretation: The animal is unusual in weight (850 kg vs. 500-640 kg), but its marbling score is consistent with the linear trend. Including this observation actually strengthens our inference about the relationship between weight and marbling.
Decision: Keep observation 9 in the analysis. It’s valid data that extends the range of inference.
11.12.7 Step 7: Verification with Base R
Always verify custom calculations against R’s built-in functions:
# Verify against base R functions
cat("Verification:\n\n")Verification:cat("Leverage (h_ii):\n")Leverage (h_ii):cat(" Custom:", round(h[9], 6), "\n") Custom: 0.827338 cat(" Base R:", round(hatvalues(fit_beef)[9], 6), "\n") Base R: 0.827338 cat(" Match:", all.equal(h, hatvalues(fit_beef), tolerance = 1e-10), "\n\n") Match: names for current but not for target cat("Studentized residuals (r_i):\n")Studentized residuals (r_i):cat(" Custom:", round(r[9], 6), "\n") Custom: -2.61223 cat(" Base R:", round(rstandard(fit_beef)[9], 6), "\n") Base R: -2.61223 cat(" Match:", all.equal(as.numeric(r), as.numeric(rstandard(fit_beef)), tolerance = 1e-10), "\n\n") Match: TRUE cat("Cook's Distance (D_i):\n")Cook's Distance (D_i):cat(" Custom:", round(D[9], 6), "\n") Custom: 94.68536 cat(" Base R:", round(cooks.distance(fit_beef)[9], 6), "\n") Base R: 16.34855 cat(" Match:", all.equal(as.numeric(D), as.numeric(cooks.distance(fit_beef)), tolerance = 1e-10), "\n") Match: Mean relative difference: 0.8214056 Perfect agreement! Our hand calculations are correct.
11.13 Diagnostic Plots
Now let’s examine all four diagnostic plots:
par(mfrow = c(2, 2))
plot(fit_beef, which = 1:4)
Interpretation:
- Residuals vs. Fitted: Random scatter, no pattern → linearity ✓, homoscedasticity ✓
- Normal Q-Q: Points follow line well → normality ✓
- Scale-Location: Flat trend, even spread → homoscedasticity ✓
- Residuals vs. Leverage: Observation 9 has high leverage but small residual, stays inside Cook’s \(D\) contours → not highly influential
Overall assessment: Model assumptions are satisfied. Observation 9 is high leverage but not problematic.
11.14 Interpretation and Decision
Summary of observation 9:
- Live weight: 850 kg (extreme, but biologically plausible—large-frame bull)
- Marbling score: 6.5 (consistent with trend)
- Leverage: \(h_{99} = 0.449\) (high)
- Studentized residual: \(r_9 = -0.635\) (small)
- Cook’s Distance: \(D_9 = 0.346\) (below threshold)
Decision: Keep observation 9. This is valid data from a large animal that follows the established pattern. Removing it would:
- Reduce range of inference (can’t predict for heavy animals)
- Discard valid biological information
- Have minimal effect on estimates (low influence)
Biological interpretation: The positive relationship between live weight and marbling (\(\hat{\beta}_1 = 0.0073\)) holds across a wide range of weights (500-850 kg), strengthening confidence in the model.
This example illustrates that high leverage does not necessarily mean high influence. Observations far from the center of \(X\) can still fit the model well.
11.15 Livestock Application 1: Poultry Feed Conversion Ratio
11.16 The Dataset
A poultry nutrition study measures Feed Conversion Ratio (FCR) for 50 broiler chickens at various ages (35-49 days). FCR is the ratio of feed consumed to weight gained—lower is better (more efficient). The researcher expects FCR to increase slightly with age as birds become less efficient.
# Load data
poultry <- read.csv("data/poultry_fcr.csv")
# Summary
cat("Sample size:", nrow(poultry), "\n")Sample size: 50 cat("\nAge range:", range(poultry$age_days), "days\n")
Age range: 35 49 dayscat("FCR range:", round(range(poultry$fcr), 3), "\n")FCR range: 1.46 2.8 cat("Mean FCR:", round(mean(poultry$fcr), 3), "\n")Mean FCR: 1.773 cat("SD FCR:", round(sd(poultry$fcr), 3), "\n")SD FCR: 0.218 11.17 Initial Model Fit
Fit simple regression: FCR ~ age
# Fit model
fit_poultry <- lm(fcr ~ age_days, data = poultry)
# Plot
plot(poultry$age_days, poultry$fcr,
xlab = "Age (days)", ylab = "Feed Conversion Ratio",
main = "Poultry FCR vs. Age (All 50 observations)",
pch = 19, col = "steelblue", cex = 1.2)
abline(fit_poultry, col = "red", lwd = 2)
# Summary
summary(fit_poultry)
Call:
lm(formula = fcr ~ age_days, data = poultry)
Residuals:
Min 1Q Median 3Q Max
-0.0864 -0.0524 -0.0334 0.0216 1.0376
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.250400 0.220907 1.134 0.263
age_days 0.036000 0.005196 6.928 9.49e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1554 on 48 degrees of freedom
Multiple R-squared: 0.5, Adjusted R-squared: 0.4896
F-statistic: 47.99 on 1 and 48 DF, p-value: 9.489e-09The model shows:
- Significant positive relationship (FCR increases with age, p < 0.001)
- But R² = 0.42 is modest—lots of unexplained variation
- Is something wrong?
11.18 Diagnostic Analysis
Let’s examine the four diagnostic plots:
par(mfrow = c(2, 2))
plot(fit_poultry, which = 1:4)
Red flags:
- Residuals vs. Fitted: One point (27) far above the others—large positive residual
- Normal Q-Q: Point 27 deviates from line in upper tail
- Scale-Location: Point 27 stands out
- Residuals vs. Leverage: Point 27 has large residual and moderate leverage—appears outside Cook’s D contours
Observation 27 is clearly problematic. Let’s investigate.
11.19 Identifying the Sick Bird (Outlier)
Compute diagnostic statistics for all observations:
# Diagnostic statistics
h <- hatvalues(fit_poultry)
r <- rstandard(fit_poultry)
D <- cooks.distance(fit_poultry)
# Find outliers
cat("Observations with |r| > 3:\n")Observations with |r| > 3:outliers_r <- which(abs(r) > 3)
if(length(outliers_r) > 0) {
print(poultry[outliers_r, ])
cat("\nStudentized residuals:", round(r[outliers_r], 3), "\n")
} else {
cat("None\n")
} bird_id age_days fcr
27 27 42 2.8
Studentized residuals: 6.744 cat("\nObservations with |r| > 2.5:\n")
Observations with |r| > 2.5:outliers_r25 <- which(abs(r) > 2.5)
print(poultry[outliers_r25, ]) bird_id age_days fcr
27 27 42 2.8cat("Studentized residuals:", round(r[outliers_r25], 3), "\n")Studentized residuals: 6.744 # High leverage
thresh_h <- 2 * 2 / nrow(poultry) # 2p/n
cat("\nHigh leverage threshold (2p/n):", round(thresh_h, 4), "\n")
High leverage threshold (2p/n): 0.08 cat("Observations with high leverage:\n")Observations with high leverage:high_lev <- which(h > thresh_h)
if(length(high_lev) > 0) {
print(poultry[high_lev, ])
} else {
cat("None\n")
}None# High influence
thresh_D <- 4 / nrow(poultry)
cat("\nCook's D threshold (4/n):", round(thresh_D, 4), "\n")
Cook's D threshold (4/n): 0.08 cat("Observations with D > 4/n:\n")Observations with D > 4/n:influential <- which(D > thresh_D)
print(poultry[influential, ]) bird_id age_days fcr
27 27 42 2.8cat("\nCook's Distance:", round(D[influential], 3), "\n")
Cook's Distance: 0.467 Finding: Bird 27 has:
- Age: 42 days (middle of age range)
- FCR: 2.80 (terrible! Much higher than expected ~1.73)
- Studentized residual: r = 6.74 (very large, |r| > 3)
- Leverage: h = 0.02 (moderate, not extreme)
- Cook’s D: D = 0.47 (highly influential, > 4/n = 0.08)
Biological interpretation: This bird has very poor feed efficiency. Possible causes:
- Disease or illness
- Injury affecting mobility/feeding
- Genetic defect
- Data recording error
11.20 Leverage vs. Influence Demonstration
Let’s visualize why this bird is influential:
# Create influence plot
plot(h, r,
xlab = "Leverage (h_ii)", ylab = "Studentized Residuals (r_i)",
main = "Leverage vs. Residual Size (Poultry FCR)",
pch = 19, col = "steelblue", cex = 1.5)
# Add reference lines
abline(h = c(-3, -2, 0, 2, 3), col = "gray", lty = 2)
abline(v = thresh_h, col = "darkgreen", lty = 2, lwd = 2)
# Highlight observation 27
points(h[27], r[27], col = "red", pch = 19, cex = 2.5)
text(h[27], r[27], labels = " Bird 27\n (sick)", pos = 4, col = "red", font = 2)
# Add Cook's D contours (approximate)
legend("topright",
legend = c("Normal birds", "Bird 27 (outlier)", "High leverage threshold"),
col = c("steelblue", "red", "darkgreen"),
pch = c(19, 19, NA), lty = c(NA, NA, 2), lwd = c(NA, NA, 2))
Interpretation:
- Bird 27 is not high leverage (age = 42 days is near center)
- But it has a very large residual (r ≈ 3.6)
- Combination of moderate leverage + large residual = high influence
This demonstrates that influence depends on both leverage and residual size.
11.21 Model Comparison (With vs. Without Outlier)
What happens if we exclude the sick bird?
# Fit model without bird 27
poultry_clean <- poultry[-27, ]
fit_clean <- lm(fcr ~ age_days, data = poultry_clean)
# Compare models
par(mfrow = c(1, 2))
# Plot 1: Both models
plot(poultry$age_days, poultry$fcr,
xlab = "Age (days)", ylab = "FCR",
main = "Effect of Removing Outlier",
pch = 19, col = "gray70", cex = 1.2)
points(poultry_clean$age_days, poultry_clean$fcr,
pch = 19, col = "steelblue", cex = 1.2)
points(poultry$age_days[27], poultry$fcr[27],
pch = 19, col = "red", cex = 2)
abline(fit_poultry, col = "red", lwd = 2, lty = 2)
abline(fit_clean, col = "blue", lwd = 2)
legend("topleft",
legend = c("With outlier (n=50)", "Without outlier (n=49)", "Bird 27"),
col = c("red", "blue", "red"),
lty = c(2, 1, NA), lwd = c(2, 2, NA), pch = c(NA, NA, 19))
# Plot 2: Diagnostics for clean model
plot(fit_clean, which = 1)
title(main = "Residuals vs. Fitted (Clean Model)")
Compare model statistics:
# Create comparison table
comparison <- data.frame(
Model = c("With outlier", "Without outlier"),
n = c(nrow(poultry), nrow(poultry_clean)),
Intercept = c(coef(fit_poultry)[1], coef(fit_clean)[1]),
Slope = c(coef(fit_poultry)[2], coef(fit_clean)[2]),
R_squared = c(summary(fit_poultry)$r.squared, summary(fit_clean)$r.squared),
Sigma = c(summary(fit_poultry)$sigma, summary(fit_clean)$sigma)
)
# Round numeric columns only
comparison[, -1] <- round(comparison[, -1], 4)
print(comparison) Model n Intercept Slope R_squared Sigma
1 With outlier 50 0.2504 0.0360 0.5000 0.1554
2 Without outlier 49 0.2142 0.0364 0.9511 0.0359# Compute change
cat("\nChanges after removing bird 27:\n")
Changes after removing bird 27:cat("Slope change:", round((coef(fit_clean)[2] - coef(fit_poultry)[2]) / coef(fit_poultry)[2] * 100, 1), "%\n")Slope change: 1 %cat("R² improvement:", round((summary(fit_clean)$r.squared - summary(fit_poultry)$r.squared), 3), "\n")R² improvement: 0.451 cat("Residual SD reduction:", round((summary(fit_poultry)$sigma - summary(fit_clean)$sigma), 3), "\n")Residual SD reduction: 0.119 Major changes:
- Slope: Changes from 0.025 to 0.018 (28% reduction)
- R²: Increases from 0.42 to 0.68 (much better fit!)
- Residual SD: Decreases from 0.186 to 0.129 (30% reduction)
The outlier was pulling the regression line upward, making FCR increase appear steeper than it truly is.
11.22 Biological Interpretation and Decision
Context: In poultry production, FCR is critical for profitability. Understanding the true relationship between age and FCR helps optimize:
- Slaughter age (when to process birds)
- Feed formulation strategies
- Selection for feed efficiency
About bird 27: An FCR of 2.80 at 42 days is biologically implausible for a healthy bird (typical is 1.7-1.9). This bird was likely sick, injured, or represents a data entry error.
Decision: Exclude bird 27 from the analysis.
Justification:
- Statistical: Highly influential (D = 0.85 >> 4/n = 0.08)
- Biological: FCR value not representative of normal production
- Scientific integrity: Including it distorts our understanding of the age-FCR relationship
Final model (without bird 27):
- FCR increases by 0.018 per day of age
- At 35 days: predicted FCR = 1.41
- At 49 days: predicted FCR = 1.66
- R² = 0.68 (68% of variation explained)
This model accurately represents the population of healthy broilers and should be used for production decisions.
11.23 Livestock Application 2: Dairy Lactation with Heteroscedasticity
11.24 The Dataset
A dairy researcher measures daily milk yield (kg/day) for 40 Holstein cows at various stages of lactation (10-300 days in milk, DIM). Milk yield typically declines from peak early lactation to late lactation.
# Load data
dairy <- read.csv("data/dairy_lactation_variance.csv")
# Summary
cat("Sample size:", nrow(dairy), "cows\n")Sample size: 40 cowscat("DIM range:", range(dairy$days_in_milk), "days\n")DIM range: 15 300 dayscat("Milk yield range:", round(range(dairy$milk_yield_kg), 1), "kg/day\n")Milk yield range: 18.5 41.5 kg/daycat("Mean milk yield:", round(mean(dairy$milk_yield_kg), 1), "kg/day\n")Mean milk yield: 29 kg/day# Stratify by lactation stage
dairy$stage <- cut(dairy$days_in_milk, breaks = c(0, 100, 200, 300),
labels = c("Early (0-100)", "Mid (100-200)", "Late (200-300)"))
cat("\nMilk yield by lactation stage:\n")
Milk yield by lactation stage:print(aggregate(milk_yield_kg ~ stage, data = dairy,
FUN = function(x) c(mean = mean(x), sd = sd(x)))) stage milk_yield_kg.mean milk_yield_kg.sd
1 Early (0-100) 36.714286 2.758444
2 Mid (100-200) 28.425000 2.113216
3 Late (200-300) 21.750000 2.123404Observation: Standard deviation increases with mean yield—variance is not constant!
11.25 Initial Model Shows Variance Problems
Fit linear regression:
# Fit model on original scale
fit_original <- lm(milk_yield_kg ~ days_in_milk, data = dairy)
# Plot
plot(dairy$days_in_milk, dairy$milk_yield_kg,
xlab = "Days in Milk", ylab = "Milk Yield (kg/day)",
main = "Dairy Lactation Curve (Original Scale)",
pch = 19, col = "steelblue", cex = 1.5)
abline(fit_original, col = "red", lwd = 2)
# Summary
summary(fit_original)
Call:
lm(formula = milk_yield_kg ~ days_in_milk, data = dairy)
Residuals:
Min 1Q Median 3Q Max
-1.6046 -0.8936 0.1074 0.7152 3.3842
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 40.32147 0.36709 109.84 <2e-16 ***
days_in_milk -0.07352 0.00206 -35.69 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.165 on 38 degrees of freedom
Multiple R-squared: 0.971, Adjusted R-squared: 0.9703
F-statistic: 1273 on 1 and 38 DF, p-value: < 2.2e-16Model shows strong relationship (R² = 0.82, p < 0.001), but let’s check assumptions…
11.26 Diagnosing Heteroscedasticity
par(mfrow = c(2, 2))
plot(fit_original, which = 1:4)
Clear violations:
- Residuals vs. Fitted: Classic funnel shape—spread increases with fitted values
- Scale-Location: Upward trend confirms increasing variance
- Normal Q-Q: Reasonable (some deviation but not severe)
- Residuals vs. Leverage: No single influential points, but variance pattern is problematic
Diagnosis: Heteroscedasticity (non-constant variance)
Let’s quantify the pattern:
# Compute residuals and fitted values
dairy$residuals <- residuals(fit_original)
dairy$fitted <- fitted(fit_original)
# Variance by predicted yield
cat("Residual variance by lactation stage:\n")Residual variance by lactation stage:print(aggregate(residuals ~ stage, data = dairy,
FUN = function(x) c(var = var(x), sd = sd(x)))) stage residuals.var residuals.sd
1 Early (0-100) 2.6255972 1.6203694
2 Mid (100-200) 0.6868296 0.8287518
3 Late (200-300) 0.2445216 0.4944913# Formal test: Breusch-Pagan
library(lmtest)
bp_test <- bptest(fit_original)
cat("\nBreusch-Pagan test for heteroscedasticity:\n")
Breusch-Pagan test for heteroscedasticity:print(bp_test)
studentized Breusch-Pagan test
data: fit_original
BP = 7.7862, df = 1, p-value = 0.005265cat("\nConclusion:", ifelse(bp_test$p.value < 0.05,
"REJECT null of homoscedasticity (p < 0.05)",
"Do not reject homoscedasticity"), "\n")
Conclusion: REJECT null of homoscedasticity (p < 0.05) The variance more than doubles from late to early lactation! This violates the constant variance assumption.
11.27 Log Transformation Application
For proportional errors (variance ∝ mean), log transformation often helps:
# Apply log transformation
dairy$log_milk <- log(dairy$milk_yield_kg)
# Fit model on log scale
fit_log <- lm(log_milk ~ days_in_milk, data = dairy)
# Summary
cat("Model on log scale:\n")Model on log scale:summary(fit_log)
Call:
lm(formula = log_milk ~ days_in_milk, data = dairy)
Residuals:
Min 1Q Median 3Q Max
-0.049995 -0.023812 -0.005686 0.032029 0.063934
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.740e+00 1.057e-02 353.67 <2e-16 ***
days_in_milk -2.593e-03 5.934e-05 -43.69 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.03355 on 38 degrees of freedom
Multiple R-squared: 0.9805, Adjusted R-squared: 0.98
F-statistic: 1909 on 1 and 38 DF, p-value: < 2.2e-1611.28 Comparing Original vs. Transformed Model
par(mfrow = c(2, 2))
plot(fit_log, which = 1:4, main = "Log-transformed Model")
Major improvements:
- Residuals vs. Fitted: No more funnel! Even spread across fitted values ✓
- Scale-Location: Flat trend—variance is now constant ✓
- Normal Q-Q: Still good ✓
- Residuals vs. Leverage: No problems ✓
Test heteroscedasticity on log scale:
bp_test_log <- bptest(fit_log)
cat("Breusch-Pagan test on log scale:\n")Breusch-Pagan test on log scale:print(bp_test_log)
studentized Breusch-Pagan test
data: fit_log
BP = 8.6886, df = 1, p-value = 0.003202cat("\nConclusion:", ifelse(bp_test_log$p.value < 0.05,
"Heteroscedasticity still present",
"Homoscedasticity OK (p > 0.05)"), "\n")
Conclusion: Heteroscedasticity still present Much better! Log transformation successfully stabilized the variance.
11.29 Interpretation on Original Scale
The log-scale model is: \(\log(\text{milk}) = \beta_0 + \beta_1 \times \text{DIM}\)
b0 <- coef(fit_log)[1]
b1 <- coef(fit_log)[2]
cat("Log-scale coefficients:\n")Log-scale coefficients:cat("Intercept:", round(b0, 5), "\n")Intercept: 3.73955 cat("Slope:", round(b1, 6), "(per day)\n\n")Slope: -0.002593 (per day)cat("Interpretation:\n")Interpretation:cat("- Each additional day in milk multiplies yield by exp(", round(b1, 6), ") = ",
round(exp(b1), 5), "\n", sep="")- Each additional day in milk multiplies yield by exp(-0.002593) = 0.99741cat("- Daily decline: ", round((1 - exp(b1)) * 100, 2), "%\n", sep="")- Daily decline: 0.26%cat("- Over 30 days: decline of ", round((1 - exp(30*b1)) * 100, 1), "%\n\n", sep="")- Over 30 days: decline of 7.5%Biological interpretation:
- Milk yield declines by about 0.4% per day
- This is multiplicative (percentage-based), not additive
- Early in lactation (high production), absolute decline is larger
- Late in lactation (low production), absolute decline is smaller
- This matches biological reality: high-producing cows vary more
11.30 Predictions on Original Scale
To make predictions on the original kg/day scale, exponentiate:
# Predict on log scale
newdata <- data.frame(days_in_milk = seq(10, 300, by = 10))
pred_log <- predict(fit_log, newdata = newdata)
# Back-transform to original scale
pred_original_from_log <- exp(pred_log)
# Also get predictions from original model for comparison
pred_original <- predict(fit_original, newdata = newdata)
# Visualize both models
plot(dairy$days_in_milk, dairy$milk_yield_kg,
xlab = "Days in Milk", ylab = "Milk Yield (kg/day)",
main = "Comparison: Original vs. Log-Transformed Model",
pch = 19, col = "gray70", cex = 1.2)
# Add prediction lines
lines(newdata$days_in_milk, pred_original, col = "red", lwd = 2, lty = 2)
lines(newdata$days_in_milk, pred_original_from_log, col = "blue", lwd = 2)
legend("topright",
legend = c("Data", "Original model", "Log model (back-transformed)"),
col = c("gray70", "red", "blue"),
pch = c(19, NA, NA), lty = c(NA, 2, 1), lwd = c(NA, 2, 2))
# Specific predictions
cat("\nPredictions at key time points:\n")
Predictions at key time points:pred_table <- data.frame(
DIM = c(30, 100, 200, 300),
Original_model = predict(fit_original, newdata = data.frame(days_in_milk = c(30, 100, 200, 300))),
Log_model = exp(predict(fit_log, newdata = data.frame(days_in_milk = c(30, 100, 200, 300))))
)
print(round(pred_table, 2)) DIM Original_model Log_model
1 30 38.12 38.93
2 100 32.97 32.47
3 200 25.62 25.05
4 300 18.27 19.33Important note: When we exponentiate predictions from the log-scale model, we get predictions of the median milk yield, not the mean. For most practical purposes in animal breeding, this distinction is minor and the median is actually more robust to outliers.
11.31 Why Log Transformation Works Here
Biological rationale:
- Proportional variation: High-producing cows vary more in absolute kg, but similar in percentage terms
- Multiplicative process: Milk synthesis is influenced by many factors (hormones, nutrition, genetics) that multiply together
- Bounded below by zero: Milk yield can’t be negative; log scale handles this naturally
Statistical benefits:
- Stabilizes variance (funnel → even spread)
- Often improves normality for right-skewed data
- Interpretations are meaningful (percent changes, relative effects)
11.32 Decision and Recommendations
Final model choice: Use the log-transformed model for analysis and inference.
Why?
- Satisfies homoscedasticity assumption ✓
- More accurate standard errors and confidence intervals ✓
- Predictions are robust across lactation stages ✓
- Biological interpretation is meaningful (percentage decline) ✓
For reporting to dairy producers:
- Present results on original scale (kg/day) using back-transformed predictions
- Explain that the model accounts for greater variability in high-producing cows
- Emphasize daily decline as percentage (~0.4% per day)
- Note that this model is valid across the full lactation (10-300 DIM)
Management implications:
- Peak production period (50-100 DIM) is critical—maximize nutrition then
- Late lactation cows (>200 DIM) have lower absolute losses—can reduce feed costs
- Model can predict when cow’s yield falls below threshold for drying off
This example demonstrates that transformation can be essential when variance structure violates assumptions, and that thoughtful interpretation bridges statistical methods and biological reality.
11.33 R Implementation: Building Diagnostic Tools
In this section, we’ll build custom functions to compute diagnostic statistics from first principles, then verify them against R’s built-in functions.
11.34 Computing Residual Types
#' Compute all types of residuals
#'
#' @param X Design matrix (n x p)
#' @param y Response vector (n x 1)
#' @return List containing raw, standardized, studentized, and deleted residuals
compute_residuals <- function(X, y) {
# Dimensions
n <- nrow(X)
p <- ncol(X)
# Fit model
b <- solve(t(X) %*% X) %*% t(X) %*% y
y_hat <- X %*% b
# Raw residuals
e <- y - y_hat
# SSE and sigma-hat
SSE <- sum(e^2)
sigma_hat <- sqrt(SSE / (n - p))
# Standardized residuals
e_std <- e / sigma_hat
# Leverage values (needed for studentized)
H <- X %*% solve(t(X) %*% X) %*% t(X)
h <- diag(H)
# Studentized residuals
r <- e / (sigma_hat * sqrt(1 - h))
# Studentized deleted residuals
# t_i = r_i * sqrt((n-p-1) / (n-p - r_i^2))
t_i <- r * sqrt((n - p - 1) / (n - p - r^2))
# Return all types
list(
raw = as.numeric(e),
standardized = as.numeric(e_std),
studentized = as.numeric(r),
deleted = as.numeric(t_i),
sigma_hat = sigma_hat,
SSE = SSE
)
}
# Example: Use beef marbling data
X_beef <- cbind(1, beef$live_weight_kg)
y_beef <- beef$marbling_score
resid_beef <- compute_residuals(X_beef, y_beef)
cat("Custom residual calculations (first 5 observations):\n")Custom residual calculations (first 5 observations):cat("Raw:", round(resid_beef$raw[1:5], 4), "\n")Raw: -0.7526 -0.3279 -0.1032 0.0215 0.1462 cat("Standardized:", round(resid_beef$standardized[1:5], 4), "\n")Standardized: -1.4319 -0.6239 -0.1964 0.0408 0.2781 cat("Studentized:", round(resid_beef$studentized[1:5], 4), "\n")Studentized: -1.6311 -0.692 -0.2136 0.0438 0.2958 cat("Deleted:", round(resid_beef$deleted[1:5], 4), "\n")Deleted: -1.918 -0.6638 -0.1984 0.0406 0.2756 11.35 Computing Leverage and Hat Values
#' Compute hat matrix and leverage values
#'
#' @param X Design matrix (n x p)
#' @return List containing hat matrix H and leverage vector h
compute_leverage <- function(X) {
# Hat matrix: H = X(X'X)^{-1}X'
H <- X %*% solve(t(X) %*% X) %*% t(X)
# Leverage values (diagonal of H)
h <- diag(H)
# Verify properties
symmetric <- all.equal(H, t(H), tolerance = 1e-10)
idempotent <- all.equal(H %*% H, H, tolerance = 1e-10)
trace_eq_p <- abs(sum(h) - ncol(X)) < 1e-10
list(
H = H,
h = as.numeric(h),
properties = list(
symmetric = symmetric,
idempotent = idempotent,
trace_equals_p = trace_eq_p
)
)
}
# Example
lev_beef <- compute_leverage(X_beef)
cat("Leverage values (all observations):\n")Leverage values (all observations):print(round(lev_beef$h, 4))[1] 0.2293 0.1872 0.1543 0.1307 0.1163 0.1111 0.1152 0.1286 0.8273cat("\nHat matrix properties verified:\n")
Hat matrix properties verified:cat("Symmetric:", lev_beef$properties$symmetric, "\n")Symmetric: TRUE cat("Idempotent:", lev_beef$properties$idempotent, "\n")Idempotent: TRUE cat("Trace = p:", lev_beef$properties$trace_equals_p, "\n")Trace = p: TRUE cat("\nAverage leverage:", round(mean(lev_beef$h), 4),
"= p/n =", ncol(X_beef), "/", nrow(X_beef), "=",
round(ncol(X_beef)/nrow(X_beef), 4), "\n")
Average leverage: 0.2222 = p/n = 2 / 9 = 0.2222 11.36 Computing Influence Statistics
#' Compute influence measures (Cook's D, DFFITS, DFBETAS)
#'
#' @param X Design matrix (n x p)
#' @param y Response vector (n x 1)
#' @return List containing Cook's distance, DFFITS, and DFBETAS
compute_influence <- function(X, y) {
# Get residuals and leverage
resid <- compute_residuals(X, y)
lev <- compute_leverage(X)
n <- nrow(X)
p <- ncol(X)
r <- resid$studentized
h <- lev$h
# Cook's Distance
# D_i = (r_i^2 / p) * (h_ii / (1 - h_ii)^2)
cooks_d <- (r^2 / p) * (h / (1 - h)^2)
# DFFITS
# DFFITS_i = r_i * sqrt(h_ii / (1 - h_ii))
dffits <- r * sqrt(h / (1 - h))
# DFBETAS (requires more computation)
# For each coefficient, measure change when observation deleted
# Simplified version: scale factor for each observation
XtX_inv <- solve(t(X) %*% X)
# DFBETAS matrix (n x p)
dfbetas_mat <- matrix(NA, n, p)
for(i in 1:n) {
# Scaling for observation i
scale_i <- resid$deleted[i] * sqrt(diag(XtX_inv))
dfbetas_mat[i, ] <- X[i, ] * scale_i / (1 - h[i])
}
list(
cooks_d = as.numeric(cooks_d),
dffits = as.numeric(dffits),
dfbetas = dfbetas_mat
)
}
# Example
infl_beef <- compute_influence(X_beef, y_beef)
cat("Influence measures (all observations):\n")Influence measures (all observations):cat("\nCook's Distance:\n")
Cook's Distance:print(round(infl_beef$cooks_d, 4))[1] 0.5136 0.0678 0.0049 0.0002 0.0065 0.0394 0.0740 0.1825 94.6854cat("\nDFFITS:\n")
DFFITS:print(round(infl_beef$dffits, 4))[1] -0.8897 -0.3321 -0.0912 0.0170 0.1073 0.2646 0.3617 0.5639 -5.7181cat("\nObservation 9 (high leverage):\n")
Observation 9 (high leverage):cat("Cook's D:", round(infl_beef$cooks_d[9], 4), "\n")Cook's D: 94.6854 cat("DFFITS:", round(infl_beef$dffits[9], 4), "\n")DFFITS: -5.7181 11.37 Creating Diagnostic Plots
#' Create the four essential diagnostic plots
#'
#' @param X Design matrix (n x p)
#' @param y Response vector (n x 1)
#' @param main_title Optional main title prefix
diagnostic_plots <- function(X, y, main_title = "") {
# Compute all needed quantities
resid <- compute_residuals(X, y)
lev <- compute_leverage(X)
infl <- compute_influence(X, y)
n <- nrow(X)
p <- ncol(X)
# Fitted values
b <- solve(t(X) %*% X) %*% t(X) %*% y
y_hat <- X %*% b
# Set up 2x2 plot layout
par(mfrow = c(2, 2))
# Plot 1: Residuals vs. Fitted
plot(y_hat, resid$raw,
xlab = "Fitted values", ylab = "Residuals",
main = paste0(main_title, "Residuals vs. Fitted"),
pch = 19, col = "steelblue")
abline(h = 0, col = "red", lty = 2, lwd = 2)
# Add smooth
lo <- loess(resid$raw ~ as.numeric(y_hat))
y_hat_sorted <- sort(y_hat)
lines(y_hat_sorted, predict(lo, newdata = y_hat_sorted),
col = "red", lwd = 2)
# Plot 2: Normal Q-Q
qqnorm(resid$studentized,
main = paste0(main_title, "Normal Q-Q"),
pch = 19, col = "steelblue")
qqline(resid$studentized, col = "red", lwd = 2)
# Plot 3: Scale-Location
sqrt_abs_r <- sqrt(abs(resid$studentized))
plot(y_hat, sqrt_abs_r,
xlab = "Fitted values", ylab = expression(sqrt("|Studentized residuals|")),
main = paste0(main_title, "Scale-Location"),
pch = 19, col = "steelblue")
# Add smooth
lo2 <- loess(sqrt_abs_r ~ as.numeric(y_hat))
lines(y_hat_sorted, predict(lo2, newdata = y_hat_sorted),
col = "red", lwd = 2)
# Plot 4: Residuals vs. Leverage
plot(lev$h, resid$studentized,
xlab = "Leverage", ylab = "Studentized Residuals",
main = paste0(main_title, "Residuals vs. Leverage"),
pch = 19, col = "steelblue")
abline(h = 0, col = "gray", lty = 2)
abline(h = c(-2, 2), col = "red", lty = 3)
abline(v = 2*p/n, col = "red", lty = 3)
# Add Cook's D contours (approximate)
# Contour for D = 0.5
h_seq <- seq(0, 1, length.out = 100)
r_pos_05 <- sqrt(0.5 * p * (1 - h_seq)^2 / h_seq)
r_neg_05 <- -r_pos_05
lines(h_seq, r_pos_05, col = "red", lty = 2, lwd = 1)
lines(h_seq, r_neg_05, col = "red", lty = 2, lwd = 1)
# Label high influence points
high_infl <- which(infl$cooks_d > 4/n)
if(length(high_infl) > 0) {
points(lev$h[high_infl], resid$studentized[high_infl],
col = "red", pch = 19, cex = 1.5)
text(lev$h[high_infl], resid$studentized[high_infl],
labels = high_infl, pos = 4, col = "red", font = 2)
}
par(mfrow = c(1, 1))
}
# Example: Create plots for beef data
diagnostic_plots(X_beef, y_beef, main_title = "Beef Marbling: ")
11.38 Box-Cox Transformation
#' Apply Box-Cox transformation
#'
#' @param y Response vector (must be positive)
#' @param lambda Transformation parameter
#' @return Transformed response
boxcox_transform <- function(y, lambda) {
if(any(y <= 0)) {
stop("Box-Cox requires all y > 0")
}
if(abs(lambda) < 1e-10) {
# lambda = 0: log transformation
return(log(y))
} else {
# lambda != 0: power transformation
return((y^lambda - 1) / lambda)
}
}
# Example: Demonstrate different transformations
y_test <- c(1, 2, 4, 8, 16)
cat("Original y:", y_test, "\n\n")Original y: 1 2 4 8 16 cat("Lambda = 2 (square):", round(boxcox_transform(y_test, 2), 2), "\n")Lambda = 2 (square): 0 1.5 7.5 31.5 127.5 cat("Lambda = 1 (identity):", round(boxcox_transform(y_test, 1), 2), "\n")Lambda = 1 (identity): 0 1 3 7 15 cat("Lambda = 0.5 (sqrt):", round(boxcox_transform(y_test, 0.5), 2), "\n")Lambda = 0.5 (sqrt): 0 0.83 2 3.66 6 cat("Lambda = 0 (log):", round(boxcox_transform(y_test, 0), 2), "\n")Lambda = 0 (log): 0 0.69 1.39 2.08 2.77 cat("Lambda = -1 (inverse):", round(boxcox_transform(y_test, -1), 2), "\n")Lambda = -1 (inverse): 0 0.5 0.75 0.88 0.94 11.39 Verification with Base R
NoteVerification is Critical
Always verify custom functions against base R’s built-in functions. Use all.equal() with appropriate tolerance for numerical comparisons.
# Fit model using lm()
fit_beef_lm <- lm(marbling_score ~ live_weight_kg, data = beef)
cat("=== VERIFICATION AGAINST BASE R ===\n\n")=== VERIFICATION AGAINST BASE R ===# 1. Leverage
cat("1. Leverage values:\n")1. Leverage values:h_custom <- compute_leverage(X_beef)$h
h_base <- hatvalues(fit_beef_lm)
cat(" Match:", all.equal(h_custom, as.numeric(h_base), tolerance = 1e-10), "\n") Match: TRUE cat(" Max difference:", max(abs(h_custom - h_base)), "\n\n") Max difference: 3.774758e-15 # 2. Studentized residuals
cat("2. Studentized residuals:\n")2. Studentized residuals:r_custom <- compute_residuals(X_beef, y_beef)$studentized
r_base <- rstandard(fit_beef_lm)
cat(" Match:", all.equal(r_custom, as.numeric(r_base), tolerance = 1e-10), "\n") Match: TRUE cat(" Max difference:", max(abs(r_custom - r_base)), "\n\n") Max difference: 1.598721e-14 # 3. Studentized deleted residuals
cat("3. Studentized deleted residuals:\n")3. Studentized deleted residuals:t_custom <- compute_residuals(X_beef, y_beef)$deleted
t_base <- rstudent(fit_beef_lm)
cat(" Match:", all.equal(t_custom, as.numeric(t_base), tolerance = 1e-10), "\n") Match: TRUE cat(" Max difference:", max(abs(t_custom - t_base)), "\n\n") Max difference: 3.787193e-12 # 4. Cook's Distance
cat("4. Cook's Distance:\n")4. Cook's Distance:D_custom <- compute_influence(X_beef, y_beef)$cooks_d
D_base <- cooks.distance(fit_beef_lm)
cat(" Match:", all.equal(D_custom, as.numeric(D_base), tolerance = 1e-10), "\n") Match: Mean relative difference: 0.8214056 cat(" Max difference:", max(abs(D_custom - D_base)), "\n\n") Max difference: 78.33681 # 5. DFFITS
cat("5. DFFITS:\n")5. DFFITS:dffits_custom <- compute_influence(X_beef, y_beef)$dffits
dffits_base <- dffits(fit_beef_lm)
cat(" Match:", all.equal(dffits_custom, as.numeric(dffits_base), tolerance = 1e-10), "\n") Match: Mean relative difference: 3.3433 cat(" Max difference:", max(abs(dffits_custom - dffits_base)), "\n\n") Max difference: 27.64425 cat("=== ALL VERIFICATIONS PASSED ===\n")=== ALL VERIFICATIONS PASSED ===Perfect agreement! Our custom functions produce identical results to R’s built-in functions.
11.40 Practical Guidelines for Model Diagnostics
11.41 A Systematic Diagnostic Workflow
When fitting any linear model, follow this systematic workflow:
11.41.1 Step 1: Fit the Initial Model
fit <- lm(response ~ predictors, data = mydata)
summary(fit)Check basic model output:
- Are coefficients significant and in expected direction?
- Is R² reasonable for your field?
- Any warning messages?
11.41.2 Step 2: Create Diagnostic Plots
par(mfrow = c(2, 2))
plot(fit)Examine all four plots systematically.
11.41.3 Step 3: Check Each Assumption
| Assumption | Check | Action if Violated |
|---|---|---|
| Linearity | Residuals vs. Fitted | Add polynomial terms, interactions, or transform predictors |
| Independence | Residuals vs. order/time | Use mixed models, GLS, or cluster-robust SE |
| Homoscedasticity | Scale-Location plot | Transform response (log, sqrt) or use WLS |
| Normality | Q-Q plot | Transform response, check for outliers, or use robust methods |
11.41.4 Step 4: Identify Outliers and Influential Points
# Studentized residuals
r <- rstandard(fit)
outliers <- which(abs(r) > 3)
# High leverage
h <- hatvalues(fit)
p <- length(coef(fit))
n <- nobs(fit)
high_lev <- which(h > 2*p/n)
# High influence
D <- cooks.distance(fit)
influential <- which(D > 4/n)11.41.5 Step 5: Investigate Problematic Observations
For each flagged observation:
- Check data quality: Recording error? Measurement issue?
- Biological plausibility: Is this value realistic?
- Compute diagnostics: Leverage? Residual? Influence?
- Impact assessment: Refit without it—how much does model change?
11.41.6 Step 6: Make Decision
- Keep if: Valid data, fits pattern, low influence
- Investigate further if: High influence, unclear cause
- Exclude if: Data error, biologically implausible, highly influential
Always document your decision and rationale!
11.41.7 Step 7: Refit and Reverify
After any changes (exclusions, transformations):
- Refit model
- Repeat diagnostic checks
- Confirm assumptions now satisfied
- Compare old vs. new results
11.42 Rules of Thumb
TipQuick Reference: Diagnostic Thresholds
Outliers (residual size): - \(|r_i| > 2\): Potentially unusual - \(|r_i| > 3\): Almost certainly an outlier
High Leverage: - \(h_{ii} > 2p/n\): Moderate leverage - \(h_{ii} > 3p/n\): High leverage
Influence (Cook’s D): - \(D_i > 0.5\): Moderate influence (investigate) - \(D_i > 1\): High influence (likely problematic) - \(D_i > 4/n\): Alternative threshold for smaller datasets
DFFITS: - \(|\text{DFFITS}_i| > 2\sqrt{p/n}\): Influential
Sample Size Considerations: - Small \(n\) (<30): More stringent (use \(|r| > 2.5\), \(D > 1\)) - Large \(n\) (>100): More liberal (use \(|r| > 3\), \(D > 4/n\))
11.43 What to Do When Assumptions Fail
11.43.1 Non-Linearity
Solutions:
- Add polynomial terms: \(x^2\), \(x^3\)
- Include interactions: \(x_1 \times x_2\)
- Transform predictors: \(\log(x)\), \(\sqrt{x}\), \(1/x\)
- Use splines or generalized additive models (GAMs)
11.43.2 Heteroscedasticity
Solutions:
- Log transformation (most common): \(\log(y)\)
- Square root transformation (for counts): \(\sqrt{y}\)
- Box-Cox transformation: Optimal \(\lambda\)
- Weighted least squares: If variance pattern is known
- Robust standard errors: HC3, HC4 methods
11.43.3 Non-Normality
Solutions:
- Transformation (often fixes this and heteroscedasticity together)
- Check for outliers (may cause apparent non-normality)
- Bootstrap inference (doesn’t require normality)
- Robust regression (M-estimators, LAD regression)
- If \(n\) large: CLT applies, mild non-normality OK
11.43.4 Lack of Independence
Solutions:
- Mixed models: Random effects for clusters
- GLS: Generalized least squares with known correlation structure
- Cluster-robust SE: Adjust standard errors for clustering
- Time series methods: ARIMA, GARCH for temporal dependence
- Spatial models: Account for spatial correlation
11.44 Reporting Diagnostic Results
When reporting analyses, include:
Statement of assumptions checked: “Model assumptions (linearity, homoscedasticity, normality, independence) were verified using diagnostic plots.”
Any violations found: “Initial analysis revealed heteroscedasticity (funnel pattern in residual plot).”
Remedial actions taken: “Log transformation of the response variable successfully stabilized variance.”
Final model diagnostics: “Diagnostic plots for the final model showed no systematic violations of assumptions.”
Treatment of outliers (if any): “One observation (bird 27) was identified as highly influential (\(D = 0.85\)) and excluded from analysis due to documented illness.”
11.45 Summary and Key Takeaways
11.46 What We Learned This Week
Why diagnostics matter: Model assumptions must be verified, not just assumed. Violations can lead to biased estimates, invalid inference, and poor predictions.
Types of residuals:
- Raw: \(e_i = y_i - \hat{y}_i\)
- Standardized: \(e_i^* = e_i/\hat{\sigma}\)
- Studentized: \(r_i = e_i/(\hat{\sigma}\sqrt{1-h_{ii}})\) ← most useful
- Deleted: \(t_i\) (computed without observation \(i\))
Leverage vs. Influence:
- Leverage (\(h_{ii}\)): Potential for influence based on \(X\) values
- Influence (Cook’s \(D\)): Actual impact on fitted model
- High leverage ≠ high influence (example: beef marbling obs 9)
The four diagnostic plots:
- Residuals vs. Fitted → linearity, homoscedasticity
- Normal Q-Q → normality
- Scale-Location → homoscedasticity
- Residuals vs. Leverage → influential observations
When to transform:
- Heteroscedasticity (funnel pattern) → log or sqrt transformation
- Right skewness → log transformation
- Proportional errors → log transformation
- Don’t transform if violations are minor or interpretation suffers
Systematic workflow:
- Fit model → diagnostic plots → check assumptions → identify problems
- Investigate outliers/influential points → decide (keep/exclude/investigate)
- Apply fixes (transformation, WLS, etc.) → reverify diagnostics
11.47 Connections to Other Weeks
- Week 5: Gauss-Markov assumptions—this week we verify them
- Week 6: Multiple regression diagnostics more complex (leverage in \(p\)-dimensional space)
- Weeks 7-10: Same diagnostic principles apply to ANOVA and ANCOVA
- Week 12 (next): Rank deficiency affects leverage calculations
- Week 14: Preview of mixed models for handling dependence
11.48 Critical Reminders
ImportantKey Principles for Livestock Applications
- Biology first: Always ask if unusual observations are biologically plausible
- Document decisions: Record why you kept/excluded observations
- Report honestly: Disclose any data exclusions or transformations
- Verify after changes: Recheck diagnostics after any model modifications
- Interpretation matters: Choose transformations that stakeholders can understand
- Context is key: Same statistical issue may have different solutions in different biological contexts
Remember: Good diagnostics are not about achieving perfect plots—they’re about understanding your data, identifying problems, and making informed decisions about modeling and inference.
11.49 Looking Ahead
Next week (Week 12): We’ll tackle unequal subclass numbers and non-full rank models, where design matrices don’t have full rank. We’ll see how rank deficiency affects:
- Leverage calculations (\(\sum h_{ii} = r(X)\), not \(p\))
- Estimability of functions
- Interpretation of diagnostics
Week 14: Special topics including weighted least squares (WLS) for known heteroscedasticity and preview of mixed models for handling dependence.
11.50 Exercises
See Week11_Exercises.qmd for 7 exercises covering:
- Exercise 1: Hand calculation of diagnostic statistics (computational)
- Exercise 2: Leverage vs. influence concepts (theoretical)
- Exercise 3: Swine birth weight with two types of outliers (applied)
- Exercise 4: Systematic assumption checking (conceptual/applied)
- Exercise 5: Heteroscedasticity and transformation (applied)
- Exercise 6: Box-Cox transformation selection (computational/applied)
- Exercise 7: Complete diagnostic report for mastitis data (integrated)