# 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