1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
R version 3.2.3 (2015-12-10) -- "Wooden Christmas-Tree"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
>
> library("multcomp")
Loading required package: mvtnorm
Loading required package: survival
Loading required package: TH.data
Loading required package: MASS
Attaching package: 'TH.data'
The following object is masked from 'package:MASS':
geyser
> tol <- sqrt(.Machine$double.eps)
> set.seed(29081975)
>
> df <- data.frame(y = rnorm(100),
+ x = runif(100),
+ z = runif(100))
>
> ### linear model
> fam <- gaussian()
> lm0 <- glm(y ~ 1, data = df, family = fam)
> lm1 <- glm(y ~ x, data = df, family = fam)
> lm2 <- glm(y ~ x + z, data = df, family = fam)
>
> gh <- glht(lm2, linfct = c("x = 0", "z = 0"))
> stopifnot(abs(anova(lm0, lm2, test = "F")[2, 6] -
+ summary(gh, test = Ftest())$test$pvalue) < tol)
> stopifnot(abs(anova(lm0, lm2, test = "Chisq")[2, 5] -
+ summary(gh, test = Chisqtest())$test$pvalue) < tol)
>
> gh <- glht(lm2, linfct = "z = 0")
> stopifnot(abs(anova(lm1, lm2, test = "F")[2, 6] -
+ summary(gh, test = Ftest())$test$pvalue) < tol)
> stopifnot(abs(anova(lm1, lm2, test = "Chisq")[2, 5] -
+ summary(gh, test = Chisqtest())$test$pvalue) < tol)
>
> ### logistic regression
> df$y <- factor(df$y < 0)
> fam <- binomial()
> lm0 <- glm(y ~ 1, data = df, family = fam)
> lm1 <- glm(y ~ x, data = df, family = fam)
> lm2 <- glm(y ~ x + z, data = df, family = fam)
>
> if (require("lmtest")) {
+
+ gh <- glht(lm2, linfct = c("x = 0", "z = 0"))
+ stopifnot(abs(waldtest(lm0, lm2, test = "Chisq")[2, 4] -
+ summary(gh, test = Chisqtest())$test$pvalue) < tol)
+
+ gh <- glht(lm2, linfct = "z = 0")
+ stopifnot(abs(waldtest(lm1, lm2, test = "Chisq")[2, 4] -
+ summary(gh, test = Chisqtest())$test$pvalue) < tol)
+ }
Loading required package: lmtest
Loading required package: zoo
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
>
>
> proc.time()
user system elapsed
0.284 0.016 0.298
|
