На основании кода ОП
df1 %>%
group_by(Treatment, Rep) %>%
summarise_at(vars(-group_cols()), list(mean = ~mean(Weight, na.rm = TRUE),
sd = ~sd(Weight, na.rm = TRUE),
se= ~sd(Weight, na.rm = TRUE)/sqrt(n()))) %>%
summarise_at(vars(mean:se), mean, na.rm = TRUE)
данные
df1 <- structure(list(Treatment = c("Line 1", "Line 1", "Line 1", "Line 1",
"Line 2", "Line 2", "Line 2", "Line 2", "Line 1", "Line 1", "Line 1",
"Line 1", "Line 2", "Line 2", "Line 2", "Line 2"), Rep = c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L),
Weight = c(NA, NA, NA, NA, 26L, 26L, 26L, 27L, 26L, 28L,
26L, 25L, 24L, 26L, 25L, NA)), class = "data.frame", row.names = c(NA,
-16L))