When to trust it, and how to fix it when it’s broken
The past several chapters/lessons have focused on linear models. Here we explore options for analysis when linear model assumptions are not met. In doing so, we review several different approaches, many of which could be (and are) the focal topic of courses and books. The goal here is to
help you figure out when you may need to employ these methods by having you understand what they do/fix
give you a start on doing these analyses in R
Citations to relevant articles papers and books are noted throughout for further exploration.
Linear models are useful for a number of reasons. They are a great way to unify most/many tests from classical statistics. In fact, most of the ranked tests we’ve developed can actually be run as linear models when n >15. For example, we can go back to our Wilcox-Mann Whitney U tests (for 2 populations) and Kruskal-Wallis (for 3+) from the comparing means among groups chapter and note the outcome from a wilcox.test
two_species_subset <- iris[iris$Species!="setosa",]
wilcox.test(Sepal.Length ~ Species, two_species_subset)
Wilcoxon rank sum test with continuity correction
data: Sepal.Length by Species
W = 526, p-value = 5.869e-07
alternative hypothesis: true location shift is not equal to 0is very close to a linear model predicting the signed rank of the data
library(car)
signed_rank = function(x) sign(x) * rank(abs(x))
Anova(lm(signed_rank(Sepal.Length) ~ Species, two_species_subset), type="III")Anova Table (Type III tests)
Response: signed_rank(Sepal.Length)
Sum Sq Df F value Pr(>F)
(Intercept) 64872 1 102.378 < 2.2e-16 ***
Species 20967 1 33.089 1.003e-07 ***
Residuals 62098 98
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In fact, we could run simulations and show that p values from these 2 approaches are highly correlated (now you know what that means) with a \(\beta\) of almost 1 (from Lindelov (n.d.)).
Linear models are also extremely robust. Consider the basic assumptions of a linear model
\[ \epsilon \approx i.i.d.\ N(\mu,\sigma) \]
Although the residuals are meant to be homoscedastic (equal or constant across all groups), it turns out the model is robust of when the largest group variance is 4-10x larger than the smallest group variance and sample sizes are approximately equal (Blanca et al. 2018; Fox 2015; A. F. Zuur, Ieno, and Elphick 2010), though highly uneven group sizes begin to cause issues (Blanca et al. 2018).
Similarly, non-normal data is not an issue in an of itself (Knief and Forstmeier 2021). This is partially because the assumptions are focused on residuals, but also because the procedure is highly robust (Blanca et al. 2017). This finding further supports the graphical consideration of assumptions , especially since many tests for normality are conservative (samples are almost never perfectly normal, and slight deviations are easier to pick up with larger sample sizes despite the fact the central limit theorem suggests this is when the issue is least important) (A. F. Zuur, Ieno, and Elphick 2010; Shatz 2023).
This is a common question. Sometimes the outcome variable is not normally distributed due to the impact of various factors. Consider this example of sparrow weight from A. F. Zuur, Ieno, and Elphick (2010). The raw data may look like this…
Sparrows <- read.table(file="data/SparrowsElphick.txt", header=TRUE)
Sparrows$fMonth<-factor(Sparrows$Month,
levels = c(5, 6, 7, 8, 9, 10),
labels = c("May", "June", "July", "August",
"Sept.", "Oct."))
Sparrows$I1 <- Sparrows$fMonth =="June" |
Sparrows$fMonth =="July" |
Sparrows$fMonth =="August"
library(ggplot2)
ggplot(Sparrows[Sparrows$I1== T, ], aes(x=wt))+
geom_histogram()+
labs( x = "Weight (g)", main = "Weight of sparrows in study", y = "Frequency")
which is skewed (albeit slightly). However, if we consider the impact of month…
library(ggplot2)
ggplot(Sparrows[Sparrows$I1== T, ], aes(x=wt, fill=fMonth))+
geom_histogram()+
labs( x = "Weight (g)", main = "Weight of sparrows in study", y = "Frequency",
fill= "Month")+
facet_col(~fMonth)
We would find the residuals to be much better behaved.
Similary, highly skewed data may arise due to unequal sample sizes (which may pose their own problems, but not outright), but models using this data may meet assumptions.
set.seed(8)
example_skewed <- data.frame(Population= c(rep("a",40),
rep("b", 30),
rep(letters[3:4], each=15)),
Growth = c(rnorm(40,1),
rnorm(30, 4),
rnorm(15, 7),
rnorm(15,10)))
ggplot(example_skewed,aes(x=Growth))+
geom_histogram()+
labs(y="Some value",
title="Example of skewed data appearing due to unequal sample size")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(example_skewed,aes(x=Growth, fill=Population))+
geom_histogram()+
labs(y="Some value",
title="Example of skewed data appearing due to unequal sample size")+
facet_col(~Population)`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
plot(lm(Growth~Population, example_skewed))
Anova(lm(Growth~Population, example_skewed), type="III") Anova Table (Type III tests)
Response: Growth
Sum Sq Df F value Pr(>F)
(Intercept) 38.75 1 32.436 1.344e-07 ***
Population 1007.71 3 281.175 < 2.2e-16 ***
Residuals 114.69 96
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1These issues have led some authors to argue violations of the linear model assumptions are less dangerous than trying new, and often less-tested, techniques that may inflate type I error rates (Knief and Forstmeier 2021; D. I. Warton et al. 2016). However, new techniques (or sometimes old techniques that are new to a given field) may be more appropriate when assumptions are clearly broken (D. Warton and Hui 2010; Reitan and Nielsen 2016; Geissinger et al. 2022). In this section we explore common-ish approaches for analyzing data when the assumptions of a linear model are broken. Our goal here is to introduce these methods. Each could be the focus of their own class, book, or much larger study. Fortunately most can be viewed as extensions to our existing knowledge, although many of the assumptions and techniques for them are less developed/still being argued about.
The various assumptions related to linear models may be prioritized on their relative importance. One such order is provided (roughly) by Gelman and Hill (2006).
Model validity
Linear relationship and additive effects of predictor variables
Errors are
independently distributed
identical (homoscedastic)
follow a normal distribution
Many datasets will violate multiple of these assumptions simultaneously, so addressing issues is often best resolved by understanding why this is happening.
A central (and often overlooked assumption) of linear models is that the relationship between the predictors (remember, not the actual data!) and the outcome variable is linear and additive. When the relationship is not linear, the resulting residuals are often not appropriately distributed as well.
To be clear, the linear model only focuses on the linear and additive relationship between the predictors and the outcome variable (this will become more important/obvious later in this section!). The model doesn’t know what the variables are. For that reason, we can add predictor variables to a model that are squares or cubes of predictors, or we can transform the outcome variable. We just need the \(Y = X\beta\) relationship to be additive and linear.
When non-linearity occurs, several options to address it exist. The best approach may depend on why the relationship is non-linear, as relationships among variables may be non-linear for a number of reasons. For example, the outcome may not actually be continuous (e.g. counts. proportions), and thus true linear relationships are not possible; outcomes may also be related to the predictors in different ways (e.g., logarithmic).
One way to address non-linearity is to transform the data (typically focusing on the dependent variable) so that the resulting variable meets the linear model assumptions (and thus uses the strengths of linear models that we noted above). As shown above, our rank-based approaches are using a similar method (not technically the same, but it works for larger sample sizes). This approach was often used in the past (e.g., arc-sin transforms of proportion data, D. Warton and Hui 2010) and supported by various approaches. For example, the Box-Cox transformation helped researchers find the best way to transform data to reduce issues with the distribution of residuals; this method also tended to impact linearity and differences in variances.
Two things should be noted regarding this approach and transformations in general:
The Box-Cox approach requires a model - it still focused on transforming data so that residuals meet assumptions. Datashould not be transformed before model assumptions are analyzed to ensure it is necessary.
Similarly (and already noted!). non-normal data may lead to models with normal residuals. However, normal data (or data transformed to be more normal) does typically lead to normal residuals, so if residuals are not normal transformations may help.
The transformed variable is now linear in respect to the predictors. This highlights the actual assumption of the model. Similarly, higher-terms (squares, cubes, etc) may be added to a linear model. The model does not care what the data represent - it only focuses on if linear relationships exist among them.
Transformations can make model interpretation and prediction difficult.
If the decision is made to transform the data, several approaches exist. Some are driven by the distribution of the data, and all depend on it. For example, log and related root transformations are useful for right-skewed data, but some can only be carried out for non-negative (e.g., square root) or strictly positive (e.g., log) values. To address this for log transformations, a small value is often added to 0 measurements.
Let’s use data (not residuals) to show how different distributions of data connect to q-q plots and consider possible fixes (always fit a model first for real analysis!). For example, we can return to our right-skewed blue jay data from the summarizing data chapter(target=“_blank”)(idea from (Hunt, n.d.)) .
set.seed(19)
blue_jays <- round(rbeta(10000,2,8)*100+runif(10000,60,80),3)
ggplot(data.frame(blue_jays),
aes(x=blue_jays)) +
geom_histogram( fill="blue", color="black") +
labs(title="Weight of Westchester blue jays",
x= "Weight (g)",
y= "Frequency")
Note the right-skewed data shows a convex curve on the Q-Q plot.
qqnorm(blue_jays)
qqline(blue_jays)
To help you understand Q-Q plots, remember they are comparing the relative spread of data (quantiles) from a target and theoretical distribution. For right-skewed data, points are shifted right (both smallest and largest observations are larger than you would expect from a normal distribution).
Conversely, we could also return to our left-skewed cardinal data
set.seed(19)
cardinals <- round(rbeta(10000,10,2)*50+runif(10000,5,10),3)
ggplot(data.frame(cardinals),
aes(x=cardinals)) +
geom_histogram( fill="red", color="black") +
labs(title="Weight of Westchester cardinals",
x= "Weight (g)",
y= "Frequency")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
and note a concave shape is seen in the Q-Q plot,as all points are shifted left.
qqnorm(cardinals)
qqline(cardinals)
For this type of data, power transformations (raising the variable to the 2, 3, or higher power) may be useful.
Meanwhile, our uniformly-distributed distributed robin data
set.seed(19)
rochester <- round(c(runif(1000,75,85)),3)
ggplot(data.frame(rochester),
aes(x=rochester)) +
geom_histogram( fill="pink", color="black") +
labs(title="Weight of Rochester robins",
x= "Weight (g)",
y= "Frequency")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
shows as a s-shape on the Q-Q plot
qqnorm(rochester)
qqline(rochester)
because it is under-dispersed (has no tails, and thus the lower points are larger than expected and the higher points are smaller than expected). Alternatively, data may be over-dispersed, like this (fake) finch data where lower points are lower than expected and higher points are higher.
set.seed(19)
library(VGAM)
finch <- rlaplace(1000,location=50, scale=4)
ggplot(data.frame(finch),
aes(x=finch)) +
geom_histogram( fill="cyan", color="black") +
labs(title="Weight of finches",
x= "Weight (g)",
y= "Frequency")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
qqnorm(finch)
qqline(finch)
Over- and under-dispersed data may mean there’s a missing factor in your analysis. For example, our bi-modal woodpecker data
set.seed(19)
woodpeckers <- round(c(rnorm(100,60,4),rnorm(100,80,4)),3)
ggplot(data.frame(woodpeckers),
aes(x=woodpeckers)) +
geom_histogram( fill="orange", color="black") +
labs(title="Weight of Westchester woodpeckers",
x= "Weight (g)",
y= "Frequency")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
is under-dispersed due the shape of the distribution, but there also signs of other issues.
qqnorm(woodpeckers)
qqline(woodpeckers)
Under-dispersion could also happen if data are bounded (e.g., by practicality or due to a measurement issue). Over-dispersion can similarly occur if the model does not account for variability (e.g., missing factors, non-linear relationship) and/or outliers, which might be related to the underlying form of the data (Payne et al. 2017) (more to come on this). In this way, over- and under-dispersion relate to both the linear relationship.
The linear relationship may be inappropriate because our data doesn’t fit it! For example, if we are modeling proportions, estimates less than 0 or below 1 do not make sense, but a linear model doesn’t account for that. The same issues arise for binary outcomes (outcome must be 0 or 1) and count data (where outcomes can’t be non-integer). Along with impacts on the assumption of linearity, data like this may also lead to unequal variances and lack of normality for residuals. While we may be able to transform the response variable to make relationships linear, another option is to translate the link function.
Although we haven’t fully explained it before, a linear model contains three components. There is always a random component that focuses on the distribution of the data. This typically requires a few parameters (or at least one!). For example, for normally-distributed data, the two parameters that describe the distribution are the mean and the variance. Note these parameters are independent (we covered all this when we introduced normality{target=“_blank”). There is also a systematic component, where a number of covariates and data points produce an estimate. Finally, there is the link function, which connects that estimate to the data distribution (this maintains the linearity of the model). Notice this means the link connects to the distribution of the data, not the data itself (which is what transforming the data is doing).
So far we have focused on the mean of the data, and the estimate is the mean, so the link has been implied. We can generalize this setup to include other distributions in the exponential family (and now others (“Applied Regression Analysis and Generalized Linear Models” 2022, ch. 15)). While we won’t develop all the math, exponential family distributions all have a canonical parameter (the mean for Gaussian, or normal, data) and a dispersion parameter (the variance for normal data). Different distributions have different relationships between these two parameters. For normal data, there is no relationship, so the variance is constant. This is why we have assumed homoscedasticity of residuals for linear models. This is not true for other families; for example, the variance in binomial data depends on the p parameter (as we noted in the chapter introducing binomial data{target=“_blank”)), and for the Poisson distribution (used for count data) the variance is the mean. Different families also use different link functions, which can impact how we interpret coefficient values and associated impacts. In general, impacts and outcomes will need to be transformed (based on the link) to the actual scale.
We can specify other data types using generalized linear models (which are different than general linear models, which is another term used to describe our “regular” linear models). These models are also known as glm or glim (we’ll use the glm jargon here). These models let us consider data that do not fit a normal approximation. Bolker (2008) has an excellent review of when various distributions should be used.
| Distribution | Typical usage, notes, and examples |
|---|---|
| Continuous | |
| Gaussian/normal | continuous variables, or variables connected to large sample sizes and not extreme values (remember the z-test approximation for \(\chi^2\) tests?) |
| Gamma | Amount of time until a certain number of events takes place (how long do you wait for a parking spot?); variance changes with mean; related to negative binomial |
| Exponential | amount of time for a single event to occur; related to geometric |
| Beta | distribution of success parameter in binomial trials with a-1 success and b-1 failures; bounded by 0 -1 (not inclusive) |
| Discrete | |
| Bernoulli/binomial | data that can be divided into successes and failures |
| Beta | Proportions that do not include 0 or 1 |
| Poisson | Fits count data where all events are independent (random); Count data where low counts and 0’s are common |
| Negative binomial | Number of failures before a set number of successes occur; Poisson-ish where variance is larger than the mean (good for overdispersion!) or where the mean (and thus variance) varies |
| Geometric | Number of trials until you get a single failure |
| Beta binomial | Binomial-ish but the success parameter (and thus variance) varies across trials |
GLM make different assumptions (though everyone does not agree on what they are!). While we still need independence of the residuals (or some extension, see below), they no longer need to be normally distributed (A. Zuur, Ieno, and Smith 2007, 86–87; A. Zuur et al. 2009, sec. 9.8.3) (though normality may hold at large sample sizes for Pearson (Agresti 2007, 87) or deviance (Montgomery, Peck, and Vining 2012, sec. 13.4.4) residuals, and some authors argue normality is required for testing (Feng, Li, and Sadeghpour 2020)) . As noted above, non-Gaussian families don’t assume a constant variance, so homogeneity assumptions would not be appropriate (Fox 2016). However, we need to ensure the correct family (and thus mean-variance link) is chosen.
In general, plotting the residuals against (A. Zuur et al. 2009)
the predicted values
each explanatory variable in model
each explanatory not retained/used in model
against time and space (if applicable)
will be useful. Any patterns may suggest missing predictors or lack of independence/correlation among measurements. Patterns in the residuals may also indicate the incorrect family has been chosen (A. Zuur et al. 2009, sec. 9.8.4); for families with fixed dispersion parameters, inspection of the estimated dispersion parameter can also be a related diagnostic.
Residuals for a Gaussian/normal glm (what we’ve been doing) are easy to calculate. We simply subtract the estimate for each data point (which is the average for similar observations!) from the observed. However, though these raw or response residuals exist for GLM, they rarely make sense. Variance may increase with mean for some distributions, for example, or outcomes may be binary (so every residual is 0 or 1!).
To account for this, we could divide the residuals by the standard deviation of the observations. This normalizes the residuals and leads to Pearson residuals. Another option is based on something called deviance. Deviance is similar to sums of squares but based on likelihood (which we shortly see is what we use to test for significance in GLM); it can also be considered the generalized form of variance or residual sums of squares. It is calculated as the difference of log-likelihoods between the focal model and the saturated model (or a model that has a coefficient for every observation and thus will have perfect fit). Deviance residuals consider how important each point is to the overall deviance. You can also compare the deviance of a model to a perfect fit (no deviance/saturated) model and an intercept-only model (worst fit, DO) and generate a pseudo-R2 measure using the formula
\[ 1-\frac{D}{D_O} \]
This is known as McFadden’s pseudo-R2 value.
Once developed, GLMs can use most of the tools we developed for basic linear models. However, a few key differences exist. Estimation of parameters and significance for these models rely on maximum-likelihood methods. It should be noted that solving for coefficients also requires iterative approaches due to non-linearity; common methods (not described here) include the Newton-Raphson or Fisher scoring methods. The key point here for practice is any noted issues with optimization are due to this process. The ratio of likelihood values from various nested tests follow a \(\chi^2\) distribution and can thus provide a p-value to compare possible models as noted in introducing G tests{target=“_blank”)).This may be represented as a Z test for single parameters (since a \(\chi^2\) distribution is a Z distribution squared - sometimes labelled a Wald test/statistic). F test reappear for distributions where we estimate the dispersion parameter (like the normal), but otherwise we can use the dispersion parameter to check that the right approach was taken.
There are numerous types of GLMs. Here we outline some of the more common ones. Each of these deserves a much fuller treatment - remember, the intent here is to get you to a place to where you can begin working with them. Harrison et al. (2018) (which also covers mixed models (coming up!)) provides an excellent review of these models, and fuller treatments are provided by A. Zuur et al. (2009) (chapters 8-11) and Fox (2016) (chapters 14 - 15). For each I provide the most relevant math, an example (or two), and some immediate extensions.
Logistic regression (D. Warton and Hui 2010) focuses on modelling yes/no or success/failure outcomes (Bernoulli or binomial data) using the logit link (thus logistic regression). The easy-to-understand issue here is that we want to focus on probability (\(\pi\)), but that can only vary from 0 to 1. Linear models are thus inappropriate, but the logit link can be used where the systematic component of the model predicts
\[ logit(\pi) = \ln(\frac{\pi}{1-\pi}), \text{where } \pi \text{ is the probability for a given set of outcomes} \]
This logit link is useful as it effectively connects the systematic component and real data (the purpose of the link). This link also allows coefficients to be interpreted as odds and odds ratios (introduced in the comparing proportions among groups- we’ll continue to see connections this chapter for logistic regression). Remember, the odds can be written as
\[ \frac{\pi}{1-\pi} \]
For a given coefficient \(\beta\), the logit link also means we can say a 1 unit increase in that single variable leads to a change in the odds ratio of \(e^{\beta}\). This also means no impact is shown as \(\beta = 0\), with positive and negative numbers implying increases or decreases in odds, which is convenient. Other link options (not discussed here) include the probit and complementary log-log approaches. For a final piece of math, note the dispersion parameter for binomial data is not estimated - it is assumed to be 1.
For example, Blais, Johnson, and Koprowski (2023) investigated how various factors impacted behavior in garter snakes (Thamnophis cyrtopsis).
.jpg)
One of their analyses considered how snake movement (labelled 0 for no movement, 1 for movement) differed among age classes (adults and juveniles).
behavior <- read.csv("data/database_Thamnophis_crytopsis_study - Copy of THCY_MH_BEH.csv",
stringsAsFactors = T)
behavior_table <- table(behavior$movement, behavior$ageclass)
behavior_table
adult juvenile
0 29 29
1 14 16You may recognize this analyses could be carried out using \(\chi^2\) tests.
chisq.test(behavior_table)
Pearson's Chi-squared test with Yates' continuity correction
data: behavior_table
X-squared = 0.0051228, df = 1, p-value = 0.9429which indicate no impact of age class on movement (p=.9429). You (hopefully) also remember we mentioned the Yate’s correction (noted in the output) was due to using a continuous distribution to fit counts. If we remove it,
chisq.test(behavior_table, correct = F)
Pearson's Chi-squared test
data: behavior_table
X-squared = 0.087924, df = 1, p-value = 0.7668we find an answer that we can replicate using glm (since the outcome is a 0 or 1 here!). We do this using the glm function in base R, which is very similar to the lm function; the only difference is we now note the distribution we are exploring using the “family” argument. Note Bernoulli data is considered a subset of binomial data and thus classified as “binomial”.
as_glm <- glm(movement~ageclass, behavior, family = "binomial")We can then use our normal approach (including summary, Anova, and glht functions) to consider the outcomes. Checking assumptions is a little trickier than for linear models (as noted above). One major check is typically ensuring the dispersion parameter is correct. We can typically check the parameter by dividing the residual deviance by its associated degrees of freedom; this information is provided by the summary function (here, 112.84/86);
summary(as_glm)
Call:
glm(formula = movement ~ ageclass, family = "binomial", data = behavior)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.7282 0.3254 -2.238 0.0252 *
ageclassjuvenile 0.1335 0.4504 0.296 0.7669
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 112.93 on 87 degrees of freedom
Residual deviance: 112.84 on 86 degrees of freedom
(15 observations deleted due to missingness)
AIC: 116.84
Number of Fisher Scoring iterations: 4For Bernoulli data like this, however, over-dispersion is not an issue (Molenberghs et al. 2012). For binomial and Poisson data, this should be approximately one. While we can still plot the model like we did for lms
plot(as_glm)
the findings may not be as useful. We only have 4 predicted values (why?) and thus only 4 residual values (why?), so we no longer expect the residuals to be normally distributed. Checking for patterns in these residuals is also complicated.
Another option (which also works for linear models) is offered by the performance package in R (another, related package not explored here, is DHARMa).
library(performance)
check_model(as_glm)
Here we see no real issues with the model. We can also check for over-dispersion
check_overdispersion(as_glm)# Overdispersion test
dispersion ratio = 1.018
p-value = 0.856No overdispersion detected.but remember for Bernoulli data (n=1) this isn’t an issue anyway.
We can note the impact of given factors using Anova
Anova(as_glm, type="III")Analysis of Deviance Table (Type III tests)
Response: movement
LR Chisq Df Pr(>Chisq)
ageclass 0.087974 1 0.7668which now uses likelihood-based tests by default. Similarly, the summary output notes details related to fitting the model (Number of Fisher scoring iterations). To interpret the coefficient from the glm, we need to consider the link function (here, the logit). If we backtransform the logit link function, we can actually say that juveniles had a \(e^\beta\) greater odds of movement than adults, which is equal to
exp(.1335)[1] 1.142821and very close to 1 (thus the high p-value given our small sample size). Note this means a coefficient of 0 means no impact, since
exp(0)[1] 1While this example simply replicates the \(\chi^2\) tests we learned earlier (deviations may arise when counts are low or high for any group, and thus \(\hat{p}\) approaches extreme values and breaks the assumptions!), the power of these models comes from the ability to consider multiple variables (which we couldn’t do before). These models are sometimes called log-linear models.
To consider this approach and extend glm to binomial data, consider work by Bueno, Machado, and Leite (2020). The authors carried out a series of experiments to determine how various factors impacted colonization of a host algae, Sargassum filipendula, by a small amphipod, Cymadusa filosa.
For one experiment, they stocked portions of aquariums with pieces of algal. Half the algae pieces had adult amphipods present(4 total; the species builds tubes to live on algae) while the other half had no adult residents. Juveniles (n=16 per patch) were initially placed on patches of natural substrate, artificial substrate, or bare areas(no substrate) in the aquarium. After 24 hours, the number of juveniles that colonized the focal algae were counted. Since the number of introduced juveniles was known, the number that did not colonize was also accounted for (raw data estimated from supplemental data given variation in sample size).
juvenile <- read.csv ("data/Supplemental_Data_Bueno_2020-S1-Presence of adults_lab.csv",
stringsAsFactors = T)since the adults column indicates the presence or absence of adults, let’s update it for clarity
library(plyr)
juvenile$adults_updated <- factor(revalue(as.character(juvenile$adults), c("0"= "absent", "4" ="present")))This is an example of a factorial-design with the following hypotheses:
The outcome, however, is a proportion. Note we can fit a linear model to the data:
juvenile_fit_lm <- lm(Colonized ~ adults_updated*source.habitat,
juvenile)
plot(juvenile_fit_lm)
The assumptions don’t even look that bad, but note the residuals are over-dispersed. We could also assess the model using
check_model(juvenile_fit_lm)
which again shows some issues.
We also know the model will make predictions that aren’t logical (outcomes outside of the 0-1 range). For this reason we will use a generalized linear model (or a log linear model) to analyze the data. We can fit a glm using the binomial family.
juvenile_fit_glm <- glm(cbind(Colonized, Not_colonized) ~ adults*source.habitat,
juvenile, family = "binomial")The summary again allows us to estimate the dispersion parameter by dividing the residual deviance by its associated degrees of freedom, and we can now consider the impact of multiple variables (including interactions, which are significant here!).
summary(juvenile_fit_glm)
Call:
glm(formula = cbind(Colonized, Not_colonized) ~ adults * source.habitat,
family = "binomial", data = juvenile)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.84583 0.31063 5.942 2.81e-09 ***
adults -0.28417 0.09507 -2.989 0.002799 **
source.habitatnatural -3.88271 0.43669 -8.891 < 2e-16 ***
source.habitatnone 1.02207 0.55484 1.842 0.065459 .
adults:source.habitatnatural 0.46383 0.13819 3.356 0.000789 ***
adults:source.habitatnone -0.38722 0.15665 -2.472 0.013441 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 347.435 on 35 degrees of freedom
Residual deviance: 97.393 on 30 degrees of freedom
AIC: 186.64
Number of Fisher Scoring iterations: 5The residual deviance/residual df (97.393/30) is ~ 3 when it is supposed to be 1. We’ll come back to this (and should), but just to show what we can do with these models, note we can use our Anova function to investigate outcomes.
Anova(juvenile_fit_glm, type="III")Analysis of Deviance Table (Type III tests)
Response: cbind(Colonized, Not_colonized)
LR Chisq Df Pr(>Chisq)
adults 9.712 1 0.001831 **
source.habitat 195.748 2 < 2.2e-16 ***
adults:source.habitat 34.738 2 2.862e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We see that there is a significant interaction between the presence of adults and source habitat on the proportion of amphipods that move to colonize algae. Given this, we can follow our typical approach to significant interactions - break the data into subsets.
However, we may want to reconsider this model given the over-dispersion issue. What should we do? One option is to use an approach that estimates the dispersion parameter. Though not fully developed here (but a little more below), this approach uses a “quasi-binomial” family.
juvenile_fit_glm_quasi <- glm(cbind(Colonized, Not_colonized) ~ adults*source.habitat,
juvenile, family = "quasibinomial")We can investigate this model using our established protocols
summary(juvenile_fit_glm_quasi)
Call:
glm(formula = cbind(Colonized, Not_colonized) ~ adults * source.habitat,
family = "quasibinomial", data = juvenile)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.8458 0.5562 3.319 0.00238 **
adults -0.2842 0.1702 -1.669 0.10544
source.habitatnatural -3.8827 0.7819 -4.966 2.56e-05 ***
source.habitatnone 1.0221 0.9934 1.029 0.31177
adults:source.habitatnatural 0.4638 0.2474 1.875 0.07060 .
adults:source.habitatnone -0.3872 0.2805 -1.381 0.17760
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasibinomial family taken to be 3.205612)
Null deviance: 347.435 on 35 degrees of freedom
Residual deviance: 97.393 on 30 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 5Anova(juvenile_fit_glm_quasi, type="III")Analysis of Deviance Table (Type III tests)
Response: cbind(Colonized, Not_colonized)
LR Chisq Df Pr(>Chisq)
adults 3.030 1 0.081751 .
source.habitat 61.064 2 5.496e-14 ***
adults:source.habitat 10.837 2 0.004435 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Since we estimated our dispersion parameter, we can also use F tests here
Anova(juvenile_fit_glm_quasi, type="III", test = "F")Analysis of Deviance Table (Type III tests)
Response: cbind(Colonized, Not_colonized)
Error estimate based on Pearson residuals
Sum Sq Df F values Pr(>F)
adults 9.712 1 3.0298 0.091998 .
source.habitat 195.748 2 30.5326 5.842e-08 ***
adults:source.habitat 34.738 2 5.4184 0.009796 **
Residuals 96.167 30
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Poisson regression focuses on count-based data. For simplicity (and to demonstrate the connections) here we note that we can almost always assess binomial data using the Poisson approach if we only focus on the number of successes. As you might imagine, this means we are ignoring data, which isn’t a great idea.
For example, we could consider the shrimp data like this
juvenile_fit_glm_poisson <- glm(Colonized ~ adults*source.habitat,
juvenile, family = "poisson")We could use our familiar protocol to assess the model.
plot(juvenile_fit_glm_poisson)
check_model(juvenile_fit_glm_poisson)
check_overdispersion(juvenile_fit_glm_poisson)# Overdispersion test
dispersion ratio = 0.840
Pearson's Chi-Squared = 25.213
p-value = 0.715No overdispersion detected.summary(juvenile_fit_glm_poisson)
Call:
glm(formula = Colonized ~ adults * source.habitat, family = "poisson",
data = juvenile)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.53897 0.11471 22.134 < 2e-16 ***
adults -0.04690 0.04260 -1.101 0.271
source.habitatnatural -1.84583 0.31063 -5.942 2.81e-09 ***
source.habitatnone 0.14660 0.15659 0.936 0.349
adults:source.habitatnatural 0.16178 0.10155 1.593 0.111
adults:source.habitatnone -0.04885 0.05972 -0.818 0.413
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 133.282 on 35 degrees of freedom
Residual deviance: 27.223 on 30 degrees of freedom
AIC: 174.77
Number of Fisher Scoring iterations: 4Noting a lack of issues, we can try
Anova(juvenile_fit_glm_poisson, type= "III")Analysis of Deviance Table (Type III tests)
Response: Colonized
LR Chisq Df Pr(>Chisq)
adults 1.218 1 0.2698
source.habitat 72.622 2 <2e-16 ***
adults:source.habitat 4.499 2 0.1054
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Similar to a binomial model, we can exponentiate parameter values to aid interpretation, but this time it means how increasing the measured value by 1 leads to a multiplicative change in the outcome. Again, a coefficient value of 0 means no impact.
So which is right? Here I might argue that “losing” data in the Poisson approach justifies the binomial (or quasi-binomial) approach.
For an example of continuous data using glms, consider the work by Li et al. (2024) on how humidity impacted responses to heat stress in parasitoids (Venturia canescens) and their hosts (Plodia interpunctella, Indian meal moths). Eggs were placed in incubators kept at ambient or increased humidity levels. Half of each group was parasitized during their 4th instar. Parasitized eggs were then exposed to no heat stress or stress for 6 or 72 hours in their 4th or 5th instar phase. Multiple traits in host and parasitoid were then measured.
One measured factor was wasp development time.
wasp_no_humidity <- read.csv ("data/Nonhumidity_hostwithparasitoids.csv",
stringsAsFactors = T)
wasp_humidity <- read.csv ("data/humidity_hostwithparasitoids.csv",
stringsAsFactors = T)
wasp <- rbind(wasp_humidity,wasp_no_humidity)Warning in `[<-.factor`(`*tmp*`, ri, value = c(NA, NA, 20230426L, 20230425L, :
invalid factor level, NA generatedwasp$Parasited_date <-as.Date(wasp$Parasited_date, "%Y%m%d")
wasp$Emerge_date <-as.Date(as.character(wasp$Emerge_date), "%Y%m%d")
wasp$develop_time <- as.numeric(wasp$Emerge_date - wasp$Parasited_date)The experimental data can be visualized as:
library(plyr)
library(Hmisc)
wasp$Duration_new <- revalue(wasp$Duration, c("C"="0", "L"="72", "S"="6"))
wasp$Duration_new <- factor(wasp$Duration_new, levels = c("0", "6", "72"))
wasp$Humidity_new <- revalue(wasp$Humidity, c("H"="high", "NH"="ambient"))
wasp$Humidity_new <- factor(wasp$Humidity_new, levels = c("ambient", "high"))
ggplot(wasp, aes(x=Duration_new, y=develop_time, color=Humidity_new))+
geom_jitter()+
labs(x="Heat Stress Duration (hours)", y= "Development time (days)",
color="Humidity")Warning: Removed 404 rows containing missing values or values outside the scale range
(`geom_point()`).
This could be analyzed using a linear model
time_model <- lm(develop_time~Humidity*Duration, wasp)but the residuals are clearly overdispersed and may be increasing with the mean.
plot(time_model)
This makes sense as the data is bounded. Negative development times are not possible, and we might assume that as the average development date increases, there is more room for variability (so variance is not constant). A gamma distribution can help these issues.
time_model_gamma <- glm(develop_time~Humidity*Duration, wasp, family = Gamma)
check_model(time_model_gamma)
and AIC suggests it is very slightly better (remember, linear models are very robust!).
AICc(time_model, time_model_gamma) df AICc
time_model 7 970.7620
time_model_gamma 7 970.0118We can then analyze the impact of various factors on development time using our standard approach.
Anova(time_model_gamma, type="III")Analysis of Deviance Table (Type III tests)
Response: develop_time
LR Chisq Df Pr(>Chisq)
Humidity 2.3234 1 0.1274
Duration 18.9323 2 7.743e-05 ***
Humidity:Duration 3.5064 2 0.1732
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Since there is not evidence for a significant interaction, we could remove it
time_model_gamma <- glm(develop_time~Humidity+Duration, wasp, family = Gamma)
Anova(time_model_gamma, type="III")Analysis of Deviance Table (Type III tests)
Response: develop_time
LR Chisq Df Pr(>Chisq)
Humidity 0.0908 1 0.7632
Duration 24.1852 2 5.601e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Results indicate heat stress, but not humidity, impact wasp development time. Since there are three levels, we could use post hoc tests to determine which treatments differ from others.
summary(glht(time_model_gamma, linfct = mcp("Duration"= "Tukey")))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: glm(formula = develop_time ~ Humidity + Duration, family = Gamma,
data = wasp)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
L - C == 0 -0.007440 0.001938 -3.839 0.000357 ***
S - C == 0 -0.001461 0.001912 -0.764 0.721799
S - L == 0 0.005979 0.001373 4.355 < 1e-04 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)This suggests the longest duration is signifcantly different than the others and has a lower value. Similary, let’s consider the coefficients.
summary(time_model_gamma)
Call:
glm(formula = develop_time ~ Humidity + Duration, family = Gamma,
data = wasp)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0505499 0.0018403 27.468 < 2e-16 ***
HumidityNH 0.0003849 0.0012771 0.301 0.763501
DurationL -0.0074400 0.0019381 -3.839 0.000174 ***
DurationS -0.0014607 0.0019124 -0.764 0.446030
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.03109605)
Null deviance: 6.2813 on 175 degrees of freedom
Residual deviance: 5.5285 on 172 degrees of freedom
(404 observations deleted due to missingness)
AIC: 968.81
Number of Fisher Scoring iterations: 4Coefficients indicate the highest predicted value for control conditions (remember, that’s the intercept condition), followed short and long stress events. This would suggest that heat stress lowers development time…but remember the graphed data (here with humidity impacts removed).
ggplot(wasp, aes(x=Duration_new, y=develop_time))+
geom_jitter()+
stat_summary(fun.data="mean_sdl", geom="crossbar", width=0.5)+
labs(x="Heat Stress Duration (hours)", y= "Development time (days)")Warning: Removed 404 rows containing non-finite outside the scale range
(`stat_summary()`).Warning: Removed 404 rows containing missing values or values outside the scale range
(`geom_point()`).
Here we see that the highest heat stress level is associated with a longer wasp development time. Why do the results seem contradictory?
Again, this is due to the link function. For gamma models, the link is the inverse, so to determine the impact of a given coefficient on the raw data we would combine the intercept and coefficients in this way.
1/(.05-.007)[1] 23.255811/(.05-.001)[1] 20.408161/(.05)[1] 20Which indicates a longer development time (23 days vs 20) for wasps exposed to more heat stress.
Beta regression (Geissinger et al. 2022) is similar to logistic/binomial regression, but it focuses on true proportion data. In other words, whereas binomial regression focuses on modeling the number of successes in a given number of trials, beta regression focuses on the distribution of the probability parameter itself. If we are exploring this, it means we know the value can’t really be 0 or 1 (or we would have no variation), so responses may not be 0 or 1 (they must be in this range). Since “accidents” do happen and may lead to zeroes or ones in the data, there are options for adjusting the data to account for these issues (p. 11, Corlatti 2021). This is useful if the outcome is between 0 and 1 (e.g. % cover, % nitrogen). Packages that can do beta analysis in R include betareg and glmmTMB.
Similar to poisson, but an additional parameter allows the mean to not be equal to the variance. Packages to build these models include MASS(glm.nb function).
All the options above still rely on a linear-link property: the systematic component of the model has a linear relationship to a predictor (in glm, we just transformed this back to the raw data as needed). However, we may instead actually want a non-linear relationship between our predictors and outcome. In fact, we may actually want to compare various linear and non-linear fits (e.g., Dornberger et al. 2016)or determine parameter values for a specified relationship.
For example, many fisheries use length of organisms as baseline data because it’s easy to collect, but we may want to use length data to estimate mass. Determining this length-weight relationship may be useful to fisheries management.
We know from basic principles the relationship between length and mass is probably not linear, because length is only capturing growth in one dimension, whereas mass is likely more closely related to volume. Given this, we might decide to fit the equation
\[ Weight =a*Length^B \] to determine this relationship.
A group (Gosnell et al. 2023) used data collected on whelks from multiple locations to consider the parameter values for this model (only 2 locations used here).
.jpg)
whelk <- read.csv("https://raw.githubusercontent.com/jsgosnell/CUNY-BioStats/master/datasets/whelk.csv")
ggplot(whelk, aes(x=Shell.Length, y = Mass, color=Location)) +
geom_point()+
labs(y= "Mass (g)", x= "Shell length (mm)", color= "Population")Warning: Removed 61 rows containing missing values or values outside the scale range
(`geom_point()`).
whelk_lm <- lm(Mass ~ Shell.Length, whelk, na.action = na.omit)
whelk_nls <- nls(Mass ~ b0 * Shell.Length^b1, whelk,
start = list(b0 = 1, b1=3), na.action = na.omit)
AICc(whelk_lm, whelk_nls) df AICc
whelk_lm 3 3947.579
whelk_nls 3 3800.106The non-linear fit is obviously better. Parameters can be found using
summary(whelk_nls)
Formula: Mass ~ b0 * Shell.Length^b1
Parameters:
Estimate Std. Error t value Pr(>|t|)
b0 2.134e-04 5.737e-05 3.719 0.000228 ***
b1 2.855e+00 5.576e-02 51.195 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 24.24 on 410 degrees of freedom
Number of iterations to convergence: 5
Achieved convergence tolerance: 7.614e-07
(61 observations deleted due to missingness)Comparing parameters among sites is slightly more difficult, as dummy variables need to be included for each site; alternatlvely, data can be transformed and a linear model considered.
Notice this approach is mechanistic - we set the relationship and found the best parameters. There are other ways to build non-linear relationships, however, and some allow algorithms to find the best relationship. While these may not be hypothesis-driven, they may prove very useful for predictions. One example of this approach is known as generalized additive modelings (or gams). These models can fit multiple functions to the data, including piece-wise functions. In this example we consider how a splines, or smooth, polynomial curves that are fit to various portions of the data and then “knotted” together, fit the data.
whelk_gam_spline <- gam(Mass ~ s(Shell.Length), data = whelk)
summary(whelk_gam_spline)
Family: gaussian
Link function: identity
Formula:
Mass ~ s(Shell.Length)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 147.321 1.198 123 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(Shell.Length) 3.98 4.993 807.4 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.908 Deviance explained = 90.8%
GCV = 598.47 Scale est. = 591.24 n = 412As noted above, this is just meant to introduce these approaches. We can compare these (and others, such as loess) using
AICc(whelk_gam_spline, whelk_nls, whelk_lm) df AICc
whelk_gam_spline 5.979836 3805.835
whelk_nls 3.000000 3800.106
whelk_lm 3.000000 3947.579A major issue for linear models is when predictors are co-linear. Mathematically speaking, perfect collinearity occurs when any column of the design (X) matrix can be derived by combining other columns. Perfect collinearity will lead to a message noting singularity issues, which is R’s way of telling you the matrix isn’t invertible (which it has to be to solve the equation).
Even partial collinearity will lead to an increase in Type II errors (A. F. Zuur, Ieno, and Elphick 2010). To put it simply, partitioning variance among related predictors is hard. For this reason, a few approaches may be used.
The first step is identifying issues. From the outset, relationships among predictor variables can be assessed using the pairs function in R. If two variables are highly correlated (r2 > .8 is a general limit), only one should be included in the model. Similarly, variance inflation factors (vif) can be assessed for the final and other models to consider this issues (all this is covered in the previous chapter that introduces multiple regression. If this occurs, some variables may need to be removed from consideration.
Covered in later chapters, but if predictors (or multiple outcomes) are correlated, methods can be used to derive statistically independent predictors (though they may be hard to interpret).
When outcome variables are linked together, a few options exist. Note this issue may be obvious from checks of assumptions, but it also may be due to experimental design.
Consider this toy example. In order to investigate impacts of climate stress on oysters, specimens are placed in individual tanks and held at normal summer (calculated using recent data) temperature or at temperatures predicted under 2 IPCC — Intergovernmental Panel on Climate Change- scenarios. Oysters were also exposed to predator cues by dropping in water from tanks with 0, low (.25/m2), or high (2/m2) densities of predators. After two months, changes in oyster shell length (growth) was measured. Twenty-five oysters were measured for each treatment combination.
You hopefully recognize this as a factorial ANOVA experiment that you know how to analyze. If you need a reminder, see the chapter on ANOVA extensions. Experiments like this are odd, however, given the space they require. It is far more common to put lots of organisms in a single container given space and costs. However, this means our measurements are connected; remember blocking and paired tests)?
There are several ways to deal with this. Here we explore each for our oyster example.
First, let’s ignore the lack of independence. This is not an option, but it let’s you see the impact.
growth_lm <- lm(growth~predator_cue*temperature, experiment)
plot(growth_lm)
library(car)
Anova(growth_lm, type="III")Anova Table (Type III tests)
Response: growth
Sum Sq Df F value Pr(>F)
(Intercept) 15.576 1 16.3673 7.265e-05 ***
predator_cue 22.085 2 11.6031 1.635e-05 ***
temperature 70.061 2 36.8090 1.753e-14 ***
predator_cue:temperature 21.376 4 5.6154 0.0002558 ***
Residuals 205.563 216
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We find significant main effects and interactions using this wrong approach.
One way of doing this focuses on the average for each unit
library(Rmisc)
averaged_experiment <- summarySE(experiment, measurevar = "growth",
groupvars = c("predator_cue", "temperature", "container"))
library(rmarkdown)
paged_table(averaged_experiment)and use that for your analysis.
average_analysis <- lm(growth~predator_cue*temperature, averaged_experiment)
Anova(average_analysis, type = "III")Error in Anova.lm(average_analysis, type = "III"): residual df = 0but that leads to an issue! Since we only get 9 average outcomes and our model requires 10 degrees of freedom (consider why), we are left with no “noise” to make the denominator for our F ratio! Even when this doesn’t happen, you have reduced your data to a much smaller number of points and are not getting credit for all your work! This is a good example of why you should analyze simulated data before you run an experiment, but there are other options.
The blocking approach we’ve already covered works might seem appropriate here.
blocking_analysis <- lm(growth~predator_cue*temperature+container, experiment)
Anova(blocking_analysis, type = "III")Error in Anova_III_lm(mod, error, singular.ok = singular.ok, ...): there are aliased coefficients in the modelbut its not? Why? This error means now our model matrix has collinearity issues. We can actually see where
alias(blocking_analysis)Model :
growth ~ predator_cue * temperature + container
Complete :
(Intercept) predator_cuenone
containerf 0 0
containerg -1 1
containerh 0 0
containeri 1 -1
predator_cuenone:temperatureelevated_scenario1 0 0
predator_cuenormal:temperatureelevated_scenario1 0 0
predator_cuenone:temperatureelevated_scenario2 -1 1
predator_cuenormal:temperatureelevated_scenario2 0 0
predator_cuenormal
containerf 0
containerg 0
containerh 1
containeri -1
predator_cuenone:temperatureelevated_scenario1 0
predator_cuenormal:temperatureelevated_scenario1 0
predator_cuenone:temperatureelevated_scenario2 0
predator_cuenormal:temperatureelevated_scenario2 1
temperatureelevated_scenario1
containerf 1
containerg 1
containerh 0
containeri -1
predator_cuenone:temperatureelevated_scenario1 0
predator_cuenormal:temperatureelevated_scenario1 0
predator_cuenone:temperatureelevated_scenario2 1
predator_cuenormal:temperatureelevated_scenario2 0
temperatureelevated_scenario2
containerf 0
containerg 1
containerh 0
containeri 0
predator_cuenone:temperatureelevated_scenario1 0
predator_cuenormal:temperatureelevated_scenario1 0
predator_cuenone:temperatureelevated_scenario2 1
predator_cuenormal:temperatureelevated_scenario2 0
containerb containerc
containerf 0 0
containerg 1 1
containerh -1 0
containeri 0 -1
predator_cuenone:temperatureelevated_scenario1 0 0
predator_cuenormal:temperatureelevated_scenario1 0 0
predator_cuenone:temperatureelevated_scenario2 1 1
predator_cuenormal:temperatureelevated_scenario2 -1 0
containerd containere
containerf -1 -1
containerg -1 0
containerh 0 -1
containeri 1 1
predator_cuenone:temperatureelevated_scenario1 1 0
predator_cuenormal:temperatureelevated_scenario1 0 1
predator_cuenone:temperatureelevated_scenario2 -1 0
predator_cuenormal:temperatureelevated_scenario2 0 -1 though the output is confusing. In general, the issue here is each unit only contributes to one level of other traits..so if we know the average impact of ambient temperatures, for example, and the impacts in two of the treatments that were held at that temperature, we can predict the other. If instead each unit contributed to multiple levels, like in our feather experiment) this isn’t an issue.
Our final option takes a new approach. It considers the units we measured as simply a sample from a larger population. Using that background, we use the information from the units to consider the distribution of sample effects we might see. The impact of unit is then considered a random-effect. For this to work, you probably want 5+ levels of the unit variable. This is because we are using the means to estimate variance (confusing?). For factors with <5 levels, random effects likely offer no benefit (Gomes 2022).
When models contain fixed (what we’ve done before) and random effects, we call them mixed-effects models. Two common packages for carrying out this analysis in R are the nlme and lme4 packages. We will focus on the lme4 package here. Random effects can be entered in the lmer (linear mixed-effects regression) function and specified as (1|Grouping Unit). One nice thing about lme4 is it will handle crossed and random effects on it’s own as long as you don’t repeat unit names. For example, we could note
mixed_analysis <- lmer(growth~predator_cue*temperature+(1|container), experiment)Once built, we need to consider assumptions. The main assumption we add here is that the random effects are normally distributed. This should be checked at each level of grouping. The check_mixed_model function (provided below) offers an automated approach for one level (also known as one-way random effects).
check_mixed_model <- function (model, model_name = NULL) {
#collection of things you might check for mixed model
par(mfrow = c(2,3))
if(length(names(ranef(model))<2)){
qqnorm(ranef(model, drop = T)[[1]], pch = 19, las = 1, cex = 1.4, main= paste(model_name,
"\n Random effects Q-Q plot"))
qqline(ranef(model, drop = T)[[1]])
}
plot(fitted(model),residuals(model), main = paste(model_name,
"\n residuals vs fitted"))
qqnorm(residuals(model), main =paste(model_name,
"\nresiduals q-q plot"))
qqline(residuals(model))
hist(residuals(model), main = paste(model_name,
"\nresidual histogram"))
}check_mixed_model(mixed_analysis)
Here we have only 9 levels of units, so the spread is not perfect. However, we also know each of these is itself an average,and averages should be normally-distributed under the central limit theorem, so we can plow ahead.
We can consider the outcome using our summary command - note the output denotes we have 225 observations and 9 grouping levels.
summary(mixed_analysis)Linear mixed model fit by REML ['lmerMod']
Formula: growth ~ predator_cue * temperature + (1 | container)
Data: experiment
REML criterion at convergence: 631.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.68605 -0.79915 -0.02646 0.70023 2.72393
Random effects:
Groups Name Variance Std.Dev.
container (Intercept) 4.6171 2.1487
Residual 0.9517 0.9755
Number of obs: 225, groups: container, 9
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.78934 2.15758 0.366
predator_cuenone 1.28982 3.05129 0.423
predator_cuenormal 0.92308 3.05129 0.303
temperatureelevated_scenario1 0.03073 3.05129 0.010
temperatureelevated_scenario2 2.06547 3.05129 0.677
predator_cuenone:temperatureelevated_scenario1 0.46256 4.31517 0.107
predator_cuenormal:temperatureelevated_scenario1 -0.24377 4.31517 -0.056
predator_cuenone:temperatureelevated_scenario2 -1.22466 4.31517 -0.284
predator_cuenormal:temperatureelevated_scenario2 -0.63903 4.31517 -0.148
Correlation of Fixed Effects:
(Intr) prdtr_cnn prdtr_cnr tmpr_1 tmpr_2 prdtr_cnn:_1 prdtr_cnr:_1
predatr_cnn -0.707
prdtr_cnrml -0.707 0.500
tmprtrlvt_1 -0.707 0.500 0.500
tmprtrlvt_2 -0.707 0.500 0.500 0.500
prdtr_cnn:_1 0.500 -0.707 -0.354 -0.707 -0.354
prdtr_cnr:_1 0.500 -0.354 -0.707 -0.707 -0.354 0.500
prdtr_cnn:_2 0.500 -0.707 -0.354 -0.354 -0.707 0.500 0.250
prdtr_cnr:_2 0.500 -0.354 -0.707 -0.354 -0.707 0.250 0.500
prdtr_cnn:_2
predatr_cnn
prdtr_cnrml
tmprtrlvt_1
tmprtrlvt_2
prdtr_cnn:_1
prdtr_cnr:_1
prdtr_cnn:_2
prdtr_cnr:_2 0.500 We can also still use Anova to get p-values. However, these are now calculated by default using likelihood-associated \(\chi^2\) tests.
Anova(mixed_analysis, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: growth
Chisq Df Pr(>Chisq)
(Intercept) 0.1338 1 0.7145
predator_cue 0.1898 2 0.9095
temperature 0.6020 2 0.7401
predator_cue:temperature 0.1837 4 0.9960You can also ask for F tests, but note the degrees of freedom associated with these tests is not clear. It’s somewhere between the “average” and “ignore” approach used above.
Anova(mixed_analysis, type = "III", test="F")Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
Response: growth
F Df Df.res Pr(>F)
(Intercept) 0.1338 1 3230137 0.7145
predator_cue 0.0949 2 3230137 0.9095
temperature 0.3010 2 3230137 0.7401
predator_cue:temperature 0.0459 4 3230137 0.9960Note this approach suggests we do not have enough data to reject the null hypothesis. Ignoring the linkages among data led to very different results. This issue (pseudopreplication) has been noted in ecology and other fields (Hurlbert 1984; Heffner, Butler, and Reilly 1996; Lazic 2010).
It should be noted that random effects match up to paired t-tests and blocking approaches for simple experiments.Remember our feather experiment (example and data provided by McDonald (2014))?
two_way_anova_example <- lm(Color_index ~ Feather_type + Bird, data=feather)
Anova(two_way_anova_example, type= "III")Anova Table (Type III tests)
Response: Color_index
Sum Sq Df F value Pr(>F)
(Intercept) 0.24133 1 26.5059 0.000119 ***
Feather_type 0.15043 1 16.5214 0.001017 **
Bird 0.21950 15 1.6072 0.184180
Residuals 0.13657 15
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Note we can get the same results using a paired t-test, but (as of early 2024) we need to use a wide dataset.
library(reshape2)
feather_wide <- dcast(feather, Bird~Feather_type, value.var="Color_index")
paged_table(feather_wide)Now we have wide data…
t.test(feather_wide$Odd, feather_wide$Typical, data=feather, paired=TRUE)
Paired t-test
data: feather_wide$Odd and feather_wide$Typical
t = -4.0647, df = 15, p-value = 0.001017
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
-0.20903152 -0.06521848
sample estimates:
mean difference
-0.137125 two_way_lme_example <- lmer(Color_index ~ Feather_type + (1|Bird), data=feather)
Anova(two_way_lme_example, type= "III", test="F")Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
Response: Color_index
F Df Df.res Pr(>F)
(Intercept) 41.816 1 28.457 4.88e-07 ***
Feather_type 16.521 1 15.000 0.001017 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1All same if we use least-squares approaches (not the default \(\chi^2\) tests that match likelihood-based methods).
Issues with errors not being equal among groups is often tied to other issues, like non-linearity or the nature of the data, which we have addressed above. Another option is to use weighted-least squares regression - this approach specifically helps when residuals are not evenly distributed among groups (or when a funnel/cone appears when you plot the model to check assumptions). For example, we could take the model on below-ground biomass that we developed in the More Anovas chapter. As a reminder, Valdez et al. (2023) wanted to consider the impact of top-down (snail grazing) and bottom- up (nutrient availability) on marsh plant (Spartina alterniflora) growth. To do this, they assigned plots to one of 3 grazer treatments and one of 2 nitrogen treatments. We focused on the impact of these treatments on standing dead mass and noted a funnel-shape when reviewing the residuals.
valdez_2023 <- read.csv("data/Spartina_alterniflora_traits.csv", stringsAsFactors = T)
sdm_model <-lm(Standing.Dead..dry..m2.~Snail.Level * Nitrogen.level, valdez_2023[valdez_2023$Snail.Level != "uncaged",])
plot(sdm_model)
To address the, we can use weighted least squares. This approach assume you built the model and then noted an issue with heteroscedasticity. To use, we can calculate a weight for each residual that is based on its variance. Doing this requires an iterative process similar to those used for estimation for generalized linear models (we’ll get to these) - below makes a value that increases with low variance.
wt_sdm <- 1 / lm(abs(sdm_model$residuals) ~ sdm_model$fitted.values)$fitted.values^2We can then add a new argument to the lm function to use these weights.
sdm_model_wls <-lm(Standing.Dead..dry..m2.~Snail.Level * Nitrogen.level, valdez_2023[valdez_2023$Snail.Level != "uncaged",], weights = wt_sdm)We can then continue on our normal route:
plot(sdm_model_wls) 
Anova(sdm_model_wls, type="III")Anova Table (Type III tests)
Response: Standing.Dead..dry..m2.
Sum Sq Df F value Pr(>F)
(Intercept) 203.228 1 76.442 1.497e-06 ***
Snail.Level 157.447 2 29.611 2.288e-05 ***
Nitrogen.level 130.416 1 49.054 1.427e-05 ***
Snail.Level:Nitrogen.level 143.013 2 26.896 3.681e-05 ***
Residuals 31.903 12
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1If you compare the two models you notice slight differences - these are minimal here due to lack of differences in variance.
summary(sdm_model)
Call:
lm(formula = Standing.Dead..dry..m2. ~ Snail.Level * Nitrogen.level,
data = valdez_2023[valdez_2023$Snail.Level != "uncaged",
])
Residuals:
Min 1Q Median 3Q Max
-3.5800 -0.8350 -0.0317 0.6600 3.7400
Coefficients:
Estimate Std. Error t value
(Intercept) 18.710 1.032 18.124
Snail.Levelremoval -14.150 1.460 -9.692
Snail.Levelsnail addition -7.567 1.460 -5.183
Nitrogen.levelwithout -15.240 1.460 -10.439
Snail.Levelremoval:Nitrogen.levelwithout 16.153 2.065 7.824
Snail.Levelsnail addition:Nitrogen.levelwithout 10.357 2.065 5.016
Pr(>|t|)
(Intercept) 4.39e-10 ***
Snail.Levelremoval 5.02e-07 ***
Snail.Levelsnail addition 0.000228 ***
Nitrogen.levelwithout 2.25e-07 ***
Snail.Levelremoval:Nitrogen.levelwithout 4.72e-06 ***
Snail.Levelsnail addition:Nitrogen.levelwithout 0.000301 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.788 on 12 degrees of freedom
Multiple R-squared: 0.9284, Adjusted R-squared: 0.8986
F-statistic: 31.14 on 5 and 12 DF, p-value: 1.792e-06summary(sdm_model_wls)
Call:
lm(formula = Standing.Dead..dry..m2. ~ Snail.Level * Nitrogen.level,
data = valdez_2023[valdez_2023$Snail.Level != "uncaged",
], weights = wt_sdm)
Weighted Residuals:
Min 1Q Median 3Q Max
-2.31784 -0.89013 -0.03702 0.92266 2.46121
Coefficients:
Estimate Std. Error t value
(Intercept) 18.710 2.140 8.743
Snail.Levelremoval -14.150 2.202 -6.426
Snail.Levelsnail addition -7.567 2.490 -3.039
Nitrogen.levelwithout -15.240 2.176 -7.004
Snail.Levelremoval:Nitrogen.levelwithout 16.153 2.322 6.956
Snail.Levelsnail addition:Nitrogen.levelwithout 10.357 2.620 3.953
Pr(>|t|)
(Intercept) 1.50e-06 ***
Snail.Levelremoval 3.27e-05 ***
Snail.Levelsnail addition 0.01030 *
Nitrogen.levelwithout 1.43e-05 ***
Snail.Levelremoval:Nitrogen.levelwithout 1.53e-05 ***
Snail.Levelsnail addition:Nitrogen.levelwithout 0.00192 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.631 on 12 degrees of freedom
Multiple R-squared: 0.8747, Adjusted R-squared: 0.8225
F-statistic: 16.75 on 5 and 12 DF, p-value: 4.789e-05Why not just always do this? Because weighted least squares implicitly assumes we know the weights. We are actually estimating them, so small datasets may lead to bad estimates and outcomes.
A very common concern regarding linear models is normality. I list it last here because this assumption is one of the least important (and the assumption is based on residuals, not data!). However, non-normal residuals are often (not always) connected to other issues, namely linearity, as noted above.
(to be added)
(to be added)