## Settings for RMarkdown http://yihui.name/knitr/options#chunk_options
opts_chunk$set(comment = "", warning = FALSE, message = FALSE, tidy = FALSE,
echo = T, fig.width = 5, fig.height = 5)
options(width = 116, scipen = 10)
setwd("~/statistics/bio232/")
## input the data
CF <- read.table("CF.dat", header = TRUE)
## Check first 6 rows
head(CF)
Age Sex Height Weight BMP FEV RV FRC TLC PEmax
1 7 0 109 13.1 68 32 258 183 137 95
2 7 1 112 12.9 65 19 449 245 134 85
3 8 0 124 14.1 64 22 441 268 147 100
4 8 1 125 16.2 67 41 234 146 124 85
5 8 0 127 21.5 93 52 202 131 104 95
6 9 0 130 17.5 68 44 308 155 118 80
## Check dimentions
dim(CF)
[1] 25 10
## Use formula for plotting
plot(PEmax ~ FEV , data = CF)
## fit a linear model
fit.fev <- lm(PEmax ~ FEV, data = CF)
summary(fit.fev)
Call:
lm(formula = PEmax ~ FEV, data = CF)
Residuals:
Min 1Q Median 3Q Max
-43.04 -24.61 -4.91 23.52 62.49
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 62.114 20.208 3.07 0.0054 **
FEV 1.354 0.555 2.44 0.0228 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 30.4 on 23 degrees of freedom
Multiple R-squared: 0.206, Adjusted R-squared: 0.171
F-statistic: 5.95 on 1 and 23 DF, p-value: 0.0228
## 95% confidence interval for those with FEV = 26
predict(fit.fev, newdata = data.frame(FEV = 26), interval = "confidence")
fit lwr upr
1 97.31 81.22 113.4
## 95% CI for everybody
predict(fit.fev, newdata = CF["FEV"], interval = "confidence")
fit lwr upr
1 105.44 92.46 118.4
2 87.84 65.83 109.8
3 91.90 72.61 111.2
4 117.62 103.11 132.1
5 132.51 109.02 156.0
6 121.68 105.19 138.2
7 100.02 85.25 114.8
8 86.48 63.52 109.4
9 94.61 77.00 112.2
10 93.25 74.82 111.7
11 114.91 101.39 128.4
12 97.31 81.22 113.4
13 123.04 105.78 140.3
14 123.04 105.78 140.3
15 104.08 90.78 117.4
16 101.38 87.17 115.6
17 128.45 107.78 149.1
18 101.38 87.17 115.6
19 113.56 100.41 126.7
20 90.54 70.38 110.7
21 112.21 99.34 125.1
22 108.15 95.52 120.8
23 139.28 110.77 167.8
24 106.79 94.04 119.5
25 132.51 109.02 156.0
## 95% prediction interval for those with FEV = 26
predict(fit.fev, newdata = data.frame(FEV = 26), interval = "prediction")
fit lwr upr
1 97.31 32.31 162.3
## Make dummy dataset for 95% prediction interval
FEV.dummy <- with(CF, seq(min(FEV), max(FEV), 0.1))
pred.limits <-
predict(fit.fev, newdata = data.frame(FEV = FEV.dummy), interval = "prediction")
pred.limits <- data.frame(FEV = FEV.dummy, pred.limits)
## Plot with ggplot2
library(ggplot2)
ggplot(CF, aes(x = FEV, y = PEmax)) +
geom_point() +
geom_smooth(method = "lm") +
geom_line(data = pred.limits, aes(x = FEV, y = lwr), linetype = 3) +
geom_line(data = pred.limits, aes(x = FEV, y = upr), linetype = 3)
