[R] trouble understanding syntax for glht methods

[Christopher W. Ryan]

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