Great question—this is such a common hurdle when moving from traditional ANOVA to more flexible statistical methods, especially with agricultural data like silage properties that often don’t fit nice normal distributions!

Let’s break this down clearly:

First, a quick clarification: LM vs. GLM #

The aov() function you’ve been using is actually a special case of a Linear Model (LM), which falls under the umbrella of Generalized Linear Models (GLMs). Traditional LMs (and ANOVA, which is LM for categorical predictors) do require normality of residuals and homogeneity of variance—this is why your Shapiro-Wilk and Levene’s tests were flagging problems.

But when people talk about using "GLM" as an alternative to ANOVA (like with R’s glm() function), we’re referring to the broader class of Generalized Linear Models—and these have very different assumptions.

Key Assumptions for GLMs (the flexible kind) #

GLMs do not require the dependent variable (your silage properties) to be normally distributed, nor do they assume constant variance (homoscedasticity). Instead, they rely on these core assumptions:

• Your dependent variable follows a distribution from the exponential family (e.g., Normal, Gamma, Poisson, Binomial, Beta—each fits different data types)
• There’s a linear relationship between the predictors (treatment, time point) and the transformed mean of the dependent variable (via a "link function")
• Observations are independent of each other (critical for all statistical models—make sure your experimental design doesn’t have clustered data without accounting for it)
• The variance of the dependent variable is related to its mean in a predictable way (e.g., Poisson: variance = mean; Gamma: variance ∝ mean²)

Applying This to Your Silage Experiment #

For your 18 physicochemical properties, here’s how to approach this:

1. Assess each variable’s distribution:
   • If a variable is continuous but right-skewed (common for things like pH, nutrient concentrations), a Gamma GLM with a log link is often a good fit.
   • If you have proportion data (e.g., dry matter percentage, 0-1 range), use a Beta regression (you’ll need the betareg package).
   • If you have count data (e.g., microbial counts), go with Poisson or Negative Binomial GLMs (Negative Binomial handles overdispersion better than Poisson).
2. Validate GLM assumptions:
   • Instead of Shapiro-Wilk, plot the model’s residuals against fitted values—you want residuals to be randomly scattered with no clear pattern.
   • Use diagnostic plots specific to your chosen distribution (e.g., for Gamma GLMs, check that Pearson residuals don’t show increasing variance with fitted values).
3. Post-hoc tests for GLMs:
   • Forget HSD.test()—the emmeans package is the gold standard for post-hoc comparisons with GLMs. It works with all exponential family distributions and lets you adjust for multiple comparisons (like Tukey) easily.

Quick R Code Example (Gamma GLM for a skewed silage variable) #

# Load required packages
library(glm2) # For stable GLM fitting (or use base glm())
library(emmeans)

# Fit Gamma GLM with treatment and time point as predictors
gamma_model <- glm(volatile_acids ~ treatment * time_point, 
                   data = silage_data, 
                   family = Gamma(link = "log"))

# Check model summary
summary(gamma_model)

# Post-hoc: Compare treatments at each time point (Tukey-adjusted)
treatment_emmeans <- emmeans(gamma_model, ~ treatment | time_point)
pairs(treatment_emmeans, adjust = "tukey")

Final Note #

Don’t force a GLM on every variable—if some of your properties do meet ANOVA assumptions, sticking with aov() + Tukey is totally fine. The goal is to match the model to your data, not the other way around!

内容的提问来源于stack exchange，提问作者Mark