Linear Model Review

\[Y_i = \beta_0 + \beta_1 X_i + \epsilon_i\]

\(Y_i\) = dependent variable (there is only 1)

\(X_i\) = independent variable(s) (there may be many)

\(B_0\) = model intercept, the overall mean of \(Y\)

\(B_1\) = how \(Y\) changes with \(X\)

\(\epsilon_i\) = model residual, the gap between the predicted value for \(Y_i\) and its observed value

The Distribution of \(Y\)

Technically in a linear model, \(\bar{y_{i.}}\), the mean value for a population, follows a normal distribution for a set of known parameters:

\[ \bar{y_{i.}} \sim 𝑁(x_i \beta, \sigma^2/r_i) \]

We often cannot observe this distribution well with our limited number of reps. But we if repeated a study an infinite number of times, a normal distribution for \(\bar{y_{i.}}\) would likely be observed.

Example Analyses 1

Example Analyses 2

Example Analyses 3

Early Lessons

• Often, the residuals are normally distributed, and hence model expectations are met!

• Linear models are considered robust to modest departures from normality; that is, the estimates will be fine and the results from ANOVA will be correct.

• The Shapiro-Wilk test is not great for asessing normality; it’s very sensitive to the size of the data set; large data sets may fail just based on sampling variation, small data sets pass because there are too few observations to adequately test.

Sampling Variation & Normality

The function rnorm() can be used to generate random values that follow a normal distribution, that we can test for normality using shapiro.test().

    Shapiro-Wilk normality test

data:  rnorm(60)
W = 0.98675, p-value = 0.7603

Sampling Variation & Normality

The Central Limit Theorum

In probability theory, the central limit theorem states that, under appropriate conditions, the distribution of a normalized version of the sample mean converges to a standard normal distribution.
–Wikipedia

What Data Are Truly Not Normal?

• Count data (unless the values are very high, e.g. hundreds)
  
• Data that lack variability; most data points are a single value (emergence and crop stand data)
• Binomial data (yes/no)
• Data with categorical outcomes (rating scales)
• Data that have many zeros
• Extremely skewed data

Count data

• Things you count: insects, calves born, trips to the grocery store (requires a generalized linear model)

No Data Variability

Binomial data

Alternatives for Non-normal data

Categorical

• Likert scales (strong disagree to agree)
• disease rating (very resistant to very tolerant)
• special models: multinomial, ordinal (very complicated)

Many zero’s

• common in count data
• common in relative abundance data from metagenomics
• zero-inflated models (complicated)

Skewed Distributions

• extreme tails on the left or more commonly, the right side of histograms
• no single solution, many other distributions that support skewed distributions: beta (for percentage), lognormal, exponential, tweedie,… (complicated)

Skewed Distributions

]

Alternatives for Non-normal Data

Non-parametric hypothesis tests (No distributional assumptions)

• Mann-Whitney \(U\) test/Wilcoxon rank-sum test for 2 samples: wilcox.test()
• Kruskall-Wallis for 3+ samples: kruskal.test()
• Wilcoxon signed-rank test for paired samples: wilcox.test(..., paired = TRUE)

Drawbacks

• only conducting a hypothesis test, no estimation

Alternatives for Non-normal Data

Permutations

1. Calculate the test statistic from your data (e.g. difference in treatment means).

2. Simulate data in which the test statistic is null (e.g. the difference in treatment means is zero).

3. Find out how often your test statistic occurs in the simulated data; the probability of this occurring is your p-value.

Permutations

Andrew Heiss’ Null Worlds

Permutation Example

Alfalfa data set

• split plot: main plot is date of harvest (4 levels), split plot is variety (3 levels)

data("Alfalfa", package = "nlme")

alf_mod <- lme(Yield ~ Variety + Date,
               random = ~1|Block/Date,
               data = Alfalfa)

Alfalfa Permuatation Example

Permutations

sim_func <- function(){
  N = nrow(Alfalfa)
  Alfalfa$yield2 <- sample(Alfalfa$Yield, size = N, replace = FALSE)
  try(perm_mod <- lme(yield2 ~ Variety + Date,
               random = ~1|Block/Date,
               data = Alfalfa))
  if(exists("perm_mod")) return(perm_mod)
}

perms <- replicate(500, sim_func(), simplify = FALSE)

Permutated ANOVA

perm_anova <- lapply(perms, anova, type = "marginal") |>  # Conduct anova
  dplyr::bind_rows() |> # some data formatting
  tibble::rownames_to_column(var = "variable") |> 
  mutate(variable = stringr::str_extract(variable, "[A-Za-z]+")) |> 
  dplyr::filter(variable != "Intercept") # remove the intercept (optional step)

  variable numDF denDF  F-value   p-value
1  Variety     2    46 1.076281 0.3492884
2     Date     3    15 1.607674 0.2294616
3  Variety     2    46 1.369906 0.2643028
4     Date     3    15 1.373157 0.2891074

(alf_aov <- anova(alf_mod, type = "marginal"))

            numDF denDF   F-value p-value
(Intercept)     1    46 209.03665  <.0001
Variety         2    46   1.91760  0.1585
Date            3    15  14.09297  0.0001

Permutated ANOVA

Additional Resources

• Zero Inflated Models and Generalized Linear Mixed Models with R (2012) by Alain Zuur, Anatoly Saveliev, and Elena Ieno

• GLMM FAQs by Ben Bolker (more resources on mixed models from Dr. Bolker)