# Reference: https://onlinecourses.science.psu.edu/stat501/node/347
readr::read_delim("http://onlinecourses.science.psu.edu/stat501/sites/onlinecourses.science.psu.edu.stat501/files/data/bloodpress.txt", delim="\t") -> univ
dplyr::glimpse(univ)
Observations: 20
Variables: 8
$ Pt     (int) 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
$ BP     (int) 105, 115, 116, 117, 112, 121, 121, 110, 110, 114, 114, 115, 114, 106, 125, 114, 106,...
$ Age    (int) 47, 49, 49, 50, 51, 48, 49, 47, 49, 48, 47, 49, 50, 45, 52, 46, 46, 46, 48, 56
$ Weight (dbl) 85.4, 94.2, 95.3, 94.7, 89.4, 99.5, 99.8, 90.9, 89.2, 92.7, 94.4, 94.1, 91.6, 87.1, ...
$ BSA    (dbl) 1.75, 2.10, 1.98, 2.01, 1.89, 2.25, 2.25, 1.90, 1.83, 2.07, 2.07, 1.98, 2.05, 1.92, ...
$ Dur    (dbl) 5.1, 3.8, 8.2, 5.8, 7.0, 9.3, 2.5, 6.2, 7.1, 5.6, 5.3, 5.6, 10.2, 5.6, 10.0, 7.4, 3....
$ Pulse  (int) 63, 70, 72, 73, 72, 71, 69, 66, 69, 64, 74, 71, 68, 67, 76, 69, 62, 70, 71, 75
$ Stress (int) 33, 14, 10, 99, 95, 10, 42, 8, 62, 35, 90, 21, 47, 80, 98, 95, 18, 12, 99, 99
mydiag <- function(vec, pch=1, bg=1, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5))
    h <- hist(vec, plot=FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col=rgb(.4,.6,.9), ...)
}
mycor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r^(4/10))
}
# las=1: Axis label style (1: easy to read)
par(adj=c(0,.5,1)[2], las=1, bg="black", col = "white", col.axis = "white", fg="white", col.sub="white", col.main="white", col.lab="white", lwd=1.5, pch=18)
pairs(univ,  diag.panel = mydiag, upper.panel = mycor)

# BSA: body surface area
lmed6 <- lm(BP ~ Weight + BSA + Age + Dur + Stress + Pulse, data = univ)
lmed4 <- lm(BP ~ Weight       + Age + Dur + Stress        , data = univ)
lmed3 <- lm(BP ~ Weight + BSA + Age                       , data = univ)
lmed2 <- lm(BP ~ Weight + BSA                             , data = univ)
lmed1 <- lm(BP ~ Weight                                   , data = univ)
lmdf <-
data_frame(
    lm6 = lmed6,
    lm4 = lmed4,
    lm3 = lmed3,
    lm2 = lmed2,
    lm1 = lmed1
) %>% t %>% data.frame
colnames(lmdf) <- names(lmed)
lmdf
lmdf$coefficients
$lm6
 (Intercept)       Weight          BSA          Age          Dur       Stress        Pulse 
-12.87047602   0.96991978   3.77649100   0.70325939   0.06838304   0.00557150  -0.08448469 

$lm4
  (Intercept)        Weight           Age           Dur        Stress 
-15.869828678   1.034128311   0.683740959   0.039888843   0.002184223 

$lm3
(Intercept)      Weight         BSA         Age 
-13.6672473   0.9058223   4.6273883   0.7016198 

$lm2
(Intercept)      Weight         BSA 
   5.653398    1.038734    5.831250 

$lm1
(Intercept)      Weight 
   2.205305    1.200931 
plot(1:20, type='n', ylim=c(-3,5))
i <- 1 # global
aa <- sapply(lmdf$residuals, function(x) { lines(x, col=colors()[25+i], pch=as.character(1:10)[i], type='b', lwd=2, lty=3); i <<- i+1 })
abline(a=0,b=0, lty=2)

Patient 20 is the eldest. Age affects Blood Pressure.

library(DAAG)
lapply(list(lmed1, lmed2, lmed3, lmed4, lmed6), vif) # lmed4 is good -- all around 1
[[1]]
Weight 
     1 

[[2]]
Weight    BSA 
4.2764 4.2764 

[[3]]
Weight    BSA    Age 
4.4036 4.2869 1.2019 

[[4]]
Weight    Age    Dur Stress 
1.2347 1.4682 1.2001 1.2411 

[[5]]
Weight    BSA    Age    Dur Stress  Pulse 
8.4170 5.3288 1.7628 1.2373 1.8348 4.4136 
lapply(list(lmed1, lmed2, lmed3, lmed4, lmed6), summary)
[[1]]

Call:
lm(formula = BP ~ Weight, data = univ)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.6933 -0.9318 -0.4935  0.7703  4.8656 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.20531    8.66333   0.255    0.802    
Weight       1.20093    0.09297  12.917 1.53e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.74 on 18 degrees of freedom
Multiple R-squared:  0.9026,    Adjusted R-squared:  0.8972 
F-statistic: 166.9 on 1 and 18 DF,  p-value: 1.528e-10


[[2]]

Call:
lm(formula = BP ~ Weight + BSA, data = univ)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.8932 -1.1961 -0.4061  1.0764  4.7524 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.6534     9.3925   0.602    0.555    
Weight        1.0387     0.1927   5.392 4.87e-05 ***
BSA           5.8313     6.0627   0.962    0.350    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.744 on 17 degrees of freedom
Multiple R-squared:  0.9077,    Adjusted R-squared:  0.8968 
F-statistic: 83.54 on 2 and 17 DF,  p-value: 1.607e-09


[[3]]

Call:
lm(formula = BP ~ Weight + BSA + Age, data = univ)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.75810 -0.24872  0.01925  0.29509  0.63030 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -13.66725    2.64664  -5.164 9.42e-05 ***
Weight        0.90582    0.04899  18.490 3.20e-12 ***
BSA           4.62739    1.52107   3.042  0.00776 ** 
Age           0.70162    0.04396  15.961 3.00e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.437 on 16 degrees of freedom
Multiple R-squared:  0.9945,    Adjusted R-squared:  0.9935 
F-statistic: 971.9 on 3 and 16 DF,  p-value: < 2.2e-16


[[4]]

Call:
lm(formula = BP ~ Weight + Age + Dur + Stress, data = univ)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.11359 -0.29586  0.01515  0.27506  0.88674 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -15.869829   3.195296  -4.967 0.000169 ***
Weight        1.034128   0.032672  31.652 3.76e-15 ***
Age           0.683741   0.061195  11.173 1.14e-08 ***
Dur           0.039889   0.064486   0.619 0.545485    
Stress        0.002184   0.003794   0.576 0.573304    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5505 on 15 degrees of freedom
Multiple R-squared:  0.9919,    Adjusted R-squared:  0.9897 
F-statistic: 458.3 on 4 and 15 DF,  p-value: 1.764e-15


[[5]]

Call:
lm(formula = BP ~ Weight + BSA + Age + Dur + Stress + Pulse, 
    data = univ)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.93213 -0.11314  0.03064  0.21834  0.48454 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -12.870476   2.556650  -5.034 0.000229 ***
Weight        0.969920   0.063108  15.369 1.02e-09 ***
BSA           3.776491   1.580151   2.390 0.032694 *  
Age           0.703259   0.049606  14.177 2.76e-09 ***
Dur           0.068383   0.048441   1.412 0.181534    
Stress        0.005572   0.003412   1.633 0.126491    
Pulse        -0.084485   0.051609  -1.637 0.125594    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4072 on 13 degrees of freedom
Multiple R-squared:  0.9962,    Adjusted R-squared:  0.9944 
F-statistic: 560.6 on 6 and 13 DF,  p-value: 6.395e-15
LS0tCnRpdGxlOiAiVmFyaWFuY2UgSW5mbGF0aW9uIEZhY3RvciIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKCmBgYHtyfQojIFJlZmVyZW5jZTogaHR0cHM6Ly9vbmxpbmVjb3Vyc2VzLnNjaWVuY2UucHN1LmVkdS9zdGF0NTAxL25vZGUvMzQ3CnJlYWRyOjpyZWFkX2RlbGltKCJodHRwOi8vb25saW5lY291cnNlcy5zY2llbmNlLnBzdS5lZHUvc3RhdDUwMS9zaXRlcy9vbmxpbmVjb3Vyc2VzLnNjaWVuY2UucHN1LmVkdS5zdGF0NTAxL2ZpbGVzL2RhdGEvYmxvb2RwcmVzcy50eHQiLCBkZWxpbT0iXHQiKSAtPiB1bml2CmRwbHlyOjpnbGltcHNlKHVuaXYpCmBgYAoKYGBge3IgZmlnLndpZHRoPTE4fQpteWRpYWcgPC0gZnVuY3Rpb24odmVjLCBwY2g9MSwgYmc9MSwgLi4uKQp7CiAgICB1c3IgPC0gcGFyKCJ1c3IiKTsgb24uZXhpdChwYXIodXNyKSkKICAgIHBhcih1c3IgPSBjKHVzclsxOjJdLCAwLCAxLjUpKQogICAgaCA8LSBoaXN0KHZlYywgcGxvdD1GQUxTRSkKICAgIGJyZWFrcyA8LSBoJGJyZWFrczsgbkIgPC0gbGVuZ3RoKGJyZWFrcykKICAgIHkgPC0gaCRjb3VudHM7IHkgPC0geS9tYXgoeSkKICAgIHJlY3QoYnJlYWtzWy1uQl0sIDAsIGJyZWFrc1stMV0sIHksIGNvbD1yZ2IoLjQsLjYsLjkpLCAuLi4pCn0KCm15Y29yIDwtIGZ1bmN0aW9uKHgsIHksIGRpZ2l0cyA9IDIsIHByZWZpeCA9ICIiLCBjZXguY29yLCAuLi4pCnsKICAgIHVzciA8LSBwYXIoInVzciIpOyBvbi5leGl0KHBhcih1c3IpKQogICAgcGFyKHVzciA9IGMoMCwgMSwgMCwgMSkpCiAgICByIDwtIGFicyhjb3IoeCwgeSkpCiAgICB0eHQgPC0gZm9ybWF0KGMociwgMC4xMjM0NTY3ODkpLCBkaWdpdHMgPSBkaWdpdHMpWzFdCiAgICB0eHQgPC0gcGFzdGUwKHByZWZpeCwgdHh0KQogICAgaWYobWlzc2luZyhjZXguY29yKSkgY2V4LmNvciA8LSAwLjgvc3Ryd2lkdGgodHh0KQogICAgdGV4dCgwLjUsIDAuNSwgdHh0LCBjZXggPSBjZXguY29yICogcl4oNC8xMCkpCn0KIyBsYXM9MTogQXhpcyBsYWJlbCBzdHlsZSAoMTogZWFzeSB0byByZWFkKQpwYXIoYWRqPWMoMCwuNSwxKVsyXSwgbGFzPTEsIGJnPSJibGFjayIsIGNvbCA9ICJ3aGl0ZSIsIGNvbC5heGlzID0gIndoaXRlIiwgZmc9IndoaXRlIiwgY29sLnN1Yj0id2hpdGUiLCBjb2wubWFpbj0id2hpdGUiLCBjb2wubGFiPSJ3aGl0ZSIsIGx3ZD0xLjUsIHBjaD0xOCkKcGFpcnModW5pdiwgIGRpYWcucGFuZWwgPSBteWRpYWcsIHVwcGVyLnBhbmVsID0gbXljb3IpCiMgQlNBOiBib2R5IHN1cmZhY2UgYXJlYQpgYGAKCmBgYHtyfQpsbWVkNiA8LSBsbShCUCB+IFdlaWdodCArIEJTQSArIEFnZSArIER1ciArIFN0cmVzcyArIFB1bHNlLCBkYXRhID0gdW5pdikKbG1lZDQgPC0gbG0oQlAgfiBXZWlnaHQgICAgICAgKyBBZ2UgKyBEdXIgKyBTdHJlc3MgICAgICAgICwgZGF0YSA9IHVuaXYpCmxtZWQzIDwtIGxtKEJQIH4gV2VpZ2h0ICsgQlNBICsgQWdlICAgICAgICAgICAgICAgICAgICAgICAsIGRhdGEgPSB1bml2KQpsbWVkMiA8LSBsbShCUCB+IFdlaWdodCArIEJTQSAgICAgICAgICAgICAgICAgICAgICAgICAgICAgLCBkYXRhID0gdW5pdikKbG1lZDEgPC0gbG0oQlAgfiBXZWlnaHQgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICwgZGF0YSA9IHVuaXYpCgpsbWRmIDwtCmRhdGFfZnJhbWUoCiAgICBsbTYgPSBsbWVkNiwKICAgIGxtNCA9IGxtZWQ0LAogICAgbG0zID0gbG1lZDMsCiAgICBsbTIgPSBsbWVkMiwKICAgIGxtMSA9IGxtZWQxCikgJT4lIHQgJT4lIGRhdGEuZnJhbWUKY29sbmFtZXMobG1kZikgPC0gbmFtZXMobG1lZCkKbG1kZgpgYGAKCmBgYHtyfQpsbWRmJGNvZWZmaWNpZW50cwpgYGAKCmBgYHtyIGZpZy53aWR0aD0xNCwgZmlnLmhlaWdodD0xMH0KcGxvdCgxOjIwLCB0eXBlPSduJywgeWxpbT1jKC0zLDUpKQppIDwtIDEgIyBnbG9iYWwKYWEgPC0gc2FwcGx5KGxtZGYkcmVzaWR1YWxzLCBmdW5jdGlvbih4KSB7IGxpbmVzKHgsIGNvbD1jb2xvcnMoKVsyNStpXSwgcGNoPWFzLmNoYXJhY3RlcigxOjEwKVtpXSwgdHlwZT0nYicsIGx3ZD0yLCBsdHk9Myk7IGkgPDwtIGkrMSB9KQphYmxpbmUoYT0wLGI9MCwgbHR5PTIpCmBgYAoKUGF0aWVudCAyMCBpcyB0aGUgZWxkZXN0LiAqQWdlKiBhZmZlY3RzIEJsb29kIFByZXNzdXJlLiAKCmBgYHtyfQpsaWJyYXJ5KERBQUcpCmxhcHBseShsaXN0KGxtZWQxLCBsbWVkMiwgbG1lZDMsIGxtZWQ0LCBsbWVkNiksIHZpZikgIyBsbWVkNCBpcyBnb29kIC0tIGFsbCBhcm91bmQgMQpgYGAKCmBgYHtyfQpsYXBwbHkobGlzdChsbWVkMSwgbG1lZDIsIGxtZWQzLCBsbWVkNCwgbG1lZDYpLCBzdW1tYXJ5KQpgYGAKCgo=