Running 'bugfix.R' [3s]
  Comparing 'bugfix.Rout' to 'bugfix.Rout.save' ... OK
  Running 'regtest-Tukey.R' [2s]
  Comparing 'regtest-Tukey.Rout' to 'regtest-Tukey.Rout.save' ... OK
  Running 'regtest-anova.R' [157s]
  Running 'regtest-fixest.R' [2s]
  Running 'regtest-interface-extended.R' [2s]
  Comparing 'regtest-interface-extended.Rout' to 'regtest-interface-extended.Rout.save' ... OK
  Running 'regtest-interface.R' [159s]
  Running 'regtest-lme.R' [166s]
  Running 'regtest-mmm.R' [168s]
  Running 'regtest-survival.R' [151s]
Running the tests in 'tests/regtest-anova.R' failed.
Complete output:
  > 
  > 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
  
  > tol <- sqrt(.Machine$double.eps)
  > set.seed(29081975)
  > 
  > df <- data.frame(y = rnorm(100), 
  +                  x = runif(100),  
  +                  z = runif(100))
  > 
  > ### linear model
  > fam <- gaussian()
  > lm0 <- glm(y ~ 1, data = df, family = fam)
  > lm1 <- glm(y ~ x, data = df, family = fam)
  > lm2 <- glm(y ~ x + z, data = df, family = fam)
  > 
  > gh <- glht(lm2, linfct = c("x = 0", "z = 0"))
  > stopifnot(abs(anova(lm0, lm2, test = "F")[2, 6] - 
  +     summary(gh, test = Ftest())$test$pvalue) < tol)
  > stopifnot(abs(anova(lm0, lm2, test = "Chisq")[2, 5] - 
  +     summary(gh, test = Chisqtest())$test$pvalue) < tol)
  > 
  > gh <- glht(lm2, linfct = "z = 0")
  > stopifnot(abs(anova(lm1, lm2, test = "F")[2, 6] - 
  +     summary(gh, test = Ftest())$test$pvalue) < tol)
  > stopifnot(abs(anova(lm1, lm2, test = "Chisq")[2, 5] - 
  +     summary(gh, test = Chisqtest())$test$pvalue) < tol)
  > 
  > ### logistic regression
  > df$y <- factor(df$y < 0)
  > fam <- binomial()
  > lm0 <- glm(y ~ 1, data = df, family = fam)
  > lm1 <- glm(y ~ x, data = df, family = fam)
  > lm2 <- glm(y ~ x + z, data = df, family = fam)
  > 
  > if (require("lmtest")) {
  + 
  +   gh <- glht(lm2, linfct = c("x = 0", "z = 0"))
  +   stopifnot(abs(waldtest(lm0, lm2, test = "Chisq")[2, 4] - 
  +       summary(gh, test = Chisqtest())$test$pvalue) < tol)
  + 
  +   gh <- glht(lm2, linfct = "z = 0")
  +   stopifnot(abs(waldtest(lm1, lm2, test = "Chisq")[2, 4] -
  +       summary(gh, test = Chisqtest())$test$pvalue) < tol)
  + }
  Loading required package: lmtest
  Loading required package: zoo
  
  Attaching package: 'zoo'
  
  The following objects are masked from 'package:base':
  
      as.Date, as.Date.numeric
  
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.23    0.20    1.36 
Running the tests in 'tests/regtest-interface.R' failed.
Complete output:
  > 
  > 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
  
  > set.seed(290875)
  > 
  > testdata <- data.frame(y = rnorm(21), 
  +                        f1 <- factor(c(rep(c("A", "B", "C"), 7))),
  +                        f2 <- factor(c(rep("D", 10), rep("E", 11))),
  +                        x <- rnorm(21))
  > 
  > # one-way ANOVA
  > coef(amod <- aov(y ~ f1, data = testdata))
  (Intercept)         f1B         f1C 
   -0.4394751   0.5151680   0.6886101 
  > glht(amod, linfct = mcp(f1 = "Dunnett"))
  
  	 General Linear Hypotheses
  
  Multiple Comparisons of Means: Dunnett Contrasts
  
  
  Linear Hypotheses:
             Estimate
  B - A == 0   0.5152
  C - A == 0   0.6886
  
  > 
  > # and a continuous covariable: ANCOVA
  > coef(lmod <- lm(y ~ f1 + x, data = testdata))
   (Intercept)          f1B          f1C            x 
  -0.434528566  0.509444592  0.686181780 -0.009491201 
  > glht(lmod, linfct = mcp(f1 = "Dunnett"))
  
  	 General Linear Hypotheses
  
  Multiple Comparisons of Means: Dunnett Contrasts
  
  
  Linear Hypotheses:
             Estimate
  B - A == 0   0.5094
  C - A == 0   0.6862
  
  > 
  > # ANCOVA with an additional factor as covariable
  > coef(lmod <- lm(y ~ f1 + f2 + x, data = testdata))
  (Intercept)         f1B         f1C         f2E           x 
  -0.40849498  0.51296437  0.69200699 -0.05266965 -0.01613183 
  > glht(lmod, linfct = mcp(f1 = "Dunnett"))
  
  	 General Linear Hypotheses
  
  Multiple Comparisons of Means: Dunnett Contrasts
  
  
  Linear Hypotheses:
             Estimate
  B - A == 0    0.513
  C - A == 0    0.692
  
  > 
  > # and with interaction terms
  > coef(lmod <- lm(y ~ f1 + f2 + f2:f1 + x, data = testdata))
  (Intercept)         f1B         f1C         f2E           x     f1B:f2E 
  -0.44532319  0.70282663  0.65613337  0.05552324 -0.03443721 -0.37862471 
      f1C:f2E 
   0.02753451 
  > glht(lmod, linfct = mcp(f1 = "Dunnett"))
  
  	 General Linear Hypotheses
  
  Multiple Comparisons of Means: Dunnett Contrasts
  
  
  Linear Hypotheses:
             Estimate
  B - A == 0   0.7028
  C - A == 0   0.6561
  
  Warning message:
  In mcp2matrix(model, linfct = linfct) :
    covariate interactions found -- default contrast might be inappropriate
  > 
  > # with contrasts as expressions
  > glht(lmod, linfct = mcp(f1 = c("B - A = 0", "C - A = 0")))
  
  	 General Linear Hypotheses
  
  Multiple Comparisons of Means: User-defined Contrasts
  
  
  Linear Hypotheses:
             Estimate
  B - A == 0   0.7028
  C - A == 0   0.6561
  
  Warning message:
  In mcp2matrix(model, linfct = linfct) :
    covariate interactions found -- default contrast might be inappropriate
  > 
  > tmp <- multcomp:::chrlinfct2matrix(c(l1 = "x1 - x2 = 2", 
  +                                       l2 = "x2 + 3 * x3 = 1"), 
  +                                       paste("x", 1:3, sep = ""))
  > 
  > stopifnot(max(abs(tmp$K - rbind(c(1, -1, 0), c(0, 1, 3)))) < sqrt(.Machine$double.eps))
  > stopifnot(max(abs(tmp$m - c(2, 1))) < sqrt(.Machine$double.eps))
  > 
  > ### coef.survreg and vcov.survreg need special tuning
  > ### thx to Z for pointing this out
  > if (require("survival")) {
  +     smod <- survreg(Surv(futime, fustat) ~ ecog.ps + rx, 
  +                     data = ovarian, dist = 'weibull')
  +     K <- diag(length(coef(smod)))
  +     rownames(K) <- names(coef(smod))
  +     glht(smod, linfct = K)
  + }
  
  	 General Linear Hypotheses
  
  Linear Hypotheses:
                   Estimate
  (Intercept) == 0   6.8967
  ecog.ps == 0      -0.3850
  rx == 0            0.5286
  
  > 
  > ### new `means' comparisons
  > amod <- aov(weight ~ dose + gesttime + number, data = litter)
  > confint(glht(amod, linfct = mcp(dose = "Means")))
  
  	 Simultaneous Confidence Intervals
  
  Multiple Comparisons of Means: Mean Contrasts
  
  
  Fit: aov(formula = weight ~ dose + gesttime + number, data = litter)
  
  Quantile = 2.5558
  95% family-wise confidence level
   
  
  Linear Hypotheses:
           Estimate lwr     upr    
  0 == 0   32.3651  30.0805 34.6498
  5 == 0   29.0127  26.6372 31.3883
  50 == 0  30.0743  27.5239 32.6246
  500 == 0 29.6899  27.1591 32.2207
  
  > 
  > proc.time()
     user  system elapsed 
     1.87    0.26    2.12 
Running the tests in 'tests/regtest-lme.R' failed.
Complete output:
  > 
  > ### mixed effects models
  > ### feature request by John Wilkinson <jnwilks@btinternet.com>
  > ### and Dieter Menne <dieter.menne@menne-biomed.de>
  > 
  > 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
  
  > 
  > lme4OK <- require("lme4")
  Loading required package: lme4
  Loading required package: Matrix
  > if (lme4OK) {
  + 
  +     data("ergoStool", package = "nlme")
  +     K <- glht(aov(effort ~ Type, data = ergoStool), mcp(Type = "Tukey"))$linfct
  + 
  +     stool.lmer <- lmer(effort ~ Type + (1 | Subject),
  +                        data = ergoStool)
  +     glme4 <- glht(stool.lmer, K)
  +     glme41 <- glht(stool.lmer, mcp(Type = "Tukey"))
  +     stopifnot(all.equal(coef(glme4), coef(glme41)))
  +     print(summary(glme41, test = Chisqtest()))
  + 
  +     nlmeOK <- require("nlme")
  +     if (nlmeOK) {
  + 
  +         stool.lme <- lme(effort ~ Type, data = ergoStool,
  +                         random = ~ 1 | Subject)
  +         gnlme <- glht(stool.lme,K)
  +         stopifnot(all.equal(coef(glme4), coef(gnlme)))
  + 
  +         gnlme2 <- glht(stool.lme, linfct = mcp(Type = "Tukey"))
  +         stopifnot(all.equal(coef(glme4), coef(gnlme2)))
  +     }
  + }
  
  	 General Linear Hypotheses
  
  Multiple Comparisons of Means: Tukey Contrasts
  
  
  Linear Hypotheses:
               Estimate
  T2 - T1 == 0   3.8889
  T3 - T1 == 0   2.2222
  T4 - T1 == 0   0.6667
  T3 - T2 == 0  -1.6667
  T4 - T2 == 0  -3.2222
  T4 - T3 == 0  -1.5556
  
  Global Test:
    Chisq DF Pr(>Chisq)
  1 67.07  3  1.812e-14
  Loading required package: nlme
  
  Attaching package: 'nlme'
  
  The following object is masked from 'package:lme4':
  
      lmList
  
  > 
  > proc.time()
     user  system elapsed 
     1.79    0.28    2.06 
Running the tests in 'tests/regtest-mmm.R' failed.
Complete output:
  > 
  > 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.46    0.26    1.70 
Running the tests in 'tests/regtest-survival.R' failed.
Complete output:
  > 
  > library("survival")
  > library("multcomp")
  Loading required package: mvtnorm
  Loading required package: TH.data
  Loading required package: MASS
  
  Attaching package: 'TH.data'
  
  The following object is masked from 'package:MASS':
  
      geyser
  
  > set.seed(29)
  > 
  > bladder1 <- bladder[bladder$enum < 5, ] 
  > bladder1$number <- as.factor(bladder1$number)
  > 
  > cm <- coxph(Surv(stop, event) ~ rx + size + number + strata(enum) + 
  +             cluster(id), data = bladder1)
  > summary(glht(cm, mcp(number = "Dunnett")))
  
  	 Simultaneous Tests for General Linear Hypotheses
  
  Multiple Comparisons of Means: Dunnett Contrasts
  
  
  Fit: coxph(formula = Surv(stop, event) ~ rx + size + number + strata(enum), 
      data = bladder1, cluster = id)
  
  Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)    
  2 - 1 == 0   0.7895     0.5124   1.541    0.494    
  3 - 1 == 0   0.3848     0.4699   0.819    0.941    
  4 - 1 == 0   0.7933     0.6487   1.223    0.727    
  5 - 1 == 0   1.3446     0.6269   2.145    0.160    
  6 - 1 == 0   2.4325     0.3440   7.071   <0.001 ***
  8 - 1 == 0   0.8000     0.3502   2.285    0.115    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  (Adjusted p values reported -- single-step method)
  
  > 
  > wm <- survreg(Surv(stop, event) ~ rx + size + number + strata(enum) + 
  +               cluster(id), data = bladder1)
  > summary(glht(wm, mcp(number = "Dunnett")))
  
  	 Simultaneous Tests for General Linear Hypotheses
  
  Multiple Comparisons of Means: Dunnett Contrasts
  
  
  Fit: survreg(formula = Surv(stop, event) ~ rx + size + number + strata(enum), 
      data = bladder1, cluster = id)
  
  Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)    
  2 - 1 == 0  -0.4609     0.3232  -1.426   0.5374    
  3 - 1 == 0  -0.2325     0.2918  -0.797   0.9328    
  4 - 1 == 0  -0.6017     0.3422  -1.758   0.3184    
  5 - 1 == 0  -1.0458     0.3844  -2.721   0.0338 *  
  6 - 1 == 0  -1.4714     0.2949  -4.989   <0.001 ***
  8 - 1 == 0  -0.3075     0.2944  -1.044   0.8089    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  (Adjusted p values reported -- single-step method)
  
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.43    0.10    1.50