Nothing
library("multcomp")
###<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 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 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))
confint(glht(simult.x12x2.x13))
## 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)
## 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)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.