Framingham example

# file location
fL<-"https://math.montana.edu/shancock/data/Framingham.txt"

#讀資料
dta <- read.table(fL, header=T)

# recode gender variable
dta$sex <- factor(dta$sex, 
                  levels=c(1, 2), 
                  labels=c("M", "F")) 

# 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)")

#定義m0模式(regression, y=sbp, x=dbp)
m0 <- lm(sbp ~ dbp, data=dta)

#定義m1模式(regression, y=sbp, x1=dbp, x2=sex)
m1 <- lm(sbp ~ dbp + sex, data=dta)

#定義m2模式(regression, y=sbp, x1=dbp, x2=sex, x3=dbp*sex)
m2 <- lm(sbp ~ dbp + sex + sex:dbp, data=dta)
multi.fit = lm(sbp ~ dbp + sex + sex:dbp, data=dta)
summary(multi.fit)

## 
## Call:
## lm(formula = sbp ~ dbp + sex + sex:dbp, data = dta)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.540  -8.990  -1.929   7.135  93.644 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  29.85394    2.14394  13.925  < 2e-16 ***
## dbp           1.22513    0.02542  48.198  < 2e-16 ***
## sexF        -22.10059    2.73844  -8.070 8.82e-16 ***
## dbp:sexF      0.30885    0.03269   9.448  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.93 on 4695 degrees of freedom
## Multiple R-squared:  0.6272, Adjusted R-squared:  0.627 
## F-statistic:  2633 on 3 and 4695 DF,  p-value: < 2.2e-16

##Findings revealed a significant 2-way interaction existed b/w sex and dbp in terms of sbp, which implied that people of different sex (M=1, F=2) demonstrated different relationship patterns b/w dbp and sbp. 

#Compute overall F-test for m0,m1,m2 models
anova(m0, m1, m2)

## Analysis of Variance Table
## 
## Model 1: sbp ~ dbp
## Model 2: sbp ~ dbp + sex
## Model 3: sbp ~ dbp + sex + sex:dbp
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1   4697 941778                                  
## 2   4696 927853  1     13925 71.800 < 2.2e-16 ***
## 3   4695 910543  1     17310 89.256 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#residual plot
plot(rstandard(m2) ~ fitted(m2), 
     cex=.5, 
     xlab="Fitted values", 
     ylab="Standardized residuals")
grid()

# The end