Introduction
Data description
3154 healthy young men aged 39-59 from the San Francisco area were assessed for their personality type. All were free from coronary heart disease at the start of the research. Eight and a half years later change in this situation was recorded.
Format
age: Age in yearsheight: Height in inchesweight: Weight in poundssdp: Systolic blood pressure in mm Hgdbp: Diastolic blood pressure in mm Hgchol: Fasting serum cholesterol in mm %behave: Behavior type which is a factor with levels A1 A2 B3 B4cigs: Number of cigarettes smoked per daydibep: Behavior type which is a factor with levels A (Agressive) B (Passive)chd: Coronary heat disease developed which is a factor with levels no yestypechd: type of coronary heart disease which is a factor with levels angina infdeath none silenttimechd: Time of CHD event or end of follow-uparcus: Arcus senilis is a factor with levels absent present
Details
The WCGS began in 1960 with 3,524 male volunteers who were employed by 11 California companies. Subjects were 39 to 59 years old and free of heart disease as determined by electrocardiogram. After the initial screening, the study population dropped to 3,154 and the number of companies to 10 because of various exclusions. The cohort comprised both blue- and white-collar employees. At baseline the following information was collected: socio-demographic including age, education, marital status, income, occupation; physical and physiological including height, weight, blood pressure, electrocardiogram, and corneal arcus; biochemical including cholesterol and lipoprotein fractions; medical and family history and use of medications; behavioral data including Type A interview, smoking, exercise, and alcohol use. Later surveys added data on anthropometry, triglycerides, Jenkins Activity Survey, and caffeine use. Average follow-up continued for 8.5 years with repeat examinations.
Source
Statistics for Epidemiology by N. Jewell (2004)
Goal
To find out the association between diseased variable and each exposure variable that is categorized. And to figure out the causal relationship and make conlusion about the inference. And fit some models to predict the outcome by those variables.
Preprocessing
Import data
# Load library
library(faraway)
data(wcgs)
# Present the structure of data
str(wcgs)Output:
## 'data.frame': 3154 obs. of 13 variables:
## $ age : int 49 42 42 41 59 44 44 40 43 42 ...
## $ height : int 73 70 69 68 70 72 72 71 72 70 ...
## $ weight : int 150 160 160 152 150 204 164 150 190 175 ...
## $ sdp : int 110 154 110 124 144 150 130 138 146 132 ...
## $ dbp : int 76 84 78 78 86 90 84 60 76 90 ...
## $ chol : int 225 177 181 132 255 182 155 140 149 325 ...
## $ behave : Factor w/ 4 levels "A1","A2","B3",..: 2 2 3 4 3 4 4 2 3 2 ...
## $ cigs : int 25 20 0 20 20 0 0 0 25 0 ...
## $ dibep : Factor w/ 2 levels "A","B": 2 2 1 1 1 1 1 2 1 2 ...
## $ chd : Factor w/ 2 levels "no","yes": 1 1 1 1 2 1 1 1 1 1 ...
## $ typechd: Factor w/ 4 levels "angina","infdeath",..: 3 3 3 3 2 3 3 3 3 3 ...
## $ timechd: int 1664 3071 3071 3064 1885 3102 3074 3071 3064 1032 ...
## $ arcus : Factor w/ 2 levels "absent","present": 1 2 1 1 2 1 1 1 1 2 ...# Create a new variable named "status"
status <- rep(1, length(wcgs$chd)) # status = 1 for all
status[which(wcgs$chd == "no")] <- 0 # status = 0 if chd == "no"
# Add status to data
wcgs <- cbind(wcgs, status)Missing value
# Check the number of NA
sum(is.na(wcgs))Output:
## [1] 14As we can see that we have 14 missing value (NA) in our raw data. Furthermore, if we observe those value type of each column, it is easy to find out that we have 2 types of value in total (factor and integer). Hence, in order to get rid of those missing value, we should replace NA of different value type by using different method. Here we use the missForest package to impute missing values. missForest runs a random forests on each variable using the observed part and predicts the na values.
# Load library
library(missForest)
# Impute missing values
wcgs.mis <- missForest(wcgs, maxiter = 5, ntree = 250)
wcgs <- wcgs.mis$ximpOutput:
## missForest iteration 1 in progress...done!
## missForest iteration 2 in progress...done!
## missForest iteration 3 in progress...done!
## missForest iteration 4 in progress...done!# Check the number of NA again
sum(is.na(wcgs))Output:
## [1] 0Split Data
In order to check the performance of our models, we need to build a training data set and a testing data set. Here I decided to split our original data into two parts by half and half. And let one of them to be the train and the other one will be the test.
train <- wcgs[ 1:1576, ] # Set row 1 to 1576 as train
test <- wcgs[1577:3152, ] # Set the rest as test
test.chd <- test[(names(test) %in% c("chd"))]
test <- test[!(names(test) %in% c("typechd", "status", "chd"))]
train <- train[!(names(train) %in% c("typechd", "status"))]Data visualization
Kaplan-Meier Plot
# Load libraries
library(survival)
library(survminer)
library(ggplot2)
# fit a model
fit <- survfit(Surv(timechd, status) ~ typechd, data = wcgs)
# KM curves
ggsurvplot(
fit ,
size = 0.3 ,
pval = TRUE ,
conf.int = TRUE , # Show confidence interval
surv.median.line = "hv" , # Specify median survival
ggtheme = theme_light() , # Change ggplot2 theme
xlab = "Days" ,
legend.title = "Type of Coronary heart disease",
legend.labs = c("Angina","Infdeath", "None", "Silent")
)
# log–rank test
data <- subset(wcgs, typechd != "none")
survdiff(Surv(timechd, status) ~ typechd, data = data)Output:
## Call:
## survdiff(formula = Surv(timechd, status) ~ typechd, data = data)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## typechd=angina 51 51 52.3 0.0304 0.0386
## typechd=infdeath 135 135 137.8 0.0573 0.1245
## typechd=silent 71 71 66.9 0.2478 0.3408
##
## Chisq= 0.3 on 2 degrees of freedom, p= 0.8From the Kaplan-Meier plot above, we can easily see that there are 3 types of Coronary heart disease. The survival probability of all of them are decreasing with time. Angina has relatively lower survival probability at the begining, and has relatively higher survival probability at the end.
We did the log-rank test with the data set which without level "none" in typechd since "none" is not a type of CHD, and our goal is to compare difference in survival between differet types of CHD.
The log-rank test resulted a p-value of p = 0.9, indicatiing that Coronary heart disease types differ insignificantly in survival. This outcome is obvious, since as we can see that although there is some difference between those three types of Coronary Heart Disease in KM curve presented above, in general, those differences are not notable.
Kernel Density Plots
# Load library
library(ggpubr)
# Kernel density plots
chd.age <- wcgs[,c("chd", "age")]
chd.height <- wcgs[,c("chd", "height")]
chd.weight <- wcgs[,c("chd", "weight")]
chd.sdp <- wcgs[,c("chd", "sdp")]
# SalePrice - KitchenQual
age.plt <- ggplot(chd.age, aes(x = age, fill = chd)) +
geom_density(alpha = 0.4) +
labs(title = "CHD distribution by Age")
# SalePrice - LotConfig
height.plt <- ggplot(chd.height, aes(x = height, fill = chd)) +
geom_density(alpha = 0.4) +
labs(title = "CHD distribution by Height")
# SalePrice - SaleType
weight.plt <- ggplot(chd.weight, aes(x = weight, fill = chd)) +
geom_density(alpha = 0.4) +
labs(title = "CHD distribution by Weight")
# SalePrice - SaleType
sdp.plt <- ggplot(chd.sdp, aes(x = sdp, fill = chd)) +
geom_density(alpha = 0.4) +
labs(title = "CHD distribution by sdp")
ggarrange(age.plt, height.plt, weight.plt, sdp.plt, ncol = 2, nrow = 2)
From these kernel density plots presented above we can see that only age has overt effects on CHD. The plot is showing that aged people tend to be more likely to have CHD. And the sdp might have slightly effects on CHD. The height and the weight does not have any association with chd due to the distribution of chd patient similarly through different height/weight.
Associations
Test for association between diseased (Response) variable and each exposure variable that is categorized.
# Behavior
behavior <- table(wcgs$chd, wcgs$behave)
behavior; chisq.test(behavior)Output:
##
## A1 A2 B3 B4
## no 234 1177 1155 331
## yes 30 148 61 18
##
## Pearson's Chi-squared test
##
## data: behavior
## X-squared = 39.916, df = 3, p-value = 1.11e-08# Behavior type
behavior_type <- table(wcgs$chd, wcgs$dibep)
behavior_type; chisq.test(behavior_type)Output:
##
## A B
## no 1486 1411
## yes 79 178
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: behavior_type
## X-squared = 39.08, df = 1, p-value = 4.069e-10# Arcus senilis
arcus_senilis <- table(wcgs$chd, wcgs$arcus)
arcus_senilis; chisq.test(arcus_senilis)Output:
##
## absent present
## no 2058 839
## yes 154 103
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: arcus_senilis
## X-squared = 13.402, df = 1, p-value = 0.0002514From the result of chi-square test above, we can clonclude that all of them have significant relationship with CHD due to the extremely low p-value.
Model Diagnosis
Logistic Regression
# Build logistic regression model
fit.glm <- glm(
chd ~ . - weight - height,
family = binomial(link = "logit") ,
data = train
)
# Check the detail of the model
summary(fit.glm)Output:
##
## Call:
## glm(formula = chd ~ . - weight - height, family = binomial(link = "logit"),
## data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1540 -0.2840 -0.1880 -0.1211 3.1893
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.4808822 1.3633689 -4.754 2.00e-06 ***
## age 0.0593623 0.0199352 2.978 0.0029 **
## sdp 0.0261042 0.0107949 2.418 0.0156 *
## dbp -0.0132608 0.0176858 -0.750 0.4534
## chol 0.0124000 0.0025709 4.823 1.41e-06 ***
## behaveA2 -0.4747303 0.3451883 -1.375 0.1690
## behaveB3 -0.8120439 0.3822421 -2.124 0.0336 *
## behaveB4 -0.5967909 0.4874605 -1.224 0.2208
## cigs 0.0209861 0.0074927 2.801 0.0051 **
## dibepB NA NA NA NA
## timechd -0.0017381 0.0001282 -13.553 < 2e-16 ***
## arcuspresent 0.3571030 0.2389478 1.494 0.1350
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 926.31 on 1575 degrees of freedom
## Residual deviance: 581.74 on 1565 degrees of freedom
## AIC: 603.74
##
## Number of Fisher Scoring iterations: 6From the model summary presented above, it is not hard to see that some variables have high significant level. height, chol, timechd have 3 stars, and Intercept, sdp, cigs have 2 stars.
Since Intercept has 2 stars significant level, which means some interception between variables still valuable, we cannot arbitrarily chose variables with high significant level. Hence, we need the variable selection part as following.
The best decision is to leave the variable selection to computer. Also, we may call it "Subset selection". And it will go both "Forward" and "Backward" directions to select a subset which has the lowest AIC to be the best subset.
# Load package
library(caret)
# Build another model
fit.glm.step <- step(fit.glm, direction = "both")Output:
## Start: AIC=603.74
## chd ~ (age + height + weight + sdp + dbp + chol + behave + cigs +
## dibep + timechd + arcus) - weight - height
9
##
##
## Step: AIC=603.74
## chd ~ age + sdp + dbp + chol + behave + cigs + timechd + arcus
##
## Df Deviance AIC
## - behave 3 586.27 602.27
## - dbp 1 582.31 602.31
## <none> 581.74 603.74
## - arcus 1 583.95 603.95
## - sdp 1 587.53 607.53
## - cigs 1 589.45 609.45
## - age 1 590.68 610.68
## - chol 1 605.32 625.32
## - timechd 1 810.54 830.54
##
## Step: AIC=602.27
## chd ~ age + sdp + dbp + chol + cigs + timechd + arcus
##
## Df Deviance AIC
## - dbp 1 586.67 600.67
## + dibep 1 583.81 601.81
## <none> 586.27 602.27
## - arcus 1 588.84 602.84
## + behave 3 581.74 603.74
## - sdp 1 591.73 605.73
## - cigs 1 595.06 609.06
## - age 1 596.66 610.66
## - chol 1 609.63 623.63
## - timechd 1 823.49 837.49
##
## Step: AIC=600.67
## chd ~ age + sdp + chol + cigs + timechd + arcus
##
## Df Deviance AIC
## + dibep 1 584.29 600.29
## <none> 586.67 600.67
## - arcus 1 589.32 601.32
## + dbp 1 586.27 602.27
## + behave 3 582.31 602.31
## - sdp 1 596.04 608.04
## - cigs 1 596.47 608.47
## - age 1 596.94 608.94
## - chol 1 609.79 621.79
## - timechd 1 823.88 835.88
##
## Step: AIC=600.29
## chd ~ age + sdp + chol + cigs + timechd + arcus + dibep
##
## Df Deviance AIC
## <none> 584.29 600.29
## - dibep 1 586.67 600.67
## - arcus 1 586.87 600.87
## + dbp 1 583.81 601.81
10
## + behave 2 582.31 602.31
## - age 1 593.31 607.31
## - cigs 1 593.36 607.36
## - sdp 1 593.45 607.45
## - chol 1 607.33 621.33
## - timechd 1 812.89 826.89# Visualize model
par(mfrow = c(2,2))
plot(fit.glm.step)
# Make prediction
test.prob <- predict(fit.glm.step, newdata = test, type = "response")
pred.glm <- rep("no", length(test.prob))
pred.glm[test.prob >= 0.5] <- "yes"
# Evaluate the model obtained on the test set
postResample(pred.glm, test.chd$chd)Output:
## Accuracy Kappa
## 0.9238579 0.3247499KNN
This study is a classification problem which reminds me that KNN could be used for this problem.
# Set to a constant result
set.seed(123)
# Set up grid and trControl
grid_knn <- expand.grid(.k = seq(2, 20, by = 1))
tr.control <- trainControl(method = "cv", number = 15)
# Build model
fit.knn <- train(
chd ~ . - weight - height ,
train ,
method = "knn" ,
trControl = tr.control ,
tuneGrid = grid_knn ,
preProcess = c("center", "scale")
)
# Present the model
fit.knnOutput:
## k-Nearest Neighbors
##
## 1576 samples
## 11 predictor
## 2 classes: 'no', 'yes'
##
## Pre-processing: centered (11), scaled (11)
## Resampling: Cross-Validated (15 fold)
## Summary of sample sizes: 1471, 1471, 1471, 1471, 1471, 1471, ...
## Resampling results across tuning parameters:
##
## k Accuracy Kappa
## 2 0.9022761 0.3564266
## 3 0.9232225 0.3955646
## 4 0.9213118 0.3804881
## 5 0.9295837 0.4169642
## 6 0.9245043 0.3505174
## 7 0.9245043 0.3545086
## 8 0.9175202 0.2629168
## 9 0.9187841 0.2634132
## 10 0.9149745 0.2241198
## 11 0.9162444 0.2189774
## 12 0.9156095 0.2173461
## 13 0.9156095 0.2035675
## 14 0.9181491 0.2223683
## 15 0.9168793 0.2001922
## 16 0.9162444 0.1829645
## 17 0.9194190 0.2195936
## 18 0.9213237 0.2395027
## 19 0.9194190 0.2199010
## 20 0.9213237 0.2309750
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 5.# Plot the results
plot(fit.knn)
# Make prediction
pred.knn <- predict(fit.knn, newdata = test)
# Evaluate the model obtained on the test set
postResample(pred.knn, test.chd$chd)Output:
## Accuracy Kappa
## 0.9276650 0.2909698Boosting
# Set to a constant result
set.seed(123)
objControl <- trainControl(
method = 'cv',
number = 20 ,
returnResamp = 'none',
summaryFunction = twoClassSummary,
classProbs = TRUE
)
# Fit a model
fit.gbm <- train(
train[,!(names(train) %in% c("weight", "height", "chd"))] ,
train[,"chd"] ,
method = "gbm" ,
trControl = objControl ,
metric = "ROC" ,
preProc = c("center", "scale")
)# Check detail of model
summary(fit.gbm)
Output:
## var rel.inf
## timechd timechd 64.1175247
## chol chol 13.1085644
## age age 8.0388713
## sdp sdp 4.6793517
## dbp dbp 3.9532382
## cigs cigs 3.9362107
## behave behave 1.2377050
## arcus arcus 0.8109386
## dibep dibep 0.1175955# Make prediction
pred.gbm <- predict(
fit.gbm,
newdata = test,
n.trees = which.min(fit.gbm$cv.error)
)
# Evaluate the model obtained on the test set
postResample(pred.gbm, test.chd$chd)Output:
## Accuracy Kappa
## 0.9244924 0.3637141Prediction
Ensemble
Since this is a classification problem so I decided to use majority voting to ensemble three models.
# The majority vote
prediction <- ifelse(pred.glm == "no" & pred.knn == "no", "no",
ifelse(pred.glm == "no" & pred.gbm == "no", "no",
ifelse(pred.knn == "no" & pred.gbm == "no", "no",
"yes")))
# Evaluate the model obtained on the test set
postResample(prediction, test.chd$chd)Output:
## Accuracy Kappa
## 0.9257614 0.3189232Conclusion
No matter what type of Coronary heart disease you get, the consquences would be similar. Although there are three types of Coronary heart disease presented in the data provided by western collaborative group study, their survival probability almost equal during the period after getting diseased till die. In another word, the difference between different types of Coronary heart disease (Angina, Infdeath, Silent) is not notable.
Furthermore, behaviors, behavior types, arcus senilis are all highly related with Coronary heart disease.
After all, we built a simple model for predicting CHD.