Christopher W. Ryan
2022-Jun-17 22:46 UTC
[R] trouble understanding syntax for glht methods
I'm running R 4.2.0 in Windows 10.
What is the proper syntax for plotting simultaneous CIs from a glht
object, choosing one's method of multiplicity correction? I seem to be
able to specify a method in summary(some_glht) that but my attempts to
do so in confint() or plot() or plot(confint()) seem to have little or
no effect on the output or on the plot produced.
library(nlme)
library(multcomp)
## subset of 20 of my observations
## call it data.temp
data.temp <- structure(list(false.date = structure(c(-7915, -9772, -10168,
-7519, -9923, -8766, -8097, -10319, -10227, -8735, -10592, -10288,
-9831, -8370, -8220, -8158, -10380, -8493, -7762, -8067), class =
"Date"),
md = c(67, 252, 368, 36, 215, 170, 225, 158, 388, 166, 239,
136, 234, 108, 27, 70, 28, 136, 128, 100), ad = c(1490, 1567,
1541, 857, 1412, 1660, 1252, 1672, 1484, 1539, 1027, 1098,
1413, 1575, 1175, 1620, 707, 1796, 1575, 1067), guv = c(1,
0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1),
season = structure(c(4L, 4L, 4L, 1L, 2L, 3L, 2L, 2L, 3L,
3L, 3L, 2L, 3L, 3L, 1L, 2L, 1L, 2L, 2L, 3L), levels = c("summer",
"fall", "winter", "spring"), class =
"factor"), md.prop c(0.0449664429530201,
0.16081684747926, 0.238805970149254, 0.0420070011668611,
0.152266288951841, 0.102409638554217, 0.179712460063898,
0.0944976076555024, 0.261455525606469, 0.107862248213125,
0.232716650438169, 0.123861566484517, 0.165605095541401,
0.0685714285714286, 0.0229787234042553, 0.0432098765432099,
0.0396039603960396, 0.0757238307349666, 0.0812698412698413,
0.0937207122774133), month = c(5, 4, 3, 6, 11, 1, 11, 10,
1, 2, 1, 11, 2, 2, 7, 9, 8, 10, 10, 12), time.index = c(89L,
28L, 15L, 102L, 23L, 61L, 83L, 10L, 13L, 62L, 1L, 11L, 26L,
74L, 79L, 81L, 8L, 70L, 94L, 84L), guv.int = c(1L, 0L, 0L,
1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L,
1L, 1L), guv.fac = structure(c(2L, 1L, 1L, 2L, 1L, 2L, 2L,
1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L), levels =
c("0",
"1"), class = "factor"), trend.temp = c(89L, 28L, 15L,
102L,
23L, 61L, 83L, 10L, 13L, 62L, 1L, 11L, 26L, 74L, 79L, 81L,
8L, 70L, 94L, 84L), my.logit = c(-3.05582997869889, -1.65216285710044,
-1.15923691048454, -3.12700417099632, -1.71693567743031,
-2.17073296188924, -1.51829680772414, -2.25991540097043,
-1.03841712788469, -2.11276561741143, -1.19303453792637,
-1.95635956492965, -1.61710078517967, -2.60884355101876,
-3.74993971087518, -3.09751496786393, -3.18841661738349,
-2.50191799561454, -2.42521746271187, -2.2690283094652)), row.names
= c(NA,
-20L), class = c("tbl_df", "tbl", "data.frame"))
## fit an AR(1) model via gls() in nlme
m.full.temp <- gls(my.logit ~ season + guv.fac + season:guv.fac +
time.index, weights = varWeights(varFunc(~guv.fac)), correlation corAR1(form = ~
time.index), data = data.temp, method = "ML")
set.seed(42)
temp <- glht(m.full.temp, linfct = c("guv.fac1 >= 0", ## summer
"guv.fac1 + seasonfall:guv.fac1 >=
0",
"guv.fac1 + seasonwinter:guv.fac1 >=
0",
"guv.fac1 + seasonspring:guv.fac1 >=
0"))
Without correction for multiplicity, versus with e.g. bonferroni
correction, the p-values from summary() behave more or less as I'd expect:
> summary(temp, adjusted(type = "none"))
[some output truncated]
Linear Hypotheses:
Estimate Std. Error z value Pr(<z)
guv.fac1 >= 0 -0.2828 1.0911 -0.259 0.3978
guv.fac1 + seasonfall:guv.fac1 >= 0 -0.4380 0.8441 -0.519 0.3019
guv.fac1 + seasonwinter:guv.fac1 >= 0 -1.0424 0.7388 -1.411 0.0791 .
guv.fac1 + seasonspring:guv.fac1 >= 0 -1.6771 0.9436 -1.777 0.0378
(Adjusted p values reported -- none method)
> summary(temp, adjusted(type = "bonferroni"))
Linear Hypotheses:
Estimate Std. Error z value Pr(<z)
guv.fac1 >= 0 -0.2828 1.0911 -0.259 1.000
guv.fac1 + seasonfall:guv.fac1 >= 0 -0.4380 0.8441 -0.519 1.000
guv.fac1 + seasonwinter:guv.fac1 >= 0 -1.0424 0.7388 -1.411 0.317
guv.fac1 + seasonspring:guv.fac1 >= 0 -1.6771 0.9436 -1.777 0.151
(Adjusted p values reported -- bonferroni method)
But confint does not:
> confint(temp, adjusted(type = "none"))
Quantile = 2.0424
95% family-wise confidence level
Linear Hypotheses:
Estimate lwr upr
guv.fac1 >= 0 -0.2828 -Inf 1.9457
guv.fac1 + seasonfall:guv.fac1 >= 0 -0.4380 -Inf 1.2858
guv.fac1 + seasonwinter:guv.fac1 >= 0 -1.0424 -Inf 0.4664
guv.fac1 + seasonspring:guv.fac1 >= 0 -1.6771 -Inf 0.2500
> confint(temp, adjusted(type = "bonferroni"))
Quantile = 2.0381
95% family-wise confidence level
Linear Hypotheses:
Estimate lwr upr
guv.fac1 >= 0 -0.2828 -Inf 1.9411
guv.fac1 + seasonfall:guv.fac1 >= 0 -0.4380 -Inf 1.2823
guv.fac1 + seasonwinter:guv.fac1 >= 0 -1.0424 -Inf 0.4633
guv.fac1 + seasonspring:guv.fac1 >= 0 -1.6771 -Inf 0.2460
Same thing with plot(confint()). Trying to use adjusted(type
"some_method_from_p.adjust") yields identical graphs regardless of the
adjustment method I type, or "none". What is the correct syntax for
specifying the multiplicity adjustment for displaying and plotting CIs?
Thanks.
--Chris Ryan