See the result for the Wind variable. In this case the t-test result is shown in summary(), and the p-value for the Wind variable is non-significant, the corresponding confidence interval is the one obtained by confint(), which uses the t-distribution. If confint.defaut(), which uses the normal distribution, is employed confidence interval does not match the t-test result. The confidence interval does not include the null value of 0 for the Wind variable. Use of confint() is preferred, but the mismatch should decrease as the sample size increases.
Agreement: summary() p AND drop1(test = “F”) p AND confint() CI
confint.default() gives normal approximation, and not good in small samples.
## Load data
data(airquality)
airquality <- airquality[complete.cases(airquality),]
## Slightly change data to get bordeline p-value
airquality <- head(airquality, 30)
airquality[c("21","23","30"),"Wind"] <- c(15,10,2.4)
## Fit linear model
lmRes <- lm(Ozone ~ Solar.R + Wind + Temp, data = airquality)
## t-test (Wald type coef/SE ~ t)
summary(lmRes)
##
## Call:
## lm(formula = Ozone ~ Solar.R + Wind + Temp, data = airquality)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.535 -12.691 -3.127 9.112 59.597
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -40.85676 30.51520 -1.339 0.19219
## Solar.R 0.01050 0.03236 0.324 0.74819
## Wind -1.75108 0.85899 -2.039 0.05179 .
## Temp 1.24203 0.41140 3.019 0.00562 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.68 on 26 degrees of freedom
## Multiple R-squared: 0.4389, Adjusted R-squared: 0.3742
## F-statistic: 6.779 on 3 and 26 DF, p-value: 0.001581
## Partial F-test (numerically equivalent)
drop1(lmRes, test = "F")
## Single term deletions
##
## Model:
## Ozone ~ Solar.R + Wind + Temp
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 8122.9 176.04
## Solar.R 1 32.89 8155.8 174.16 0.1053 0.748189
## Wind 1 1298.31 9421.2 178.49 4.1557 0.051785 .
## Temp 1 2847.59 10970.5 183.05 9.1147 0.005619 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## t-distribution based confidence interval (should match t-test p)
confint(lmRes)
## 2.5 % 97.5 %
## (Intercept) -103.58164096 21.86812582
## Solar.R -0.05602161 0.07702190
## Wind -3.51675087 0.01458832
## Temp 0.39639227 2.08767624
## Normal distribution based CI (not necessary match t-test p)
confint.default(lmRes)
## 2.5 % 97.5 %
## (Intercept) -100.66544241 18.95192727
## Solar.R -0.05292889 0.07392917
## Wind -3.43466155 -0.06750101
## Temp 0.43570776 2.04836075
See the result for the fibrinogen variable. In this case, summary() shows the Wald test using normal approximation. The corresponding confidence interval is obtained by confint.default(), which also uses normal approximation. However, confint() gives the profile likelihood based confidence interval, which is considered more appropriate. This confidence interval does not include the null value of 0 and disagrees with the summary() result. The p-value corresponding to confint() can be obtained by likelihood ratio test, i.e., drop1(model, test = “LRT”).
Agreement 1: summary() p AND confint.default() CI
Agreement 2: drop1(test = “LRT”) p AND confint() CI
1 and 2 are different things (normal approximation-based vs likelihood ratio-based). 2 is preferred.
## Load data
library(HSAUR2)
data(plasma)
## Slightly change data to get bordeline p-value (reduce effect of fibrinogen)
plasma["17","fibrinogen"] <- 3.45
## Fit logistic model
glmRes <- glm(ESR ~ fibrinogen + globulin, family = binomial, data = plasma)
## Z-test (Wald type coef/SE ~ Z)
summary(glmRes)
##
## Call:
## glm(formula = ESR ~ fibrinogen + globulin, family = binomial,
## data = plasma)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9712 -0.6156 -0.3471 -0.2096 2.2604
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -12.9135 5.8320 -2.214 0.0268 *
## fibrinogen 1.9101 0.9786 1.952 0.0510 .
## globulin 0.1593 0.1192 1.337 0.1813
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 30.885 on 31 degrees of freedom
## Residual deviance: 23.050 on 29 degrees of freedom
## AIC: 29.05
##
## Number of Fisher Scoring iterations: 5
## Likelihood ratio test
drop1(glmRes, test = "LRT")
## Single term deletions
##
## Model:
## ESR ~ fibrinogen + globulin
## Df Deviance AIC LRT Pr(>Chi)
## <none> 23.050 29.050
## fibrinogen 1 28.945 32.945 5.8957 0.01518 *
## globulin 1 25.019 29.019 1.9692 0.16053
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Profile likelihood confidence interval (should match likelihood ratio test p)
confint(glmRes)
## 2.5 % 97.5 %
## (Intercept) -27.43604938 -3.1995151
## fibrinogen 0.32994133 4.3063200
## globulin -0.06201605 0.4294997
## Normal distribution based CI (should match Z-test result)
confint.default(glmRes)
## 2.5 % 97.5 %
## (Intercept) -24.343937001 -1.4830603
## fibrinogen -0.007947785 3.8282127
## globulin -0.074269943 0.3928369
## Based on normal distribution
confint.default
## function (object, parm, level = 0.95, ...)
## {
## cf <- coef(object)
## pnames <- names(cf)
## if (missing(parm))
## parm <- pnames
## else if (is.numeric(parm))
## parm <- pnames[parm]
## a <- (1 - level)/2
## a <- c(a, 1 - a)
## pct <- format.perc(a, 3)
## fac <- qnorm(a)
## ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm,
## pct))
## ses <- sqrt(diag(vcov(object)))[parm]
## ci[] <- cf[parm] + ses %o% fac
## ci
## }
## <bytecode: 0x7fcb1d5516b0>
## <environment: namespace:stats>
## Based on t-distribution
confint.lm
## function (object, parm, level = 0.95, ...)
## {
## cf <- coef(object)
## pnames <- names(cf)
## if (missing(parm))
## parm <- pnames
## else if (is.numeric(parm))
## parm <- pnames[parm]
## a <- (1 - level)/2
## a <- c(a, 1 - a)
## fac <- qt(a, object$df.residual)
## pct <- format.perc(a, 3)
## ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm,
## pct))
## ses <- sqrt(diag(vcov(object)))[parm]
## ci[] <- cf[parm] + ses %o% fac
## ci
## }
## <bytecode: 0x7fcb1d845d10>
## <environment: namespace:stats>
## Based on profile
MASS:::confint.glm
## function (object, parm, level = 0.95, trace = FALSE, ...)
## {
## pnames <- names(coef(object))
## if (missing(parm))
## parm <- seq_along(pnames)
## else if (is.character(parm))
## parm <- match(parm, pnames, nomatch = 0L)
## message("Waiting for profiling to be done...")
## utils::flush.console()
## object <- profile(object, which = parm, alpha = (1 - level)/4,
## trace = trace)
## confint(object, parm = parm, level = level, trace = trace,
## ...)
## }
## <bytecode: 0x7fcb1d9af3f8>
## <environment: namespace:MASS>