How to handle Pseudoreplication - Aggregate observations per group for linear mixed model or take individual observations

Question

I´m trying to fit a linear mixed model to test if my dependent variable (B) is different between my treatments (A with 2 levels A&B). The study was conducted at severals sampling sites. To account for the paired data (measurements within one site) I would take sampling site as random effect into the model:

lmer(B ~ A + (1|Site))

However, to also account for the varibailty within each treamtent, serveral measurements were taken per site per treatment (i.e. field replicats, or pseudoreplicates, not the same number per treatment, as treatment B is expected to have a higher variabilty than treatment A). I now wonder if I need to aggregate my data first and take the mean value over the pseudoreplicates (for each treatment for each site) as dependent variable or every single observation? Or as a third option include some nested structure into the random effect to account for the field replicates.

My data looks something like this:

Site	Treatment A	Observation	Measurement B
1	A	1	5.3
1	A	2	5.3
1	A	3	8.1
1	B	1	7.2
1	B	2	5.2
2	A	1	2.6
2	A	2	3.3
2	A	3	1.3
2	B	1	4.1
2	B	2	6.3

Answer 1

I would not recommend aggregation as it would not allow you to account for within-treatment variability and would likely lose information about the variability of measurements, reducing statistical power.

I would rather suggest a model with nested random effects, where you account for the fact that observations are nested within treatments, which are in turn nested within sites. This model could look something like lmer(B ~ A + (1|Site/Treatment)), assuming that Treatment is nested within Site. This captures the variability at both the site level and the treatment level within each site. The notation (1|Site/Treatment) indicates that both site and treatment are random effects, with treatment nested within site.