alanocallaghan / contrast

R contrast package

I work on statistical models and their applications for biomedical data. At the moment I am working on Bayesian modelling of eQTLs using allele specific expression in celltype-specific studies of immune-mediated diseases. Previously I worked on quantifying expression heterogeneity in single cell RNAseq data using a Bayesian approach.

I maintain some R packages:

• CRAN packages:
  • bayefdr
  • contrast
• Bioc packages:

  Package	Devel	Release	Issues	Pull requests
  snifter
  densvis
  scater
  BASiCStan

I currently maintain two parts of the Bioconductor book ecosystem:

Book/Chapter	Devel	Release	Issues	Pull requests
OSCA	devel	release
OSCA.multisample	devel	release

I also contribute to:

• CRAN packages:
  • heatmaply
• Bioc:

  Package	Devel	Release	Issues	Pull requests
  BASiCS

Some other stuff:

• plotlyutils: cool plotting ideas with plotly
• shinystats: shiny apps to demonstrate statistical ideas.

contrast's Issues

Functions in formula causes error: "factor(s) not in design"

When using splines in the models the contrast the generateData throws due to the variable not being within the dataset:

library(splines)

lm_fit_3 <- lm(mpg ~ ns(wt, 2) + am, data = mtcars)
summary_3 <- summary(lm_fit_3)
anova_3 <- anova(lm_fit_3)

test_that("function in formula works as expected", {
  a_5 <- list(am = 0, wt = 2)
  b_5 <- list(am = 1, wt = 2)
  c_5 <- contrast(lm_fit_3, a_5, b_5)
  expect_equivalent(
    summary_3$coefficients["am", "Pr(>|t|)"],
    c_5$Pvalue
  )

  # bad inputs
  a_6 <- list(am = "a", wt = 2)
  b_6 <- list(am = "b", wt = 2)
  expect_error(
    contrast(lm_fit_3, a_6, b_6),
    "variable 'am' was fitted with type \"numeric\" but type \"factor\""
  )
})

Throws:

Error in generateData(fit = list(coefficients = c(`(Intercept)` = 34.6323803160965, : 
  factor(s) not in design: diet group

The problem lies in generateData's data check (https://github.com/Alanocallaghan/contrast/blob/dadc7b08cedab2c255bb230143f4c42c9f3b06fa/R/misc.R#L38-L41)

A correct handler would be:

   wh <- sapply(fnam, function(x) any(grepl(x, nam, fixed = TRUE)))
   if (!all(wh))
      stop(paste("factor(s) not in design:",
                 paste(names(factors)[!wh], collapse = " ")))

using fun doesn't work

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(contrast)

nSim <- 1000
a <- 4
b <- 4.1

simData <- rbind(
  data.frame(trt="A", y=rpois(n=nSim, lambda=a)),
  data.frame(trt="B", y=rpois(n=nSim, lambda=b))
)

the_fit <- glm(y ~ -1 + trt, family = poisson(link = log), data = simData)
the_contrast <- contrast::contrast(the_fit, list(trt="B"), list(trt="A"), fcFunc = exp)

print(the_contrast, X = FALSE)
#> glm model parameter contrast
#> 
#>     Contrast      S.E.        Lower      Upper    t   df Pr(>|t|)
#> 1 0.03523938 0.0222804 -0.008455873 0.07893463 1.58 1998   0.1139
print(the_contrast, X = FALSE, fun = exp)
#> glm model parameter contrast
#> Error in fun(w$SE): non-numeric argument to mathematical function

Contrast with interaction terms

From John Sorkin via email

Does the package automatically handle an interaction term, e.g. for the model and contrast statement below,

fit1 <- lm(y~ x + z + x*z)

contrast(x=10, z=10)

Will the contrast statement include in its calculation the x*z interaction? Some testing I have done suggests it does not.

If I am correct (that the contrast statement written above does NOT include the interaction term in its computation), is there a way to specify the interaction in the contrast statement?

I have tried contrast(x=10, z=10, xz=100), contrast(x=10, z=10, xz=100), contrast(x=10, z=10, “xz”=100) none of which run; all produce an error message.

In a model with an offset, will the contrast statement work and include the offset in its computation, for example

fit2 <- lm(y~ x + z + x*z + offset=offset)

contrast(fit2, x=10, y=10)

Some experimentation suggests the contrast function does not account for the offset.

Assuming I am correct, that the offset term is NOT honored in the contrast statement, is there a way to write a contrast statement that will work?