Problem statement
Complete the regression analysis of the blood pressure data from the Framingham Heart Study in the R script . In there a gender difference in regression slopes?
Data management
# file location
fL<-"https://math.montana.edu/shancock/data/Framingham.txt"
#
dta <- read.table(fL, header=T)
# inspect data structure
str(dta)
## 'data.frame': 4699 obs. of 10 variables:
## $ sex : int 1 1 1 1 1 2 1 1 1 1 ...
## $ sbp : int 120 130 144 92 162 212 140 174 142 115 ...
## $ dbp : int 80 78 90 66 98 118 85 102 94 70 ...
## $ scl : int 267 192 207 231 271 182 276 259 242 242 ...
## $ chdfate : int 1 1 1 1 1 1 1 1 1 1 ...
## $ followup: int 18 35 109 147 169 199 201 209 265 278 ...
## $ age : int 55 53 61 48 39 61 44 39 47 60 ...
## $ bmi : num 25 28.4 25.1 26.2 28.4 ...
## $ month : int 8 12 8 11 11 2 6 11 5 10 ...
## $ id : int 2642 4627 2568 4192 3977 659 2290 4267 2035 3587 ...
knitr::kable(head(dta))
1 |
120 |
80 |
267 |
1 |
18 |
55 |
25.0 |
8 |
2642 |
1 |
130 |
78 |
192 |
1 |
35 |
53 |
28.4 |
12 |
4627 |
1 |
144 |
90 |
207 |
1 |
109 |
61 |
25.1 |
8 |
2568 |
1 |
92 |
66 |
231 |
1 |
147 |
48 |
26.2 |
11 |
4192 |
1 |
162 |
98 |
271 |
1 |
169 |
39 |
28.4 |
11 |
3977 |
2 |
212 |
118 |
182 |
1 |
199 |
61 |
33.3 |
2 |
659 |
# recode gender variable
dta$sex <- factor(dta$sex,
levels=c(1, 2),
labels=c("M", "F"))
Visualization
# lattice plot
library(lattice)
xyplot(sbp ~ dbp | sex, #依據性別畫舒張壓和收縮壓
data=dta, cex=.5,
type=c("p","g","r"),
xlab="Diastolic pressure (mmHg)", #舒張壓
ylab="Systolic pressure (mmHg)") #收縮壓

Analysis
#Gender difference in regression slopes
#install.packages("lsmeans")
library(lsmeans)
## 載入需要的套件:emmeans
## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.
#
m0 <- lm(sbp ~ dbp, data=dta)
#
m1 <- lm(sbp ~ dbp + sex, data=dta)
#sex:dbp是交互作用
m2 <- lm(sbp ~ dbp + sex + sex:dbp, data=dta)
#
knitr::kable(anova(m0, m1, m2))
4697 |
941778.1 |
NA |
NA |
NA |
NA |
4696 |
927853.3 |
1 |
13924.79 |
71.79987 |
0 |
4695 |
910543.1 |
1 |
17310.17 |
89.25579 |
0 |
# Obtain slopes
m2$coefficients
## (Intercept) dbp sexF dbp:sexF
## 29.8539356 1.2251337 -22.1005883 0.3088544
m.lst <- lstrends(m2, "sex", var="dbp")
knitr::kable(m.lst)
M |
1.225134 |
0.0254189 |
4695 |
1.175301 |
1.274967 |
F |
1.533988 |
0.0205576 |
4695 |
1.493685 |
1.574291 |
# Compare slopes
knitr::kable(pairs(m.lst))
M - F |
-0.3088544 |
0.0326916 |
4695 |
-9.447528 |
0 |
#
plot(rstandard(m2) ~ fitted(m2),
cex=.5,
xlab="Fitted values",
ylab="Standardized residuals")
grid()

The end