Is this really the most practical way to return the p-value of a linear model object (lm) in R?
What is the most practical way to extract the global p-value of a linear model lm
? I usually get results from summary
and put statistics on F tests and degrees of freedom in pf
:
set.seed(1)
n <- 10
x <- 1:10
y <- 2*x+rnorm(n)
fit <- lm(y ~ x)
summary(fit) # global p-value: 1.324e-08
fstat <- summary(fit)$fstat
pval <- pf(fstat[1], fstat[2], fstat[3], lower.tail = FALSE)
pval
source to share
Check out the broom package:
library(broom)
set.seed(1)
n <- 10
x <- 1:10
y <- 2*x+rnorm(n)
fit <- lm(y ~ x)
glance(fit)
# r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
# 1 0.9851881 0.9833366 0.8090653 532.1048 1.324022e-08 2 -10.95491 27.90982 28.81758 5.236693 8
glance(fit)$p.value
# [1] 1.324022e-08
tidy(fit)
# term estimate std.error statistic p.value
# 1 (Intercept) -0.1688236 0.55269681 -0.3054542 7.678170e-01
# 2 x 2.0547321 0.08907516 23.0673979 1.324022e-08
source to share
Since you asked for this:
Below is a boneless implementation that skips bells and whistles (and checks) lm
. As a result, it happens faster, but you will use it at your own peril and risk, that is, apply the warnings to help("lm.fit")
. Due to laziness, the code for calculating the F-statistic has been pulled from the source summary.lm
and only slightly modified (so please consider licence()
and citation("stats")
).
fit1 <- lm.fit(cbind(1, x), y)
fstats <- function(obj) {
p <- obj$rank
rdf <- obj$df.residual
r <- obj$residuals
f <- obj$fitted.values
mss <- sum((f - mean(f))^2)
rss <- sum(r^2)
resvar <- rss/rdf
df.int <- 1L #assumes there is always an intercept
fstatistic <- c(value = (mss/(p - df.int))/resvar,
numdf = p - df.int, dendf = rdf)
fstatistic["pval"] <- pf(fstatistic[1L],
fstatistic[2L],
fstatistic[3L], lower.tail = FALSE)
fstatistic
}
fstats(fit1)
# value numdf dendf pval
#5.321048e+02 1.000000e+00 8.000000e+00 1.324022e-08
source to share