retrieving p-values in lm
On Fri, 2005-12-09 at 14:19 +0100, Patrick Kuss wrote:
Dear list, I want to retrieve the p-value of a two-polynomial regression. For a one-polynomial lm I can easily do this with: summary(lm(b~a, data=c)[[4]][[8]]. But how do I find the final p-value in the two-polynomial regression? Under $coefficients I don't find it Any suggestions? Patrick alt <-(2260,2183,2189,1930,2435, 2000,2100,2050,2020,2470, 1700,2310,2090,1560,2060, 1790,1940,2100,2250,2010) H <- c(0.2034,0.1845,0.2053,0.1788,0.2196, 0.2037,0.1655,0.2176,0.1844,0.2033, 0.1393,0.2019,0.1975,0.1490,0.1917, 0.2180,0.2064,0.1943,0.2139,0.1320) X <- data.frame(alt,H) lm.res <- summary(lm(H~alt,data=X)) lm.res p1 <- lm.res[[4]][[8]] p1 lm.res.2 <- summary(lm(H~alt+I(alt^2),data=X)) lm.res.2 str(lm.res.2) # where is p p2 <- lm.res.2[[???]][[????]]
First, you might want to review Chapter 11: Statistical Models in R in An Introduction to R, which is available with your R installation or from the main R web site under Documentation. Specifically, page 53 describes the extractor functions to be used for getting model information. In this case using coef() will extract the model coefficients in both cases:
coef(lm.res)
Estimate Std. Error t value Pr(>|t|) (Intercept) 6.245371e-02 4.713400e-02 1.325024 0.20173833 alt 6.179038e-05 2.261665e-05 2.732074 0.01368545
coef(lm.res.2)
Estimate Std. Error t value Pr(>|t|) (Intercept) -9.433748e-02 3.133627e-01 -0.3010488 0.7670283 alt 2.178857e-04 3.091330e-04 0.7048283 0.4904618 I(alt^2) -3.838002e-08 7.579576e-08 -0.5063610 0.6191070 In both models, the coefficients are present if you review the structure as you have in your code above:
names(lm.res)
[1] "call" "terms" "residuals" "coefficients" [5] "aliased" "sigma" "df" "r.squared" [9] "adj.r.squared" "fstatistic" "cov.unscaled"
names(lm.res.2)
[1] "call" "terms" "residuals" "coefficients" [5] "aliased" "sigma" "df" "r.squared" [9] "adj.r.squared" "fstatistic" "cov.unscaled" So, you can get the term p values by using:
coef(lm.res)[, 4]
(Intercept) alt 0.20173833 0.01368545
coef(lm.res.2)[, 4]
(Intercept) alt I(alt^2) 0.7670283 0.4904618 0.6191070 In terms of the overall model p value, this is actually calculated when you display (print) the model. It is not stored as part of the model object itself. If you review the code for print.summary.lm() using:
stats:::print.summary.lm
...
pf(x$fstatistic[1], x$fstatistic[2], x$fstatistic[3],
lower.tail = FALSE)
...
Where the first argument is the F statistic and the other two are the
degrees of freedom:
lm.res$fstatistic
value numdf dendf 7.464231 1.000000 18.000000
lm.res.2$fstatistic
value numdf dendf 3.706139 2.000000 17.000000 So, in the case of your two models:
x <- lm.res pf(x$fstatistic[1], x$fstatistic[2], x$fstatistic[3],
lower.tail = FALSE)
value
0.01368545
x <- lm.res.2 pf(x$fstatistic[1], x$fstatistic[2], x$fstatistic[3],
lower.tail = FALSE)
value
0.0461472
HTH,
Marc Schwartz