mmm.RdCalculation of correlation between test statistics from multiple marginal models using the score decomposition
mmm(...)
mlf(...)A names argument list containing fitted models (mmm) or
definitions of linear functions (mlf). If
only one linear function is defined for mlf,
it will be applied to all models in mmm by
glht.mlf.
Estimated correlations of the estimated parameters of interest from the multiple marginal models are obtained using a stacked version of the i.i.d. decomposition of parameter estimates by means of score components (first derivatives of the log likelihood). The method is less conservative than the Bonferroni correction. The details are provided by Pipper, Ritz and Bisgaard (2012).
The implementation assumes that the model were fitted to the same data,
i.e., the rows of the matrices returned by estfun belong to the
same observations for each model.
The reference distribution is always multivariate normal, if you want
to use the multivariate t, please specify the corresponding degrees of
freedom as an additional df argument to glht.
Observations with missing values contribute zero to the score function.
Models have to be fitted using na.exclude as na.action
argument.
An object of class mmm or mlf, basically a named list of the
arguments with a special method for glht being available for
the latter. vcov, estfun, and
bread methods are available for objects of class
mmm.
Christian Bressen Pipper, Christian Ritz and Hans Bisgaard (2011), A Versatile Method for Confirmatory Evaluation of the Effects of a Covariate in Multiple Models, Journal of the Royal Statistical Society, Series C (Applied Statistics), 61, 315–326.
### replicate analysis of Hasler & Hothorn (2011),
### A Dunnett-Type Procedure for Multiple Endpoints,
### The International Journal of Biostatistics: Vol. 7: Iss. 1, Article 3.
### DOI: 10.2202/1557-4679.1258
library("sandwich")
### see ?coagulation
if (require("SimComp")) {
data("coagulation", package = "SimComp")
### level "S" is the standard, "H" and "B" are novel procedures
coagulation$Group <- relevel(coagulation$Group, ref = "S")
### fit marginal models
(m1 <- lm(Thromb.count ~ Group, data = coagulation))
(m2 <- lm(ADP ~ Group, data = coagulation))
(m3 <- lm(TRAP ~ Group, data = coagulation))
### set-up Dunnett comparisons for H - S and B - S
### for all three models
g <- glht(mmm(Thromb = m1, ADP = m2, TRAP = m3),
mlf(mcp(Group = "Dunnett")), alternative = "greater")
### joint correlation
cov2cor(vcov(g))
### simultaneous p-values adjusted by taking the correlation
### between the score contributions into account
summary(g)
### simultaneous confidence intervals
confint(g)
### compare with
if (FALSE) { # \dontrun{
library("SimComp")
SimCiDiff(data = coagulation, grp = "Group",
resp = c("Thromb.count","ADP","TRAP"),
type = "Dunnett", alternative = "greater",
covar.equal = TRUE)
} # }
### use sandwich variance matrix
g <- glht(mmm(Thromb = m1, ADP = m2, TRAP = m3),
mlf(mcp(Group = "Dunnett")),
alternative = "greater", vcov. = sandwich)
summary(g)
confint(g)
}
#> Loading required package: SimComp
#>
#> Simultaneous Confidence Intervals
#>
#> Fit: NULL
#>
#> Quantile = -2.2694
#> 95% family-wise confidence level
#>
#>
#> Linear Hypotheses:
#> Estimate lwr upr
#> Thromb: B - S <= 0 0.12170 -0.07616 Inf
#> Thromb: H - S <= 0 0.04351 -0.17922 Inf
#> ADP: B - S <= 0 0.21211 0.03817 Inf
#> ADP: H - S <= 0 0.08422 -0.06893 Inf
#> TRAP: B - S <= 0 0.10525 -0.20317 Inf
#> TRAP: H - S <= 0 0.07109 -0.24455 Inf
#>
### attitude towards science data
data("mn6.9", package = "TH.data")
### one model for each item
mn6.9.y1 <- glm(y1 ~ group, family = binomial(),
na.action = na.omit, data = mn6.9)
mn6.9.y2 <- glm(y2 ~ group, family = binomial(),
na.action = na.omit, data = mn6.9)
mn6.9.y3 <- glm(y3 ~ group, family = binomial(),
na.action = na.omit, data = mn6.9)
mn6.9.y4 <- glm(y4 ~ group, family = binomial(),
na.action = na.omit, data = mn6.9)
### test all parameters simulaneously
summary(glht(mmm(mn6.9.y1, mn6.9.y2, mn6.9.y3, mn6.9.y4),
mlf(diag(2))))
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Linear Hypotheses:
#> Estimate Std. Error z value Pr(>|z|)
#> mn6.9.y1: 1 == 0 1.47323 0.06649 22.156 <1e-04 ***
#> mn6.9.y1: 2 == 0 0.14833 0.09638 1.539 0.5806
#> mn6.9.y2: 1 == 0 -1.38713 0.06476 -21.419 <1e-04 ***
#> mn6.9.y2: 2 == 0 -0.19598 0.09455 -2.073 0.2330
#> mn6.9.y3: 1 == 0 0.44024 0.05306 8.298 <1e-04 ***
#> mn6.9.y3: 2 == 0 -0.37449 0.07417 -5.049 <1e-04 ***
#> mn6.9.y4: 1 == 0 0.16537 0.05197 3.182 0.0107 *
#> mn6.9.y4: 2 == 0 -0.37132 0.07357 -5.047 <1e-04 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> (Adjusted p values reported -- single-step method)
#>
### group differences
summary(glht(mmm(mn6.9.y1, mn6.9.y2, mn6.9.y3, mn6.9.y4),
mlf("group2 = 0")))
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Linear Hypotheses:
#> Estimate Std. Error z value Pr(>|z|)
#> mn6.9.y1: group2 == 0 0.14833 0.09638 1.539 0.409
#> mn6.9.y2: group2 == 0 -0.19598 0.09455 -2.073 0.144
#> mn6.9.y3: group2 == 0 -0.37449 0.07417 -5.049 <1e-04 ***
#> mn6.9.y4: group2 == 0 -0.37132 0.07357 -5.047 <1e-04 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> (Adjusted p values reported -- single-step method)
#>
### alternative analysis of Klingenberg & Satopaa (2013),
### Simultaneous Confidence Intervals for Comparing Margins of
### Multivariate Binary Data, CSDA, 64, 87-98
### http://dx.doi.org/10.1016/j.csda.2013.02.016
### see supplementary material for data description
### NOTE: this is not the real data but only a subsample
influenza <- structure(list(
HEADACHE = c(1L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L,
0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L,
1L, 1L), MALAISE = c(0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L,
0L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L,
0L), PYREXIA = c(0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L,
1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L
), ARTHRALGIA = c(0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L,
0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L
), group = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L), .Label = c("pla", "trt"), class = "factor"), Freq = c(32L,
165L, 10L, 23L, 3L, 1L, 4L, 2L, 4L, 2L, 1L, 1L, 1L, 1L, 167L,
1L, 11L, 37L, 7L, 7L, 5L, 3L, 3L, 1L, 2L, 4L, 2L)), .Names = c("HEADACHE",
"MALAISE", "PYREXIA", "ARTHRALGIA", "group", "Freq"), row.names = c(1L,
2L, 3L, 5L, 9L, 36L, 43L, 50L, 74L, 83L, 139L, 175L, 183L, 205L,
251L, 254L, 255L, 259L, 279L, 281L, 282L, 286L, 302L, 322L, 323L,
366L, 382L), class = "data.frame")
influenza <- influenza[rep(1:nrow(influenza), influenza$Freq), 1:5]
### Fitting marginal logistic regression models
(head_logreg <- glm(HEADACHE ~ group, data = influenza,
family = binomial()))
#>
#> Call: glm(formula = HEADACHE ~ group, family = binomial(), data = influenza)
#>
#> Coefficients:
#> (Intercept) grouptrt
#> -0.8664 -0.1384
#>
#> Degrees of Freedom: 499 Total (i.e. Null); 498 Residual
#> Null Deviance: 594.8
#> Residual Deviance: 594.3 AIC: 598.3
(mala_logreg <- glm(MALAISE ~ group, data = influenza,
family = binomial()))
#>
#> Call: glm(formula = MALAISE ~ group, family = binomial(), data = influenza)
#>
#> Coefficients:
#> (Intercept) grouptrt
#> -1.9551 0.4114
#>
#> Degrees of Freedom: 499 Total (i.e. Null); 498 Residual
#> Null Deviance: 422.7
#> Residual Deviance: 420 AIC: 424
(pyre_logreg <- glm(PYREXIA ~ group, data = influenza,
family = binomial()))
#>
#> Call: glm(formula = PYREXIA ~ group, family = binomial(), data = influenza)
#>
#> Coefficients:
#> (Intercept) grouptrt
#> -2.389 -1.158
#>
#> Degrees of Freedom: 499 Total (i.e. Null); 498 Residual
#> Null Deviance: 215.8
#> Residual Deviance: 208.1 AIC: 212.1
(arth_logreg <- glm(ARTHRALGIA ~ group, data = influenza,
family = binomial()))
#>
#> Call: glm(formula = ARTHRALGIA ~ group, family = binomial(), data = influenza)
#>
#> Coefficients:
#> (Intercept) grouptrt
#> -2.6178 -0.1337
#>
#> Degrees of Freedom: 499 Total (i.e. Null); 498 Residual
#> Null Deviance: 237.8
#> Residual Deviance: 237.7 AIC: 241.7
### Simultaneous inference for log-odds
xy.sim <- glht(mmm(head = head_logreg,
mala = mala_logreg,
pyre = pyre_logreg,
arth = arth_logreg),
mlf("grouptrt = 0"))
summary(xy.sim)
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Linear Hypotheses:
#> Estimate Std. Error z value Pr(>|z|)
#> head: grouptrt == 0 -0.1384 0.1990 -0.695 0.9181
#> mala: grouptrt == 0 0.4114 0.2538 1.621 0.3354
#> pyre: grouptrt == 0 -1.1580 0.4460 -2.596 0.0356 *
#> arth: grouptrt == 0 -0.1337 0.3661 -0.365 0.9918
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> (Adjusted p values reported -- single-step method)
#>
confint(xy.sim)
#>
#> Simultaneous Confidence Intervals
#>
#> Fit: NULL
#>
#> Quantile = 2.4753
#> 95% family-wise confidence level
#>
#>
#> Linear Hypotheses:
#> Estimate lwr upr
#> head: grouptrt == 0 -0.13837 -0.63087 0.35412
#> mala: grouptrt == 0 0.41140 -0.21680 1.03960
#> pyre: grouptrt == 0 -1.15795 -2.26207 -0.05384
#> arth: grouptrt == 0 -0.13371 -1.03995 0.77253
#>
### Artificial examples
### Combining linear regression and logistic regression
set.seed(29)
y1 <- rnorm(100)
y2 <- factor(y1 + rnorm(100, sd = .1) > 0)
x1 <- gl(4, 25)
x2 <- runif(100, 0, 10)
m1 <- lm(y1 ~ x1 + x2)
m2 <- glm(y2 ~ x1 + x2, family = binomial())
### Note that the same explanatory variables are considered in both models
### but the resulting parameter estimates are on 2 different scales
### (original and log-odds scales)
### Simultaneous inference for the same parameter in the 2 model fits
summary(glht(mmm(m1 = m1, m2 = m2), mlf("x12 = 0")))
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Linear Hypotheses:
#> Estimate Std. Error z value Pr(>|z|)
#> m1: x12 == 0 -0.3537 0.2850 -1.241 0.321
#> m2: x12 == 0 -0.6409 0.5765 -1.112 0.392
#> (Adjusted p values reported -- single-step method)
#>
### Simultaneous inference for different parameters in the 2 model fits
summary(glht(mmm(m1 = m1, m2 = m2),
mlf(m1 = "x12 = 0", m2 = "x13 = 0")))
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Linear Hypotheses:
#> Estimate Std. Error z value Pr(>|z|)
#> m1: x12 == 0 -0.3537 0.2850 -1.241 0.362
#> m2: x13 == 0 -0.8264 0.5824 -1.419 0.270
#> (Adjusted p values reported -- single-step method)
#>
### Simultaneous inference for different and identical parameters in the 2
### model fits
summary(glht(mmm(m1 = m1, m2 = m2),
mlf(m1 = c("x12 = 0", "x13 = 0"), m2 = "x13 = 0")))
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Linear Hypotheses:
#> Estimate Std. Error z value Pr(>|z|)
#> m1: x12 == 0 -0.3537 0.2850 -1.241 0.418
#> m1: x13 == 0 -0.4220 0.2849 -1.481 0.286
#> m2: x13 == 0 -0.8264 0.5824 -1.419 0.317
#> (Adjusted p values reported -- single-step method)
#>
### Examples for binomial data
### Two independent outcomes
y1.1 <- rbinom(100, 1, 0.45)
y1.2 <- rbinom(100, 1, 0.55)
group <- factor(rep(c("A", "B"), 50))
m1 <- glm(y1.1 ~ group, family = binomial)
m2 <- glm(y1.2 ~ group, family = binomial)
summary(glht(mmm(m1 = m1, m2 = m2),
mlf("groupB = 0")))
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Linear Hypotheses:
#> Estimate Std. Error z value Pr(>|z|)
#> m1: groupB == 0 0.41502 0.40905 1.015 0.523
#> m2: groupB == 0 0.08161 0.40407 0.202 0.974
#> (Adjusted p values reported -- single-step method)
#>
### Two perfectly correlated outcomes
y2.1 <- rbinom(100, 1, 0.45)
y2.2 <- y2.1
group <- factor(rep(c("A", "B"), 50))
m1 <- glm(y2.1 ~ group, family = binomial)
m2 <- glm(y2.2 ~ group, family = binomial)
summary(glht(mmm(m1 = m1, m2 = m2),
mlf("groupB = 0")))
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Linear Hypotheses:
#> Estimate Std. Error z value Pr(>|z|)
#> m1: groupB == 0 0.2427 0.4028 0.603 0.547
#> m2: groupB == 0 0.2427 0.4028 0.603 0.547
#> (Adjusted p values reported -- single-step method)
#>
### use sandwich covariance matrix
summary(glht(mmm(m1 = m1, m2 = m2),
mlf("groupB = 0"), vcov. = sandwich))
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Linear Hypotheses:
#> Estimate Std. Error z value Pr(>|z|)
#> m1: groupB == 0 0.2427 0.4028 0.603 0.547
#> m2: groupB == 0 0.2427 0.4028 0.603 0.547
#> (Adjusted p values reported -- single-step method)
#>