Mixed model with raw data vs. linear model with aggregated data

Question

I have a designed field study in which 4 treatments (Tx = CT, RT, RL & TS) were applied in 4 blocks (Bk A-D). Measurements were taken in 4 places in each experimental unit; some measurements failed although there were at least 2 measurements in each unit. Response variable is Rv. I accept alpha = 0.1 probably of type I error. I analyzed the data in lme4 by comparing treatment + block to treatment alone (full dataset ‘ds1’, n=54)

m1 = lmer(Rv ~ Tx + (1|Bk)) 
m0 = lmer(Rv ~ (1|Bk), 
anova(m1, m0), 

The likelihood ratio test of the anova() gives p=0.3, indicating there is no chance that treatment is significant. Nevertheless the graph of the treatments encouraged me to investigate further so I tried the procedure outlined by Simon Wood (2006), of averaging by treatment and block, and using a linear model to estimate treatment effects, then estimating the variance of the random block (bk) effect from residual variances. (aggregated dataset n = 16)

ds2 <- ds1[ ,  .(Rv = mean(Rv)), by =c('Bk', 'Tx')] # aggregate by Tx & Bk 
lm1 = lm(rv ~ tx + bk)

Here, treatment (tx) is significant at p = 0.07, which is different than the answer provided by lme4. Wood (2006) says “assume that the data are balanced with respect to the model, meaning that for each factor or interaction in the model, the same number of data have been collected at each of its levels”. That appears to be untrue in my case, because of the varying number of measurements in each experimental unit. Should I just count Treatment (Tx) as non-significant, or is there any justification for using the Aggregated data / linear model analysis?

Tx=c('TS','TS','TS','TS','CT','CT','CT','CT','TS','TS','TS','TS','RL','RL','RL','RL','RT','RT','RT','RT','RL','RL','RT','RT','RT','RT','RL','RL','RL','RT','RT','RT','TS','TS','CT','CT','CT','CT','RL','RL','RL','RT','RT','RT','RT','CT','CT','CT','CT','CT','CT','TS','TS')
Bk= c('A','A','A','A','B','B','B','B','D','D','D','D','D','D','D','D','D','D','D','D','A','A','A','A','A','A','C','C','C','C','C',
'C','C','C','C','C','C','C','B','B','B','B','B','B','B','D','D','D','A','A','A','B','B')
Rv=c(2.08,2.08,2.52,3.42,2.8,5.57,2.53,3.69,1.55,1.45,3.98,3.19,2.3,2.09,2.26,2.1,3.21,2.99,2.11,2.09,1.64,1.74,1.66,6.41,1.86,2.71,1.83,0.86,2.37,1.05,1.37,2.08,1.09,1.44,0.6,1.24,3.32,1.34,1.86,4.54,2.7,2.5,4.93,2.85,3.42,2.77,2.71,4.11,5.29,2.16,3.15,4.58,2.89)

Answer 1

1. you have replication within block:Tx so m1 is a kind of pseudoreplication -- it will have > nominal Type I error. Two models that account for this are ~ Tx + (Tx|Bk) and ~ Tx + (1|Bk) + (1|Bk:Tx).
2. The random intercept model will not equal the aggregated model even with balance, because the RI model is effectively treating the replicates within block:Tx as independent.
3. The aggregated model is an alternative to the random intercept and slope model in item 1 but it has poorly behaved Type I error depending on structure of variance.
4. A model that is equivalent to Random Intercept and Slope model if balanced is the multivariate repeated measures ANOVA model in the afex package m3 <- aov_4(Rv ~ Tx + (Tx|Bk))

Some of this is explored in the largely complete chapter, including the simulation near the bottom (although I don't describe the simulation or results very thoroughly).

https://middleprofessor.github.io/benchbiostats/chapters/lmm.html