Efficiently compute basic statistical inference from regressions with multiple outcomes

This function simplifies calculating p-values from linear models in which there are many outcome variables, such as in voxel-wise regressions. To perform such an analysis in R, you can concatenate the outcome variables column-wise into an n by p matrix y, where there are n subjects and p outcomes (see Examples). Calling lm(y~x) calculates the coefficients, but statistical inference is not provided. This function provides basic statistical inference efficiently.

bigLMStats(mylm, lambda = 0, includeIntercept = FALSE)

Arguments

mylm: Object of class lm.
lambda: Value of ridge penalty for inverting ill-conditioned matrices.
includeIntercept: Whether or not to include p-values for intercept term in result.

Value

A list containing objects:

fstat: F-statistic of whole model (one value per outcome).
pval.model: p-value of model (one value per outcome).
beta: Values of coefficients (one value per predictor per outcome).
beta.std: Standard error of coefficients.
beta.t: T-statistic of coefficients.
beta.pval: p-value of coefficients.

Author

Kandel BM.

Examples



nsub <- 100
set.seed(1500)
x <- 1:nsub
y <- matrix(c(x + rnorm(nsub), sin(x)), nrow = nsub)
x <- cbind(x, x^2)
y1 <- y[, 1]
y2 <- y[, 2]
lm1 <- lm(y1 ~ x)
lm2 <- lm(y2 ~ x)
mylm <- lm(y ~ x)

myest <- bigLMStats(mylm)
print(paste(
  "R beta estimates for first outcome is", summary(lm1)$coefficients[-1, 1],
  "and for second outcome is", summary(lm2)$coefficients[-1, 1]
))
#> [1] "R beta estimates for first outcome is 0.996800064361502 and for second outcome is -0.000873326597445242"  
#> [2] "R beta estimates for first outcome is 3.96575611836239e-05 and for second outcome is -3.0419087751237e-06"
print(paste("and our estimate is", as.numeric(myest$beta[, 1]), as.numeric(myest$beta[, 2])))
#> [1] "and our estimate is 0.996800064361502 -0.000873326597445242"  
#> [2] "and our estimate is 3.96575611836239e-05 -3.0419087751237e-06"
print(paste(
  "R std error estimate for first outcome is", summary(lm1)$coefficients[-1, 2],
  "and for second outcome is", summary(lm2)$coefficients[-1, 2],
  "and our estimate is", myest$beta.std[, 1], myest$beta.std[, 2]
))
#> [1] "R std error estimate for first outcome is 0.0145452856084367 and for second outcome is 0.0100592155571058 and our estimate is 0.0145452856084366 0.0100592155571058"        
#> [2] "R std error estimate for first outcome is 0.000139527376091849 and for second outcome is 9.64942174398549e-05 and our estimate is 0.000139527376091849 9.64942174398549e-05"
print(paste(
  "R t value estimate for first outcome is", summary(lm1)$coefficients[-1, 3],
  "and for second outcome is", summary(lm2)$coefficients[-1, 3],
  "and our estimate is", myest$beta.t[, 1], myest$beta.t[, 2]
))
#> [1] "R t value estimate for first outcome is 68.5308003703503 and for second outcome is -0.0868185588118081 and our estimate is 68.5308003703503 -0.0868185588118081"  
#> [2] "R t value estimate for first outcome is 0.284227814601185 and for second outcome is -0.0315242597518316 and our estimate is 0.284227814601185 -0.0315242597518316"
print(paste(
  "R pval for first outcome is", summary(lm1)$coefficients[-1, 4],
  "and for second outcome is", summary(lm2)$coefficients[-1, 4],
  "and our estimate is", myest$beta.pval[, 1], myest$beta.pval[, 2]
))
#> [1] "R pval for first outcome is 5.73699385771355e-84 and for second outcome is 0.930994703501836 and our estimate is 5.73699385771339e-84 0.930994703501836"
#> [2] "R pval for first outcome is 0.776841657572342 and for second outcome is 0.974916219349365 and our estimate is 0.776841657572342 0.974916219349365"