MATH 6627: Understanding Output in Regression
Georges Monette
November 2014
This is an example of a .Rmd file containing R Markdown code. With R Markdown, you can
1. Combine R code, comments and mathematical formulas in one file
2. Produce output in a single .html file
3. Write reproducible analyses – rerunning the same file produces the same output.
History: The original version was written in September 2013. It was revised in September 2014 to incorporate changes in the ‘spida’ package, in particular, the ‘gpanel.fit’ function to draw estimated coefficients and error bands.
WARNING: Changes not finished – only completed up to comment line below. Some of the material has been transfered to Regression in R 2014 xx xx_vX.R
Our theme:
With complex regression models:
Most estimated coefficients are likely to be of little interest, and
the answers to important questions won’t be in your estimated coefficients
This R Markdown script is a discussion about understanding output in regression.
Regression is at the conceptual core of statistics.
Initially, we focus on linear regresion where our goal is to be able to - interpret every parameter estimated in regression outupt, focussing initially on regression coefficients - formulate which parameters and functions of parameters are relevant to a subject-matter question - estimate and test parameters and groups of parameters, by - formulating linear hypotheses for regression coefficients and applying Wald tests, or - reparametrizing the model and testing hypotheses for the reparametrized model, or - formulating and fitting null models followed by a likelihood ratio test to compare the full model with the null model
We want to understand the role and limitations of a number of tools: - estimated regression coefficients - Wald tests for linear hypotheses - ‘anova’ for a likelihood ratio test of two models - ‘Anova’ (in the ‘car’ package) sums of squares: Type I, Type II and Type III and their use and interpretation
Understanding Regression
Consider four approaches to regression: statistical theory, the mathematical formulation of models, computational expression of models in the language of a statistical package and the graphical representation of models.
To master regression you need to know how to go from one representation to another and you need to know how to work within each representation to solve problems.
Statistical
\[ Y = X \beta + \epsilon,\: \epsilon \sim N(0,\sigma^2) \] and all the theory that follows, e.g. \[ Var(\hat{\beta}) = \sigma^2 \left( X'X \right)^{-1} \]
Mathematical
\[ y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \epsilon_i \\ \frac{\partial E(y)}{\partial x} = \beta_1 + 2\,\beta_2 x \]
Computing
library(car)
fit <- lm(income ~ education * type, data = Prestige)
summary(fit)|
| Call:
| lm(formula = income ~ education * type, data = Prestige)
|
| Residuals:
| Min 1Q Median 3Q Max
| -6330.8 -1769.2 -356.8 1166.5 17326.2
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) -1865.0 3682.3 -0.506 0.6137
| education 866.0 436.4 1.984 0.0502 .
| typeprof -3068.4 7191.8 -0.427 0.6706
| typewc 3646.5 9274.0 0.393 0.6951
| education:typeprof 234.0 617.3 0.379 0.7055
| education:typewc -569.2 884.8 -0.643 0.5216
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 3333 on 92 degrees of freedom
| (4 observations deleted due to missingness)
| Multiple R-squared: 0.4105, Adjusted R-squared: 0.3785
| F-statistic: 12.81 on 5 and 92 DF, p-value: 1.856e-09Graphical: data space
library(car)
library(lattice)
xyplot(income ~ education | type, Prestige, type = c('p','r','smooth'))
### Graphical: beta space Figure: Confidence ellipse for two parameters jointly. The blue ellipse has 95% coverage in 2 dimensions and its perpendicular shadows onto the vertical and horizontal axes form Scheffe 95% confidence intervals for testing in a space of dimension 2. The similar shadows of the red ellipse provide ordinary 95% confidence intervals.
library(spida2)
library(latticeExtra)
fit <- lm(income ~ education + women, data = Prestige)
summary(fit)|
| Call:
| lm(formula = income ~ education + women, data = Prestige)
|
| Residuals:
| Min 1Q Median 3Q Max
| -7257.6 -1160.1 -238.6 681.1 16044.3
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) -1491.998 1162.299 -1.284 0.202
| education 944.881 103.731 9.109 9.60e-15 ***
| women -64.056 8.921 -7.180 1.31e-10 ***
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 2839 on 99 degrees of freedom
| Multiple R-squared: 0.5618, Adjusted R-squared: 0.5529
| F-statistic: 63.46 on 2 and 99 DF, p-value: < 2.2e-16plot(rbind(cell(fit),0),type= 'n',
xlab = expression(beta[education]),
ylab = expression(beta[women]))
lines(cell(fit,dfn=2), type = 'l', col = 'blue')
lines(cell(fit,dfn=1), type = 'l', col = 'red')
abline(h=0)
abline(v=0)
points(c(0,0), pch = 18)
Interpreting Regression Coefficients: Smoking and Life Expectancy
With complex models:
- most regression coefficients are of little interest and
- most interesting questions are not answered with the regression coefficients.
Why do we pay attention to regression output? Because it may make some sense for very simple additive models – but even then it is fraught with subtle traps most analysts do not understand.
We will illustrate these with the Smoking and Life Expectancy example.
fit.hiv2 <- lm( LifeExp ~ log(HE) * (smoke + I(smoke^2)) + hiv+special, dd ,
na.action = na.exclude)
summary(fit.hiv2)|
| Call:
| lm(formula = LifeExp ~ log(HE) * (smoke + I(smoke^2)) + hiv +
| special, data = dd, na.action = na.exclude)
|
| Residuals:
| Min 1Q Median 3Q Max
| -9.0373 -2.3005 0.2043 2.0760 9.7344
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 3.283e+01 2.674e+00 12.280 < 2e-16 ***
| log(HE) 6.091e+00 5.024e-01 12.124 < 2e-16 ***
| smoke 3.642e-02 7.520e-03 4.844 3.31e-06 ***
| I(smoke^2) -1.518e-05 3.946e-06 -3.846 0.000181 ***
| hiv -7.351e-01 7.593e-02 -9.681 < 2e-16 ***
| special -1.822e+01 2.137e+00 -8.526 2.11e-14 ***
| log(HE):smoke -4.878e-03 1.155e-03 -4.223 4.30e-05 ***
| log(HE):I(smoke^2) 2.007e-06 5.726e-07 3.504 0.000614 ***
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 3.63 on 141 degrees of freedom
| (45 observations deleted due to missingness)
| Multiple R-squared: 0.8613, Adjusted R-squared: 0.8544
| F-statistic: 125 on 7 and 141 DF, p-value: < 2.2e-16To interpret these coefficient, we consider the mathematical formula for the model:
Letting \(\eta = E(y|HE,Smoke,HIV,Special)\) \[ \begin{aligned} \eta = &\beta_0 +\beta_1 \times ln(HE) \\ & +\beta_2 \times Smoke \\ & +\beta_3 \times Smoke^2 \\ & +\beta_4 \times HIV \\ & +\beta_5 \times Special \\ & +\beta_6 \times ln(HE) \,Smoke \\ & +\beta_7 \times ln(HE) \,Smoke^2 \end{aligned} \]
To understand the interpretation of the coefficients \(\beta_i\), we differentiate \(\eta\) with respect to each of the independent variables:
\[ \begin{aligned} \frac{\partial \eta}{\partial HE} & = \beta_1 \frac{1}{HE}+\beta_6 \frac{Smoke}{HE} + \beta_7 \frac{Smoke^2}{HE} \\ \frac{\partial \eta}{\partial Smoke} & = \beta_2 + 2 \beta_3 Smoke +\beta_6 ln(HE) + 2 \beta_7 Smoke\, ln(HE) \\ \frac{\partial \eta}{\partial HIV} & = \beta_4 \\ \frac{\partial \eta}{\partial Special} & = \beta_5 \\ \frac{\partial^2 \eta}{\partial HE^2} & = \beta_1 \frac{-1}{HE^2}+\beta_6 \frac{-Smoke}{HE^2} + \beta_7 \frac{-Smoke^2}{HE^2} \\ \frac{\partial^2 \eta}{\partial Smoke^2} & = 2 \beta_3 \\ \frac{\partial^2 \eta}{\partial HE \, \partial Smoke} & = \beta_6 \frac{1}{HE} + 2 \beta_7 \frac{Smoke}{HE} \\ \end{aligned} \]
Thus \(\beta_2\) is the partial derivative of \(\eta\) with respect to \(Smoke\) when \(ln(HE) = Smoke = 0\).
When \(ln(HE) = 5\) and \(Smoke = 4\), the partial derivative of \(\eta\) with respect to \(Smoke\) is
\[ \frac{\partial \eta}{\partial Smoke} = \beta_2 + 8 \beta_3 +5 \beta_6 + 40 \beta_7 \] whose estimator is \[ \hat{\beta_2} + 8 \hat{\beta_3} +5 \hat{\beta_6} + 40 \hat{\beta_7} \] which we can express as a linear transformation of the \(\hat{\beta}\) vector. Letting \[L = \left[ {\begin{array}{*{20}{c}}0&0&1&8&0&0&5&{40}\end{array}} \right]\] we have: \[\hat{\phi} = L \hat{\beta} = \left[ {\begin{array}{*{20}{c}}0&0&1&8&0&0&5&{40}\end{array}} \right]\widehat {\left[ {\begin{array}{*{20}{c}}{{\beta _0}}\\{{\beta _1}}\\{{\beta _2}}\\{{\beta _3}}\\{{\beta _4}}\\{{\beta _5}}\\{{\beta _6}}\\{{\beta _7}}\end{array}} \right]}\]
If we wish to simultaneously estimate the ‘effect’ of Smoke and the ‘effect’ of HE given values of HE and Smoke, we can form the \(L\) matrix: \[ L = \left[ {\begin{array}{*{20}{c}} 0 & 0 & 1 & 2\, Smoke & 0 & 0 & ln(HE) & 2\, Smoke \, ln(HE) \\ 0 & 1 & 0 & 0 & 0 & 0 & \frac{Smoke}{HE} &\frac{Smoke^2}{HE} \end{array}} \right] \] and \[ \hat{\phi} = L \hat{\beta} \] In this case \(\hat{\phi}\) is a column vector of length 2.
In both cases, inference about \(\phi\) uses the fact that \[ Var(\hat{\phi}) = L Var(\hat{\beta}) L' \] and \[ Var(\hat{\beta}) = \sigma^2 (X'X)^{-1} \] With a normal linear model in which \[ Y = X \beta + \epsilon, \:\epsilon \sim N(0,\sigma^2 \,I) \] we have that \[ \left( {\hat{\phi} - \phi} \right)'\left( {s^2 L(X'X)^{-1} L'} \right)^{-1} \left( {\hat{\phi} - \phi} \right) \sim h \times F_{h,\nu} \] where \(h\) is the number of rows of \(L\) (assuming that \(L\) is of full row rank) and \(\nu = n - p\) where \(n\) and \(p\) are the number of rows and columns of \(X\) respectively, again assuming that \(X\) is of full column rank.
We can compute these quantities in R from a fitted model.
L <- evalq( rbind(
c( 0,0,1, 2*smoke, 0 ,0 , log(HE), 2*smoke*log(HE)),
c( 0,1,0, 0 , 0, 0, smoke / HE, smoke^2/HE)), envir = list( smoke = 4, HE = exp(5)))
L| [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
| [1,] 0 0 1 8 0 0 5.00000000 40.0000000
| [2,] 0 1 0 0 0 0 0.02695179 0.1078072\(\hat{\beta}\):
coef(fit.hiv2)| (Intercept) log(HE) smoke
| 3.283134e+01 6.091328e+00 3.642387e-02
| I(smoke^2) hiv special
| -1.517663e-05 -7.350980e-01 -1.822331e+01
| log(HE):smoke log(HE):I(smoke^2)
| -4.878054e-03 2.006775e-06\(\hat{\phi} = L \hat{\beta}\):
(phihat <- L %*% coef(fit.hiv2))| [,1]
| [1,] 0.01199246
| [2,] 6.09119673\(s^2 (X'X)^{-1}\):
vcov(fit.hiv2)| (Intercept) log(HE) smoke I(smoke^2)
| (Intercept) 7.148161e+00 -1.298157e+00 -1.437171e-02 5.526193e-06
| log(HE) -1.298157e+00 2.524328e-01 2.379512e-03 -8.798735e-07
| smoke -1.437171e-02 2.379512e-03 5.655121e-05 -2.706223e-08
| I(smoke^2) 5.526193e-06 -8.798735e-07 -2.706223e-08 1.556800e-11
| hiv -9.426054e-03 -2.319091e-03 1.603252e-05 -6.423962e-09
| special 4.335683e-01 -1.122056e-01 -1.123946e-03 4.625775e-07
| log(HE):smoke 2.438561e-03 -4.363636e-04 -8.441015e-06 3.949943e-09
| log(HE):I(smoke^2) -9.000438e-07 1.540858e-07 3.944264e-09 -2.226900e-12
| hiv special log(HE):smoke
| (Intercept) -9.426054e-03 4.335683e-01 2.438561e-03
| log(HE) -2.319091e-03 -1.122056e-01 -4.363636e-04
| smoke 1.603252e-05 -1.123946e-03 -8.441015e-06
| I(smoke^2) -6.423962e-09 4.625775e-07 3.949943e-09
| hiv 5.765616e-03 -1.218220e-03 2.000199e-06
| special -1.218220e-03 4.568101e+00 2.166940e-04
| log(HE):smoke 2.000199e-06 2.166940e-04 1.334160e-06
| log(HE):I(smoke^2) -2.425728e-10 -8.002814e-08 -6.010943e-10
| log(HE):I(smoke^2)
| (Intercept) -9.000438e-07
| log(HE) 1.540858e-07
| smoke 3.944264e-09
| I(smoke^2) -2.226900e-12
| hiv -2.425728e-10
| special -8.002814e-08
| log(HE):smoke -6.010943e-10
| log(HE):I(smoke^2) 3.279108e-13\(\hat{Var}(\hat{\phi}) = L (s^2 (X'X)^{-1}) L'\):
(Vphihat <- L %*% vcov(fit.hiv2) %*% t(L))| [,1] [,2]
| [1,] 5.453266e-06 0.0001967712
| [2,] 1.967712e-04 0.2524093523To test the hypothesis that \(\phi = 0\), we have
\(F = \hat{\phi}'\left({ \hat{Var}(\hat{\phi})}\right) ^{-1}\hat{\phi}/h\)
(Ftest <- (t(phihat) %*% solve(Vphihat) %*% phihat)/2)| [,1]
| [1,] 78.447591-pf(Ftest,2, fit.hiv2$df.residual)| [,1]
| [1,] 0pf(Ftest,2, fit.hiv2$df.residual, lower.tail = FALSE)| [,1]
| [1,] 1.254508e-23Functions to test linear hypotheses
The functions ‘lht’ in the ‘car’ package and ‘wald’ in the ‘spida’ package can be used to test General Linear Hypotheses.
require(car)
lht(fit.hiv2,L)| Linear hypothesis test
|
| Hypothesis:
| smoke + 8 I(smoke^2) + 5 log(HE):smoke + 40 log(HE):I(smoke^2) = 0
| log(HE) + 0.0269517879963419 log(HE):smoke + 0.107807151985367 log(HE):I(smoke^2) = 0
|
| Model 1: restricted model
| Model 2: LifeExp ~ log(HE) * (smoke + I(smoke^2)) + hiv + special
|
| Res.Df RSS Df Sum of Sq F Pr(>F)
| 1 143 3925.4
| 2 141 1858.0 2 2067.4 78.448 < 2.2e-16 ***
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1wald(fit.hiv2, L)| numDF denDF F.value p.value
| 1 2 141 78.44759 <.00001
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| 1 0.011992 0.002335 141 5.135465 <.00001 0.007376 0.016609
| 2 6.091197 0.502404 141 12.124111 <.00001 5.097979 7.084414The ‘lht’ function can take a right-hand side to test hypotheses of the form \(H_0 : \phi = \phi_0\). ‘wald’ can only test \(H_0 : \phi = 0\). The ‘wald’ function can handle \(L\) matrices that are row rank deficient. For example, in some uses of ‘wald’, the ‘L’ matrix is the whole design matrix.
The second argument of the ‘wald’ matrix can be a regular expression that is matched against the names of terms in the model. All terms matched by the regular expression are simultaneously tested. Thus one can test the ‘overall’ siginificance of an independent variable by testing whether all terms containing that variable are equal to 0. One can also use this approach to test higher-order interactions.
wald(fit.hiv2, "smoke") # is there statistical evidence that 'smoke' improves prediction?| numDF denDF F.value p.value
| smoke 4 141 8.024538 1e-05
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| smoke 0.036424 0.007520 141 4.843564 <.00001 0.021557
| I(smoke^2) -0.000015 0.000004 141 -3.846440 0.00018 -0.000023
| log(HE):smoke -0.004878 0.001155 141 -4.223209 0.00004 -0.007162
| log(HE):I(smoke^2) 0.000002 0.000001 141 3.504458 0.00061 0.000001
|
| Coefficients Upper 0.95
| smoke 0.051290
| I(smoke^2) -0.000007
| log(HE):smoke -0.002595
| log(HE):I(smoke^2) 0.000003wald(fit.hiv2, "HE") # ditto for HE| numDF denDF F.value p.value
| HE 3 141 87.10337 <.00001
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| log(HE) 6.091328 0.502427 141 12.123808 <.00001 5.098064
| log(HE):smoke -0.004878 0.001155 141 -4.223209 0.00004 -0.007162
| log(HE):I(smoke^2) 0.000002 0.000001 141 3.504458 0.00061 0.000001
|
| Coefficients Upper 0.95
| log(HE) 7.084592
| log(HE):smoke -0.002595
| log(HE):I(smoke^2) 0.000003wald(fit.hiv2, ":") # ditto for interactions?| numDF denDF F.value p.value
| : 2 141 9.2372 0.00017
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| log(HE):smoke -0.004878 0.001155 141 -4.223209 0.00004 -0.007162
| log(HE):I(smoke^2) 0.000002 0.000001 141 3.504458 0.00061 0.000001
|
| Coefficients Upper 0.95
| log(HE):smoke -0.002595
| log(HE):I(smoke^2) 0.000003wald(fit.hiv2, "2)" ) # ditto for quadratic terms?| numDF denDF F.value p.value
| 2) 2 141 8.835643 0.00024
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| I(smoke^2) -1.5e-05 4e-06 141 -3.846440 0.00018 -2.3e-05
| log(HE):I(smoke^2) 2.0e-06 1e-06 141 3.504458 0.00061 1.0e-06
|
| Coefficients Upper 0.95
| I(smoke^2) -7e-06
| log(HE):I(smoke^2) 3e-06There are many strategies for potentially simplifying large models. One is to attack higher-order interactions and simplify the model by dropping groups of interactions that are not significant. Results may depend on the precise strategy. Another approach is to drop all terms for selected independent variables if they are not sufficiently significant in an overall test. The end result will typically very much depend on the original approach selected. The choice of approach should be guided by many factors: which null hypotheses are likely to be reasonable, the interpretive value of having a simple additive model versus the added validity of estimating conditional effects that are not averaged over levels of variables that may be important, etc. There’s a good discussion of these problems in Snijders & Bosker.
Estimating effects over a grid
In a model with interactions and non-linear functions of some independent variables, it is often necessary to characterize how effects (partial derivatives) and inferences about effects vary over a range of predictors.
The ‘effects’ package by John Fox can help with this task. It can also be done in a more laborious but possibly more flexible way with the ‘Lfx’ function in the ‘spida’ package. The ‘Lfx’ function generates an expression which can then be edited to generate large \(L\) matrices. The result of the wald test applied to this \(L\) matrix can be transformed into a data frame for plotting.
Lfx(fit.hiv2)| list( 1,
| 1 * M(log(HE)),
| 1 * smoke,
| 1 * M(I(smoke^2)),
| 1 * hiv,
| 1 * special,
| 1 * M(log(HE)) * smoke,
| 1 * M(log(HE)) * M(I(smoke^2))
| )The expression generated by ‘Lfx’ can be edited to generate desired effects. The result is then fed back to ‘Lfx’ along with a data frame on which to evaluate the edited expression. Note that the ‘M’ functions preserved the shape of multi-term blocks in the design matrix so that multiplying them by 0 is a way of generating a block of 0s of the right dimension. In the following, we edit the expression to estimate the effect of smoking by differentiating with respect to ‘smoke’:
pred <- expand.grid(HE = c(50,150,500, 1000, 1500, 5000), smoke = seq(10,2000,20), hiv = 0, special = 0)
head(pred) # first 6 lines of 'pred'| HE smoke hiv special
| 1 50 10 0 0
| 2 150 10 0 0
| 3 500 10 0 0
| 4 1000 10 0 0
| 5 1500 10 0 0
| 6 5000 10 0 0L <- Lfx(fit.hiv2,
list( 0,
0 * M(log(HE)),
1 ,
1 * M(I(2*smoke)),
0 * hiv,
0 * special,
1 * M(log(HE)) * 1,
1 * M(log(HE)) * M(I(2*smoke))
), pred)
dim(L)| [1] 600 8head(L)| (Intercept) log(HE) smoke I(smoke^2) hiv special log(HE):smoke
| 1 0 0 1 20 0 0 3.912023
| 2 0 0 1 20 0 0 5.010635
| 3 0 0 1 20 0 0 6.214608
| 4 0 0 1 20 0 0 6.907755
| 5 0 0 1 20 0 0 7.313220
| 6 0 0 1 20 0 0 8.517193
| log(HE):I(smoke^2)
| 1 78.24046
| 2 100.21271
| 3 124.29216
| 4 138.15511
| 5 146.26441
| 6 170.34386ww <- wald(fit.hiv2, L)| Warning in wald(fit.hiv2, L): Poorly conditioned L matrix, calculated numDF
| may be incorrectww <- as.data.frame(ww)
head(ww)| coef se U2 L2 HE smoke hiv
| 1 0.0171942856 0.003272967 0.023740220 0.010648352 50 10 0
| 2 0.0118792893 0.002313918 0.016507124 0.007251454 150 10 0
| 3 0.0060545675 0.001764820 0.009584207 0.002524927 500 10 0
| 4 0.0027011781 0.001883655 0.006468487 -0.001066131 1000 10 0
| 5 0.0007395711 0.002094140 0.004927851 -0.003448709 1500 10 0
| 6 -0.0050851507 0.003067648 0.001050145 -0.011220447 5000 10 0
| special
| 1 0
| 2 0
| 3 0
| 4 0
| 5 0
| 6 0xyplot( coef ~ smoke, ww, groups = HE, auto.key = list(columns = 3, lines = T, points = F),type = 'l')
With labels that are more informative:
Figure 1: Change in Life Expectancy associated with a increase in cigarette consumption of 1 cigarette per day per capita for different levels of health expenditures per capita per year (US$). 
| (Intercept) log(HE) smoke I(smoke^2) hiv special log(HE):smoke
| 1 0 0 1 8 0 0 5
| log(HE):I(smoke^2)
| 1 40
| attr(,"data")
| smoke HE hiv special
| 1 4 148.4132 0 0| numDF denDF F.value p.value
| At smoke = 4, HE = 148.4 1 141 26.373 <.00001
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| 1 0.011992 0.002335 141 5.135465 <.00001 0.007376 0.016609Exercises
- Carry out a similar process to estimate the ‘effect’ of health expenditures per capita.
- Study the relative contribution of private versus public health expenditures on life expectancy.
- Explore the ‘effects’ package and compare its functionality with ‘Lfx’
Wald tests vs Likelihood Ratio Tests (LRT)
Let’s consider a test for the need for a quadratic term in ‘smoke’. There are two terms in the model that contain the quadratic term and a test to remove it involves more than one parameter. We need a test of \[ H_0: \beta_3 = \beta_7 = 0 \] We cannot simply test each hypothesis \(H_0: \beta_3 = 0\) and \(H_0: \beta_7 = 0\) separately. We will see many examples where individual hypotheses are not significant, yet a joint hypothesis is highly significant. This is not a example of this phenomenon since the p-value for each hypothesis is small. Nevertheless, a test of a joint hypothesis needs to be carried out correctly. We consider two ways: a Wald test and a Likelihood Ratio Test executed with the ‘anova’ function, a clear misnomer.
## Wald test using indices of coefficients
wald(fit.hiv2, c(4,8))| numDF denDF F.value p.value
| 1 2 141 8.835643 0.00024
|
| Estimate Std.Error DF t-value p-value Lower 0.95
| I(smoke^2) -1.5e-05 4e-06 141 -3.846440 0.00018 -2.3e-05
| log(HE):I(smoke^2) 2.0e-06 1e-06 141 3.504458 0.00061 1.0e-06
|
| Upper 0.95
| I(smoke^2) -7e-06
| log(HE):I(smoke^2) 3e-06wald(fit.hiv2, list("Quadratic in smoke" =c(4,8)))| numDF denDF F.value p.value
| Quadratic in smoke 2 141 8.835643 0.00024
|
| Estimate Std.Error DF t-value p-value Lower 0.95
| I(smoke^2) -1.5e-05 4e-06 141 -3.846440 0.00018 -2.3e-05
| log(HE):I(smoke^2) 2.0e-06 1e-06 141 3.504458 0.00061 1.0e-06
|
| Upper 0.95
| I(smoke^2) -7e-06
| log(HE):I(smoke^2) 3e-06## Wald test using regular expression
wald(fit.hiv2, "2")| numDF denDF F.value p.value
| 2 2 141 8.835643 0.00024
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| I(smoke^2) -1.5e-05 4e-06 141 -3.846440 0.00018 -2.3e-05
| log(HE):I(smoke^2) 2.0e-06 1e-06 141 3.504458 0.00061 1.0e-06
|
| Coefficients Upper 0.95
| I(smoke^2) -7e-06
| log(HE):I(smoke^2) 3e-06## Likelihood ratio test
# We need to fit the 'null' model
fit0 <- update(fit.hiv2, .~ log(HE)*smoke + hiv + special)
summary(fit0)|
| Call:
| lm(formula = LifeExp ~ log(HE) + smoke + hiv + special + log(HE):smoke,
| data = dd, na.action = na.exclude)
|
| Residuals:
| Min 1Q Median 3Q Max
| -9.5894 -2.2654 0.0006 2.4404 9.6020
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 3.630e+01 2.082e+00 17.438 < 2e-16 ***
| log(HE) 5.728e+00 3.613e-01 15.854 < 2e-16 ***
| smoke 1.132e-02 3.149e-03 3.596 0.000444 ***
| hiv -7.617e-01 7.896e-02 -9.647 < 2e-16 ***
| special -1.802e+01 2.243e+00 -8.032 3.23e-13 ***
| log(HE):smoke -1.660e-03 4.628e-04 -3.586 0.000459 ***
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 3.824 on 143 degrees of freedom
| (45 observations deleted due to missingness)
| Multiple R-squared: 0.8439, Adjusted R-squared: 0.8384
| F-statistic: 154.6 on 5 and 143 DF, p-value: < 2.2e-16# Then compare the null model with the 'full' model:
# By default, 'anova' uses an F distribution for the LRT
# taking advantage of the linear gaussian model
anova(fit0, fit.hiv2) | Analysis of Variance Table
|
| Model 1: LifeExp ~ log(HE) + smoke + hiv + special + log(HE):smoke
| Model 2: LifeExp ~ log(HE) * (smoke + I(smoke^2)) + hiv + special
| Res.Df RSS Df Sum of Sq F Pr(>F)
| 1 143 2090.8
| 2 141 1858.0 2 232.86 8.8356 0.0002426 ***
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Using the general asymptotic distribution, chi-square, for the LRT gives a slightly
# different but very close result
anova(fit0, fit.hiv2, test="LRT")| Analysis of Variance Table
|
| Model 1: LifeExp ~ log(HE) + smoke + hiv + special + log(HE):smoke
| Model 2: LifeExp ~ log(HE) * (smoke + I(smoke^2)) + hiv + special
| Res.Df RSS Df Sum of Sq Pr(>Chi)
| 1 143 2090.8
| 2 141 1858.0 2 232.86 0.0001455 ***
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The Wald test and the LRT using the F statistics give identical results. This is the luxury of working with a Gaussian homoskedastic independent linear model. Exercises ———– 1. Explore the pros and cons of Wald tests versus Likelihood Ratio Tests.
Interpreting sequential tests
## Type I: sequential tests
anova(fit.hiv2) # sequential - Type 1 tests| Analysis of Variance Table
|
| Response: LifeExp
| Df Sum Sq Mean Sq F value Pr(>F)
| log(HE) 1 8631.5 8631.5 655.0400 < 2.2e-16 ***
| smoke 1 152.3 152.3 11.5572 0.0008779 ***
| I(smoke^2) 1 399.7 399.7 30.3307 1.675e-07 ***
| hiv 1 1228.0 1228.0 93.1880 < 2.2e-16 ***
| special 1 878.4 878.4 66.6576 1.637e-13 ***
| log(HE):smoke 1 81.6 81.6 6.1932 0.0139881 *
| log(HE):I(smoke^2) 1 161.8 161.8 12.2812 0.0006136 ***
| Residuals 141 1858.0 13.2
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Type 2: Each term added last except for higher-order interactions
require(car)
Anova(fit.hiv2) # Type 2 is default| Anova Table (Type II tests)
|
| Response: LifeExp
| Sum Sq Df F value Pr(>F)
| log(HE) 3199.9 1 242.8357 < 2.2e-16 ***
| smoke 129.5 1 9.8297 0.0020900 **
| I(smoke^2) 71.0 1 5.3901 0.0216845 *
| hiv 1235.0 1 93.7227 < 2.2e-16 ***
| special 957.9 1 72.6974 2.109e-14 ***
| log(HE):smoke 235.0 1 17.8355 4.297e-05 ***
| log(HE):I(smoke^2) 161.8 1 12.2812 0.0006136 ***
| Residuals 1858.0 141
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Type 3: Each term added last
require(car)
Anova(fit.hiv2, type = 3) | Anova Table (Type III tests)
|
| Response: LifeExp
| Sum Sq Df F value Pr(>F)
| (Intercept) 1987.03 1 150.794 < 2.2e-16 ***
| log(HE) 1936.86 1 146.987 < 2.2e-16 ***
| smoke 309.14 1 23.460 3.315e-06 ***
| I(smoke^2) 194.96 1 14.795 0.0001809 ***
| hiv 1235.00 1 93.723 < 2.2e-16 ***
| special 957.94 1 72.697 2.109e-14 ***
| log(HE):smoke 235.02 1 17.835 4.297e-05 ***
| log(HE):I(smoke^2) 161.83 1 12.281 0.0006136 ***
| Residuals 1857.97 141
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Type 3 is identical to regression output except that it uses equivalent F tests
## and a single test for terms with multiple degrees of freedom
summary(fit.hiv2)|
| Call:
| lm(formula = LifeExp ~ log(HE) * (smoke + I(smoke^2)) + hiv +
| special, data = dd, na.action = na.exclude)
|
| Residuals:
| Min 1Q Median 3Q Max
| -9.0373 -2.3005 0.2043 2.0760 9.7344
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 3.283e+01 2.674e+00 12.280 < 2e-16 ***
| log(HE) 6.091e+00 5.024e-01 12.124 < 2e-16 ***
| smoke 3.642e-02 7.520e-03 4.844 3.31e-06 ***
| I(smoke^2) -1.518e-05 3.946e-06 -3.846 0.000181 ***
| hiv -7.351e-01 7.593e-02 -9.681 < 2e-16 ***
| special -1.822e+01 2.137e+00 -8.526 2.11e-14 ***
| log(HE):smoke -4.878e-03 1.155e-03 -4.223 4.30e-05 ***
| log(HE):I(smoke^2) 2.007e-06 5.726e-07 3.504 0.000614 ***
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 3.63 on 141 degrees of freedom
| (45 observations deleted due to missingness)
| Multiple R-squared: 0.8613, Adjusted R-squared: 0.8544
| F-statistic: 125 on 7 and 141 DF, p-value: < 2.2e-16Working with factors
Controlling for WHO regions provides a non-trivial example of the use of factors in regression.
When you create a data frame in R, non-numeric variables are automatically turned into factors. Factors are both a strength of R and a frequent source of annoyance and confusion. See traps and pitfalls with factors.
Let’s create a small data frame to illustrate how factors work:
set.seed(147)
sdf <- data.frame( x = c(1:7,6:10),
g = rep(c('a','b','c'),c(2,5,5)))
sdf| x g
| 1 1 a
| 2 2 a
| 3 3 b
| 4 4 b
| 5 5 b
| 6 6 b
| 7 7 b
| 8 6 c
| 9 7 c
| 10 8 c
| 11 9 c
| 12 10 csdf$y <- with(sdf, x + c(1,0,2)[g]+.5* rnorm(12))
sdf| x g y
| 1 1 a 2.120108
| 2 2 a 2.853852
| 3 3 b 2.722852
| 4 4 b 4.825593
| 5 5 b 4.623128
| 6 6 b 5.732132
| 7 7 b 7.221543
| 8 6 c 7.060536
| 9 7 c 9.465820
| 10 8 c 9.748093
| 11 9 c 11.580972
| 12 10 c 11.661232sdf$g| [1] a a b b b b b c c c c c
| Levels: a b cunclass(sdf$g) # the innards of g | [1] 1 1 2 2 2 2 2 3 3 3 3 3
| attr(,"levels")
| [1] "a" "b" "c"# g is actually a numeric variable consisting of indices into a
# vector of 'levels'.
sfit <- lm( y ~ x + g, sdf, na.action = na.exclude)
summary(sfit)|
| Call:
| lm(formula = y ~ x + g, data = sdf, na.action = na.exclude)
|
| Residuals:
| Min 1Q Median 3Q Max
| -0.7367 -0.3465 -0.1574 0.2736 0.8536
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 0.9074 0.4421 2.052 0.0742 .
| x 1.0530 0.1250 8.421 3.01e-05 ***
| gb -1.1476 0.6449 -1.779 0.1131
| gc 0.5716 0.9408 0.608 0.5603
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 0.5662 on 8 degrees of freedom
| Multiple R-squared: 0.9797, Adjusted R-squared: 0.972
| F-statistic: 128.4 on 3 and 8 DF, p-value: 4.178e-07- Note that there is no term called ‘ga’ although there are 3 levels: ‘a’, ‘b’ and ‘c’
- The ‘missing’ level, ‘a’, is called the reference level
- Each term shows a comparison with the reference level To work out what the coefficients for ‘gb’ and ‘gc’ mean, you need to look at the X matrix:
model.matrix(sfit, na.action = na.exclude)| (Intercept) x gb gc
| 1 1 1 0 0
| 2 1 2 0 0
| 3 1 3 1 0
| 4 1 4 1 0
| 5 1 5 1 0
| 6 1 6 1 0
| 7 1 7 1 0
| 8 1 6 0 1
| 9 1 7 0 1
| 10 1 8 0 1
| 11 1 9 0 1
| 12 1 10 0 1
| attr(,"assign")
| [1] 0 1 2 2
| attr(,"contrasts")
| attr(,"contrasts")$g
| [1] "contr.treatment"model.matrix(~ x + g, sdf)| (Intercept) x gb gc
| 1 1 1 0 0
| 2 1 2 0 0
| 3 1 3 1 0
| 4 1 4 1 0
| 5 1 5 1 0
| 6 1 6 1 0
| 7 1 7 1 0
| 8 1 6 0 1
| 9 1 7 0 1
| 10 1 8 0 1
| 11 1 9 0 1
| 12 1 10 0 1
| attr(,"assign")
| [1] 0 1 2 2
| attr(,"contrasts")
| attr(,"contrasts")$g
| [1] "contr.treatment"model.matrix(~ g, sdf)| (Intercept) gb gc
| 1 1 0 0
| 2 1 0 0
| 3 1 1 0
| 4 1 1 0
| 5 1 1 0
| 6 1 1 0
| 7 1 1 0
| 8 1 0 1
| 9 1 0 1
| 10 1 0 1
| 11 1 0 1
| 12 1 0 1
| attr(,"assign")
| [1] 0 1 1
| attr(,"contrasts")
| attr(,"contrasts")$g
| [1] "contr.treatment"If you work through the model: \[ E(y|x,g) = \beta_0 + \beta_x x + \beta_{gb} gb + \beta_{gc} gc \] where \(gb = 1\) if \(g=b\) and 0 otherwise, and \(gc = 1\) if \(g=c\) and 0 otherwise, you see that
- \(\beta_{gb}\) is the difference between the expected level for group ‘b’ versus the reference group ‘a’ keeping x constant and
- \(\beta_{gc}\) is the same comparison for group ‘c’ compared with the reference group ‘a’.
# Plotting fits within groups and panels using latticeExtra:
# See more elegant but perhaps less flexible approaches in
# the 'car' and in the 'effects' package by John Fox
pred <- expand.grid( x = 0:13, g = levels( sdf$g))
# the values over which we want to see predicted lines
# every combination of x and g
pred <- merge( sdf, pred, all = T) # merge with data
pred$y1 <- predict(sfit, newdata = pred) # the predicted value
pred <- pred [ order(pred$x),] # order so lines won't be interrupted
require(latticeExtra)
require(spida)| Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
| logical.return = TRUE, : there is no package called 'spida'td(cex=2, lty=1:3, # this gives you optional control over line styles, colour, etc.
pch = 16:18,
col = c('red','blue','magenta'),
lwd =2) # from spida
xyplot( y ~ x , pred, groups = g , auto.key = list(columns = 3,lines = T),
y1 = pred$y1, subscripts = T,
sub = "compare adjusted and unadjusted differences between groups") +
glayer( panel.lines( x, y1[subscripts],...,type = 'l')) +
layer( panel.abline( v = c(0,3), col = 'grey')) +
glayer( panel.abline( h = mean(y,na.rm=T),...))
## an alternative ... but
xyplot( y ~ x, sdf, groups = g, type = c('p','r'))
### Exercise: 1. Draw by hand the values of estimated coefficients in the plot. 2. How would you estimate the differences between horizontal lines? Does ‘g’ matter? ——————
summary(sfit) # p-values not significant|
| Call:
| lm(formula = y ~ x + g, data = sdf, na.action = na.exclude)
|
| Residuals:
| Min 1Q Median 3Q Max
| -0.7367 -0.3465 -0.1574 0.2736 0.8536
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 0.9074 0.4421 2.052 0.0742 .
| x 1.0530 0.1250 8.421 3.01e-05 ***
| gb -1.1476 0.6449 -1.779 0.1131
| gc 0.5716 0.9408 0.608 0.5603
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 0.5662 on 8 degrees of freedom
| Multiple R-squared: 0.9797, Adjusted R-squared: 0.972
| F-statistic: 128.4 on 3 and 8 DF, p-value: 4.178e-07## But
wald(sfit, "g") # simultaneous test that both are 0: different answer!| numDF denDF F.value p.value
| g 2 8 9.711573 0.00724
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| gb -1.147573 0.644944 8 -1.779337 0.11306 -2.634817
| gc 0.571585 0.940783 8 0.607563 0.56032 -1.597865
|
| Coefficients Upper 0.95
| gb 0.339671
| gc 2.741035Using the GLH (General Linear Hypothesis)
Lmu.6 <-list( "at x = 6" = rbind(
'g = a' = c(1,6,0,0),
'g = b' = c(1,6,1,0),
'g = c' = c(1,6,0,1)))
Lmu.6| $`at x = 6`
| [,1] [,2] [,3] [,4]
| g = a 1 6 0 0
| g = b 1 6 1 0
| g = c 1 6 0 1wald(sfit, Lmu.6)| numDF denDF F.value p.value
| at x = 6 3 8 613.0644 <.00001
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| g = a 7.225664 0.690607 8 10.46277 1e-05 5.633120 8.818207
| g = b 6.078091 0.282401 8 21.52290 <.00001 5.426872 6.729309
| g = c 7.797249 0.355897 8 21.90875 <.00001 6.976550 8.617948Ldiff <- rbind(
'b - a' = c(0,0,1,0),
'c - a' = c(0,0,0,1),
'c - b' = c(0,0,-1,1))
Ldiff| [,1] [,2] [,3] [,4]
| b - a 0 0 1 0
| c - a 0 0 0 1
| c - b 0 0 -1 1wald(sfit,Ldiff)| numDF denDF F.value p.value
| 1 2 8 9.711573 0.00724
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| b - a -1.147573 0.644944 8 -1.779337 0.11306 -2.634817 0.339671
| c - a 0.571585 0.940783 8 0.607563 0.56032 -1.597865 2.741035
| c - b 1.719159 0.518616 8 3.314900 0.01062 0.523229 2.915088wald(sfit, 'g')| numDF denDF F.value p.value
| g 2 8 9.711573 0.00724
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| gb -1.147573 0.644944 8 -1.779337 0.11306 -2.634817
| gc 0.571585 0.940783 8 0.607563 0.56032 -1.597865
|
| Coefficients Upper 0.95
| gb 0.339671
| gc 2.741035This illustrated the crucial point that separate tests of \[ H_0: \beta_1 = 0 \] and \[ H_0: \beta_2 = 0 \] can yield very different ‘conclusions’ that a test of the joint hypothesis: \[ H_0: \beta_1 = \beta_2 = 0 \] Later, we will see how the relationship between confidence ellipses(oids) and tests makes this clear. ### Reparametrization to answer different questions
sdf$g2 <- relevel(sdf$g, 'b') # makes 'b' the reference level
fitr <- lm( y ~ x + g2, sdf)
summary(fitr)|
| Call:
| lm(formula = y ~ x + g2, data = sdf)
|
| Residuals:
| Min 1Q Median 3Q Max
| -0.7367 -0.3465 -0.1574 0.2736 0.8536
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) -0.2402 0.6746 -0.356 0.7310
| x 1.0530 0.1250 8.421 3.01e-05 ***
| g2a 1.1476 0.6449 1.779 0.1131
| g2c 1.7192 0.5186 3.315 0.0106 *
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 0.5662 on 8 degrees of freedom
| Multiple R-squared: 0.9797, Adjusted R-squared: 0.972
| F-statistic: 128.4 on 3 and 8 DF, p-value: 4.178e-07fitr2 <- lm( y ~ x + g2 -1 , sdf) # dropping the intercept
summary(fitr2)|
| Call:
| lm(formula = y ~ x + g2 - 1, data = sdf)
|
| Residuals:
| Min 1Q Median 3Q Max
| -0.7367 -0.3465 -0.1574 0.2736 0.8536
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| x 1.0530 0.1250 8.421 3.01e-05 ***
| g2b -0.2402 0.6746 -0.356 0.7310
| g2a 0.9074 0.4421 2.052 0.0742 .
| g2c 1.4790 1.0319 1.433 0.1897
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 0.5662 on 8 degrees of freedom
| Multiple R-squared: 0.9961, Adjusted R-squared: 0.9941
| F-statistic: 508.3 on 4 and 8 DF, p-value: 1.176e-09fitr3 <- lm( y ~ I(x-6) + g2 -1 , sdf) # recentering x
summary(fitr3) # compare with earlier|
| Call:
| lm(formula = y ~ I(x - 6) + g2 - 1, data = sdf)
|
| Residuals:
| Min 1Q Median 3Q Max
| -0.7367 -0.3465 -0.1574 0.2736 0.8536
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| I(x - 6) 1.0530 0.1250 8.421 3.01e-05 ***
| g2b 6.0781 0.2824 21.523 2.29e-08 ***
| g2a 7.2257 0.6906 10.463 6.05e-06 ***
| g2c 7.7972 0.3559 21.909 1.99e-08 ***
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 0.5662 on 8 degrees of freedom
| Multiple R-squared: 0.9961, Adjusted R-squared: 0.9941
| F-statistic: 508.3 on 4 and 8 DF, p-value: 1.176e-09wald(sfit, Lmu.6)| numDF denDF F.value p.value
| at x = 6 3 8 613.0644 <.00001
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| g = a 7.225664 0.690607 8 10.46277 1e-05 5.633120 8.818207
| g = b 6.078091 0.282401 8 21.52290 <.00001 5.426872 6.729309
| g = c 7.797249 0.355897 8 21.90875 <.00001 6.976550 8.617948Equivalent models
What makes the last three models equivalent?
summary(lm( model.matrix(fitr) ~ model.matrix(sfit)-1)) # note the Resid. SE| Warning in summary.lm(object, ...): essentially perfect fit: summary may be
| unreliable
| Warning in summary.lm(object, ...): essentially perfect fit: summary may be
| unreliable
| Warning in summary.lm(object, ...): essentially perfect fit: summary may be
| unreliable
| Warning in summary.lm(object, ...): essentially perfect fit: summary may be
| unreliable| Response (Intercept) :
|
| Call:
| lm(formula = `(Intercept)` ~ model.matrix(sfit) - 1)
|
| Residuals:
| Min 1Q Median 3Q Max
| -7.504e-16 -4.690e-17 0.000e+00 4.690e-17 7.504e-16
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| model.matrix(sfit)(Intercept) 1.000e+00 2.966e-16 3.371e+15 <2e-16
| model.matrix(sfit)x -3.752e-17 8.390e-17 -4.470e-01 0.667
| model.matrix(sfit)gb -6.379e-16 4.327e-16 -1.474e+00 0.179
| model.matrix(sfit)gc -5.253e-16 6.312e-16 -8.320e-01 0.429
|
| model.matrix(sfit)(Intercept) ***
| model.matrix(sfit)x
| model.matrix(sfit)gb
| model.matrix(sfit)gc
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 3.799e-16 on 8 degrees of freedom
| Multiple R-squared: 1, Adjusted R-squared: 1
| F-statistic: 2.079e+31 on 4 and 8 DF, p-value: < 2.2e-16
|
|
| Response x :
|
| Call:
| lm(formula = x ~ model.matrix(sfit) - 1)
|
| Residuals:
| Min 1Q Median 3Q Max
| -8.742e-16 -8.917e-17 -2.706e-17 1.417e-16 8.742e-16
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| model.matrix(sfit)(Intercept) -1.026e-15 3.623e-16 -2.831e+00 0.0221
| model.matrix(sfit)x 1.000e+00 1.025e-16 9.759e+15 <2e-16
| model.matrix(sfit)gb -8.052e-16 5.285e-16 -1.523e+00 0.1661
| model.matrix(sfit)gc -9.918e-16 7.709e-16 -1.287e+00 0.2342
|
| model.matrix(sfit)(Intercept) *
| model.matrix(sfit)x ***
| model.matrix(sfit)gb
| model.matrix(sfit)gc
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 4.64e-16 on 8 degrees of freedom
| Multiple R-squared: 1, Adjusted R-squared: 1
| F-statistic: 5.459e+32 on 4 and 8 DF, p-value: < 2.2e-16
|
|
| Response g2a :
|
| Call:
| lm(formula = g2a ~ model.matrix(sfit) - 1)
|
| Residuals:
| Min 1Q Median 3Q Max
| -5.639e-17 -2.701e-17 -8.327e-18 3.385e-17 6.024e-17
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| model.matrix(sfit)(Intercept) 1.000e+00 3.675e-17 2.721e+16 < 2e-16
| model.matrix(sfit)x -4.223e-17 1.039e-17 -4.063e+00 0.00362
| model.matrix(sfit)gb -1.000e+00 5.361e-17 -1.865e+16 < 2e-16
| model.matrix(sfit)gc -1.000e+00 7.820e-17 -1.279e+16 < 2e-16
|
| model.matrix(sfit)(Intercept) ***
| model.matrix(sfit)x **
| model.matrix(sfit)gb ***
| model.matrix(sfit)gc ***
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 4.706e-17 on 8 degrees of freedom
| Multiple R-squared: 1, Adjusted R-squared: 1
| F-statistic: 2.258e+32 on 4 and 8 DF, p-value: < 2.2e-16
|
|
| Response g2c :
|
| Call:
| lm(formula = g2c ~ model.matrix(sfit) - 1)
|
| Residuals:
| Min 1Q Median 3Q Max
| -1.054e-17 -6.564e-18 -8.950e-20 4.320e-18 1.981e-17
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| model.matrix(sfit)(Intercept) 0.000e+00 7.899e-18 0.00e+00 1
| model.matrix(sfit)x 0.000e+00 2.234e-18 0.00e+00 1
| model.matrix(sfit)gb 0.000e+00 1.152e-17 0.00e+00 1
| model.matrix(sfit)gc 1.000e+00 1.681e-17 5.95e+16 <2e-16 ***
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 1.012e-17 on 8 degrees of freedom
| Multiple R-squared: 1, Adjusted R-squared: 1
| F-statistic: 1.222e+34 on 4 and 8 DF, p-value: < 2.2e-16Each model matrix spans exactly the same linear space. Thus their columns are just different bases for the same space and the \(\beta\)s for one model are just a linear transformation of the \(\betas\)s for the other model. ### Exercises:
## What do the coefficients estimate in each of the following?:
summary( fit1 <- lm( y ~ g - 1, sdf)) # an example where '=' would not work|
| Call:
| lm(formula = y ~ g - 1, data = sdf)
|
| Residuals:
| Min 1Q Median 3Q Max
| -2.8428 -0.4108 -0.1774 0.9497 2.1965
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| ga 2.4870 1.1855 2.098 0.0653 .
| gb 5.0250 0.7498 6.702 8.83e-05 ***
| gc 9.9033 0.7498 13.209 3.39e-07 ***
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 1.677 on 9 degrees of freedom
| Multiple R-squared: 0.9613, Adjusted R-squared: 0.9485
| F-statistic: 74.59 on 3 and 9 DF, p-value: 1.118e-06model.matrix(fit1) | ga gb gc
| 1 1 0 0
| 2 1 0 0
| 3 0 1 0
| 4 0 1 0
| 5 0 1 0
| 6 0 1 0
| 7 0 1 0
| 8 0 0 1
| 9 0 0 1
| 10 0 0 1
| 11 0 0 1
| 12 0 0 1
| attr(,"assign")
| [1] 1 1 1
| attr(,"contrasts")
| attr(,"contrasts")$g
| [1] "contr.treatment"summary( fit2 <- lm( y ~ x + g - 1, sdf)) # an example where '=' would not work|
| Call:
| lm(formula = y ~ x + g - 1, data = sdf)
|
| Residuals:
| Min 1Q Median 3Q Max
| -0.7367 -0.3465 -0.1574 0.2736 0.8536
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| x 1.0530 0.1250 8.421 3.01e-05 ***
| ga 0.9074 0.4421 2.052 0.0742 .
| gb -0.2402 0.6746 -0.356 0.7310
| gc 1.4790 1.0319 1.433 0.1897
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 0.5662 on 8 degrees of freedom
| Multiple R-squared: 0.9961, Adjusted R-squared: 0.9941
| F-statistic: 508.3 on 4 and 8 DF, p-value: 1.176e-09model.matrix(fit2)| x ga gb gc
| 1 1 1 0 0
| 2 2 1 0 0
| 3 3 0 1 0
| 4 4 0 1 0
| 5 5 0 1 0
| 6 6 0 1 0
| 7 7 0 1 0
| 8 6 0 0 1
| 9 7 0 0 1
| 10 8 0 0 1
| 11 9 0 0 1
| 12 10 0 0 1
| attr(,"assign")
| [1] 1 2 2 2
| attr(,"contrasts")
| attr(,"contrasts")$g
| [1] "contr.treatment"Factor with interaction:
sfit2 <- lm( y ~ x * g, sdf)
summary(sfit2)|
| Call:
| lm(formula = y ~ x * g, data = sdf)
|
| Residuals:
| Min 1Q Median 3Q Max
| -0.57949 -0.34154 -0.07762 0.29828 0.79094
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 1.3864 1.4179 0.978 0.366
| x 0.7337 0.8968 0.818 0.445
| gb -1.3133 1.7596 -0.746 0.484
| gc -0.5363 2.1597 -0.248 0.812
| x:gb 0.2566 0.9189 0.279 0.789
| x:gc 0.3979 0.9189 0.433 0.680
|
| Residual standard error: 0.6341 on 6 degrees of freedom
| Multiple R-squared: 0.9809, Adjusted R-squared: 0.9649
| F-statistic: 61.51 on 5 and 6 DF, p-value: 4.499e-05From which you might conclude that ‘nothing is significant’!
This illustrates that it is often wrong wrong wrong to form conclusions on the basis of scanning p-values in regression output.
Plotting the fitted model:
pred$y2 <- predict( sfit2, newdata = pred)
xyplot( y ~ x, pred, groups = g, type = c('p','r'))
# or
xyplot( y ~ x, pred, groups = g,
subscripts = T, y2 = pred$y2, y1 = pred$y1) +
glayer( panel.xyplot( x, y2[subscripts], ..., type = 'l'))
# or
xyplot( y ~ x, pred, groups = g,
ylim = c(0,13), auto.key = T,
subscripts = T, y2 = pred$y2, y1 = pred$y1) +
glayer( panel.xyplot( x, y2[subscripts], ..., type = 'l')) +
layer( panel.abline( v = c(0,6)))
The model is: \[
\begin{aligned}
E(y|x,g) = & \beta_0 + \beta_x x + \beta_{gb} gb + \beta_{gc} gc\\
& + \beta_{x:gb} x \times gb + \beta_{x:gc} x \times gc \\
\end{aligned}
\] Taking partial derivatives, we see that \(\beta_{gb}\) is the difference between between group ‘b’ minus group ‘a’ when \(x = 0\). There might not be strong evidence of differences between groups outside the range of the data.
What would happen if we were to explore the difference between group ‘b’ and group ‘c’ when \(x = 6\):
L.bc.6 <- rbind( 'c - b|x=6'=
c(0,0,-1,1,-6,6))
wald( sfit2, L.bc.6)| numDF denDF F.value p.value
| 1 1 6 7.293301 0.03555
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| c - b|x=6 1.62458 0.60156 6 2.700611 0.03555 0.152615 3.096545We could also do this by reparametrizing:
sfit2.x6 <- lm( y ~ I(x-6) * relevel(g,'b'), sdf)
summary(sfit2.x6)|
| Call:
| lm(formula = y ~ I(x - 6) * relevel(g, "b"), data = sdf)
|
| Residuals:
| Min 1Q Median 3Q Max
| -0.57949 -0.34154 -0.07762 0.29828 0.79094
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 6.0154 0.3473 17.320 2.37e-06 ***
| I(x - 6) 0.9904 0.2005 4.939 0.00261 **
| relevel(g, "b")a -0.2266 4.0750 -0.056 0.95746
| relevel(g, "b")c 1.6246 0.6016 2.701 0.03555 *
| I(x - 6):relevel(g, "b")a -0.2566 0.9189 -0.279 0.78939
| I(x - 6):relevel(g, "b")c 0.1413 0.2836 0.498 0.63611
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 0.6341 on 6 degrees of freedom
| Multiple R-squared: 0.9809, Adjusted R-squared: 0.9649
| F-statistic: 61.51 on 5 and 6 DF, p-value: 4.499e-05Here are some other ways of exploring the model:
wald(sfit2, ":")| numDF denDF F.value p.value
| : 2 6 0.1890466 0.83249
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| x:gb 0.256648 0.918898 6 0.27930 0.78939 -1.991815 2.505111
| x:gc 0.397910 0.918898 6 0.43303 0.68013 -1.850553 2.646374wald(sfit2, "g")| numDF denDF F.value p.value
| g 4 6 3.965856 0.06566
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| gb -1.313275 1.759556 6 -0.746367 0.48365 -5.618754
| gc -0.536269 2.159667 6 -0.248311 0.81217 -5.820784
| x:gb 0.256648 0.918898 6 0.279300 0.78939 -1.991815
| x:gc 0.397910 0.918898 6 0.433030 0.68013 -1.850553
|
| Coefficients Upper 0.95
| gb 2.992205
| gc 4.748246
| x:gb 2.505111
| x:gc 2.646374wald(sfit2, "x")| numDF denDF F.value p.value
| x 3 6 18.97152 0.00183
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| x 0.733744 0.896753 6 0.818223 0.44450 -1.460531 2.928020
| x:gb 0.256648 0.918898 6 0.279300 0.78939 -1.991815 2.505111
| x:gc 0.397910 0.918898 6 0.433030 0.68013 -1.850553 2.646374Using type 2 Anova gives you tests for ‘g’ and ‘x’ that assume that higher-order interactions involving ‘g’ and ‘x’ are all 0. Note that the error term used is the error term for the full model including interactions. This can lead to inconsistencies with tests based on a model in which interactions has been dropped.
Anova(sfit2) # type 2 anova| Anova Table (Type II tests)
|
| Response: y
| Sum Sq Df F value Pr(>F)
| x 22.7323 1 56.5365 0.0002865 ***
| g 6.2264 2 7.7427 0.0217785 *
| x:g 0.1520 2 0.1890 0.8324941
| Residuals 2.4125 6
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(sfit) # here the gain in degrees of freedom outweighs the increase in SSE| Anova Table (Type II tests)
|
| Response: y
| Sum Sq Df F value Pr(>F)
| x 22.7323 1 70.9133 3.013e-05 ***
| g 6.2264 2 9.7116 0.007243 **
| Residuals 2.5645 8
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Using Lfx with factors
The ‘M’ function associated with ‘Lfx’ can generate code to test for differences between factor levels
Lfx(sfit2)| list( 1,
| 1 * x,
| 1 * M(g),
| 1 * x * M(g)
| )The idea is to use the ‘Lfx’ expression to difference and then to apply to a data frame.
dpred <- expand.grid( x= 0:12, g = levels(sdf$g) , g0 = levels(sdf$g))
dim(dpred)| [1] 117 3some(dpred)| x g g0
| 10 9 a a
| 18 4 b a
| 23 9 b a
| 24 10 b a
| 48 8 a b
| 52 12 a b
| 70 4 c b
| 79 0 a c
| 84 5 a c
| 106 1 c c# we don't need to compare g with g0 with the same levels
# and we only need comparisons in one direction (perhaps not)
dpred <- subset( dpred, g0 < g)
dim(dpred)| [1] 39 3some(dpred)| x g g0
| 16 2 b a
| 18 4 b a
| 20 6 b a
| 31 4 c a
| 33 6 c a
| 36 9 c a
| 39 12 c a
| 71 5 c b
| 75 9 c b
| 76 10 c bLfx(sfit2)| list( 1,
| 1 * x,
| 1 * M(g),
| 1 * x * M(g)
| )# 'difference' g - g0, just like differentiating wrt to g
# except that 'M' function generates differences
Lmat <- Lfx( sfit2,
list( 0,
0 * x,
1 * M(g,g0),
1 * x * M(g,g0)
), dpred)
Lmat| (Intercept) x gb gc x:gb x:gc
| 14 0 0 1 0 0 0
| 15 0 0 1 0 1 0
| 16 0 0 1 0 2 0
| 17 0 0 1 0 3 0
| 18 0 0 1 0 4 0
| 19 0 0 1 0 5 0
| 20 0 0 1 0 6 0
| 21 0 0 1 0 7 0
| 22 0 0 1 0 8 0
| 23 0 0 1 0 9 0
| 24 0 0 1 0 10 0
| 25 0 0 1 0 11 0
| 26 0 0 1 0 12 0
| 27 0 0 0 1 0 0
| 28 0 0 0 1 0 1
| 29 0 0 0 1 0 2
| 30 0 0 0 1 0 3
| 31 0 0 0 1 0 4
| 32 0 0 0 1 0 5
| 33 0 0 0 1 0 6
| 34 0 0 0 1 0 7
| 35 0 0 0 1 0 8
| 36 0 0 0 1 0 9
| 37 0 0 0 1 0 10
| 38 0 0 0 1 0 11
| 39 0 0 0 1 0 12
| 66 0 0 -1 1 0 0
| 67 0 0 -1 1 -1 1
| 68 0 0 -1 1 -2 2
| 69 0 0 -1 1 -3 3
| 70 0 0 -1 1 -4 4
| 71 0 0 -1 1 -5 5
| 72 0 0 -1 1 -6 6
| 73 0 0 -1 1 -7 7
| 74 0 0 -1 1 -8 8
| 75 0 0 -1 1 -9 9
| 76 0 0 -1 1 -10 10
| 77 0 0 -1 1 -11 11
| 78 0 0 -1 1 -12 12
| attr(,"data")
| x g g0
| 14 0 b a
| 15 1 b a
| 16 2 b a
| 17 3 b a
| 18 4 b a
| 19 5 b a
| 20 6 b a
| 21 7 b a
| 22 8 b a
| 23 9 b a
| 24 10 b a
| 25 11 b a
| 26 12 b a
| 27 0 c a
| 28 1 c a
| 29 2 c a
| 30 3 c a
| 31 4 c a
| 32 5 c a
| 33 6 c a
| 34 7 c a
| 35 8 c a
| 36 9 c a
| 37 10 c a
| 38 11 c a
| 39 12 c a
| 66 0 c b
| 67 1 c b
| 68 2 c b
| 69 3 c b
| 70 4 c b
| 71 5 c b
| 72 6 c b
| 73 7 c b
| 74 8 c b
| 75 9 c b
| 76 10 c b
| 77 11 c b
| 78 12 c bwald(sfit2, Lmat)| numDF denDF F.value p.value
| 1 4 6 3.965856 0.06566
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| 14 -1.313275 1.759556 6 -0.746367 0.48365 -5.618754 2.992205
| 15 -1.056627 1.061052 6 -0.995829 0.35778 -3.652928 1.539675
| 16 -0.799979 0.918898 6 -0.870585 0.41745 -3.048442 1.448485
| 17 -0.543331 1.500555 6 -0.362087 0.72970 -4.215056 3.128394
| 18 -0.286683 2.312510 6 -0.123970 0.90539 -5.945191 5.371826
| 19 -0.030035 3.183157 6 -0.009436 0.99278 -7.818940 7.758870
| 20 0.226613 4.075049 6 0.055610 0.95746 -9.744673 10.197899
| 21 0.483261 4.976777 6 0.097103 0.92581 -11.694473 12.660995
| 22 0.739909 5.883820 6 0.125753 0.90404 -13.657280 15.137098
| 23 0.996557 6.794050 6 0.146681 0.88819 -15.627885 17.620999
| 24 1.253205 7.706338 6 0.162620 0.87616 -17.603525 20.109935
| 25 1.509853 8.620030 6 0.175156 0.86672 -19.582602 22.602308
| 26 1.766501 9.534723 6 0.185270 0.85912 -21.564127 25.097129
| 27 -0.536269 2.159667 6 -0.248311 0.81217 -5.820784 4.748246
| 28 -0.138358 1.566112 6 -0.088345 0.93248 -3.970495 3.693779
| 29 0.259552 1.389244 6 0.186830 0.85795 -3.139805 3.658909
| 30 0.657462 1.759556 6 0.373652 0.72151 -3.648017 4.962942
| 31 1.055373 2.439432 6 0.432631 0.68040 -4.913702 7.024447
| 32 1.453283 3.239501 6 0.448613 0.66946 -6.473490 9.380056
| 33 1.851193 4.089823 6 0.452634 0.66672 -8.156242 11.858629
| 34 2.249104 4.964643 6 0.453024 0.66645 -9.898941 14.397148
| 35 2.647014 5.852988 6 0.452250 0.66698 -11.674731 16.968759
| 36 3.044924 6.749518 6 0.451132 0.66774 -13.470552 19.560400
| 37 3.442835 7.651358 6 0.449964 0.66854 -15.279363 22.165032
| 38 3.840745 8.556828 6 0.448851 0.66929 -17.097058 24.778548
| 39 4.238655 9.464886 6 0.447829 0.66999 -18.921088 27.398398
| 66 0.777006 1.933745 6 0.401814 0.70174 -3.954698 5.508711
| 67 0.918268 1.665645 6 0.551299 0.60134 -3.157417 4.993954
| 68 1.059531 1.403640 6 0.754845 0.47891 -2.375054 4.494115
| 69 1.200793 1.151900 6 1.042446 0.33738 -1.617805 4.019391
| 70 1.342055 0.918898 6 1.460505 0.19445 -0.906408 3.590519
| 71 1.483318 0.722985 6 2.051657 0.08604 -0.285764 3.252399
| 72 1.624580 0.601560 6 2.700611 0.03555 0.152615 3.096545
| 73 1.765842 0.601560 6 2.935438 0.02610 0.293878 3.237807
| 74 1.907105 0.722985 6 2.637820 0.03865 0.138023 3.676186
| 75 2.048367 0.918898 6 2.229155 0.06734 -0.200096 4.296830
| 76 2.189630 1.151900 6 1.900885 0.10604 -0.628968 5.008227
| 77 2.330892 1.403640 6 1.660605 0.14786 -1.103692 5.765476
| 78 2.472154 1.665645 6 1.484203 0.18829 -1.603531 6.547840ww <- as.data.frame(wald(sfit2, Lmat))
head(ww)| coef se U2 L2 x g g0
| 14 -1.31327459 1.7595563 2.205838 -4.832387 0 b a
| 15 -1.05662663 1.0610524 1.065478 -3.178731 1 b a
| 16 -0.79997867 0.9188983 1.037818 -2.637775 2 b a
| 17 -0.54333072 1.5005547 2.457779 -3.544440 3 b a
| 18 -0.28668276 2.3125101 4.338337 -4.911703 4 b a
| 19 -0.03003481 3.1831572 6.336280 -6.396349 5 b aww$gap <- with(ww, paste( g, '-', g0))
xyplot( coef ~ x | gap, ww)
td( col = c('black', 'blue','blue'), lty = 1, lwd = 2)
xyplot( coef +I(coef+2*se) + I(coef-2*se) ~ x | gap, ww, type = 'l') +
layer_( panel.abline( h = 0))
# or
xyplot( coef +I(coef+2*se) + I(coef-2*se) ~ x | gap, ww, type = 'l',
ylab = "Estimated difference plus or minus 2 SEs",
xlim = c(0,12))+
layer_( panel.abline( h = 0))
# or
xyplot( coef +I(coef+2*se) + I(coef-2*se) ~ x | gap, ww, type = 'l',
ylab = list("Estimated difference plus or minus 2 SEs", cex = 1.3, font = 2),
xlim = c(0,12))+
layer_( panel.abline( h = 0))
Note that if the significant gap between ‘c’ and ‘b’ around \(x=6\) is a question inspired by the data and not a ‘prior’ hypothesis, then some adjustment should be made for multiplicity. ### Exercises
- Explore how to modify the appearance of the ‘strips’, i.e. where it says ‘c - b’
- Plot approximate ‘95% confidence bands’ for each group with +/- 2 SEs
- Plot approximate ‘95% prediction bands’ for each group.
# reset the random seed
set.seed(NULL) Using WHO regions as predictors of Life Expectancy
fitr <- lm( LifeExp ~ (smoke + I(smoke^2)) * region + hiv + special, dd)
summary(fitr)|
| Call:
| lm(formula = LifeExp ~ (smoke + I(smoke^2)) * region + hiv +
| special, data = dd)
|
| Residuals:
| Min 1Q Median 3Q Max
| -10.3209 -3.1528 0.1424 3.0536 10.7776
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 5.645e+01 1.512e+00 37.324 < 2e-16 ***
| smoke 9.687e-03 1.343e-02 0.721 0.472152
| I(smoke^2) 1.503e-05 1.811e-05 0.830 0.408190
| regionAMR 1.079e+01 2.846e+00 3.792 0.000226 ***
| regionEMR -1.088e+00 3.106e+00 -0.350 0.726834
| regionEUR 2.211e+01 3.982e+00 5.551 1.5e-07 ***
| regionSEAR 1.146e+01 3.977e+00 2.881 0.004634 **
| regionWPR 8.297e+00 7.772e+00 1.068 0.287671
| hiv -2.405e-01 1.112e-01 -2.162 0.032401 *
| special -1.006e+01 2.958e+00 -3.402 0.000885 ***
| smoke:regionAMR 1.876e-02 1.693e-02 1.108 0.269968
| smoke:regionEMR 1.857e-02 1.532e-02 1.212 0.227560
| smoke:regionEUR -9.910e-03 1.430e-02 -0.693 0.489629
| smoke:regionSEAR 2.485e-03 2.294e-02 0.108 0.913918
| smoke:regionWPR 4.084e-03 2.122e-02 0.192 0.847700
| I(smoke^2):regionAMR -3.295e-05 1.988e-05 -1.658 0.099715 .
| I(smoke^2):regionEMR -2.432e-05 1.844e-05 -1.319 0.189551
| I(smoke^2):regionEUR -1.559e-05 1.817e-05 -0.858 0.392392
| I(smoke^2):regionSEAR -2.559e-05 2.434e-05 -1.051 0.295049
| I(smoke^2):regionWPR -1.788e-05 1.938e-05 -0.923 0.357846
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 4.556 on 132 degrees of freedom
| (42 observations deleted due to missingness)
| Multiple R-squared: 0.8045, Adjusted R-squared: 0.7764
| F-statistic: 28.6 on 19 and 132 DF, p-value: < 2.2e-16wald(fitr,":")| numDF denDF F.value p.value
| : 10 132 3.198898 0.00104
|
| Coefficients Estimate Std.Error DF t-value p-value
| smoke:regionAMR 0.018757 0.016932 132 1.107787 0.26997
| smoke:regionEMR 0.018572 0.015320 132 1.212307 0.22756
| smoke:regionEUR -0.009910 0.014304 132 -0.692837 0.48963
| smoke:regionSEAR 0.002485 0.022944 132 0.108305 0.91392
| smoke:regionWPR 0.004084 0.021225 132 0.192431 0.84770
| I(smoke^2):regionAMR -0.000033 0.000020 132 -1.657885 0.09972
| I(smoke^2):regionEMR -0.000024 0.000018 132 -1.318706 0.18955
| I(smoke^2):regionEUR -0.000016 0.000018 132 -0.858101 0.39239
| I(smoke^2):regionSEAR -0.000026 0.000024 132 -1.051283 0.29505
| I(smoke^2):regionWPR -0.000018 0.000019 132 -0.922703 0.35785
|
| Coefficients Lower 0.95 Upper 0.95
| smoke:regionAMR -0.014736 0.052251
| smoke:regionEMR -0.011732 0.048876
| smoke:regionEUR -0.038205 0.018384
| smoke:regionSEAR -0.042901 0.047871
| smoke:regionWPR -0.037900 0.046069
| I(smoke^2):regionAMR -0.000072 0.000006
| I(smoke^2):regionEMR -0.000061 0.000012
| I(smoke^2):regionEUR -0.000052 0.000020
| I(smoke^2):regionSEAR -0.000074 0.000023
| I(smoke^2):regionWPR -0.000056 0.000020wald(fitr,"2):")| numDF denDF F.value p.value
| 2): 5 132 2.093536 0.07012
|
| Coefficients Estimate Std.Error DF t-value p-value
| I(smoke^2):regionAMR -3.3e-05 2.0e-05 132 -1.657885 0.09972
| I(smoke^2):regionEMR -2.4e-05 1.8e-05 132 -1.318706 0.18955
| I(smoke^2):regionEUR -1.6e-05 1.8e-05 132 -0.858101 0.39239
| I(smoke^2):regionSEAR -2.6e-05 2.4e-05 132 -1.051283 0.29505
| I(smoke^2):regionWPR -1.8e-05 1.9e-05 132 -0.922703 0.35785
|
| Coefficients Lower 0.95 Upper 0.95
| I(smoke^2):regionAMR -7.2e-05 6.0e-06
| I(smoke^2):regionEMR -6.1e-05 1.2e-05
| I(smoke^2):regionEUR -5.2e-05 2.0e-05
| I(smoke^2):regionSEAR -7.4e-05 2.3e-05
| I(smoke^2):regionWPR -5.6e-05 2.0e-05## Using data values for prediction instead of creating a separate prediction data frame
# This can work with curvilinear models if the data is sufficiently dense
dd$yq <- predict( fitr, dd)
xyplot( LifeExp ~ smoke | region, dd)
# reorder for nice lines:
dd <- dd[order(dd$region,dd$smoke),]
xyplot( LifeExp ~ smoke | region, dd, subscripts = T, yq = dd$yq) +
layer( panel.xyplot( x, yq[subscripts],..., type = 'b', col = 'red'))
# presents a problem because of missing hiv values
# try again keeping non-missing data together to avoid interrupting lines
dd <- dd[order(is.na(dd$yq),dd$region,dd$smoke),]
xyplot( LifeExp ~ smoke | region, dd, subscripts = T, yq = dd$yq) +
layer( panel.xyplot( x, yq[subscripts],..., type = 'b', col = 'red'))
## Including Health Expenditure
fitrhe <- lm( LifeExp ~ (smoke + I(smoke^2)) * region * log(HE) + hiv + special, dd)
summary(fitrhe)|
| Call:
| lm(formula = LifeExp ~ (smoke + I(smoke^2)) * region * log(HE) +
| hiv + special, data = dd)
|
| Residuals:
| Min 1Q Median 3Q Max
| -6.9890 -1.4067 0.1123 1.2544 9.7261
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 5.088e+01 7.646e+00 6.654 1.12e-09 ***
| smoke -1.832e-02 5.188e-02 -0.353 0.724587
| I(smoke^2) 4.173e-05 7.885e-05 0.529 0.597687
| regionAMR -1.375e+01 1.544e+01 -0.890 0.375292
| regionEMR 6.606e+00 1.782e+01 0.371 0.711655
| regionEUR -8.512e+00 1.948e+01 -0.437 0.662904
| regionSEAR 1.698e+02 8.098e+01 2.097 0.038272 *
| regionWPR -6.176e+00 2.585e+01 -0.239 0.811624
| log(HE) 1.527e+00 1.679e+00 0.910 0.364981
| hiv -3.930e-01 1.003e-01 -3.917 0.000155 ***
| special -1.296e+01 2.343e+00 -5.532 2.13e-07 ***
| smoke:regionAMR 6.425e-02 8.185e-02 0.785 0.434105
| smoke:regionEMR 2.158e-02 6.486e-02 0.333 0.740033
| smoke:regionEUR 2.975e-02 5.868e-02 0.507 0.613198
| smoke:regionSEAR -1.262e+00 6.146e-01 -2.053 0.042407 *
| smoke:regionWPR 1.699e-02 7.631e-02 0.223 0.824199
| I(smoke^2):regionAMR -5.165e-05 9.443e-05 -0.547 0.585509
| I(smoke^2):regionEMR -4.252e-05 8.567e-05 -0.496 0.620633
| I(smoke^2):regionEUR -4.762e-05 7.938e-05 -0.600 0.549800
| I(smoke^2):regionSEAR 1.839e-03 8.954e-04 2.054 0.042279 *
| I(smoke^2):regionWPR -3.620e-05 8.313e-05 -0.435 0.664067
| smoke:log(HE) 5.178e-03 9.762e-03 0.530 0.596867
| I(smoke^2):log(HE) -6.228e-06 1.370e-05 -0.454 0.650385
| regionAMR:log(HE) 3.852e+00 2.780e+00 1.386 0.168672
| regionEMR:log(HE) -1.179e+00 3.844e+00 -0.307 0.759613
| regionEUR:log(HE) 3.195e+00 2.932e+00 1.090 0.278285
| regionSEAR:log(HE) -3.376e+01 1.688e+01 -2.000 0.047953 *
| regionWPR:log(HE) 2.831e+00 4.291e+00 0.660 0.510830
| smoke:regionAMR:log(HE) -1.044e-02 1.392e-02 -0.750 0.454840
| smoke:regionEMR:log(HE) -2.242e-03 1.204e-02 -0.186 0.852601
| smoke:regionEUR:log(HE) -6.670e-03 1.044e-02 -0.639 0.524198
| smoke:regionSEAR:log(HE) 2.685e-01 1.292e-01 2.077 0.040092 *
| smoke:regionWPR:log(HE) -4.551e-03 1.300e-02 -0.350 0.726878
| I(smoke^2):regionAMR:log(HE) 6.664e-06 1.592e-05 0.419 0.676298
| I(smoke^2):regionEMR:log(HE) 5.200e-06 1.448e-05 0.359 0.720256
| I(smoke^2):regionEUR:log(HE) 6.979e-06 1.376e-05 0.507 0.613012
| I(smoke^2):regionSEAR:log(HE) -3.900e-04 1.869e-04 -2.087 0.039160 *
| I(smoke^2):regionWPR:log(HE) 5.384e-06 1.424e-05 0.378 0.705969
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 3.194 on 111 degrees of freedom
| (45 observations deleted due to missingness)
| Multiple R-squared: 0.9154, Adjusted R-squared: 0.8873
| F-statistic: 32.48 on 37 and 111 DF, p-value: < 2.2e-16length(coef(fitrhe))| [1] 38# should have > 380 observations using Harrell's rules of thumb for valid regression
wald(fitrhe, ":") | numDF denDF F.value p.value
| : 27 111 1.37843 0.12542
|
| Coefficients Estimate Std.Error DF t-value p-value
| smoke:regionAMR 0.064254 0.081849 111 0.785035 0.43411
| smoke:regionEMR 0.021575 0.064860 111 0.332641 0.74003
| smoke:regionEUR 0.029747 0.058678 111 0.506947 0.61320
| smoke:regionSEAR -1.261784 0.614556 111 -2.053164 0.04241
| smoke:regionWPR 0.016993 0.076313 111 0.222673 0.82420
| I(smoke^2):regionAMR -0.000052 0.000094 111 -0.546953 0.58551
| I(smoke^2):regionEMR -0.000043 0.000086 111 -0.496344 0.62063
| I(smoke^2):regionEUR -0.000048 0.000079 111 -0.599894 0.54980
| I(smoke^2):regionSEAR 0.001839 0.000895 111 2.054470 0.04228
| I(smoke^2):regionWPR -0.000036 0.000083 111 -0.435470 0.66407
| smoke:log(HE) 0.005178 0.009762 111 0.530439 0.59687
| I(smoke^2):log(HE) -0.000006 0.000014 111 -0.454460 0.65039
| regionAMR:log(HE) 3.851782 2.780032 111 1.385517 0.16867
| regionEMR:log(HE) -1.179181 3.844186 111 -0.306744 0.75961
| regionEUR:log(HE) 3.194526 2.932055 111 1.089518 0.27829
| regionSEAR:log(HE) -33.758269 16.879974 111 -1.999901 0.04795
| regionWPR:log(HE) 2.830507 4.290773 111 0.659673 0.51083
| smoke:regionAMR:log(HE) -0.010442 0.013922 111 -0.750002 0.45484
| smoke:regionEMR:log(HE) -0.002242 0.012041 111 -0.186234 0.85260
| smoke:regionEUR:log(HE) -0.006670 0.010440 111 -0.638910 0.52420
| smoke:regionSEAR:log(HE) 0.268460 0.129242 111 2.077190 0.04009
| smoke:regionWPR:log(HE) -0.004551 0.012996 111 -0.350166 0.72688
| I(smoke^2):regionAMR:log(HE) 0.000007 0.000016 111 0.418627 0.67630
| I(smoke^2):regionEMR:log(HE) 0.000005 0.000014 111 0.359028 0.72026
| I(smoke^2):regionEUR:log(HE) 0.000007 0.000014 111 0.507212 0.61301
| I(smoke^2):regionSEAR:log(HE) -0.000390 0.000187 111 -2.087196 0.03916
| I(smoke^2):regionWPR:log(HE) 0.000005 0.000014 111 0.378247 0.70597
|
| Coefficients Lower 0.95 Upper 0.95
| smoke:regionAMR -0.097934 0.226443
| smoke:regionEMR -0.106950 0.150100
| smoke:regionEUR -0.086528 0.146022
| smoke:regionSEAR -2.479567 -0.044001
| smoke:regionWPR -0.134226 0.168211
| I(smoke^2):regionAMR -0.000239 0.000135
| I(smoke^2):regionEMR -0.000212 0.000127
| I(smoke^2):regionEUR -0.000205 0.000110
| I(smoke^2):regionSEAR 0.000065 0.003614
| I(smoke^2):regionWPR -0.000201 0.000129
| smoke:log(HE) -0.014165 0.024521
| I(smoke^2):log(HE) -0.000033 0.000021
| regionAMR:log(HE) -1.657038 9.360602
| regionEMR:log(HE) -8.796691 6.438330
| regionEUR:log(HE) -2.615536 9.004588
| regionSEAR:log(HE) -67.207064 -0.309474
| regionWPR:log(HE) -5.671945 11.332960
| smoke:regionAMR:log(HE) -0.038029 0.017146
| smoke:regionEMR:log(HE) -0.026102 0.021617
| smoke:regionEUR:log(HE) -0.027359 0.014018
| smoke:regionSEAR:log(HE) 0.012358 0.524561
| smoke:regionWPR:log(HE) -0.030303 0.021202
| I(smoke^2):regionAMR:log(HE) -0.000025 0.000038
| I(smoke^2):regionEMR:log(HE) -0.000024 0.000034
| I(smoke^2):regionEUR:log(HE) -0.000020 0.000034
| I(smoke^2):regionSEAR:log(HE) -0.000760 -0.000020
| I(smoke^2):regionWPR:log(HE) -0.000023 0.000034wald(fitrhe, "2") | numDF denDF F.value p.value
| 2 12 111 0.7935878 0.65615
|
| Coefficients Estimate Std.Error DF t-value p-value
| I(smoke^2) 0.000042 0.000079 111 0.529253 0.59769
| I(smoke^2):regionAMR -0.000052 0.000094 111 -0.546953 0.58551
| I(smoke^2):regionEMR -0.000043 0.000086 111 -0.496344 0.62063
| I(smoke^2):regionEUR -0.000048 0.000079 111 -0.599894 0.54980
| I(smoke^2):regionSEAR 0.001839 0.000895 111 2.054470 0.04228
| I(smoke^2):regionWPR -0.000036 0.000083 111 -0.435470 0.66407
| I(smoke^2):log(HE) -0.000006 0.000014 111 -0.454460 0.65039
| I(smoke^2):regionAMR:log(HE) 0.000007 0.000016 111 0.418627 0.67630
| I(smoke^2):regionEMR:log(HE) 0.000005 0.000014 111 0.359028 0.72026
| I(smoke^2):regionEUR:log(HE) 0.000007 0.000014 111 0.507212 0.61301
| I(smoke^2):regionSEAR:log(HE) -0.000390 0.000187 111 -2.087196 0.03916
| I(smoke^2):regionWPR:log(HE) 0.000005 0.000014 111 0.378247 0.70597
|
| Coefficients Lower 0.95 Upper 0.95
| I(smoke^2) -0.000115 0.000198
| I(smoke^2):regionAMR -0.000239 0.000135
| I(smoke^2):regionEMR -0.000212 0.000127
| I(smoke^2):regionEUR -0.000205 0.000110
| I(smoke^2):regionSEAR 0.000065 0.003614
| I(smoke^2):regionWPR -0.000201 0.000129
| I(smoke^2):log(HE) -0.000033 0.000021
| I(smoke^2):regionAMR:log(HE) -0.000025 0.000038
| I(smoke^2):regionEMR:log(HE) -0.000024 0.000034
| I(smoke^2):regionEUR:log(HE) -0.000020 0.000034
| I(smoke^2):regionSEAR:log(HE) -0.000760 -0.000020
| I(smoke^2):regionWPR:log(HE) -0.000023 0.000034wald(fitrhe, "2|:.*:") # quadratic terms and 3 and higher way interaction| numDF denDF F.value p.value
| 2|:.*: 17 111 1.286303 0.21443
|
| Coefficients Estimate Std.Error DF t-value p-value
| I(smoke^2) 0.000042 0.000079 111 0.529253 0.59769
| I(smoke^2):regionAMR -0.000052 0.000094 111 -0.546953 0.58551
| I(smoke^2):regionEMR -0.000043 0.000086 111 -0.496344 0.62063
| I(smoke^2):regionEUR -0.000048 0.000079 111 -0.599894 0.54980
| I(smoke^2):regionSEAR 0.001839 0.000895 111 2.054470 0.04228
| I(smoke^2):regionWPR -0.000036 0.000083 111 -0.435470 0.66407
| I(smoke^2):log(HE) -0.000006 0.000014 111 -0.454460 0.65039
| smoke:regionAMR:log(HE) -0.010442 0.013922 111 -0.750002 0.45484
| smoke:regionEMR:log(HE) -0.002242 0.012041 111 -0.186234 0.85260
| smoke:regionEUR:log(HE) -0.006670 0.010440 111 -0.638910 0.52420
| smoke:regionSEAR:log(HE) 0.268460 0.129242 111 2.077190 0.04009
| smoke:regionWPR:log(HE) -0.004551 0.012996 111 -0.350166 0.72688
| I(smoke^2):regionAMR:log(HE) 0.000007 0.000016 111 0.418627 0.67630
| I(smoke^2):regionEMR:log(HE) 0.000005 0.000014 111 0.359028 0.72026
| I(smoke^2):regionEUR:log(HE) 0.000007 0.000014 111 0.507212 0.61301
| I(smoke^2):regionSEAR:log(HE) -0.000390 0.000187 111 -2.087196 0.03916
| I(smoke^2):regionWPR:log(HE) 0.000005 0.000014 111 0.378247 0.70597
|
| Coefficients Lower 0.95 Upper 0.95
| I(smoke^2) -0.000115 0.000198
| I(smoke^2):regionAMR -0.000239 0.000135
| I(smoke^2):regionEMR -0.000212 0.000127
| I(smoke^2):regionEUR -0.000205 0.000110
| I(smoke^2):regionSEAR 0.000065 0.003614
| I(smoke^2):regionWPR -0.000201 0.000129
| I(smoke^2):log(HE) -0.000033 0.000021
| smoke:regionAMR:log(HE) -0.038029 0.017146
| smoke:regionEMR:log(HE) -0.026102 0.021617
| smoke:regionEUR:log(HE) -0.027359 0.014018
| smoke:regionSEAR:log(HE) 0.012358 0.524561
| smoke:regionWPR:log(HE) -0.030303 0.021202
| I(smoke^2):regionAMR:log(HE) -0.000025 0.000038
| I(smoke^2):regionEMR:log(HE) -0.000024 0.000034
| I(smoke^2):regionEUR:log(HE) -0.000020 0.000034
| I(smoke^2):regionSEAR:log(HE) -0.000760 -0.000020
| I(smoke^2):regionWPR:log(HE) -0.000023 0.000034fitr2 <- lm( LifeExp ~ (smoke + log(HE) + region)^2 + hiv + special, dd)
summary(fitr2)|
| Call:
| lm(formula = LifeExp ~ (smoke + log(HE) + region)^2 + hiv + special,
| data = dd)
|
| Residuals:
| Min 1Q Median 3Q Max
| -6.9873 -1.8511 0.1576 1.7496 8.8413
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 4.534e+01 3.590e+00 12.628 < 2e-16 ***
| smoke 1.726e-02 4.899e-03 3.522 0.000594 ***
| log(HE) 2.623e+00 8.903e-01 2.947 0.003816 **
| regionAMR 6.004e+00 6.061e+00 0.991 0.323748
| regionEMR 3.025e+00 8.003e+00 0.378 0.706085
| regionEUR -4.070e+00 7.195e+00 -0.566 0.572584
| regionSEAR 7.659e+00 6.956e+00 1.101 0.272966
| regionWPR -2.520e+00 5.886e+00 -0.428 0.669332
| hiv -4.235e-01 9.424e-02 -4.494 1.55e-05 ***
| special -1.352e+01 2.183e+00 -6.195 7.35e-09 ***
| smoke:log(HE) -6.207e-04 6.668e-04 -0.931 0.353725
| smoke:regionAMR -1.142e-02 4.176e-03 -2.735 0.007128 **
| smoke:regionEMR -8.776e-03 4.180e-03 -2.099 0.037745 *
| smoke:regionEUR -1.372e-02 3.758e-03 -3.651 0.000379 ***
| smoke:regionSEAR -1.321e-02 4.980e-03 -2.652 0.009019 **
| smoke:regionWPR -1.100e-02 4.111e-03 -2.675 0.008450 **
| log(HE):regionAMR 8.114e-01 1.180e+00 0.687 0.493061
| log(HE):regionEMR 3.281e-01 1.700e+00 0.193 0.847248
| log(HE):regionEUR 2.389e+00 1.156e+00 2.067 0.040730 *
| log(HE):regionSEAR 7.592e-01 1.455e+00 0.522 0.602617
| log(HE):regionWPR 2.073e+00 1.159e+00 1.789 0.075940 .
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 3.254 on 128 degrees of freedom
| (45 observations deleted due to missingness)
| Multiple R-squared: 0.8988, Adjusted R-squared: 0.883
| F-statistic: 56.83 on 20 and 128 DF, p-value: < 2.2e-16wald(fitr2, ':')| numDF denDF F.value p.value
| : 11 128 2.478905 0.00752
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| smoke:log(HE) -0.000621 0.000667 128 -0.930769 0.35373 -0.001940
| smoke:regionAMR -0.011419 0.004176 128 -2.734757 0.00713 -0.019681
| smoke:regionEMR -0.008776 0.004180 128 -2.099411 0.03774 -0.017047
| smoke:regionEUR -0.013722 0.003758 128 -3.651365 0.00038 -0.021158
| smoke:regionSEAR -0.013207 0.004980 128 -2.651798 0.00902 -0.023061
| smoke:regionWPR -0.010997 0.004111 128 -2.674936 0.00845 -0.019132
| log(HE):regionAMR 0.811401 1.180355 128 0.687421 0.49306 -1.524133
| log(HE):regionEMR 0.328123 1.699935 128 0.193021 0.84725 -3.035488
| log(HE):regionEUR 2.389164 1.155748 128 2.067201 0.04073 0.102319
| log(HE):regionSEAR 0.759205 1.454596 128 0.521935 0.60262 -2.118962
| log(HE):regionWPR 2.073361 1.158788 128 1.789249 0.07594 -0.219500
|
| Coefficients Upper 0.95
| smoke:log(HE) 0.000699
| smoke:regionAMR -0.003157
| smoke:regionEMR -0.000505
| smoke:regionEUR -0.006286
| smoke:regionSEAR -0.003352
| smoke:regionWPR -0.002862
| log(HE):regionAMR 3.146934
| log(HE):regionEMR 3.691734
| log(HE):regionEUR 4.676008
| log(HE):regionSEAR 3.637371
| log(HE):regionWPR 4.366221wald(fitr2, 'HE):|:log')| numDF denDF F.value p.value
| HE):|:log 6 128 0.9991023 0.42891
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| smoke:log(HE) -0.000621 0.000667 128 -0.930769 0.35373 -0.001940
| log(HE):regionAMR 0.811401 1.180355 128 0.687421 0.49306 -1.524133
| log(HE):regionEMR 0.328123 1.699935 128 0.193021 0.84725 -3.035488
| log(HE):regionEUR 2.389164 1.155748 128 2.067201 0.04073 0.102319
| log(HE):regionSEAR 0.759205 1.454596 128 0.521935 0.60262 -2.118962
| log(HE):regionWPR 2.073361 1.158788 128 1.789249 0.07594 -0.219500
|
| Coefficients Upper 0.95
| smoke:log(HE) 0.000699
| log(HE):regionAMR 3.146934
| log(HE):regionEMR 3.691734
| log(HE):regionEUR 4.676008
| log(HE):regionSEAR 3.637371
| log(HE):regionWPR 4.366221fitr3 <- lm( LifeExp ~ region* smoke + log(HE)+ hiv + special, dd)
summary(fitr3)|
| Call:
| lm(formula = LifeExp ~ region * smoke + log(HE) + hiv + special,
| data = dd)
|
| Residuals:
| Min 1Q Median 3Q Max
| -7.016 -1.863 0.165 1.627 8.700
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 41.364150 1.472822 28.085 < 2e-16 ***
| regionAMR 8.712382 1.438295 6.057 1.31e-08 ***
| regionEMR 3.799335 1.759763 2.159 0.032631 *
| regionEUR 9.842506 1.887832 5.214 6.84e-07 ***
| regionSEAR 10.122815 1.839269 5.504 1.83e-07 ***
| regionWPR 7.460923 1.969657 3.788 0.000229 ***
| smoke 0.010666 0.002712 3.933 0.000134 ***
| log(HE) 3.688028 0.314332 11.733 < 2e-16 ***
| hiv -0.501073 0.079030 -6.340 3.25e-09 ***
| special -14.766685 2.013362 -7.334 1.91e-11 ***
| regionAMR:smoke -0.009904 0.003175 -3.120 0.002216 **
| regionEMR:smoke -0.007235 0.003068 -2.358 0.019822 *
| regionEUR:smoke -0.011760 0.002835 -4.148 5.93e-05 ***
| regionSEAR:smoke -0.009663 0.004494 -2.150 0.033335 *
| regionWPR:smoke -0.008164 0.003110 -2.625 0.009666 **
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 3.254 on 134 degrees of freedom
| (45 observations deleted due to missingness)
| Multiple R-squared: 0.894, Adjusted R-squared: 0.883
| F-statistic: 80.76 on 14 and 134 DF, p-value: < 2.2e-16wald(fitr3, ":")| numDF denDF F.value p.value
| : 5 134 4.25484 0.00127
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| regionAMR:smoke -0.009904 0.003175 134 -3.119787 0.00222 -0.016183
| regionEMR:smoke -0.007235 0.003068 134 -2.357952 0.01982 -0.013303
| regionEUR:smoke -0.011760 0.002835 134 -4.147578 0.00006 -0.017368
| regionSEAR:smoke -0.009663 0.004494 134 -2.150191 0.03334 -0.018551
| regionWPR:smoke -0.008164 0.003110 134 -2.625207 0.00967 -0.014314
|
| Coefficients Upper 0.95
| regionAMR:smoke -0.003625
| regionEMR:smoke -0.001166
| regionEUR:smoke -0.006152
| regionSEAR:smoke -0.000775
| regionWPR:smoke -0.002013wald(fitr3, 'region')| numDF denDF F.value p.value
| region 10 134 8.121058 <.00001
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| regionAMR 8.712382 1.438295 134 6.057439 <.00001 5.867686
| regionEMR 3.799335 1.759763 134 2.159003 0.03263 0.318829
| regionEUR 9.842506 1.887832 134 5.213655 <.00001 6.108703
| regionSEAR 10.122815 1.839269 134 5.503717 <.00001 6.485062
| regionWPR 7.460923 1.969657 134 3.787931 0.00023 3.565285
| regionAMR:smoke -0.009904 0.003175 134 -3.119787 0.00222 -0.016183
| regionEMR:smoke -0.007235 0.003068 134 -2.357952 0.01982 -0.013303
| regionEUR:smoke -0.011760 0.002835 134 -4.147578 0.00006 -0.017368
| regionSEAR:smoke -0.009663 0.004494 134 -2.150191 0.03334 -0.018551
| regionWPR:smoke -0.008164 0.003110 134 -2.625207 0.00967 -0.014314
|
| Coefficients Upper 0.95
| regionAMR 11.557078
| regionEMR 7.279840
| regionEUR 13.576310
| regionSEAR 13.760568
| regionWPR 11.356561
| regionAMR:smoke -0.003625
| regionEMR:smoke -0.001166
| regionEUR:smoke -0.006152
| regionSEAR:smoke -0.000775
| regionWPR:smoke -0.002013wald(fitr3, 'HE')| numDF denDF F.value p.value
| HE 1 134 137.6607 <.00001
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| log(HE) 3.688028 0.314332 134 11.73289 <.00001 3.066333 4.309723Anova(fitr3)| Anova Table (Type II tests)
|
| Response: LifeExp
| Sum Sq Df F value Pr(>F)
| region 634.67 5 11.9873 1.339e-09 ***
| smoke 6.95 1 0.6561 0.41936
| log(HE) 1457.70 1 137.6607 < 2.2e-16 ***
| hiv 425.67 1 40.1992 3.250e-09 ***
| special 569.62 1 53.7926 1.910e-11 ***
| region:smoke 225.27 5 4.2548 0.00127 **
| Residuals 1418.94 134
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fitr3)| Analysis of Variance Table
|
| Response: LifeExp
| Df Sum Sq Mean Sq F value Pr(>F)
| region 5 8942.0 1788.40 168.8905 < 2.2e-16 ***
| smoke 1 141.2 141.20 13.3348 0.0003725 ***
| log(HE) 1 1788.2 1788.22 168.8740 < 2.2e-16 ***
| hiv 1 337.6 337.64 31.8856 9.375e-08 ***
| special 1 537.9 537.94 50.8014 5.711e-11 ***
| region:smoke 5 225.3 45.05 4.2548 0.0012701 **
| Residuals 134 1418.9 10.59
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1wald(fitr3, 'smoke')| numDF denDF F.value p.value
| smoke 6 134 3.655058 0.00216
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| smoke 0.010666 0.002712 134 3.933292 0.00013 0.005303
| regionAMR:smoke -0.009904 0.003175 134 -3.119787 0.00222 -0.016183
| regionEMR:smoke -0.007235 0.003068 134 -2.357952 0.01982 -0.013303
| regionEUR:smoke -0.011760 0.002835 134 -4.147578 0.00006 -0.017368
| regionSEAR:smoke -0.009663 0.004494 134 -2.150191 0.03334 -0.018551
| regionWPR:smoke -0.008164 0.003110 134 -2.625207 0.00967 -0.014314
|
| Coefficients Upper 0.95
| smoke 0.016029
| regionAMR:smoke -0.003625
| regionEMR:smoke -0.001166
| regionEUR:smoke -0.006152
| regionSEAR:smoke -0.000775
| regionWPR:smoke -0.002013library(p3d)
Plot3d(LifeExp ~ smoke + HE | region, dd)| region col
| 1 AFR blue
| 2 AMR green
| 3 EMR orange
| 4 EUR magenta
| 5 SEAR cyan
| 6 WPR red
|
| Use left mouse to rotate, middle mouse (or scroll) to zoom, right mouse to change perspectiveFit3d( fitr3, other.vars=list(hiv=0,special=0))| Warning in log(HE): NaNs produced#Id3d(par=2, labels = dd$country)
par3d(windowRect=c(10,10,700,700))
rgl.snapshot('regions.png')
### Exercise Explore this data further producing informative graphs.