# R Help 10: Model Building

R Help 10: Model Building

## R

### Martians (underspecified model)

• Fit a multiple linear regression model of weight vs height + water.
• Fit a simple linear regression model of weight vs height.
• Create a scatterplot of weight vs height with points marked by water level and regression lines for each model.

attach(martian)

model.1 <- lm(weight ~ height + water)
summary(model.1)
#              Estimate Std. Error t value Pr(>|t|)
# (Intercept) -1.220194   0.320978  -3.801  0.00421 **
# height       0.283436   0.009142  31.003 1.85e-10 ***
# water        0.111212   0.005748  19.348 1.22e-08 ***
# ---
# Residual standard error: 0.1305 on 9 degrees of freedom

model.2 <- lm(weight ~ height)
summary(model.2)
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept) -4.14335    1.75340  -2.363   0.0397 *
# height       0.38893    0.04543   8.561 6.48e-06 ***
# ---
# Residual standard error: 0.808 on 10 degrees of freedom

plot(x=height, y=weight, col=water/10+1,
panel.last = c(lines(sort(height[water==0]),
fitted(model.1)[water==0][order(height[water==0])],
col=1),
lines(sort(height[water==10]),
fitted(model.1)[water==10][order(height[water==10])],
col=2),
lines(sort(height[water==20]),
fitted(model.1)[water==20][order(height[water==20])],
col=3),
lines(sort(height), fitted(model.2)[order(height)], col=4)))
legend("topleft", col=1:4, pch=1, lty=1, inset=0.02,
legend=c("Model 1, water=0", "Model 1, water=10",
"Model 1, water=20", "Model 2"))

detach(martian)

### Cement hardening (variable selection using stepwise regression)

• Create a scatterplot matrix of the data.
• Use the add1 and drop1 functions to conduct stepwise regression.

attach(cement)

pairs(cement)

model.0 <- lm(y ~ 1)
add1(model.0, ~ x1 + x2 + x3 + x4, test="F")
# Model:
# y ~ 1
#        Df Sum of Sq     RSS    AIC F value    Pr(>F)
#               2715.76 71.444
# x1      1   1450.08 1265.69 63.519 12.6025 0.0045520 **
# x2      1   1809.43  906.34 59.178 21.9606 0.0006648 ***
# x3      1    776.36 1939.40 69.067  4.4034 0.0597623 .
# x4      1   1831.90  883.87 58.852 22.7985 0.0005762 ***

model.4 <- lm(y ~ x4)
add1(model.4, ~ . + x1 + x2 + x3, test="F")
# Model:
# y ~ x4
#        Df Sum of Sq    RSS    AIC  F value    Pr(>F)
#               883.87 58.852
# x1      1    809.10  74.76 28.742 108.2239 1.105e-06 ***
# x2      1     14.99 868.88 60.629   0.1725    0.6867
# x3      1    708.13 175.74 39.853  40.2946 8.375e-05 ***

model.14 <- lm(y ~ x1 + x4)
drop1(model.14, ~ ., test="F")
# Model:
# y ~ x1 + x4
#        Df Sum of Sq     RSS    AIC F value    Pr(>F)
#                 74.76 28.742
# x1      1     809.1  883.87 58.852  108.22 1.105e-06 ***
# x4      1    1190.9 1265.69 63.519  159.30 1.815e-07 ***

add1(model.14, ~ . + x2 + x3, test="F")
# Model:
# y ~ x1 + x4
#        Df Sum of Sq    RSS    AIC F value  Pr(>F)
#               74.762 28.742
# x2      1    26.789 47.973 24.974  5.0259 0.05169 .
# x3      1    23.926 50.836 25.728  4.2358 0.06969 .

model.124 <- lm(y ~ x1 + x2 + x4)
drop1(model.124, ~ ., test="F")
# Model:
# y ~ x4 + x1 + x2
#        Df Sum of Sq    RSS    AIC  F value    Pr(>F)
#                47.97 24.974
# x1      1    820.91 868.88 60.629 154.0076 5.781e-07 ***
# x2      1     26.79  74.76 28.742   5.0259   0.05169 .
# x4      1      9.93  57.90 25.420   1.8633   0.20540

model.12 <- lm(y ~ x1 + x2)
add1(model.12, ~ . + x3 + x4, test="F")
# Model:
# y ~ x1 + x2
#        Df Sum of Sq    RSS    AIC F value Pr(>F)
#               57.904 25.420
# x3      1    9.7939 48.111 25.011  1.8321 0.2089
# x4      1    9.9318 47.973 24.974  1.8633 0.2054

detach(cement)

### IQ and body size (variable selection using stepwise regression)

• Create a scatterplot matrix of the data.
• Use the add1 and drop1 functions to conduct stepwise regression.

attach(iqsize)

pairs(iqsize)

model.0 <- lm(PIQ ~ 1)
add1(model.0, ~ Brain + Height + Weight, test="F")
# Model:
# PIQ ~ 1
#        Df Sum of Sq   RSS    AIC F value  Pr(>F)
#               18895 237.94
# Brain   1   2697.09 16198 234.09  5.9945 0.01935 *
# Height  1    163.97 18731 239.61  0.3151 0.57802
# Weight  1      0.12 18894 239.94  0.0002 0.98806

model.1 <- lm(PIQ ~ Brain)
add1(model.1, ~ . + Height + Weight, test="F")
# Model:
# PIQ ~ Brain
#        Df Sum of Sq   RSS    AIC F value   Pr(>F)
#               16198 234.09
# Height  1   2875.65 13322 228.66  7.5551 0.009399 **
# Weight  1    940.94 15256 233.82  2.1586 0.150705

model.12 <- lm(PIQ ~ Brain + Height)
drop1(model.12, ~ ., test="F")
# Model:
# PIQ ~ Brain + Height
#        Df Sum of Sq   RSS    AIC F value    Pr(>F)
#               13322 228.66
# Brain   1    5408.8 18731 239.61 14.2103 0.0006045 ***
# Height  1    2875.6 16198 234.09  7.5551 0.0093991 **

add1(model.12, ~ . + Weight, test="F")
# Model:
# PIQ ~ Brain + Height
#        Df Sum of Sq   RSS    AIC F value Pr(>F)
#               13322 228.66
# Weight  1 0.0031633 13322 230.66       0 0.9977

detach(iqsize)

### Blood pressure (variable selection using stepwise regression)

• Create scatterplot matrices of the data.
• Use the add1 and drop1 functions to conduct stepwise regression.

attach(bloodpress)

pairs(bloodpress[,c(2:5)])
pairs(bloodpress[,c(2,6:8)])

model.0 <- lm(BP  ~ 1)
add1(model.0, ~ Age + Weight + BSA + Dur + Pulse + Stress, test="F")
# Model:
# BP ~ 1
#        Df Sum of Sq    RSS    AIC  F value    Pr(>F)
#               560.00 68.644
# Age     1    243.27 316.73 59.247  13.8248 0.0015737 **
# Weight  1    505.47  54.53 24.060 166.8591 1.528e-10 ***
# BSA     1    419.86 140.14 42.938  53.9270 8.114e-07 ***
# Dur     1     48.02 511.98 68.851   1.6883 0.2102216
# Pulse   1    291.44 268.56 55.946  19.5342 0.0003307 ***
# Stress  1     15.04 544.96 70.099   0.4969 0.4898895

model.2 <- lm(BP ~ Weight)
add1(model.2, ~ . + Age + BSA + Dur + Pulse + Stress, test="F")
# Model:
# BP ~ Weight
#        Df Sum of Sq    RSS     AIC  F value    Pr(>F)
#               54.528  24.060
# Age     1    49.704  4.824 -22.443 175.1622 2.218e-10 ***
# BSA     1     2.814 51.714  25.000   0.9251   0.34962
# Dur     1     6.095 48.433  23.689   2.1393   0.16181
# Pulse   1     8.940 45.588  22.478   3.3338   0.08549 .
# Stress  1     9.660 44.868  22.160   3.6601   0.07273 .

model.12 <- lm(BP ~ Age + Weight)
drop1(model.12, ~ ., test="F")
# Model:
# BP ~ Age + Weight
#        Df Sum of Sq    RSS     AIC F value    Pr(>F)
#                 4.82 -22.443
# Age     1    49.704  54.53  24.060  175.16 2.218e-10 ***
# Weight  1   311.910 316.73  59.247 1099.20 < 2.2e-16 ***

add1(model.12, ~ . + BSA + Dur + Pulse + Stress, test="F")
# Model:
# BP ~ Age + Weight
#        Df Sum of Sq    RSS     AIC F value   Pr(>F)
#               4.8239 -22.443
# BSA     1   1.76778 3.0561 -29.572  9.2550 0.007764 **
# Dur     1   0.17835 4.6456 -21.196  0.6143 0.444639
# Pulse   1   0.49557 4.3284 -22.611  1.8319 0.194719
# Stress  1   0.16286 4.6611 -21.130  0.5591 0.465486

model.123 <- lm(BP ~ Age + Weight + BSA)
drop1(model.123, ~ ., test="F")
# Model:
# BP ~ Age + Weight + BSA
#        Df Sum of Sq    RSS     AIC F value    Pr(>F)
#                3.056 -29.572
# Age     1    48.658 51.714  25.000 254.740 3.002e-11 ***
# Weight  1    65.303 68.359  30.581 341.886 3.198e-12 ***
# BSA     1     1.768  4.824 -22.443   9.255  0.007764 **

add1(model.123, ~ . + Dur + Pulse + Stress, test="F")
# Model:
# BP ~ Age + Weight + BSA
#        Df Sum of Sq    RSS     AIC F value Pr(>F)
#               3.0561 -29.572
# Dur     1   0.33510 2.7210 -29.894  1.8473 0.1942
# Pulse   1   0.04111 3.0150 -27.842  0.2045 0.6576
# Stress  1   0.21774 2.8384 -29.050  1.1507 0.3004

detach(bloodpress)

### Cement hardening (variable selection using best subsets regression)

• Use the regsubsets function in the leaps package to conduct variable selection using exhaustive search (i.e., best subsets regression). Note that the nbest=2 argument returns the best two models with 1, 2, ..., k predictors.
• Fit models with all four predictors (assumed unbiased) and just two predictors to retrieve the information needed to calculate $$C_p$$ for the model with just two predictors by hand.
• Fit model with $$x_1$$, $$x_2$$, and $$x_4$$ and note the variance inflation factors for $$x_2$$ and $$x_4$$ are very high.
• Fit model with $$x_1$$, $$x_2$$, and $$x_3$$ and note the variance inflation factors are acceptable.
• Fit model with $$x_1$$ and $$x_2$$ and note the variance inflation factors are acceptable, adjusted $$R^2$$ is high, and a residual analysis and normality test yields no concerns.

attach(cement)

library(leaps)

subset <- regsubsets(y ~ x1 + x2 + x3 + x4, method="exhaustive", nbest=2, data=cement)
cbind(summary(subset)$outmat, round(summary(subset)$adjr2, 3), round(summary(subset)$cp, 1)) # x1 x2 x3 x4 # 1 ( 1 ) " " " " " " "*" "0.645" "138.7" # 1 ( 2 ) " " "*" " " " " "0.636" "142.5" # 2 ( 1 ) "*" "*" " " " " "0.974" "2.7" # 2 ( 2 ) "*" " " " " "*" "0.967" "5.5" # 3 ( 1 ) "*" "*" " " "*" "0.976" "3" # 3 ( 2 ) "*" "*" "*" " " "0.976" "3" # 4 ( 1 ) "*" "*" "*" "*" "0.974" "5" model.1234 <- lm(y ~ x1 + x2 + x3 + x4) model.12 <- lm(y ~ x1 + x2) SSE.k <- sum(residuals(model.12)^2) # SSE_k = 57.90448 MSE.all <- summary(model.1234)$sigma^2 # MSE_all = 5.982955
params <- summary(model.12)$df[1] # k+1 = 3 n <- sum(summary(model.1234)$df[1:2]) # n = 13
SSE.k/MSE.all + 2*params - n # Cp = 2.678242

model.14 <- lm(y ~ x1 + x4)

SSE.k <- sum(residuals(model.14)^2) # SSE_k = 74.76211
params <- summary(model.14)$df[1] # k+1 = 3 SSE.k/MSE.all + 2*params - n # Cp = 5.495851 model.124 <- lm(y ~ x1 + x2 + x4) library(car) vif(model.124) # x1 x2 x4 # 1.06633 18.78031 18.94008 model.123 <- lm(y ~ x1 + x2 + x3) vif(model.123) # x1 x2 x3 # 3.251068 1.063575 3.142125 summary(model.12) # Estimate Std. Error t value Pr(>|t|) # (Intercept) 52.57735 2.28617 23.00 5.46e-10 *** # x1 1.46831 0.12130 12.11 2.69e-07 *** # x2 0.66225 0.04585 14.44 5.03e-08 *** # --- # Residual standard error: 2.406 on 10 degrees of freedom # Multiple R-squared: 0.9787, Adjusted R-squared: 0.9744 # F-statistic: 229.5 on 2 and 10 DF, p-value: 4.407e-09 vif(model.12) # x1 x2 # 1.055129 1.055129 plot(x=fitted(model.12), y=rstandard(model.12), panel.last = abline(h=0, lty=2)) qqnorm(rstandard(model.12), main="", datax=TRUE) qqline(rstandard(model.12), datax=TRUE) library(nortest) ad.test(rstandard(model.12)) # A = 0.6136, p-value = 0.08628 detach(cement) ### IQ and body size (variable selection using best subsets regression) • Load the iqsize data. • Use the regsubsets function in the leaps package to conduct variable selection using exhaustive search (i.e., best subsets regression). • Fit model with Brain and Height and note the variance inflation factors are acceptable, adjusted $$R^2$$ is as good as it gets with this dataset, and a residual analysis and normality test yields no concerns. iqsize <- read.table("~/path-to-data/iqsize.txt", header=T) attach(iqsize) subset <- regsubsets(PIQ ~ Brain + Height + Weight, method="exhaustive", nbest=2, data=iqsize) cbind(summary(subset)$outmat, round(summary(subset)$adjr2, 3), round(summary(subset)$cp, 1))
#          Brain Height Weight
# 1  ( 1 ) "*"   " "    " "    "0.119"  "7.3"
# 1  ( 2 ) " "   "*"    " "    "-0.019" "13.8"
# 2  ( 1 ) "*"   "*"    " "    "0.255"  "2"
# 2  ( 2 ) "*"   " "    "*"    "0.146"  "6.9"
# 3  ( 1 ) "*"   "*"    "*"    "0.233"  "4"

model.12 <- lm(PIQ ~ Brain + Height)
summary(model.12)
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept) 111.2757    55.8673   1.992 0.054243 .
# Brain         2.0606     0.5466   3.770 0.000604 ***
# Height       -2.7299     0.9932  -2.749 0.009399 **
# ---
# Residual standard error: 19.51 on 35 degrees of freedom
# Multiple R-squared:  0.2949,  Adjusted R-squared:  0.2546
# F-statistic: 7.321 on 2 and 35 DF,  p-value: 0.002208

vif(model.12)
#    Brain   Height
# 1.529463 1.529463

plot(x=fitted(model.12), y=rstandard(model.12),
panel.last = abline(h=0, lty=2))

qqnorm(rstandard(model.12), main="", datax=TRUE)
qqline(rstandard(model.12), datax=TRUE)

ad.test(rstandard(model.12)) # A = 0.2629, p-value = 0.6829

detach(iqsize)

### Blood pressure (variable selection using best subsets regression)

• Use the regsubsets function in the leaps package to conduct variable selection using exhaustive search (i.e., best subsets regression).
• Fit model with Age and Weight and note the variance inflation factors are acceptable, adjusted $$R^2$$ can't get much better, and a residual analysis and normality test yields no concerns.

attach(bloodpress)

subset <- regsubsets(BP ~ Age + Weight + BSA + Dur + Pulse + Stress,
method="exhaustive", nbest=2, data=bloodpress)
cbind(summary(subset)$outmat, round(summary(subset)$adjr2, 3),
round(summary(subset)$cp, 1)) # Age Weight BSA Dur Pulse Stress # 1 ( 1 ) " " "*" " " " " " " " " "0.897" "312.8" # 1 ( 2 ) " " " " "*" " " " " " " "0.736" "829.1" # 2 ( 1 ) "*" "*" " " " " " " " " "0.99" "15.1" # 2 ( 2 ) " " "*" " " " " " " "*" "0.91" "256.6" # 3 ( 1 ) "*" "*" "*" " " " " " " "0.994" "6.4" # 3 ( 2 ) "*" "*" " " " " "*" " " "0.991" "14.1" # 4 ( 1 ) "*" "*" "*" "*" " " " " "0.994" "6.4" # 4 ( 2 ) "*" "*" "*" " " " " "*" "0.994" "7.1" # 5 ( 1 ) "*" "*" "*" " " "*" "*" "0.994" "7" # 5 ( 2 ) "*" "*" "*" "*" "*" " " "0.994" "7.7" # 6 ( 1 ) "*" "*" "*" "*" "*" "*" "0.994" "7" model.12 <- lm(BP ~ Age + Weight) summary(model.12) # Estimate Std. Error t value Pr(>|t|) # (Intercept) -16.57937 3.00746 -5.513 3.80e-05 *** # Age 0.70825 0.05351 13.235 2.22e-10 *** # Weight 1.03296 0.03116 33.154 < 2e-16 *** # --- # Residual standard error: 0.5327 on 17 degrees of freedom # Multiple R-squared: 0.9914, Adjusted R-squared: 0.9904 # F-statistic: 978.2 on 2 and 17 DF, p-value: < 2.2e-16 vif(model.12) # Age Weight # 1.198945 1.198945 plot(x=fitted(model.12), y=rstandard(model.12), panel.last = abline(h=0, lty=2)) qqnorm(rstandard(model.12), main="", datax=TRUE) qqline(rstandard(model.12), datax=TRUE) ad.test(rstandard(model.12)) # A = 0.275, p-value = 0.6225 detach(bloodpress) ### Peruvian blood pressure (variable selection using best subsets regression) • Load the peru data. • Use the regsubsets function in the leaps package to conduct variable selection using exhaustive search (i.e., best subsets regression). • Fit the best 5-predictor and 4-predictor models. • Calculate AIC and BIC by hand. • Use the stepAIC function in the MASS package to conduct variable selection using a stepwise algorithm based on AIC or BIC. peru <- read.table("~/path-to-data/peru.txt", header=T) attach(peru) fraclife <- Years/Age n <- length(Systol) # 39 subset <- regsubsets(Systol ~ Age + Years + fraclife + Weight + Height + Chin + Forearm + Pulse, method="exhaustive", nbest=2, data=peru) cbind(summary(subset)$outmat, round(summary(subset)$rsq, 3), round(summary(subset)$adjr2, 3), round(summary(subset)$cp, 1), round(sqrt(summary(subset)$rss/(n-c(rep(1:7,rep(2,7)),8)-1)), 4))
#          Age Years fraclife Weight Height Chin Forearm Pulse
# 1  ( 1 ) " " " "   " "      "*"    " "    " "  " "     " "   "0.272" "0.252" "30.5" "11.3376"
# 1  ( 2 ) " " " "   "*"      " "    " "    " "  " "     " "   "0.076" "0.051" "48.1" "12.7697"
# 2  ( 1 ) " " " "   "*"      "*"    " "    " "  " "     " "   "0.473" "0.444" "14.4" "9.7772"
# 2  ( 2 ) " " "*"   " "      "*"    " "    " "  " "     " "   "0.421" "0.389" "19.1" "10.2512"
# 3  ( 1 ) " " " "   "*"      "*"    " "    "*"  " "     " "   "0.503" "0.461" "13.7" "9.6273"
# 3  ( 2 ) " " "*"   "*"      "*"    " "    " "  " "     " "   "0.49"  "0.447" "14.8" "9.7509"
# 4  ( 1 ) "*" "*"   "*"      "*"    " "    " "  " "     " "   "0.597" "0.55"  "7.2"  "8.7946"
# 4  ( 2 ) "*" "*"   "*"      " "    "*"    " "  " "     " "   "0.525" "0.469" "13.7" "9.5502"
# 5  ( 1 ) "*" "*"   "*"      "*"    " "    "*"  " "     " "   "0.639" "0.584" "5.5"  "8.4571"
# 5  ( 2 ) "*" "*"   "*"      "*"    " "    " "  "*"     " "   "0.631" "0.576" "6.1"  "8.5417"
# 6  ( 1 ) "*" "*"   "*"      "*"    " "    "*"  "*"     " "   "0.649" "0.583" "6.6"  "8.4663"
# 6  ( 2 ) "*" "*"   "*"      "*"    "*"    "*"  " "     " "   "0.643" "0.576" "7.1"  "8.5337"
# 7  ( 1 ) "*" "*"   "*"      "*"    "*"    "*"  "*"     " "   "0.661" "0.584" "7.5"  "8.4556"
# 7  ( 2 ) "*" "*"   "*"      "*"    " "    "*"  "*"     "*"   "0.655" "0.577" "8"    "8.522"
# 8  ( 1 ) "*" "*"   "*"      "*"    "*"    "*"  "*"     "*"   "0.666" "0.577" "9"    "8.5228"

model.5 <- lm(Systol ~ Age + Years + fraclife + Weight + Chin)
summary(model.5)
#              Estimate Std. Error t value Pr(>|t|)
# (Intercept)  109.3590    21.4843   5.090 1.41e-05 ***
# Age           -1.0120     0.3059  -3.308 0.002277 **
# Years          2.4067     0.7426   3.241 0.002723 **
# fraclife    -110.8112    27.2795  -4.062 0.000282 ***
# Weight         1.0976     0.2980   3.683 0.000819 ***
# Chin          -1.1918     0.6140  -1.941 0.060830 .
# ---
# Residual standard error: 8.457 on 33 degrees of freedom
# Multiple R-squared:  0.6386,  Adjusted R-squared:  0.5839
# F-statistic: 11.66 on 5 and 33 DF,  p-value: 1.531e-06

k <- 5
n*log(sum(residuals(model.5)^2))-n*log(n)+2*(k+1) # AIC = 172.0151
n*log(sum(residuals(model.5)^2))-n*log(n)+log(n)*(k+1) # BIC = 181.9965

model.4 <- lm(Systol ~ Age + Years + fraclife + Weight)
summary(model.4)
#              Estimate Std. Error t value Pr(>|t|)
# (Intercept)  116.8354    21.9797   5.316 6.69e-06 ***
# Age           -0.9507     0.3164  -3.004 0.004971 **
# Years          2.3393     0.7714   3.032 0.004621 **
# fraclife    -108.0728    28.3302  -3.815 0.000549 ***
# Weight         0.8324     0.2754   3.022 0.004742 **
# ---
# Residual standard error: 8.795 on 34 degrees of freedom
# Multiple R-squared:  0.5974,  Adjusted R-squared:   0.55
# F-statistic: 12.61 on 4 and 34 DF,  p-value: 2.142e-06

k <- 4
n*log(sum(residuals(model.4)^2))-n*log(n)+2*(k+1) # AIC = 174.2316
n*log(sum(residuals(model.4)^2))-n*log(n)+log(n)*(k+1) # BIC = 182.5494

library(MASS)
subset.aic <- stepAIC(lm(Systol ~ Age + Years + fraclife + Weight + Height +
Chin + Forearm + Pulse), direction="both", k=2)
# Step:  AIC=172.02
# Systol ~ Age + Years + fraclife + Weight + Chin

subset.bic <- stepAIC(lm(Systol ~ Age + Years + fraclife + Weight + Height +
Chin + Forearm + Pulse), direction="both", k=log(n))
# Step:  AIC=182
# Systol ~ Age + Years + fraclife + Weight + Chin

detach(peru)

### Measurements of college students (variable selection using stepwise regression)

• Use the add1 and drop1 functions to conduct stepwise regression.
• Use the regsubsets function to conduct variable selection using backward elimination.
• Use the regsubsets function to conduct variable selection using forward selection.

attach(physical)

gender <- ifelse(Sex=="Female",1,0)

model.0 <- lm(Height ~ 1)
add1(model.0, ~ LeftArm + LeftFoot + LeftHand + HeadCirc + nose + gender, test="F")
# Model:
# Height ~ 1
#          Df Sum of Sq     RSS    AIC  F value    Pr(>F)
#                 1054.75 164.46
# LeftArm   1    590.21  464.53 121.35  67.3396 5.252e-11 ***
# LeftFoot  1    707.42  347.33 105.36 107.9484 2.172e-14 ***
# LeftHand  1    143.59  911.15 158.41   8.3525  0.005570 **
# HeadCirc  1    189.24  865.51 155.58  11.5880  0.001272 **
# nose      1     85.25  969.49 161.82   4.6605  0.035412 *
# gender    1    533.24  521.51 127.72  54.1923 1.181e-09 ***

model.2 <- lm(Height ~ LeftFoot)
add1(model.2, ~ . + LeftArm + LeftHand + HeadCirc + nose + gender, test="F")
# Model:
# Height ~ LeftFoot
#         Df Sum of Sq    RSS     AIC F value    Pr(>F)
#                 347.33 105.361
# LeftArm   1   107.143 240.18  87.074 23.1967 1.305e-05 ***
# LeftHand  1    15.359 331.97 104.874  2.4059    0.1269
# HeadCirc  1     2.313 345.01 106.994  0.3486    0.5575
# nose      1     1.449 345.88 107.131  0.2178    0.6427
# gender    1    15.973 331.35 104.772  2.5066    0.1194

model.12 <- lm(Height ~ LeftArm + LeftFoot)
drop1(model.12, ~ ., test="F")
# Model:
# Height ~ LeftArm + LeftFoot
#          Df Sum of Sq    RSS     AIC F value    Pr(>F)
#                 240.18  87.074
# LeftArm   1    107.14 347.33 105.361  23.197 1.305e-05 ***
# LeftFoot  1    224.35 464.53 121.353  48.572 5.538e-09 ***

add1(model.12, ~ . + LeftHand + HeadCirc + nose + gender, test="F")
# Model:
# Height ~ LeftArm + LeftFoot
#          Df Sum of Sq    RSS    AIC F value Pr(>F)
#                 240.18 87.074
# LeftHand  1    3.7854 236.40 88.200  0.8167 0.3704
# HeadCirc  1    1.4016 238.78 88.752  0.2994 0.5867
# nose      1    0.4463 239.74 88.971  0.0950 0.7592
# gender    1    3.7530 236.43 88.207  0.8096 0.3725

subset <- regsubsets(Height ~ LeftArm + LeftFoot + LeftHand + HeadCirc + nose + gender,
method="backward", data=physical)

subset <- regsubsets(Height ~ LeftArm + LeftFoot + LeftHand + HeadCirc + nose + gender,
method="forward", data=physical)

detach(physical)

 [1] Link ↥ Has Tooltip/Popover Toggleable Visibility