env <- new.env()
load("../data.rda", envir = env)
df <- get("CH08TA01", envir = env)
names(df) <- c("Y", "X1", "X2")
df <- transform(df,
x1 = scale(X1, scale = 0.4),
x2 = scale(X2, scale = 10),
x1sq = scale(X1, scale = 0.4)^2,
x2sq = scale(X2, scale = 10)^2,
x1x2 = scale(X1,, 0.4) * scale(X2,, 10))
fit <- lm(Y ~ ., df)
model.frame(df)
## Y X1 X2 x1 x2 x1sq x2sq x1x2
## 1 150 0.6 10 -1 -1 1 1 1
## 2 86 1.0 10 0 -1 0 1 0
## 3 49 1.4 10 1 -1 1 1 -1
## 4 288 0.6 20 -1 0 1 0 0
## 5 157 1.0 20 0 0 0 0 0
## 6 131 1.0 20 0 0 0 0 0
## 7 184 1.0 20 0 0 0 0 0
## 8 109 1.4 20 1 0 1 0 0
## 9 279 0.6 30 -1 1 1 1 -1
## 10 235 1.0 30 0 1 0 1 0
## 11 224 1.4 30 1 1 1 1 1
par(mfrow = c(2, 2), pch = 19)
plot(resid(fit) ~ fitted(fit), xlab = "Fitted", ylab = "Residual")
title("(a) Residual Plot Against Fitted Values")
plot(resid(fit) ~ x1, df, xlab = "x1", ylab = "Residual")
title("(b) Residual Plot against x1")
plot(resid(fit) ~ x2, df, xlab = "x2", ylab = "Residual")
title("(c) Residual Plot against x2")
qqnorm(resid(fit), xlab = "Expected", ylab = "Residual", main = "")
qqline(resid(fit))
title("(d) Normal Probability Plot")
A lot of analysis was performed over these pages, both shown and not shown. For that reason I will follow a similar process. However, results that can be easily obtained by summary functions will be ignored. The F-test that is performed is not that in the summary results (either here or in the SAS output the authors provide). The summary results show the F-test as to whether or not all the predictor variables are relevant. This statistic is 10.57, which is still less than the 19.2. Recall from Chapter 3 that the alr3 library contains pureErrorAnova for an extended ANOVA table.
library(alr3)
fit <- lm(Y ~ ., df) # Restate the model
summary(fit)
##
## Call:
## lm(formula = Y ~ ., data = df)
##
## Residuals:
## 1 2 3 4 5 6 7 8 9 10
## -21.46 9.26 12.20 41.93 -5.84 -31.84 21.16 -25.40 -20.46 7.26
## 11
## 13.20
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 151.43 45.46 3.33 0.0208 *
## X1 -139.58 33.04 -4.22 0.0083 **
## X2 7.55 1.32 5.71 0.0023 **
## x1 NA NA NA NA
## x2 NA NA NA NA
## x1sq 27.39 20.34 1.35 0.2359
## x2sq -10.61 20.34 -0.52 0.6244
## x1x2 11.50 16.19 0.71 0.5092
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 32.4 on 5 degrees of freedom
## Multiple R-squared: 0.914, Adjusted R-squared: 0.827
## F-statistic: 10.6 on 5 and 5 DF, p-value: 0.0109
pureErrorAnova(fit) # alr3 function
## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 18704 18704 26.63 0.036 *
## X2 1 34201 34201 48.70 0.020 *
## x1sq 1 1646 1646 2.34 0.265
## x2sq 1 285 285 0.41 0.589
## x1x2 1 529 529 0.75 0.477
## Residuals 5 5240 1048
## Lack of fit 3 3836 1279 1.82 0.374
## Pure Error 2 1405 702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(Y ~ 1, df), fit) # Explicit F-Test
## Analysis of Variance Table
##
## Model 1: Y ~ 1
## Model 2: Y ~ X1 + X2 + x1 + x2 + x1sq + x2sq + x1x2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 10 60606
## 2 5 5240 5 55366 10.6 0.011 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate F* (p 302)
pureErrorAnova(fit)[" Lack of fit", "Mean Sq"] / pureErrorAnova(fit)[" Pure Error", "Mean Sq"] # (6.68b)
## [1] 1.82
# Correlation (p 301)
with(df, rbind(
"X and Xsq" = round(c("1" = cor(X1, X1^2), "2" = cor(X2, X2^2)), 3),
"x and xsq" = round(c("1" = cor(x1, x1^2), "2" = cor(x2, x2^2)), 3)))
## 1 2
## X and Xsq 0.991 0.986
## x and xsq 0.000 0.000
anova(lm(Y ~ x1 + x2, df), fit) # partial F-test
## Analysis of Variance Table
##
## Model 1: Y ~ x1 + x2
## Model 2: Y ~ X1 + X2 + x1 + x2 + x1sq + x2sq + x1x2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 8 7700
## 2 5 5240 3 2460 0.78 0.55
fit <- lm(Y ~ x1 + x2, df) # refit model (p 304)
# Diagnostic Plots
par(mfrow = c(2, 2), pch = 19) # Same as FIGURE 8.5
plot(resid(fit) ~ fitted(fit), xlab = "Fitted", ylab = "Residual")
title("(a) Residual Plot Against Fitted Values")
plot(resid(fit) ~ x1, df, xlab = "x1", ylab = "Residual")
title("(b) Residual Plot against x1")
plot(resid(fit) ~ x2, df, xlab = "x2", ylab = "Residual")
title("(c) Residual Plot against x2")
qqnorm(resid(fit), xlab = "Expected", ylab = "Residual", main = "")
title("(d) Normal Probability Plot")
qqline(resid(fit))
confint(lm(Y ~ X1 + X2, df)) # Bonferroni procedure on original coefs. (p 305)
## 2.5 % 97.5 %
## (Intercept) 64.618 256.55
## X1 -212.602 -66.56
## X2 4.629 10.47
Here we make use of the function I defined in chapter 6 for FIGURE 6.7.
lm3d <- function(model, res, ...) {
ivs <- labels(terms(model))
x <- model.frame(model)[, ivs[1]] # 1st independent variable
y <- model.frame(model)[, ivs[2]] # 2nd independent variable
xr <- seq(min(x), max(x), len = res) # equidistant sequence from range of 1st iv
yr <- seq(min(y), max(y), len = res) # equidistant sequence from range of 2nd iv
# Create a list object of these sequences to be turned into a grid
newdata <- list(xr, yr)
names(newdata) <- ivs
newdata <- do.call('expand.grid', newdata)
zr <- matrix(predict(model, newdata), res, res)
persp(xr, yr, zr, ...)
}
lm3d(lm(Y ~ X1 + X2, df), res = 30, ticktype = "detailed", shade = 0.5, expand = 2/3,
xlab = "Charge Rate", ylab = "Temp", zlab = "Cycles",
theta = 320, phi = 30)
df <- get("CH07TA01", envir = env)
names(df) <- c("X1", "X2", "X3", "Y")
# Add centered variables
df <- transform(df,
x1 = scale(X1, T, F),
x2 = scale(X2, T, F),
x3 = scale(X3, T, F))
fit <- lm(Y ~ (X1 + X2 + X3)^2, df)
cor(model.matrix(fit)[, -1]) # Interactions, high correlation (p 312)
## X1 X2 X3 X1:X2 X1:X3 X2:X3
## X1 1.0000 0.92384 0.45778 0.9888 0.9003 0.8907
## X2 0.9238 1.00000 0.08467 0.9663 0.6720 0.6536
## X3 0.4578 0.08467 1.00000 0.3324 0.7877 0.8064
## X1:X2 0.9888 0.96634 0.33239 1.0000 0.8345 0.8219
## X1:X3 0.9003 0.67197 0.78770 0.8345 1.0000 0.9984
## X2:X3 0.8907 0.65361 0.80641 0.8219 0.9984 1.0000
fit <- lm(Y ~ (x1 + x2 + x3)^2, df) # fit centered values instead
cor(model.matrix(fit)[, -1]) # less correlation, but still present
## x1 x2 x3 x1:x2 x1:x3 x2:x3
## x1 1.0000 0.92384 0.45778 -0.4770 -0.17342 -0.2216
## x2 0.9238 1.00000 0.08467 -0.4298 -0.17254 -0.1437
## x3 0.4578 0.08467 1.00000 -0.2159 -0.03041 -0.2354
## x1:x2 -0.4770 -0.42979 -0.21589 1.0000 0.23283 0.2919
## x1:x3 -0.1734 -0.17254 -0.03041 0.2328 1.00000 0.8905
## x2:x3 -0.2216 -0.14366 -0.23537 0.2919 0.89051 1.0000
summary(fit)
##
## Call:
## lm(formula = Y ~ (x1 + x2 + x3)^2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.927 -0.873 0.168 1.293 3.331
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.52689 1.07363 19.12 6.7e-11 ***
## x1 3.43781 3.57867 0.96 0.35
## x2 -2.09472 3.03677 -0.69 0.50
## x3 -1.61634 1.90721 -0.85 0.41
## x1:x2 0.00888 0.03085 0.29 0.78
## x1:x3 -0.08479 0.07342 -1.15 0.27
## x2:x3 0.09042 0.09200 0.98 0.34
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.6 on 13 degrees of freedom
## Multiple R-squared: 0.823, Adjusted R-squared: 0.741
## F-statistic: 10.1 on 6 and 13 DF, p-value: 0.000296
anova(fit)
## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 352 352 52.22 6.7e-06 ***
## x2 1 33 33 4.92 0.045 *
## x3 1 12 12 1.71 0.213
## x1:x2 1 1 1 0.22 0.646
## x1:x3 1 3 3 0.40 0.538
## x2:x3 1 7 7 0.97 0.344
## Residuals 13 88 7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(Y ~ x1 + x2 + x3, df), fit) # partial F-test: F = 0.53 < qf(.95, 3,13)
## Analysis of Variance Table
##
## Model 1: Y ~ x1 + x2 + x3
## Model 2: Y ~ (x1 + x2 + x3)^2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 16 98.4
## 2 13 87.7 3 10.7 0.53 0.67
df <- get("CH08TA02", envir = env)
names(df) <- c("Y", "X1", "X2")
df <- transform(df, X2 = factor(X2, labels = c("Mutual", "Stock")))
fit <- lm(Y ~ X1 + X2, df)
When using abline on a fitted model, it is really just supplying the slope and intercept coefficients to the relevant parameters. We can use this fact to adjust those parameters for the two situations. In particular, the default parameters to abline are of the form
abline(intercept, slope)
Under a qualitative predictor, when it is not the reference level, we see that coefficient addition to the intercept term. This is what we mimic below given how abline works.
plot(Y ~ X1, df, pch = 19, ylim = c(0, 40), xlim = c(0, 310),
xlab = "Size of Firm", ylab = "Months Elapsed")
abline(coef(fit)[[1]], coef(fit)[[2]],) # Mutual Firm
abline(coef(fit)[[1]] + coef(fit)[[3]], coef(fit)[[2]], col = "red") # Stock Firm
If you simply tried to cbind a factor value, it will return the way R stores it as 1 and 2. There are 2 approaches used below to access the information we want. First, we can use the model.matrix for the fitted model. This returns the X matrix used in the linear model which contains the binary encoding for a dummy variable (multiple columns for multiple category predictor). The other approach is to realize that if you remove the class (unclass) from that variable, it returns as the storage 1 and 2 values, numerically. We can then offset that by 1 to get a binary encoding as desired. This is performed below to obtain the interaction terms.
data.frame(
"Months" = df$Y,
"Firm Size" = df$X1,
"Firm Type" = df$X2,
"Indicator" = model.matrix(fit)[, 3],
"Interaction" = df$X1 * (unclass(df$X2) - 1)
)
## Months Firm.Size Firm.Type Indicator Interaction
## 1 17 151 Mutual 0 0
## 2 26 92 Mutual 0 0
## 3 21 175 Mutual 0 0
## 4 30 31 Mutual 0 0
## 5 22 104 Mutual 0 0
## 6 0 277 Mutual 0 0
## 7 12 210 Mutual 0 0
## 8 19 120 Mutual 0 0
## 9 4 290 Mutual 0 0
## 10 16 238 Mutual 0 0
## 11 28 164 Stock 1 164
## 12 15 272 Stock 1 272
## 13 11 295 Stock 1 295
## 14 38 68 Stock 1 68
## 15 31 85 Stock 1 85
## 16 21 224 Stock 1 224
## 17 20 166 Stock 1 166
## 18 13 305 Stock 1 305
## 19 30 124 Stock 1 124
## 20 14 246 Stock 1 246
Also included are the results for the confidence intervals of the parameters. The formal test is inherent in summary as the p-values and printed stars demonstrate these coeffficients are significant at an alpha = 0.01 level. This is confirmed by the respective confidence interval not including zero.
summary(fit)
##
## Call:
## lm(formula = Y ~ X1 + X2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.692 -1.704 -0.439 1.921 6.341
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.87407 1.81386 18.68 9.1e-13 ***
## X1 -0.10174 0.00889 -11.44 2.1e-09 ***
## X2Stock 8.05547 1.45911 5.52 3.7e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.22 on 17 degrees of freedom
## Multiple R-squared: 0.895, Adjusted R-squared: 0.883
## F-statistic: 72.5 on 2 and 17 DF, p-value: 4.77e-09
anova(fit)
## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 1188 1188 114.5 5.7e-09 ***
## X2 1 316 316 30.5 3.7e-05 ***
## Residuals 17 176 10
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(fit)
## 2.5 % 97.5 %
## (Intercept) 30.0472 37.70098
## X1 -0.1205 -0.08298
## X2Stock 4.9770 11.13391
fit <- lm(Y ~ X1 + X2, df)
plot(Y ~ X1, df, col = X2, pch = 19, ylim = c(0, 40), xlim = c(0, 310),
xlab = "Size of Firm", ylab = "Months Elapsed")
abline(coef(fit)[[1]], coef(fit)[[2]]) # Mutual Firm
abline(coef(fit)[[1]] + coef(fit)[[3]], coef(fit)[[2]], col = "red") # Stock Firm
Note that the author does an illustration that exaggerates the change in slope. As is evidenced from this example, the slope is relatively unchanged so that the two lines appear parallel. The difference between the two is the coefficient X1X2 = -0.000417
Here I set a new model that nests the X2 levels. By excluding the intercept term, the coefficients 1 and 3 apply to the first X2 level (mutual firm) and the coefficients 2 and 4 apply to the second X2 level (stock firm).
fit <- lm(Y ~ X2/X1 - 1, df) # Using nested model for ease of presentation
plot(Y ~ X1, df, col = X2, pch = 19, ylim = c(0, 40), xlim = c(0, 310),
xlab = "Size of Firm", ylab = "Months Elapsed")
abline(coef(fit)[[1]], coef(fit)[[3]]) # Mutual Firm
abline(coef(fit)[[2]], coef(fit)[[4]], col = "red") # Stock Firm
fit <- lm(Y ~ X1 * X2, df) # Going back to intended model
summary(fit)
##
## Call:
## lm(formula = Y ~ X1 * X2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.714 -1.706 -0.456 1.931 6.326
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.838369 2.440650 13.86 2.5e-10 ***
## X1 -0.101531 0.013053 -7.78 8.0e-07 ***
## X2Stock 8.131250 3.654052 2.23 0.041 *
## X1:X2Stock -0.000417 0.018331 -0.02 0.982
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.32 on 16 degrees of freedom
## Multiple R-squared: 0.895, Adjusted R-squared: 0.875
## F-statistic: 45.5 on 3 and 16 DF, p-value: 4.67e-08
anova(fit)
## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 1188 1188 107.8 1.6e-08 ***
## X2 1 316 316 28.7 6.4e-05 ***
## X1:X2 1 0 0 0.0 0.98
## Residuals 16 176 11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(Y ~ X1 + X2, df), fit)
## Analysis of Variance Table
##
## Model 1: Y ~ X1 + X2
## Model 2: Y ~ X1 * X2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 17 176
## 2 16 176 1 0.00571 0 0.98
df <- get("CH08TA05", envir = env)
names(df) <- c("Y", "X1", "X2")
df <- transform(df, X2 = factor(X2, labels = c("Line1", "Line2")))
data.frame(
"Scrap" = df$Y,
"Line_Speed" = df$X1,
"Indicator" = df$X2)
## Scrap Line_Speed Indicator
## 1 218 100 Line2
## 2 248 125 Line2
## 3 360 220 Line2
## 4 351 205 Line2
## 5 470 300 Line2
## 6 394 255 Line2
## 7 332 225 Line2
## 8 321 175 Line2
## 9 410 270 Line2
## 10 260 170 Line2
## 11 241 155 Line2
## 12 331 190 Line2
## 13 275 140 Line2
## 14 425 290 Line2
## 15 367 265 Line2
## 16 140 105 Line1
## 17 277 215 Line1
## 18 384 270 Line1
## 19 341 255 Line1
## 20 215 175 Line1
## 21 180 135 Line1
## 22 260 200 Line1
## 23 361 275 Line1
## 24 252 155 Line1
## 25 422 320 Line1
## 26 273 190 Line1
## 27 410 295 Line1
fit <- lm(Y ~ X2/X1 - 1, df) # Using nested model for ease of presentation
plot(Y ~ X1, df, col = X2, ylim = c(100, 500), xlim = c(100, 350), pch = 19,
xlab = "Line Speed", ylab = "Amount of Scrap")
abline(coef(fit)[[1]], coef(fit)[[3]]) # Production Line 1
abline(coef(fit)[[2]], coef(fit)[[4]], col = "red") # Production Line 2
The author mentions a residual plot against X2 and a normal probability plot (page 331), but did not include them. I have included them below. Note that 'fit' is still referring to the fitted model, but the results are the same as lm(Y ~ X1*X2, df).
plot(fitted(fit), resid(fit), col = df$X2, pch = 19,
xlim = c(100, 500), ylim = c(-40,40),
ylab = "Residual", xlab = "Fitted")
title("(a) Production Line 1")
abline(0, 0)
plot(jitter(as.numeric(df$X2), 0.1), resid(fit), pch = 19,
xlab = "Production Line", ylab = "Residual", sub = "jittered added to X2")
qqnorm(resid(fit), xlab = "Expected", ylab = "Residual", main = "")
qqline(resid(fit))
title("Normal Probability Plot")
fit <- lm(Y ~ X1 * X2, df) # No longer using the nested model
summary(fit)
##
## Call:
## lm(formula = Y ~ X1 * X2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.50 -11.06 2.78 14.82 39.51
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.5745 20.8697 0.36 0.7200
## X1 1.3220 0.0926 14.27 6.4e-13 ***
## X2Line2 90.3909 28.3457 3.19 0.0041 **
## X1:X2Line2 -0.1767 0.1288 -1.37 0.1835
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.8 on 23 degrees of freedom
## Multiple R-squared: 0.945, Adjusted R-squared: 0.937
## F-statistic: 131 on 3 and 23 DF, p-value: 1.34e-14
anova(fit)
## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 149661 149661 347.55 2.2e-15 ***
## X2 1 18694 18694 43.41 1.0e-06 ***
## X1:X2 1 810 810 1.88 0.18
## Residuals 23 9904 431
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1