# R lsmeans in a model with two residuals (split-plot) [closed]

Let's assume the model:

lm(VAR ~ A * B + (A : R), data)


which produces this ANOVA:

Analysis of Variance Table

Response: VAR
Df  Sum Sq Mean Sq F value    Pr(>F)
A          2 2444.07 1222.04 71.4330 1.086e-14 ***
B          3 2370.92  790.31 46.1966 8.675e-14 ***
A:B        6 1376.40  229.40 13.4094 1.189e-08 ***
A:R       15  379.18   25.28  1.4776    0.1548
Residuals 45  769.84   17.11
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


This is an analysis of a split-plot design, with two residuals (A:R and Residuals). We know that, for this type of design, the A means should be tested with the 'A:R' residuals and B means with 'Residuals' (which also makes the F value of A wrong, but this is not the point).

I am trying to perform a post-hoc Tukey test using the lsmeans package, but I could not manage to get my confidence intervals right, considering the mentioned before. It is very important for me to use the package because some of my variables are unbalanced.

For clarification, if I try confint(pairs(lsmeans(model, pairwise ~ A))), I get:

NOTE: Results may be misleading due to involvement in interactions
contrast      estimate       SE df  lower.CL  upper.CL
control - T1 13.999669 1.193994 45 11.105889 16.893450
control - T2  9.400151 1.193994 45  6.506371 12.293932
T1 - T2      -4.599518 1.193994 45 -7.493299 -1.705738

Results are averaged over the levels of: B, R
Confidence level used: 0.95
Conf-level adjustment: tukey method for comparing a family of 3 estimates


Where, for the first line, the upper.CL value should be 17.76957.

How can I use the package for this situation?

Sample data:

structure(list(A = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L
), .Label = c("control", "T1", "T2"), class = "factor"), B = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L,
3L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L,
3L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("placebo", "D1", "D2",
"D3"), class = c("ordered", "factor")), VAR = c(49.4461686554634,
52.8401943062405, 53.1493049522644, 51.4218622201876, 49.5739625365273,
44.7711140039805, 45.5892748670425, 49.2934242984433, 48.1124355048776,
44.373148667141, 52.8742228937389, 44.751429590018, 47.6410725365848,
51.0997809793751, 34.6423104049863, 44.863129361936, 52.5070664501272,
45.246006557, 42.4742980179221, 50.921529132397, 46.673788967046,
47.9566067375011, 46.8931285459289, 52.4795249482172, 50.9005883761804,
48.4487973095543, 48.0595406895989, 43.7080988978334, 50.6670099107763,
57.613741932035, 42.3081168342102, 35.2048805044803, 40.4268626953135,
30.5542203074074, 36.4722288755433, 44.2375951451213, 28.9947194917908,
29.9990191696814, 32.8560294911122, 29.3787723397584, 17.8108088004077,
31.9639530468079, 15.0297564986717, 11.6103644309025, 23.2769424869608,
23.9774437540106, 18.7609811936338, 21.3422470375469, 49.4044826321325,
40.2495469338607, 47.1503186948931, 50.8081648047913, 45.2199530578737,
41.3630619250218, 38.6704818185838, 40.4887847543725, 35.8724077935333,
38.3050685669546, 42.9633536530607, 38.2957146954417, 36.0365213440119,
37.1530824135316, 37.1754729697022, 40.5131592687021, 38.7368678093011,
29.3870360968715, 31.0906971673814, 34.4100275541074, 40.1061102694397,
28.107451337622, 32.8203640321119, 29.6630255846913), R = structure(c(1L,
2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L,
6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L,
4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L,
2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L,
6L, 1L, 2L, 3L, 4L, 5L, 6L), .Label = c("1", "2", "3", "4", "5",
"6"), class = "factor")), .Names = c("A", "B", "VAR", "R"), row.names = c(NA,
-72L), class = "data.frame")


You should fit a mixed model to these data. I suggest using the lme4 package, and then if I understand it right, you'd do this:

library(lme4)
mod = lmer(VAR ~ A * B + (1 | R / A), data)


This model has A and B as fixed effects, R as random blocks, and A:R as whole plots. The lsmeans package does support these models, and lots of others. You should also install the pbkrtest package, which will be used to obtain the Kenward-Rogers degrees of freedom.

PS - You ask for pairwise comparisons twice in the command you show. It should be:

confint(pairs(lsmeans(model, ~ A)))

• Fine, thanks for the answer. Unfortunately I can't change the model, I need to stick with this nested design and OLS. Do you have any suggestion for manually fixing this CL values to fit my needs? Jul 31, 2015 at 13:50
• Given your requirements I guess your best option is to do it by hand. Jul 31, 2015 at 14:48