With multiple categorical predictor of erosion
, flowspeed
, and species
, this is similar to what's called "repeated measures analysis of variance." That's a linear model in which you account for the repeated measurements on the same individuals. The problem is that classic repeated measures analysis of variance assumes equal numbers of observations in each treatment/species group, which you don't have.
One way to deal with the different numbers of observations is a linear mixed model. You use erosion
, flowspeed
, and species
as fixed-effect predictors, the angle
as the outcome, and treat the individual plants as providing a random effect.
You set up 1 row of data for each observation, with annotations of angle
, erosion
, flowspeed
, species
, and an ID
indicating the individual plant. It's best with only 6 combinations of erosion
with flowspeed
to code them as categorical predictors, not as numeric. Include height
and weight
also on each line if you want to include those variables in the analysis. Use a set of ID
values from 1 to 153 instead of re-numbering from 1 within each species. Otherwise the software will think that the plants with ID=1
are all the same individual, and a member of all 3 species!
With the lme4
package in R, you could start with something like:
myModel <- lmer(angle ~ erosion*flowspeed*species + (1|ID), data = myData)
That allows for different associations with angle
depending on the combinations of erosion
and flowspeed
and species
. It takes the repeated measurements into account by estimating different intercepts (estimated angle
at reference levels of erosion
and flowspeed
and species
) for the 153 individuals (ID). You don't need to worry about terminology like "nested." The software will correctly interpret the distribution of ID
values among the treatment/species combinations.
That will return a large number of fixed-effect regression coefficients: by my quick count, 2 for erosion
, 1 for flowspeed
, 2 for species
, 2 for erosion:flowspeed
interactions, 4 for erosion:species
interactions, 2 for flowspeed:species
interactions, and 4 for erosion:flowspeed:species
interactions. Do NOT spend much time trying to figure those coefficients out individually. They describe the model in a way that subsequent analysis with other software will make clearer. You will also get an estimate of the variance among the ID
-specific intercept values.
I'd recommend using the Anova()
function in the R car
package to evaluate the overall associations of each of erosion
and flowspeed
and species
, and their sets of interactions, with the angle
outcome. The "Type II" default analysis provided by that function handles different numbers of observations properly, while the standard anova()
or aov()
functions in R don't.
You then can use post-modeling software like that in the emmeans
package to evaluate and compare predicted angle
values among combinations of fixed-effect predictors.
You do have to check whether the assumptions of the linear model are reasonably well met. The main issue with categorical predictors is whether the ranges of residuals (differences between observed and predicted angle
values) are similar over the range of predicted values. If that's not the case, you might have to consider some pre-transformation of the angle
values. A reasonably normal distribution of the residuals is a plus, but not so critical when you have a large number of observations.
The above doesn't incorporate height
and weight
explicitly in the model. It includes them implicitly in the ID
values and the corresponding random intercepts. You could add them as explicit predictors in the model. If you do, think carefully about the form in which to include them. If you just include them as linear terms, you are implicitly assuming that the angle
is linearly and additively associated with each of height
and weight
on top of all of the other effects associated with erosion
and flowspeed
and species
. Is that reasonable?
Finally, there is one limitation to the study design that you need to address in discussing your results. As all plants received the same treatments in the same order, you can't rule out the possibility of some time- or exposure-dependence of the results. That is, the results of later treatment combinations might not just depend on erosion
and flowspeed
and species
but also on the time elapsed or treatments previously experienced.