BIOSCI220 Notes: Regression (Lecture 9) and Multivariate Analysis (Lecture 10)

Regression and Multivariate Analysis Notes (BIOSCI220)

Regression in R (Lecture 9)

Purpose: Build and interpret linear regression models using R, including how to handle continuous and categorical predictors, and how to assess model fit and uncertainty.

7.1 Fit linear regression models with different predictor structures

Single continuous explanatory variable
- Example: fit a model predicting billdepthmm from billlengthmm
- R syntax (example dataset penguins_nafree):
```
slm <- lm(bill_depth_mm ~ bill_length_mm, data = penguins_nafree)
summary(slm)$coef
```
- Interpretation: intercept and slope describe the mean billdepthmm as a linear function of billlengthmm.
Null model (no explanatory variable)
- Model: Yi = α + εi
- R syntax:
```
slm_null <- lm(bill_depth_mm ~ 1, data = penguins_nafree)
summary(slm_null)$coef
```
- Interpretation: intercept equals the average response; no predictor effect.
Single categorical explanatory variable (using dummy variables)
- Example: using species (Adelie, Chinstrap, Gentoo) as predictor
- R syntax (factor levels: Adelie as reference):
- Model: billdepthmm ~ species
- In the output, the intercept corresponds to the baseline (reference) species (Adelie), and the coefficients for other species reflect differences from the baseline.
A continuous and a categorical factor explanatory variable
- Model: billdepthmm ~ billlengthmm + species
- R syntax:
```
slm.sp <- lm(bill_depth_mm ~ bill_length_mm + species, data = penguins_nafree)
summary(slm.sp)$coef
```
- Interpretation: main effect of billlengthmm plus differences by species relative to the reference species, with the baseline intercept adjusted accordingly.
An interaction term in the explanatory variables
- Model: billdepthmm ~ billlengthmm * species
- R syntax:
```
slm.int <- lm(bill_depth_mm ~ bill_length_mm * species, data = penguins_nafree)
summary(slm.int)$coef
```
- Alternative specification (equivalent): billdepthmm ~ billlengthmm + species + billlengthmm:species
- Interpretation: interaction terms allow the slope of billdepthmm on billlengthmm to differ by species.
Data context
- penguinsnafree is the penguins dataset with missing values removed for the variables of interest (billdepthmm, billlength_mm, species).
- Key outputs used for interpretation: coefficients, standard errors, t-values, p-values, and overall model fit metrics.

7.2 Why include interaction terms?

Interactions test whether the effect of one predictor on the response depends on the level of another predictor.
In ecology/biology contexts, the relationship between a continuous measure (e.g., billlengthmm) and a response (billdepthmm) may differ across species or groups.
Including interactions helps uncover heterogeneous effects that a purely additive model (no interactions) would miss.

7.3 Differences between the operators : and * in an R model-fitting formula

: (colon) denotes an interaction term only, with no main effects
- Example: billlengthmm:species adds an interaction between billlengthmm and species but does not include billlengthmm or species as main effects by itself.
* (asterisk) expands to both main effects and the interaction
- Example: billlengthmm * species expands to billlengthmm + species + billlengthmm:species
Practical takeaway: use * when you want main effects plus their interaction; use : when you want only the interaction term(s).

7.4 Point predictions from linear regression models

Predictions from models with a single continuous predictor
- Example: predict for a new billlengthmm value using the slm model
- R: predict(slm, newdata = data.frame(bill_length_mm = 50))
Predictions with a continuous and a categorical predictor
- Example: predict from slm.sp (billdepthmm ~ billlengthmm + species) for a given species (e.g., Chinstrap) and billlengthmm
- R: newdata <- data.frame(bill_length_mm = 50, species = factor("Chinstrap", levels = c("Adelie","Chinstrap","Gentoo")))
 predict(slm.sp, newdata = newdata)
Predictions with an interaction term
- Example: predict from slm.int for a given billlengthmm and species, which includes the interaction term
- R: use the same predict function with an appropriate newdata frame that includes billlengthmm and species (levels must be valid for the model)
Example interpretation (from slides)
- For a model with both billlengthmm and species (no interaction): intercept and coefficients represent baseline and group differences.
- For a model with interaction (billlengthmm * species): predicted billdepthmm depends on both billlengthmm and species, with separate slopes by species.

7.5 Model selection using anova() and AIC()

Use anova() to compare nested models (models with the same response and a subset of predictors)
- Example: anova(slm_null, slm) compares a null model to a model with billlengthmm as a predictor.
Use AIC() to compare non-nested models or to compare multiple models at once; lower AIC indicates a better trade-off between fit and complexity
- Example: AIC(slm_null, slm, slm_sp, slm_int)
Practical guidance:
- Adding more predictors or interactions does not guarantee a better model; check AIC and diagnostic plots.
- Use model diagnostics to assess assumptions (normality, homoscedasticity, etc.).

7.6 Confidence intervals for population parameters

CI concept: a range of values around a population parameter (e.g., a regression coefficient) that is believed to contain the true value with a given probability (e.g., 95%).
Interpretation (typical phrasing):
- 95% CI for a parameter means: There is a 95% probability that the true parameter lies within the interval, or we are 95% confident that the interval contains the true parameter.
Computing CIs in R:
- For coefficients: confint(model) or confint(model, level = 0.95)
- You can specify a different confidence level, e.g., confint(slm_sp, level = 0.99)
Interpreting intercepts and coefficients in context:
- The intercept represents the mean of the response at the reference level of the categorical predictor (or the mean of the response when the continuous predictor is zero, depending on centering).
- Coefficients for categorical predictors represent differences from the reference category.
- In a model with multiple predictors, you can interpret predicted means for a given combination of predictor values and, with CIs, understand the uncertainty around those estimates.
Example interpretation from slides (using the penguin data):
- We are 95% confident that the average bill depth of an Adelie penguin is between 9.2 and 11.9 mm, given the other variables in the model.
- The average bill depth for Chinstrap penguins is between 1.5 and 2.4 mm shallower than that for Adelie penguins, given the other variables in the model.
Summary takeaway for Section 7:
- You can fit different regression structures, compare them with anova() and AIC(), and obtain predictions and confidence intervals for interpretation. The choice of model should balance explanatory power with parsimony and ensure assumptions are reasonable.

Multivariate Data Analysis (Lecture 10)

Focus: Techniques for analyzing data with many variables per observation, reducing dimensionality, and exploring relationships among samples.

8.1 What is multivariate data?

Description: Each observation has many variables; data are arranged with one row per observation and many columns for variables.
Challenge: Interpretation and inference can be more complex due to interdependencies among variables.
Purpose: Ordinate data to provide a useful approximation of relationships in a high-dimensional space.

8.2 Dimension reduction: why and how

Goal: Represent data in fewer dimensions while preserving as much information about relationships among observations as possible.
Why dimension reduction?
- Improves interpretability and visualization (2-D or 3-D plots).
- Facilitates communication of patterns and relationships.
Common approaches (examples): PCA, PCoA, cluster analysis, non-metric multidimensional scaling (nMDS).
Decision factors: data type, research questions, and desired interpretation.

8.3 Non-metric multidimensional scaling (nMDS)

What it is: A dimension-reduction tool that uses a distance matrix and rank-order information to place samples in a low-dimensional space (usually 2-D or 3-D).
Key features:
- Works with essentially any MV data.
- Aims to preserve the rank order of pairwise distances among samples when embedded in a few dimensions.
- Iterative algorithm that seeks a configuration minimizing Kruskal’s stress.
Core elements:
- Distance matrix: computes pairwise dissimilarities between samples (based on the original data).
- Rank-based preservation: only the order of distances matters, not the exact distances.
Output interpretation:
- Points close together in the 2-D plot indicate similar samples; farther apart indicates greater dissimilarity.

8.4 Distance measures: Euclidean vs Bray-Curtis

Euclidean distance (L2):
- Formula for two samples j and k across p variables:
- d_{ ext{Euc}}(j,k) = igg(rac{1}{p} igg)
  ight?
- Practical standard form (unweighted): d{ ext{Euc}}(j,k) = igg(igg(igl(x{1j}-x{1k}igr)^2 + igl(x{2j}-x{2k}igr)^2 + \n \, \cdots + igl(x{pj}-x_{pk}igr)^2igg)igg)^{1/2}
- Use: when data are continuous and on similar scales; sensitive to variable scales and zeros; standardization may be needed.
Bray-Curtis distance (BC):
- Formula (for samples j and k across variables i):
- $d{ ext{BC}}(j,k) = \frac{\bigg(\bigl|x{1j}-x{1k}\bigr| + \bigl|x{2j}-x{2k}\bigr| + \ \, \bigl|x{pj}-x{pk}\bigr|\bigg)}{\bigl(x{1j}+x{1k} + x{2j}+x{2k} + \ \, \bigl(x{pj}+x_{pk}\bigr)\bigr)}$
- Also written as: BC(j,k) = sumi |Yi,j - Yi,k| / sumi (Yi,j + Yi,k)
- Common in ecology/metagenomics where data are non-negative and often sparse (many zeros).
- Weighting: BC gives different emphasis to small vs large differences; double zeros are treated as uninformative (see later notes).
Practical guidance: choose distance measure based on data type, scale, sparsity, and research question.

8.5 Implementing nMDS in R and visualizing results

General workflow:
- Compute a distance matrix between samples using an appropriate measure (e.g., Euclidean with standardization or Bray-Curtis for ecological counts).
- Run nMDS on the distance matrix to obtain a 2-D (or 3-D) solution.
- Plot the 2-D coordinates (ordination space) and optionally add sample labels/colors by group.
- Assess goodness-of-fit via Kruskal’s stress statistic.
Example code outline (using vegan package in R):

  library(vegan)
  # data_matrix: observations x variables (MV data)
  # 1) distance matrix (e.g., Bray-Curtis)
  dist_mat <- vegdist(data_matrix, method = 'bray')
  # 2) run nMDS
  nmds_res <- metaMDS(dist_mat, k = 2)
  # 3) plot
  plot(nmds_res $points[,1], nmds_res$ points[,2], type='n')
  text(nmds_res $points[,1], nmds_res$ points[,2], labels = rownames(data_matrix))
  # 4) examine stress
  nmds_res$stress

Interpretation: The two NMDS axes (NMDS1, NMDS2) are non-metric representations of dissimilarities between samples; distances on the plot reflect ranked relationships (not exact values).

8.6 How to interpret and diagnose nMDS outputs

Stress (Kruskal’s stress): a goodness-of-fit measure for the NMDS solution.
- Lower stress indicates a better representation of the rank-order distances in the reduced dimensions.
- Guidelines (typical interpretation):
- Stress < 0.05: excellent
- 0.05 ≤ Stress < 0.10: good
- 0.10 ≤ Stress < 0.20: fair to moderate
- > 0.20: poor fit (caution in interpretation)
Global vs local minima:
- NMDS solutions may converge to local minima; multiple random starts can help locate a global optimum.
- Visual inspection and multiple runs help assess stability.
Kruskal’s stress concept (under the hood):
- The algorithm iteratively adjusts sample positions to minimize the discrepancy between observed distances and ordination distances, quantified by stress.
Practical implications for analysis:
- The chosen distance metric and data preprocessing (scaling, standardization, zero handling) affect the NMDS results.
- Always report stress and consider biological/interpretive relevance of the ordination pattern rather than absolute distances.

8.7 Distance measures in practice and taking notes from hands-on slides

Euclidean vs Bray-Curtis in practice:
- Euclidean distance is intuitive for continuous, standardized data and is sensitive to scale.
- Bray-Curtis is widely used for ecological data with many zeros and non-negative counts; it emphasizes proportional differences and is robust to some zeros.
Distances and data transformations:
- For Euclidean, standardizing or normalizing variables can prevent dominated dimensions by highly variable variables.
- For Bray-Curtis, data typically remain non-negative; total counts can influence BC values, so preprocessing may be needed depending on study design.

8.8 Practical takeaways from the module

You can create and compare many models and distance-based representations; choose methods that align with data types and research questions.
Model selection balances fit and complexity; model diagnostics are essential to validate assumptions.
nMDS provides a flexible, rank-based visualization for multivariate relationships but requires careful interpretation and reporting of stress values.
The concepts of distance, rank, and ordination underpin many exploratory data analysis approaches in biology and ecology.
Administrative and reflection notes:
- The module includes a practical H5P tutorial on nMDS with interactive quizzes to reinforce understanding.
- The emphasis is on interpreting multivariate relationships rather than on exact numerical equivalences.