1 Purpose

This lesson follows the data-exploration protocol of Zuur, Ieno, and Elphick (2010).

The goal is to understand the data and identify issues that could affect a later model. This is exploratory analysis, not a confirmatory test of whether vegetation causes changes in sparrow abundance. Keep biological hypotheses and the sampling design in view when deciding what to model.

1.1 Learning goals

  1. Load and inspect a dataset provided by an R package.
  2. Check row structure, missingness, ranges, and unusual observations.
  3. Visualize the response and predictor distributions with ggplot2.
  4. Screen vegetation predictors for collinearity without using a mechanical deletion rule.
  5. Inspect response–predictor patterns and consider whether a count model may be appropriate.

1.2 Packages

The car package provides the VIF diagnostic used below.

library(ggversa)
## Warning: package 'ggversa' was built under R version 4.5.3
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.5.3
## Warning: package 'ggplot2' was built under R version 4.5.3
## Warning: package 'tidyr' was built under R version 4.5.3
## Warning: package 'dplyr' was built under R version 4.5.3
## Warning: package 'lubridate' was built under R version 4.5.3
library(car)

2 Step 1: Load the data

VegSamplesV1 contains vegetation measurements and the number of saltmarsh sparrows banded at sampled plots. Load the dataset directly from ggversa; this avoids a machine-specific file path and makes the report easier to reproduce.

data("VegSamplesV1", package = "ggversa")
veg <- as_tibble(VegSamplesV1)

dim(veg)
## [1] 60 19
names(veg)
##  [1] "Year"                "Site"                "UniversalPlotName"  
##  [4] "Banded"              "PtCountsum"          "Avgmaxht"           
##  [7] "Avgdens"             "ht.thatch"           "S.patens"           
## [10] "Distichlis"          "S.alternifloraShort" "S.alternifloraTall" 
## [13] "Juncus"              "Bare"                "Other"              
## [16] "Phragmites"          "Shrub"               "Tallsedge"          
## [19] "Water"

There should be 60 rows and 19 columns. Each row represents a sampled plot. Variables include year and site identifiers, the response Banded, vegetation structure, and percent-cover measurements.

3 Step 2: Check the data structure

Start by looking at a few rows and the column types. Names such as S.patens contain periods; these are valid R names and can be used directly in tidyselect or formulas.

slice_head(veg, n = 6)
glimpse(veg)
## Rows: 60
## Columns: 19
## $ Year                <int> 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 20…
## $ Site                <fct> "EastRiverMarsh", "EastRiverMarsh", "EastRiverMars…
## $ UniversalPlotName   <fct> 02PLOT012, 02PLOT013, 02PLOT014, 02PLOT015, 02PLOT…
## $ Banded              <int> 31, 40, 15, 27, 20, 17, 17, 29, 13, 37, 35, 47, 27…
## $ PtCountsum          <int> 53, 40, 15, 33, 17, 30, 49, 39, 22, 23, 41, 26, 23…
## $ Avgmaxht            <dbl> 33.19, 48.22, 33.81, 28.34, 28.75, 24.19, 39.34, 3…
## $ Avgdens             <dbl> 23.30, 27.63, 26.40, 45.63, 38.88, 53.60, 18.20, 4…
## $ ht.thatch           <dbl> 4.67, 6.56, 3.33, 7.33, 3.44, 3.44, 2.11, 5.89, 5.…
## $ S.patens            <dbl> 21.31, 11.67, 29.09, 52.22, 47.84, 60.42, 17.35, 5…
## $ Distichlis          <dbl> 41.47, 52.78, 17.80, 19.76, 12.78, 0.01, 2.67, 10.…
## $ S.alternifloraShort <dbl> 32.22, 23.61, 35.89, 5.56, 15.49, 26.44, 42.56, 20…
## $ S.alternifloraTall  <dbl> 2.78, 9.72, 0.00, 3.33, 10.00, 0.00, 25.03, 6.44, …
## $ Juncus              <dbl> 0.00, 0.00, 2.22, 7.78, 0.00, 0.00, 0.56, 5.56, 0.…
## $ Bare                <dbl> 1.67, 0.00, 8.33, 2.22, 10.14, 2.01, 11.29, 1.78, …
## $ Other               <dbl> 0.56, 0.00, 0.00, 8.89, 3.75, 0.00, 0.56, 0.00, 1.…
## $ Phragmites          <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.…
## $ Shrub               <dbl> 0.00, 0.00, 6.67, 0.00, 0.00, 0.00, 0.00, 0.00, 0.…
## $ Tallsedge           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ Water               <dbl> 0.00, 2.22, 0.00, 0.25, 0.00, 11.11, 0.00, 0.00, 4…

Check whether the identifiers and numeric variables have sensible types. For example, Site and UniversalPlotName should be text, while Banded and vegetation measurements should be numeric.

veg %>%
  count(Year, Site, name = "n_plots") %>%
  arrange(Year, Site)
veg %>%
  summarise(
    n_rows = n(),
    n_unique_plot_names = n_distinct(UniversalPlotName),
    duplicated_plot_names = n_rows - n_unique_plot_names
  )

If a plot name is not unique, determine whether it is expected to recur across years before treating it as a duplicate. Repeated plots or sites can also mean observations are not independent; keep that in mind for later modelling.

4 Step 3: Check missing values and plausible ranges

Summarize missing values by column. Missingness is not automatically an error, but it changes which observations are available for each calculation.

missing_by_column <- veg %>%
  summarise(across(everything(), ~ sum(is.na(.x)))) %>%
  pivot_longer(
    cols = everything(),
    names_to = "variable",
    values_to = "n_missing"
  ) %>%
  filter(n_missing > 0) %>%
  arrange(desc(n_missing))

missing_by_column

Check numeric ranges. Percent-cover variables should be between 0 and 100. Do not assume all cover columns for a plot must sum to exactly 100: survey definitions and rounding may affect the total.

cover_vars <- c(
  "S.patens", "Distichlis", "S.alternifloraShort",
  "S.alternifloraTall", "Juncus", "Bare", "Other",
  "Phragmites", "Shrub", "Tallsedge", "Water"
)

veg %>%
  summarise(across(
    all_of(cover_vars),
    list(min = ~ min(.x, na.rm = TRUE), max = ~ max(.x, na.rm = TRUE))
  )) %>%
  pivot_longer(
    cols = everything(),
    names_to = c("variable", ".value"),
    names_sep = "_"
  )

If a value is outside its expected range, check the source record and measurement protocol. Do not silently clip or remove it.

5 Step 4: Explore the response (Banded)

First inspect the range, centre, spread, and number of zeros. Banded is a count, so its distribution does not need to be normal. A histogram and a Cleveland-style dot plot help identify unusual values without declaring them errors.

veg %>%
  summarise(
    n = n(),
    n_missing = sum(is.na(Banded)),
    minimum = min(Banded, na.rm = TRUE),
    median = median(Banded, na.rm = TRUE),
    mean = mean(Banded, na.rm = TRUE),
    maximum = max(Banded, na.rm = TRUE),
    n_zero = sum(Banded == 0, na.rm = TRUE),
    proportion_zero = mean(Banded == 0, na.rm = TRUE)
  )
ggplot(veg, aes(x = Banded)) +
  geom_histogram(binwidth = 5, boundary = 0, color = "white", fill = "#367A8A") +
  labs(
    title = "Distribution of banded saltmarsh sparrows",
    x = "Number banded per plot",
    y = "Number of plots"
  ) +
  theme_minimal(base_size = 12)

veg %>%
  mutate(plot_order = row_number()) %>%
  ggplot(aes(x = plot_order, y = Banded)) +
  geom_point(alpha = 0.8, color = "#B8573C", size = 2) +
  labs(
    title = "Banded sparrows by row in the dataset",
    x = "Row order in the data",
    y = "Number banded per plot"
  ) +
  theme_minimal(base_size = 12)

The second graph uses row order only to locate observations for follow-up. It does not imply a trend over time. If a value looks unusual, inspect its plot, year, site, and source record. A point outside a boxplot whisker—or far from the other points—is a prompt to investigate, not a reason to delete it.

Check how often each count occurs, including zero:

veg %>%
  count(Banded, name = "n_plots") %>%
  ggplot(aes(x = Banded, y = n_plots)) +
  geom_col(fill = "#367A8A") +
  scale_x_continuous(breaks = scales::breaks_pretty()) +
  labs(
    title = "Frequency of each banded-sparrow count",
    x = "Number banded per plot",
    y = "Number of plots"
  ) +
  theme_minimal(base_size = 12)

6 Step 5: Explore predictor distributions

These plots show the shape and scale of vegetation measurements. Free scales make each panel readable, but axes differ between panels; compare shapes within each panel, not bar heights across variables.

predictor_vars <- c("Avgmaxht", "Avgdens", "ht.thatch", cover_vars)

veg %>%
  select(all_of(predictor_vars)) %>%
  pivot_longer(
    cols = everything(),
    names_to = "variable",
    values_to = "value"
  ) %>%
  ggplot(aes(x = value)) +
  geom_histogram(bins = 20, color = "white", fill = "#6E9F76") +
  facet_wrap(vars(variable), scales = "free", ncol = 4) +
  labs(
    title = "Distributions of vegetation predictors",
    x = "Observed value",
    y = "Number of plots"
  ) +
  theme_minimal(base_size = 10)

Several cover measurements may have many zeros or be strongly right-skewed. That can be ecologically meaningful (for example, a plant may be absent from many plots). Do not transform variables solely to make their raw histograms look normal; model assumptions concern the response distribution and, for a Gaussian linear model, conditional residuals.

7 Step 6: Screen for collinearity among predictors

Collinearity means that predictors contain overlapping information. Start with a correlation heat map. This is a screening tool: pairwise correlation does not capture every multivariable relationship, and correlations can be sensitive to extreme values or many zeros.

cor_mat <- veg %>%
  select(all_of(predictor_vars)) %>%
  cor(use = "pairwise.complete.obs")

cor_long <- as.data.frame(as.table(cor_mat)) %>%
  as_tibble() %>%
  rename(var_x = Var1, var_y = Var2, r = Freq) %>%
  mutate(
    var_x = factor(var_x, levels = predictor_vars),
    var_y = factor(var_y, levels = predictor_vars)
  ) %>%
  filter(as.integer(var_x) < as.integer(var_y))

ggplot(cor_long, aes(x = var_x, y = var_y, fill = r)) +
  geom_tile(color = "white") +
  geom_text(aes(label = sprintf("%.2f", r)), size = 2.6) +
  scale_fill_gradient2(
    low = "#3B78A1", mid = "white", high = "#B8573C",
    midpoint = 0, limits = c(-1, 1), name = "Pearson r"
  ) +
  coord_equal() +
  labs(
    title = "Pairwise correlations among vegetation predictors",
    x = NULL,
    y = NULL
  ) +
  theme_minimal(base_size = 10) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid = element_blank()
  )

Percent-cover variables are constrained measurements: increasing one cover category can reduce the relative cover available to others. Strong negative correlations may therefore reflect the composition of the sample, not a direct biological interaction.

7.1 Calculate variance inflation factors (VIF)

The supplied Appendix S1 script uses corvif() from HighstatLib.r to examine the vegetation predictors before simplifying its regression. Here we use car::vif() for the same diagnostic. With numeric predictors in an ordinary linear model, VIF measures how strongly each predictor can be explained by the other predictors together. It can reveal redundancy that is not obvious from pairwise correlations.

For predictor j,

\[ \mathrm{VIF}_j = \frac{1}{1-R_j^2}, \]

where R-squared comes from regressing that predictor on all the other predictors. A VIF of 1 means no linear redundancy. A VIF of 5 means the estimated coefficient variance is five times what it would be for an orthogonal predictor with the same spread; its standard error is inflated by the square root of 5. Values above 3–5 deserve attention in this exercise, but these are discussion points rather than universal pass/fail rules.

7.1.1 Fit the full model and inspect VIF

Use the same 14 vegetation predictors as the supplied example. Keep the same complete observations for all the comparisons below so a change in VIF reflects a change in predictors rather than a different set of rows.

veg_model <- veg %>%
  select(Banded, all_of(predictor_vars)) %>%
  drop_na()

m_full <- lm(
  Banded ~ Avgmaxht + Avgdens + ht.thatch +
    S.patens + Distichlis + S.alternifloraShort +
    S.alternifloraTall + Juncus + Bare + Other +
    Phragmites + Shrub + Tallsedge + Water,
  data = veg_model
)

vif_full <- car::vif(m_full)

vif_full_table <- tibble(
  variable = names(vif_full),
  VIF = as.numeric(vif_full)
) %>%
  arrange(desc(VIF))

vif_full_table

In the uploaded data, S.alternifloraTall and S.patens have VIFs of approximately 159, followed by S.alternifloraShort at 121. These variables contain substantial overlapping information in this predictor set. A high VIF does not mean a variable is ecologically unimportant or has a weak association with sparrow counts.

ggplot(vif_full_table, aes(x = VIF, y = reorder(variable, VIF))) +
  geom_point(size = 3, color = "#367A8A") +
  geom_vline(xintercept = 5, linetype = "dashed", color = "#B8573C") +
  labs(
    title = "Collinearity in the full vegetation model",
    subtitle = "Dashed line: VIF = 5, a discussion guide for this exercise",
    x = "Variance inflation factor",
    y = NULL
  ) +
  theme_minimal(base_size = 11)

7.1.2 Remove one candidate and reassess

To begin the discussion, try removing S.alternifloraTall. This choice follows the reduced predictor set in the supplied Appendix example. Its VIF is almost tied with S.patens, so the tiny difference between them is not a scientific reason to prefer one. In a real analysis, decide which variable to retain using the ecological question, measurement quality, and the habitat feature it represents.

# update() keeps the existing response and predictors except the named term.
m_no_tall <- update(m_full, . ~ . - S.alternifloraTall)

vif_no_tall <- car::vif(m_no_tall)

vif_comparison <- vif_full_table %>%
  rename(full_model = VIF) %>%
  left_join(
    tibble(
      variable = names(vif_no_tall),
      without_tall = as.numeric(vif_no_tall)
    ),
    by = "variable"
  ) %>%
  arrange(desc(full_model))

vif_comparison

After this single removal, the VIF for S.patens falls from about 159.35 to 4.43 and the VIF for Juncus falls from 45.00 to 1.63. This is why we reassess the whole predictor set after each change. Do not remove every variable that had a large VIF in the original model.

Avgmaxht now has the largest VIF, about 6.11. It is a candidate for discussion, not an automatic deletion: if vegetation height is central to the research question, consider retaining it and reassessing which overlapping measurements are needed instead.

7.1.3 A brief look at the reduced set in the original example

The supplied script also omits Avgmaxht and Avgdens from its reduced vegetation model. The following optional teaching comparison reproduces that predictor set. It shows the result of those choices without introducing a full variable-selection procedure.

m_reduced_example <- update(m_no_tall, . ~ . - Avgmaxht )

vif_reduced <- car::vif(m_reduced_example)

vif_comparison <- vif_comparison %>%
  left_join(
    tibble(
      variable = names(vif_reduced),
      original_reduced_set = as.numeric(vif_reduced)
    ),
    by = "variable"
  )

vif_comparison %>%
  mutate(across(where(is.numeric), ~ round(.x, 2)))

An NA in this table means the predictor was omitted from that model. In the uploaded data, the largest VIF in this reduced set is about 2.13. This indicates less collinearity; it does not establish better prediction, a better ecological explanation, or a final model. The linear models here reproduce the original collinearity example. The response distribution and sampling design still need attention before making inferences.

7.2 From collinearity to variable selection: the next conversation

We now have an initial reason to discuss which variables to retain. Ask:

  • Which habitat features are essential to the biological question?
  • Which predictors measure similar information, and which measurement should represent it?
  • Are some predictors too sparsely observed to estimate their effects reliably?
  • Is the proposed model too complex for 60 sampled plots?

For this exercise, identify one candidate to remove, give a biological or measurement-based reason, and recalculate VIF. Keep a short record of the decision. Formal model selection can be introduced later; VIF alone does not rank ecological importance, and repeatedly removing predictors until p-values become small would make the resulting evidence misleading.

8 Step 7: Inspect relationships with the response

Plot a small set of biologically motivated predictors against Banded. The smooth curves describe patterns in these 60 plots; they are not evidence of causation and do not account for repeated sampling or other covariates.

relationship_vars <- c("Juncus", "Shrub", "S.patens", "Avgmaxht")

veg %>%
  select(Banded, Site, all_of(relationship_vars)) %>%
  pivot_longer(
    cols = all_of(relationship_vars),
    names_to = "predictor",
    values_to = "value"
  ) %>%
  ggplot(aes(x = value, y = Banded)) +
  geom_point(aes(color = Site), alpha = 0.8, size = 2) +
  geom_smooth(
    aes(group = 1),
    method = "loess",
    formula = y ~ x,
    se = TRUE,
    color = "#222222",
    linewidth = 0.7
  ) +
  facet_wrap(vars(predictor), scales = "free_x") +
  labs(
    title = "Exploratory relationships with banded sparrows",
    x = "Predictor value",
    y = "Number banded per plot",
    color = "Site"
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom")
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at -0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 0.0022278
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at -0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger
## Warning: Failed to fit group 1.
## Caused by error in `predLoess()`:
## ! NA/NaN/Inf in foreign function call (arg 5)

Look for curvature, groups, uneven spread, and observations that dominate the pattern. If plots were sampled repeatedly at the same site, a site-level structure may need to be considered in the final model.

9 Step 8: Check count-model considerations

Because Banded is a non-negative count, a Poisson or negative-binomial model may be more suitable than ordinary least squares. Before choosing a model, examine the sampling process, zeros, independence, and whether counts represent comparable effort. A high variance-to-mean ratio can flag overdispersion, but it is not a complete model diagnostic by itself.

veg %>%
  summarise(
    n = sum(!is.na(Banded)),
    mean_count = mean(Banded, na.rm = TRUE),
    variance_count = var(Banded, na.rm = TRUE),
    variance_to_mean = variance_count / mean_count,
    proportion_zero = mean(Banded == 0, na.rm = TRUE)
  )

If you later fit a model, begin with a small, pre-specified set of predictors rather than all 14 vegetation variables in a 60-row dataset. For example, this Poisson model illustrates syntax; it is not a final recommended analysis:

count_model <- glm(
  Banded ~ Juncus + Shrub,
  family = poisson(link = "log"),
  data = veg
)

summary(count_model)

# A rough overdispersion check; values well above 1 deserve investigation.
sum(residuals(count_model, type = "pearson")^2) /
  df.residual(count_model)

Check the final model’s residuals and design assumptions, and consider a negative-binomial model if overdispersion remains. If sampling effort differs among plots, include an appropriate effort term or offset. If sites or plots are revisited, account for that dependence rather than treating all rows as independent.

10 Step 9: Summarize the exploration

Before moving on, write down short answers to these questions:

  1. Are there missing values, duplicate identifiers, or values outside plausible ranges? How did you verify them?
  2. Is Banded sparse or strongly skewed? How many zeros are present?
  3. Which predictors are strongly correlated, and could the pattern arise because cover values are compositional?
  4. Do plots from the same site or year appear to form groups?
  5. Which small set of predictors is justified by your ecological question before looking at model p-values?
  6. What model family and dependence structure fit the response and sampling process?

Keep a record of any exclusions or transformations and the reason for each. Exploration should guide the analysis and reveal data problems; it should not be used to search repeatedly for a model that produces a desired p-value.

11 References and data source