Rafael, se entendi corretamente, você fará um experimento em DBC com (sub)-repetições dentro da parcela, certo? É o mesmo caso de quando se mede altura (da planta, da primeira espiga/vagen) em várias plantas na mesma parcela, ou conta-se o número de vagens, enfim. Parece que a sua dúvida é quanto a análise usando valores individuais (4 x 4 x 10 = 160 registros) ou com os valores agregados (4 x 4 = 16 registros). Se você usar uma agregação de soma ou média, para comparar os tratamentos o resultado é o mesmo (F e p-valor), como demonstra o código abaixo a partir de dados gerados.

> # Tabela de dados artificiais (em inclusão de efeitos, só ruído branco).
> da <- expand.grid(bloc = gl(4, 1),
+                   trt = gl(4, 1),
+                   rept = gl(10, 1))
> da$parc <- with(da, interaction(bloc, trt))
> set.seed(2017)
> da$y <- rnorm(nrow(da))
>
> # Modelo com observações dentro de parcela (dois estratos).
> m0 <- aov(y ~ bloc + trt + Error(parc/rept), data = da)
> summary(m0)

Error: parc
          Df Sum Sq Mean Sq F value Pr(>F)
bloc       3  0.945  0.3150   0.442  0.728
trt        3  1.433  0.4776   0.671  0.591
Residuals  9  6.408  0.7120              

Error: parc:rept
           Df Sum Sq Mean Sq F value Pr(>F)
Residuals 144  150.2   1.043              
>
> # Redução com média por parcela.
> db <- aggregate(y ~ bloc + trt, data = da, FUN = mean)
>
> # Modelo feito sobre as médias das parcelas.
> m1 <- lm(y ~ bloc + trt, data = db)
> anova(m1)
Analysis of Variance Table

Response: y
          Df  Sum Sq  Mean Sq F value Pr(>F)
bloc       3 0.09449 0.031496  0.4423 0.7285
trt        3 0.14329 0.047764  0.6708 0.5911
Residuals  9 0.64081 0.071201              
>

Usando a média, o quadro de anova do estrato `parc` é igual ao quadro de anova utilizando a média como resposta. Portanto, em termos de inferência sobre os tratamentos, é a mesma coisa. A análise com estrato é mais informativa porque te permite compreender a maginitude das variância entre parcelas e dentro de parcelas quando você obtém os componentes de variância, o que é útil para decisões de planejamento e dimensionamento amostral. Por outro lado, é mais fácil trabalhar com objetos de classe lm (m1) do que aov.list (m0), então, depois de ver o resultado do m0 eu seguiria a análise com o objeto m1, onde posso facilmente aplicar uma série de métodos e outras funções.

Se o número de observações for diferente em cada parcela, é necessário representar isso com peso na análise (lm(..., weigths)). Os pesos foram educadamente ignorados nessa análise que fiz porque o número de observações foi o mesmo. Veja alguns exemplos em: http://leg.ufpr.br/~walmes/cursoR/geneticaEsalq/script07.html.

À diposição.

Walmes.