R Under development (unstable) (2023-06-12 r84534 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

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

> 
> ###<FIXME> compare results of mmod and glht.mlf </FIXME>
> 
> ### code by Christian Ritz
> "mmod" <- function(modelList, varName, seType = "san")
+ {
+     require(multcomp, quietly = TRUE)
+     require(sandwich, quietly = TRUE)
+ 
+     if (length(seType) == 1) {seType <- rep(seType, length(modelList))}
+     if (length(varName) == 1) {varName <- rep(varName, length(modelList))}    
+ 
+     ## Extracting score contributions from the individual model fits    
+     makeIIDdecomp <- function(modelObject, varName) 
+     {
+         numObsUsed <- ifelse(inherits(modelObject, "coxph"), modelObject$n, nrow(modelObject$model))
+         iidVec0 <- bread(modelObject)[varName, , drop = FALSE] %*% t(estfun(modelObject))
+         moNAac <- modelObject$na.action
+         numObs <- numObsUsed + length(moNAac)
+         iidVec <- rep(0, numObs)
+         if (!is.null(moNAac))
+         {
+             iidVec[-moNAac] <- sqrt(numObs/numObsUsed) * iidVec0
+         }
+         else {
+             iidVec <- iidVec0
+         }
+         list(iidVec = iidVec, numObsUsed = numObsUsed, numObs = numObs)
+     }
+     numModels <- length(modelList)
+     if (identical(length(varName), 1))
+     {
+         varName <- rep(varName, numModels)
+     }
+     iidList <- mapply(makeIIDdecomp, modelList, varName, SIMPLIFY = FALSE)
+     iidresp <- matrix(as.vector(unlist(lapply(iidList, function(listElt) {listElt[[1]]}))), nrow = numModels, byrow = TRUE)
+     pickFct <- function(modelObject, varName, matchStrings) 
+     {
+         as.vector(na.omit((coef(summary(modelObject))[varName, ])[matchStrings]))
+     }
+           
+     ## Retrieving parameter estimates from the individual fits
+     estVec <- as.vector(unlist(mapply(pickFct, modelList, varName, MoreArgs = list(matchStrings = c("Estimate", "coef")))))
+     # "Estimate" or "coef" used in glm(), lm() and coxph() summary output, respectively
+     
+     ## Calculating the estimated variance-covariance matrix of the parameter estimates
+     numObs <- iidList[[1]]$numObs
+     covar <- (iidresp %*% t(iidresp)) / numObs
+     vcMat <- covar / numObs  # Defining the finite-sample variance-covariance matrix
+ 
+     ## Replacing sandwich estimates by model-based standard errors
+     modbas <- seType == "mod"
+     if (any(modbas))
+     {  
+         corMat <- cov2cor(vcMat)
+         ## Retrieving standard errors for the specified estimate from the individual fits
+         modSE <- as.vector(unlist(mapply(pickFct, modelList, varName, MoreArgs = list(matchStrings = c("Std. Error", "se(coef)")))))
+         
+         sanSE <- sqrt(diag(vcMat))
+         sanSE[modbas] <- modSE[modbas]
+         vcMat <- diag(sanSE) %*% corMat %*% diag(sanSE)
+     } 
+ 
+     ## Naming the parameter vector (easier way to extract the names of the model fits provided as a list in the first argument?)
+     names1 <- sub("list", "", deparse(substitute(modelList)), fixed = TRUE)
+     names2 <- sub("(", "", names1, fixed = TRUE)
+     names3 <- sub(")", "", names2, fixed = TRUE)
+     names4 <- sub(" ", "", names3, fixed = TRUE)
+     names(estVec) <- unlist(strsplit(names4, ","))
+              
+     return(parm(coef = estVec, vcov = vcMat, df = 0))
+ }
>     
> 
> 
> 
> 
> set.seed(29)
> ## Combining linear regression and logistic regression
> 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
> simult.x12 <- mmod(list(m1, m2), c("x12", "x12"))
> summary(glht(simult.x12))

	 Simultaneous Tests for General Linear Hypotheses

Linear Hypotheses:
        Estimate Std. Error z value Pr(>|z|)
m1 == 0  -0.3537     0.2808  -1.260    0.312
m2 == 0  -0.6409     0.5681  -1.128    0.382
(Adjusted p values reported -- single-step method)

> 
> ## Simultaneous inference for different parameters in the 2 model fits
> simult.x12.x13 <- mmod(list(m1, m2), c("x12", "x13"))
> summary(glht(simult.x12.x13))

	 Simultaneous Tests for General Linear Hypotheses

Linear Hypotheses:
        Estimate Std. Error z value Pr(>|z|)
m1 == 0  -0.3537     0.2808  -1.260    0.351
m2 == 0  -0.8264     0.5874  -1.407    0.276
(Adjusted p values reported -- single-step method)

> 
> ## Simultaneous inference for different and identical parameters in the 2 model fits
> simult.x12x2.x13 <- mmod(list(m1, m1, m2), c("x12", "x13", "x13"))
> summary(glht(simult.x12x2.x13))

	 Simultaneous Tests for General Linear Hypotheses

Linear Hypotheses:
         Estimate Std. Error z value Pr(>|z|)
m1 == 0   -0.3537     0.2808  -1.260    0.407
m1 == 0   -0.4220     0.2801  -1.507    0.274
 m2 == 0  -0.8264     0.5874  -1.407    0.323
(Adjusted p values reported -- single-step method)

> confint(glht(simult.x12x2.x13))

	 Simultaneous Confidence Intervals

Fit: NULL

Quantile = 2.3087
95% family-wise confidence level
 

Linear Hypotheses:
         Estimate lwr     upr    
m1 == 0  -0.3537  -1.0019  0.2945
m1 == 0  -0.4220  -1.0687  0.2247
 m2 == 0 -0.8264  -2.1825  0.5297

> 
> 
> ## Examples for binomial data
> 
> ## Two independent outcomes
> y1.1 <- rbinom(100, 1, 0.5)
> y1.2 <- rbinom(100, 1, 0.5)
> group <- factor(rep(c("A", "B"), 50))
> 
> modely1.1 <- glm(y1.1 ~ group, family = binomial)
> modely1.2 <- glm(y1.2 ~ group, family = binomial)
> 
> mmObj.y1 <- mmod(list(modely1.1, modely1.2), "groupB")
> simult.y1 <- glht(mmObj.y1)
> summary(simult.y1)

	 Simultaneous Tests for General Linear Hypotheses

Linear Hypotheses:
               Estimate Std. Error z value Pr(>|z|)  
modely1.1 == 0   0.8473     0.4186   2.024    0.084 .
modely1.2 == 0   0.2404     0.4008   0.600    0.796  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

> 
> ## Two perfectly correlated outcomes
> y2.1 <- rbinom(100, 1, 0.5)
> y2.2 <- y2.1
> group <- factor(rep(c("A", "B"), 50))
> 
> modely2.1 <- glm(y2.1 ~ group, family = binomial)
> modely2.2 <- glm(y2.2 ~ group, family = binomial)
> 
> mmObj.y2 <- mmod(list(modely2.1, modely2.2), "groupB")
> simult.y2 <- glht(mmObj.y2)
> summary(simult.y2)

	 Simultaneous Tests for General Linear Hypotheses

Linear Hypotheses:
               Estimate Std. Error z value Pr(>|z|)
modely2.1 == 0   0.2412     0.4015   0.601    0.548
modely2.2 == 0   0.2412     0.4015   0.601    0.548
(Adjusted p values reported -- single-step method)

> 
> 
> proc.time()
   user  system elapsed 
   1.09    0.14    1.21