1 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?

2 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))
sex sbp dbp scl chdfate followup age bmi month id
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")) 

3 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)")  #收縮壓

4 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))
Res.Df RSS Df Sum of Sq F Pr(>F)
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)
sex dbp.trend SE df lower.CL upper.CL
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))
contrast estimate SE df t.ratio p.value
M - F -0.3088544 0.0326916 4695 -9.447528 0
#
plot(rstandard(m2) ~ fitted(m2), 
     cex=.5, 
     xlab="Fitted values", 
     ylab="Standardized residuals")
grid()

5 The end