Multivariate analysis for repeated measures with dependent treatment groups and independent control groups

I conducted experiments to test the effect of a treatment on multiple response variables. Measurements were taken over consecutive rounds (i.e., time steps) of an experiment, with each round containing a treatment group and a control group. The treatment groups were dependent across rounds. The control groups were independent across rounds (e.g., controls accounted for environmental conditions, etc. during the round).

I want to determine if there are significant differences in response variables between the control and treatment groups, and whether differences between the control and treatment groups change across rounds (i.e., over time). What is the appropriate strategy to answer these questions?

Simulated example: I want to test the effect of soil reuse on plant growth variables (e.g., growth rate, stem height, etc). In round 1 of the experiment, the control group has 3 plots of “new” soil, and the treatment group has 3 plots of “reused” soil that have already been used once to grow plants. In round 2, the control group has 3 plots of “new” soil again, and the treatment group has the same 3 plots of “reused” soil that have now been used twice, etc for round 3. Rounds occur at equal time intervals.

R code for simulated data:

set.seed(1000)
data <- data.frame(round = c(rep(1,6), rep(2,6), rep(3,6)),
treatment = rep( c(rep("new",3), rep("reused",3)), 3),
soil_plot = rep( c("A", "B", "C", "D", "E", "F"), 3),
var1 = rnorm(18, 700, 200),
var2 = rnorm(18, 1, 0.08),
var3 = rnorm(18, 7.5, 0.9),
var4 = rnorm(18, 0.5, 0.02))


My attempted strategies have included MANOVA with and without an Error term in the model to specify repeated subjects. The response variables do not have high Pearson correlations (< 0.6), so I don’t think I could assume some response variables behave similarly and therefore model a reduced number of response variables.

# without Error term
manova1 <- manova(cbind(var1, var2, var3, var4) ~ treatment*as.ordered(round), data = data)
summary(manova1)

# with Error term
manova2 <- manova(cbind(var1, var2, var3, var4) ~ treatment*as.ordered(round) + Error(soil_plot), data = data)
summary(manova2)


If there was only one response variable, the best solution I’ve come across is a linear mixed effects model with the soil plots as a random factor (using package nlme). Putting multiple response variables in this function, grouped by cbind(), seems to only use the first response variable listed.

install.packages("nlme")
library(nlme)
install.packages("car")
library(car) # for Anova

lme_var1 <- lme(var1 ~ as.ordered(round)*treatment, random = ~1 | soil_plot, correlation = corAR1(form = ~ round | soil_plot), data = data, method = "REML")
Anova(lme_var1)

lme_all <- lme(cbind(var1, var2, var3, var4) ~ as.ordered(round)*treatment, random = ~1 | soil_plot, correlation = corAR1(form = ~ round | soil_plot), data = data, method = "REML")
Anova(lme_all)


However, none of these strategies takes into account that the control groups are independent across rounds AND that treatment groups are dependent. In the simulated data, I could rename the control group plots so that they are each unique (instead of A,B,C repeated in each round), however this would create an unbalanced model.