# Measuring the effect of treatment on variable over time

I am trying to establish whether a treatment had a significant effect on a variable over time. I have raspberry plants treated with $$5$$ concentrations of compost ($$0$$%, $$5$$%, $$10$$%, $$15$$% and $$20$$%) and am interested in the effect of concentration on various plant attributes over time (each attribute tested separately). For example, the effect of treatment (concentration of compost) on plant height over a period of time.

I have measured each value over several weeks, so they have many time stamps. I think running a simple ANOVA will compare the treatments overall based on each treatment's values over time and a repeated measures ANOVA will compare treatments in each time stamp and overall but to the progress.

I'm also considering comparing regressions but not sure what is the right way to go about it. What would be the best way to compare progress as well as overall result? P.S. I'm running my analysis in R.

Edit - I need to consider 2 things:

1. The effect of concentration might not be linear. Meaning the effect might only show from a certain concentration (say, $$10$$% or above) and might have a bell shaped effect (showing optimal concentration).

2. I might need to treat each concentration as a separate group and not as a linear effect.

A very easy way to accomplish this would be through a mixed model, with the treatment as a fixed effect and random slope, and the period of measurement as the random intercept. As an example below with some simulated data in R:

#### Load Library and Set Seed ####
library(lmerTest)
set.seed(123)

#### Create Grid ####
n <- 100 # plants
t <- 3 # times
df <- expand.grid(plant_id = 1:n, time = 1:t)

#### Simulate Fixed Effects ####
df$$concentration <- sample(c(0, 5, 10, 15, 20), n * t, replace = TRUE) df$$fixed_effect <- 2 * df$concentration #### Simulate Random Effects #### random_intercept <- rnorm(n) df$$random_effect <- random_intercept[df$$plant_id] #### Simulate RSE #### df$error <- rnorm(n * t)

#### Simulate Response Variable (Height) ####
df$$height <- df$$fixed_effect + df$$random_effect + df$$error


We can fit a very simple random intercept model where we fit concentration and time as main effects and an interaction between time and concentration, modeling the cluster of each individual plant to capture repeated measures.

#### Fit Model ####
model <- lmer(height
~ concentration
* time
+ (1|plant_id),
data = df)
summary(model)


The summary of the model can be shown below, which shows the fixed and random effects. Plant ID, the individual plants, vary by about 1.3921 inches in height. The concentration of chemicals models the effect we wanted, about a 2 inch increase in height per unit increase of concentration. Time and the interaction between concentration and time have no effect because I didn't model that into the simulation, and so the effect is near zero:

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: height ~ concentration * time + (1 | plant_id)
Data: df

REML criterion at convergence: 1016.5

Scaled residuals:
Min       1Q   Median       3Q      Max
-2.45043 -0.55583  0.02086  0.55367  2.17557

Random effects:
Groups   Name        Variance Std.Dev.
plant_id (Intercept) 1.3921   1.1799
Residual             0.9303   0.9645
Number of obs: 300, groups:  plant_id, 100

Fixed effects:
Estimate Std. Error         df t value Pr(>|t|)
(Intercept)         -0.026806   0.302805 284.436465  -0.089    0.930
concentration        1.974084   0.024410 232.164396  80.872   <2e-16 ***
time                -0.027776   0.128062 221.194446  -0.217    0.828
concentration:time   0.009915   0.011094 231.049136   0.894    0.372
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) cncntr time
concentratn -0.782
time        -0.854  0.783
cncntrtn:tm  0.729 -0.929 -0.846


The summary shows that there may be a relationship between concentration and height, but plotting will give us a much better understanding of if that is actually the case (in terms of magnitude of the effect).

#### Plot Regression Lines ####
library(tidyverse)
df %>%
mutate(plant_id = factor(plant_id)) %>%
ggplot(aes(x=concentration,
y=height))+
geom_point(position = "jitter")+
geom_line(stat = "smooth",
method = "lm",
se=F,
color="gray",
alpha = .4,
aes(group=plant_id))+
geom_smooth(color="steelblue")+
theme_bw()+
labs(x="Concentration",
y="Height",
title = "Concentration x Height")


The gray lines are the random intercepts of each plant (which have multiple measurements due to time) and the overall effect is in blue. We can see that the intercepts vary by about the amount we suggested earlier, whereas the slopes are fairly similar. Grouping them by time gives us a lot more simple relationship where the overall effect doesn't have any substantial change by time, indicating that the effect is constant across time points.

df %>%
mutate(plant_id = factor(plant_id)) %>%
ggplot(aes(x=concentration,
y=height))+
geom_point(position = "jitter")+
geom_line(stat = "smooth",
method = "lm",
se=F,
color="gray",
alpha = .4,
aes(group=plant_id))+
geom_smooth(color="steelblue")+
theme_bw()+
labs(x="Concentration",
y="Height",
title = "Concentration x Height")+
facet_wrap(~time)


This is a simple example and may not model in whether or not changes differ in magnitude across plants. For that, you would need a more complex approach (such as entering random slopes into the model). In any case, a few good references on how to use mixed models can be found below.

### Edit

You asked about nonlinear effects in the comments. There are a couple things you could consider. The most flexible would be a generalized additive mixed model (GAMM), where the effects are treated as nonlinear and do not require any special transformations of the terms (indeed, you often need to scale the predictors in interactions for linear mixed models for stability reasons). Another easier but less flexible option is to include polynomials into the regression equation, but these can sometimes have issues with interpolation at the boundaries of the distribution or can lead to some overfitting. This will depend a lot on the nature of curvilinearity. I provide a useful reference text from PeerJ on modeling GAMMs.

#### References

• Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1). https://doi.org/10.18637/jss.v067.i01
• Bolker, B. M., Brooks, M. E., Clark, C. J., Geange, S. W., Poulsen, J. R., Stevens, M. H. H., & White, J.-S. S. (2009). Generalized linear mixed models: A practical guide for ecology and evolution. Trends in Ecology & Evolution, 24(3), 127–135. https://doi.org/10.1016/j.tree.2008.10.008
• Brown, V. A. (2021). An introduction to linear mixed-effects modeling in R. Advances in Methods and Practices in Psychological Science, 4(1), 1–19. https://doi.org/10.1177/2515245920960351
• Harrison, X. A., Donaldson, L., Correa-Cano, M. E., Evans, J., Fisher, D. N., Goodwin, C. E. D., Robinson, B. S., Hodgson, D. J., & Inger, R. (2018). A brief introduction to mixed effects modelling and multi-model inference in ecology. PeerJ, 6, e4794. https://doi.org/10.7717/peerj.4794
• Meteyard, L., & Davies, R. A. I. (2020). Best practice guidance for linear mixed-effects models in psychological science. Journal of Memory and Language, 112, 104092. https://doi.org/10.1016/j.jml.2020.104092
• Pedersen, E. J., Miller, D. L., Simpson, G. L., & Ross, N. (2019). Hierarchical generalized additive models in ecology: An introduction with mgcv. PeerJ, 7, e6876. https://doi.org/10.7717/peerj.6876
• Thanks @Shawn Hemelstrand, this is very thorough and I'll definitely try it. However, as I've added in the edit, what happens if the effect is not linear? like a bell shaped or with a threshold level from which it starts climbing? Should I treat each concentration as an independent group? Oct 30, 2023 at 3:41
• I have edited my answer accordingly. Oct 30, 2023 at 3:55