Should I model the effect of a treatment or the difference? I have some questions about how to answer one of my hypothesis for my master thesis. If there is a difference/effect of treatment in the normoxic O2 uptake and critical O2 values between treatments in an aquatic invertebrate. I have done 3 toxicity experiments, lasting 22 days. And another 3 lasting 5 days with consecutive respirometry experiments lasting 48h. An aquatic invertebrate species was exposed to 2 different test media and one control treatment. The treatments in question is acidic Al-rich water to simulate acid rain, acidic Al-poor water and reference/control water (untreated freshwater).
Acidic Al-rich water is assumed to be very toxic at first (because of polymerization), then the toxicity decreases with time. So a gradient/level is created inside the experiment trays (See attatched figure at the bottom).
For the toxicity experiment lasting 22 days I have analysed the survival data with 'coxph' in the survivalAnalysis package in R.
And for the toxicity experiment with respirometry analysis, I have used the respR package to extract normoxic O2 and critical O2 values.
Now to my "problem", in the literature I have seen some using ANOVA, Tukey post-hoc analysis, Linear Mixed Models or General Additive Models for comparing normoxic O2 uptake and critical O2 values. But, since I have a gradient in my data, where "Level 1" (see figure) is assumed to be most toxic for the acidic Al-rich treatment, I want to test for an effect of polymerizing aluminium on the normoxic O2 uptake and the critical O2 values.
Since 3 animals at a time is exposed to one of the treatments, inside the same chamber, in the same tray (see figure) then moved to the respirometry chambers for 48h. And they can't travel between chambers, trays or treatments I was thinking maybe I could use linear mixed models because of this statement In Harrison et al. 2018 he writes this about LMM “…we might measure several chicks from the same clutch, and several clutches from different females…”. This is somewhat similar to my experiment and animals.
Also, I think my data is nested, Harrison continues to write: “Here, female ID is said to be nested within woodland: each woodland contains multiple females unique to that woodland (that never move among woodlands).”
Female ID is similar to my ‘animal ID’ and woodland is similar to my ‘treatment’ or ‘level’, since animals never travel/move from one level or treatment to another.
My model with 'lme4' in R would look like this:
Predictor: normoxic O2 and O2crit
Fixed effect: treatment
Random effect: animal id nested in Level
lmer(normO2 ~ treatment + (1|Level/animalID) 
lmer(O2crit ~ treatment + (1|Level/animalID) 
I have also added a sample of my data normoxic O2 uptake and critical O2 values in a separate link (bottom).
Am I just overcomplicating and can use the ANOVA or Tukey post-hoc test? Or should I continue with LMM or move to GAM because of non-linearity? Or is all this wrong and I should use another analysis? Feel free to comment, and please tell me if something is unclear and I will refrase and specify more clearly.  Thank you for taking the time to help me with this.
Edit: level is toxicity level, where Level 1 in theory is assumed to be most toxic for acidic Al-rich treatment'. For other acidic Al-poor and control the toxicity is assumed to constant, since there is no polymerization of aluminium.
For the acidic Al-rich treatment, each level receives the same experimental condition, Level 1, most toxic, Level 2 less toxic and Level 3 even less. There are no differences 'inside' the level.
For the treatments acidic Al-poor & reference/control treatment. The toxicity is similar between all levels.
At the moment there is no animal-individual-covariate I want to control for.
The clarify my drawing, it is hatchery trays and troughs with customized chambers inside the troughs.

Sample from datasets

Based on what you've said so far, I think you can radically simplify:
agg_df <- aggregate(cbind(normO2, O2crit) ~ Level + treatment, data = orig_data, FUN = mean)
lm(normO2 ~ treatment*Level, data = agg_df)
lm(O2crit ~ treatment*Level, data = agg_df)


*

*The rationale for aggregating (taking the mean across animals within Level/treatment combinations) is that, as explained by Murtaugh (2007), when the responses are nested (no predictor variables vary within clusters), any inferences at the level of the fixed effects are unaffected by aggregating to the cluster level. (Things get a little more complicated with an unbalanced design and much more complicated with a non-Gaussian response variable ...)

*You should (I think) treat Level as a fixed effect because (1) the levels are not exchangeable (i.e., Level 1 is qualitatively different from levels 2 and 3) and (2) three levels aren't really enough to fit a random effect reliably anyway.

*I'm not sure what you mean by "nonlinearity" - I think you should treat both level and treatment as categorical predictors (i.e. make them factors); the fact that level is ordered doesn't make that much difference to the analysis (you could treat level as an ordered factor, but that also makes the output a bit more confusing).


Murtaugh, Paul A. “Simplicity and Complexity in Ecological Data Analysis.” Ecology 88, no. 1 (2007): 56–62.
