FAQ/rHSD

(There is also example in the post-hoc talk demo file available at the Graduate statistics talk page from here)

install.packages(c("MASS"))
install.packages(c("multcomp"))
library(multcomp)
library(MASS)
library(car)
library(nlme)

> oats

	B	V	N	Y
1	I	Victory	0.0cwt	111
2	I	Victory	0.2cwt	130
3	I	Victory	0.4cwt	157
4	I	Victory	0.6cwt	174
5	I	Golden.rain	0.0cwt	117
6	I	Golden.rain	0.2cwt	114
7	I	Golden.rain	0.4cwt	161
8	I	Golden.rain	0.6cwt	141
9	I	Marvellous	0.0cwt	105
10	I	Marvellous	0.2cwt	140
11	I	Marvellous	0.4cwt	118
12	I	Marvellous	0.6cwt	156
13	II	Victory	0.0cwt	61
14	II	Victory	0.2cwt	91
15	II	Victory	0.4cwt	97
16	II	Victory	0.6cwt	100
17	II	Golden.rain	0.0cwt	70
18	II	Golden.rain	0.2cwt	108
19	II	Golden.rain	0.4cwt	126
20	II	Golden.rain	0.6cwt	149
21	II	Marvellous	0.0cwt	96
22	II	Marvellous	0.2cwt	124
23	II	Marvellous	0.4cwt	121
24	II	Marvellous	0.6cwt	144
25	III	Victory	0.0cwt	68
26	III	Victory	0.2cwt	64
27	III	Victory	0.4cwt	112
28	III	Victory	0.6cwt	86
29	III	Golden.rain	0.0cwt	60
30	III	Golden.rain	0.2cwt	102
31	III	Golden.rain	0.4cwt	89
32	III	Golden.rain	0.6cwt	96
33	III	Marvellous	0.0cwt	89
34	III	Marvellous	0.2cwt	129
35	III	Marvellous	0.4cwt	132
36	III	Marvellous	0.6cwt	124
37	IV	Victory	0.0cwt	74
38	IV	Victory	0.2cwt	89
39	IV	Victory	0.4cwt	81
40	IV	Victory	0.6cwt	122
41	IV	Golden.rain	0.0cwt	64
42	IV	Golden.rain	0.2cwt	103
43	IV	Golden.rain	0.4cwt	132
44	IV	Golden.rain	0.6cwt	133
45	IV	Marvellous	0.0cwt	70
46	IV	Marvellous	0.2cwt	89
47	IV	Marvellous	0.4cwt	104
48	IV	Marvellous	0.6cwt	117
49	V	Victory	0.0cwt	62
50	V	Victory	0.2cwt	90
51	V	Victory	0.4cwt	100
52	V	Victory	0.6cwt	116
53	V	Golden.rain	0.0cwt	80
54	V	Golden.rain	0.2cwt	82
55	V	Golden.rain	0.4cwt	94
56	V	Golden.rain	0.6cwt	126
57	V	Marvellous	0.0cwt	63
58	V	Marvellous	0.2cwt	70
59	V	Marvellous	0.4cwt	109
60	V	Marvellous	0.6cwt	99
61	VI	Victory	0.0cwt	53
62	VI	Victory	0.2cwt	74
63	VI	Victory	0.4cwt	118
64	VI	Victory	0.6cwt	113
65	VI	Golden.rain	0.0cwt	89
66	VI	Golden.rain	0.2cwt	82
67	VI	Golden.rain	0.4cwt	86
68	VI	Golden.rain	0.6cwt	104
69	VI	Marvellous	0.0cwt	97
70	VI	Marvellous	0.2cwt	99
71	VI	Marvellous	0.4cwt	119
72	VI	Marvellous	0.6cwt	121

Aov.mod <- aov(Y ~ N + V + Error(B/V), data = oats) 
Lme.mod <- lme(Y ~ N + V, random = ~1 | B/V, data = oats) 
summary(Aov.mod)

Error = B

				Df	Sum Sq	Mean Sq	F
Residuals				5	15875	3175.1

Error = B:V

	Df	Sum Sq	Mean Sq	F
V	2	1786.4	893.18	1.4853
Residuals	10	6013.3	601.33

Error = Within

	Df	Sum Sq	Mean Sq	F
N	3	20020.5	6673.5	41.053
Residuals	51	8290.5	162.6

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(Lme.mod)

	numDF	denDF	F-value	p-value
(Intercept)	1	51	245.14566	<.0001
N	3	51	41.05276	<.0001
V	2	10	1.48534	0.2724

summary(Lme.mod)

Linear mixed-effects model fit by REML

Data: oats

AIC			BIC	logLik
586.0688			605.7756	-284.0344

Random effects:

Formula: ~1 | B
- (Intercept)

StdDev: 14.64486

Formula: ~1 | V %in% B
- (Intercept) Residual

StdDev: 10.47344 12.74987

Fixed effects: Y ~ N + V

	Value	Std.Error	DF	t-value	p-value
(Intercept)	79.91667	8.220345	51	9.721814	0.0000
N0.2cwt	19.50000	4.249956	51	4.588283	0.0000
N0.4cwt	34.83333	4.249956	51	8.196164	0.0000
N0.6cwt	44.00000	4.249956	51	10.353049	0.0000
VMarvellous	5.29167	7.078905	10	0.747526	0.4720
VVictory	-6.87500	7.078905	10	-0.971195	0.3544

Correlation:

(Intr)	N0.2cw	N0.4cw	N06.cw	VMrv11
-0.259	-	-	-	-
-0.259	0.500	-	-	-
-0.259	0.500	0.500	-	-
-0.431	0.000	0.000	0.000	-
-0.431	0.000	0.000	0.000	0.500

Standardized Within-Group Residuals:

Min				Q1	Med	Q3	Max
-1.84134117				-0.66279719	-0.06694235	0.63822499	1.66 066729

Number of Observations: 72

Number of Groups:

|||||| B || V %in% B ||

summary(glht(Lme.mod, linfct=mcp(V="Tukey")))

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: lme.formula(fixed = Y ~ N + V, data = oats, random = ~1 | B/V)

Linear Hypotheses:

	Estimate	Std. Error	z value	Pr(>\|z\|)
Marvellous - Golden.rain == 0	5.292	7.079	0.748	0.735
Victory - Golden.rain == 0	-6.875	7.079	-0.971	0.595
Victory - Marvellous == 0	-12.167	7.079	-1.719	0.198

(Adjusted p values reported -- single-step method)

Similarly for the within subjects effect

summary(glht(Lme.mod, linfct=mcp(N="Tukey")))

Using Tukey's test for main effects to correct for multiplicity in a repeated measures ANOVA using R