tests/testthat/test-prettify.R

In hofnerb/papeR: A Toolbox for Writing Pretty Papers and Reports

library("papeR")
context("prettify functions")

set.seed(1234)

################################################################################
## Test computation of CIs when data is part of the call
## (i.e. not only a link is passed but really the data)
################################################################################

test_that("computation of CIs when data is part of the call", {
    model_fit <- function(my_data, model_class) {
        do.call(model_class, list(y ~ x, data = my_data))
    }

    for (model_class in c("lm", "glm")) {
        x <- rnorm(100)
        y <- rnorm(100, mean = 2 * x)
        data <- data.frame(y = y, x = x)

        ## fit model with data argument
        mod <- do.call(model_class, list(y ~ x, data = data))
        psm1 <- prettify(summary(mod))
        rm(data)
        psm1a <- prettify(summary(mod))

        ## fit model without data argument
        mod2 <- do.call(model_class, list(y ~ x))
        psm2 <- prettify(summary(mod2))

        ## fit model in different environment
        mod3 <- model_fit(data.frame(y = y, x = x), model_class)
        psm3 <- prettify(summary(mod3))

        ## change data and compute summary
        x <- rnorm(100)
        y <- rnorm(100, mean = 2 * x)
        data <- data.frame(y = y, x = x)

        psm4 <- prettify(summary(mod))

        expect_equal(psm1, psm1a)
        expect_equal(psm1, psm2)
        expect_equal(psm1, psm3)
        expect_equal(psm1, psm4)
    }
})

################################################################################
## Test computation of CIs when data is *not* part of the call
################################################################################

rm(list = ls())
model_fit2 <- function(my_data, model_class) {
    switch(model_class,
           lm = lm(y ~ x, data = my_data),
           glm = glm(y ~ x, data = my_data),
           coxph = coxph(Surv(y, event) ~ x, data = my_data),
           lme = lme(y ~ x, random = ~ 1 | group, data = my_data),
           lmer = lmer(y ~ x + (1 | group), data = my_data))
}

fit_model <- function(model_class =  c("lm", "glm", "coxph", "lme", "lmer")) {
    x <- rnorm(100)
    y <- rnorm(100, mean = 2 * x)
    data <- data.frame(y = y, x = x)

    if (model_class %in% c("lme", "lmer")) {
        group <- sample(1:2, 100, replace = TRUE)
        data$group <- group
    }

    if (model_class %in% "coxph") {
        event <- as.logical(sample(0:1, 100, replace = TRUE))
        data$event <- event
        data$y <- exp(y)
    }

    ## fit model with data argument
    mod <- switch(model_class,
                  lm = lm(y ~ x, data = data),
                  glm = glm(y ~ x, data = data),
                  coxph = coxph(Surv(y, event) ~ x, data = data),
                  lme = lme(y ~ x, random = ~ 1 | group, data = data),
                  lmer = lmer(y ~ x + (1 | group), data = data))
    psm1 <- prettify(summary(mod))
    data_tmp <- data
    rm(data)
    psm1a <- try(prettify(summary(mod)), silent = TRUE)

    ## fit model without data argument
    mod2 <- switch(model_class,
                   lm = lm(y ~ x),
                   glm = glm(y ~ x),
                   coxph = coxph(Surv(y, event) ~ x),
                   lme = lme(y ~ x, random = ~ 1 | group),
                   lmer = lmer(y ~ x + (1 | group)))
    psm2 <- prettify(summary(mod2))

    ## fit model in different environment
    mod3 <- model_fit2(data_tmp, model_class)
    psm3 <- try(prettify(summary(mod3)), silent = TRUE)

    ## change data and compute summary
    x <- rnorm(100)
    y <- rnorm(100, mean = 2 * x)
    data <- data.frame(y = y, x = x)

    if (model_class %in% c("lme", "lmer")) {
        group <- sample(1:2, 100, replace = TRUE)
        data$group <- group
    }

    if (model_class %in% "coxph") {
        event <- as.logical(sample(0:1, 100, replace = TRUE))
        data$event <- event
        data$y <- exp(y)
    }

    psm4 <- prettify(summary(mod))

    ret <- list(psm1, psm1a, psm2, psm3, psm4)
    names(ret) <- c("standard", "data removed from global environment",
                    "no data argument in call", "local environment",
                    "data in global environment changed")
    return(ret)
}


## check lm interface
test_that("lm interface works", {
    expect_warning(res <- fit_model("lm"))
    expect_equivalent(res[[1]], res[[3]])
    expect_equivalent(res[[1]], res[[4]])
    ## differences in CIs as different data is used
    expect_match(all.equal(res[[1]], res[[5]], check.attributes = FALSE)[1], "lower")
    expect_match(all.equal(res[[1]], res[[5]], check.attributes = FALSE)[2], "upper")
    ## CI dropped. Other values equal
    expect_equivalent(res[[1]][, -(3:4)], res[[2]])
})

## check glm interface
test_that("glm interface works", {
    expect_warning(res <- fit_model("glm"))
    expect_equivalent(res[[1]], res[[3]])
    expect_equivalent(res[[1]], res[[4]])
    ## differences in CIs as different data is used
    expect_match(all.equal(res[[1]], res[[5]], check.attributes = FALSE)[1], "lower")
    expect_match(all.equal(res[[1]], res[[5]], check.attributes = FALSE)[2], "upper")
    ## CI dropped. Other values equal
    expect_equivalent(res[[1]][, -(3:4)], res[[2]])
})

### check lme interface
if (require("nlme", quietly = TRUE)) {
    test_that("nlme interface works", {
        expect_warning(res <- fit_model("lme"))
        expect_equivalent(res[[1]], res[[3]])
        expect_equivalent(res[[1]], res[[4]])
        ## differences in CIs as different data is used
        expect_match(all.equal(res[[1]], res[[5]], check.attributes = FALSE)[1], "lower")
        expect_match(all.equal(res[[1]], res[[5]], check.attributes = FALSE)[2], "upper")
        ## CI dropped. Other values equal
        stopifnot(all.equal(res[[1]][, -(3:4)], res[[2]], check.attributes = FALSE))
    })
}

## check coxph interfaces
if (require("survival", quietly = TRUE)) {
    test_that("survival interface works", {
        expect_warning(res <- fit_model("coxph"))
        expect_equivalent(res[[1]], res[[3]])
        expect_is(res[[2]], "try-error")
        expect_is(res[[4]], "try-error")
        ## differences in CIs as different data is used
        expect_match(all.equal(res[[1]], res[[5]], check.attributes = FALSE)[1], "lower")
        expect_match(all.equal(res[[1]], res[[5]], check.attributes = FALSE)[2], "upper")
    })
}

## check lmer interface
if (require("lme4", quietly = TRUE, warn.conflicts = FALSE)) {
    test_that("lme4 interface works", {
        expect_warning(res <- fit_model("lmer"))
        if (packageDescription("lme4")$Version >= 1) {
            expect_equivalent(res[[1]], res[[3]])
            expect_is(res[[2]], "try-error")
            expect_is(res[[4]], "try-error")
            ## differences in CIs as different data is used
            expect_match(all.equal(res[[1]], res[[5]], check.attributes = FALSE)[1], "lower")
            expect_match(all.equal(res[[1]], res[[5]], check.attributes = FALSE)[2], "upper")
        }
    })
}

################################################################################
## Test OR
################################################################################

test_that("OR are included", {
    x <- rnorm(100)
    y <- rbinom(100, 1, make.link("logit")$linkinv(x * 2))
    data <- data.frame(x, y)
    mod <- glm(y ~ x, data = data, family = binomial)
    ps <- prettify(summary(mod))
    expect_true("Odds Ratio" %in% names(ps))
})

################################################################################
## Test specification of CIs
################################################################################

fit_model <- function(model_class =  c("lm", "glm", "coxph", "lme", "lmer")) {
    x <- rnorm(100)
    y <- rnorm(100, mean = 2 * x)
    data <- data.frame(y = y, x = x)

    if (model_class %in% c("lme", "lmer")) {
        group <- sample(1:2, 100, replace = TRUE)
        data$group <- group
    }

    if (model_class %in% "coxph") {
        event <- as.logical(sample(0:1, 100, replace = TRUE))
        data$event <- event
        data$y <- exp(y)
    }

    ## fit model with data argument
    mod <- switch(model_class,
                  lm = lm(y ~ x, data = data),
                  glm = glm(y ~ x, data = data),
                  coxph = coxph(Surv(y, event) ~ x, data = data),
                  lme = lme(y ~ x, random = ~ 1 | group, data = data),
                  lmer = lmer(y ~ x + (1 | group), data = data))
    return(list(data = data, model = mod))
}

if (require("nlme", quietly = TRUE) && require("lme4", quietly = TRUE) && packageDescription("lme4")$Version >= 1) {
    test_that("CIs can be hand specified", {
        for (model_class in c("lm", "glm", "lme", "lmer")) {
            ## fit model
            RES <- fit_model(model_class)
            mod <- RES$mod
            data <- RES$data
            CI <- matrix(c(1, 2, 3, 4), ncol = 2)
            ps <- prettify(summary(mod), confint = CI)
            expect_equivalent(as.matrix(ps[,3:4]), CI, info = model_class)
        }
    })
}

################################################################################
## Test anova
################################################################################

if (require("nlme", quietly = TRUE) && require("survival", quietly = TRUE) && require("lme4", quietly = TRUE) && packageDescription("lme4")$Version >= 1) {
    test_that("prettify.anova works", {
        for (model_class in c("lm", "glm", "lme", "lmer", "coxph")) {
            ## fit model
            RES <- fit_model(model_class)
            mod <- RES$mod
            data <- RES$data
            if (model_class == "lm") {
                ps_anova <- prettify(anova(mod))
                nc <- ncol(ps_anova)
                expect_match(ps_anova[, nc - 1], "<0.001")
                expect_match(as.character(ps_anova[, nc]),
                             "\\*\\*\\*")
            }
            if (require("car", quietly = TRUE)) {
                ps_anova <- prettify(Anova(mod))
                ps_anova_lbl <- prettify(Anova(mod),
                                         labels = c(x = "Predictor x",
                                                    z = "nonsense"))
                nc <- ncol(ps_anova)
                if (model_class %in% c("lm", "glm", "lmer", "coxph")) {
                    expect_match(ps_anova[, nc - 1],
                                 "<0.001", info = model_class)
                    expect_match(as.character(ps_anova[, nc]),
                                 "\\*\\*\\*", info = model_class)
                    expect_match(ps_anova_lbl[, 1],
                                 "Predictor x", info = model_class)
                }
                if (model_class == "lme") {
                    expect_match(ps_anova[2, nc - 1],
                                 "<0.001", info = model_class)
                    expect_match(as.character(ps_anova[2, nc]),
                                 "\\*\\*\\*", info = model_class)
                    expect_match(ps_anova_lbl[2, 1],
                                 "Predictor x", info = model_class)
                }
            }
        }
    })
}

################################################################################
## Test HR
################################################################################

if (require("survival", quietly = TRUE)) {
    test_that("survival works", {
        RES <- fit_model("coxph")
        mod <- RES$mod
        data <- RES$data
        CI <- matrix(c(1, 2), ncol = 2)

        ps1 <- prettify(summary(mod), HR = TRUE)
        ps2 <- prettify(summary(mod), HR = FALSE)
        ps3 <- prettify(summary(mod), confint = CI, HR = TRUE)
        ps4 <- prettify(summary(mod), confint = CI, HR = FALSE)

        expect_equivalent(ps1[,2], coef(mod))
        expect_equivalent(ps1[,3], exp(coef(mod)))
        expect_equivalent(ps2[,2], coef(mod))
        expect_equivalent(ps1[,4:5], exp(ps2[, 3:4]))
        expect_equivalent(as.matrix(ps3[,4:5]), exp(CI))
        expect_equivalent(as.matrix(ps4[,3:4]), CI)
    })
}


################################################################################
## Extra column for value labels
################################################################################

fit_anova <- function(model_class =  c("lm", "glm", "coxph", "lme", "lmer")) {
    x <- as.factor(sample(c("Gr. A", "Gr. B", "Gr. C"), 100, replace = TRUE))
    y <- rnorm(100, mean = 2 * (x == "Gr. B") - 1 * (x == "Gr. C"))
    data <- data.frame(y = y, x = x)

    if (model_class %in% c("lme", "lmer")) {
        group <- sample(1:2, 100, replace = TRUE)
        data$group <- group
    }

    if (model_class %in% "coxph") {
        event <- as.logical(sample(0:1, 100, replace = TRUE))
        data$event <- event
        data$y <- exp(y)
    }

    ## fit model with data argument
    mod <- switch(model_class,
                  lm = lm(y ~ x, data = data),
                  glm = glm(y ~ x, data = data),
                  coxph = coxph(Surv(y, event) ~ x, data = data),
                  lme = lme(y ~ x, random = ~ 1 | group, data = data),
                  lmer = lmer(y ~ x + (1 | group), data = data))
    return(list(data = data, model = mod))
}


if (require("nlme", quietly = TRUE) && require("survival", quietly = TRUE) && require("lme4", quietly = TRUE) && packageDescription("lme4")$Version >= 1) {
    test_that("prettify.anova works", {
        for (model_class in c("lm", "glm", "lme", "lmer", "coxph")) {
            ## fit model
            RES <- fit_anova(model_class)
            mod <<- RES$mod
            data <- RES$data

            sum1 <- prettify(summary(mod), confint = FALSE)
            sum2 <- prettify(summary(mod), confint = FALSE, extra.column = TRUE)
            expect_equal(colnames(sum2)[2], "Factor Level", info = model_class)
            if (model_class == "coxph") {
                expect_equal(sum1[1, 1], "x: Gr. B", info = model_class)
                expect_equal(sum2[1, 2], "Gr. B", info = model_class)
                expect_true(all(sum2[,1] == c("x", "")))
            } else {
                expect_equal(sum1[2, 1], "x: Gr. B", info = model_class)
                expect_equal(sum2[2, 2], "Gr. B", info = model_class)
                expect_true(all(sum2[2:3,1] == c("x", "")))
            }
        }
    })
}

################################################################################
## Test that trailing zeros are not dropped
################################################################################

if (require("survival", quietly = TRUE)) {
    test_that("trailing zeros are not dropped", {
        data(cancer, package = "survival")
        mod <- coxph(Surv(futime, fustat) ~ age, data = ovarian)
        x <- prettify(summary(mod), digit = 3, env = .GlobalEnv)
        expect_match(x[1,5], "1.30")
    })
}