Statistical Models and Methods (MATH2038 UNNC) Spring Coursework
Zhongyu FAN - 20412336
2024-05-06
1 Part(a)
1.1 Introduction
The purpose of this report is to investigate associations between the birthweight of children and various potential explanatory variables and to establish a predictive model for the birthweight of children. One main purpose of this report is to identify variables related to birth weight, which can then be further investigated by medical professionals to understand any possible causal relationships.
In our training dataset, we have 327 records containing information on four numerical variables (Age of mother - age, Gestation period - gest, Pre-pregnancy weight of mother - weight, Rate of growth of child in the first trimester - rate), two categorical variables (Sex of child, with levels ‘Female’ and ‘Male’ - sex, Whether the mother smoked during pregnancy, with levels ‘No’, ‘Light’, and ‘Heavy’ - smokes), and the response variable (birth weight of children - bwt).
Here we can see the structure of the data as follows.
## 'data.frame': 327 obs. of 7 variables:
## $ age : int 25 34 29 23 26 30 27 26 40 19 ...
## $ gest : int 294 252 280 280 266 245 273 273 287 287 ...
## $ sex : Factor w/ 2 levels "Female","Male": 1 1 2 1 1 2 1 1 2 2 ...
## $ smokes: Factor w/ 3 levels "Heavy","Light",..: 3 3 1 2 2 1 3 3 3 3 ...
## $ weight: num 65.2 58.6 57.3 65.9 59.1 49.1 69.5 75 82 58.3 ...
## $ rate : num -0.01618 0.00678 0.03293 -0.02179 0.06346 ...
## $ bwt : int 3870 2440 3400 2820 2940 2400 3080 3600 3600 3930 ...1.2 Exploratory Data Analysis
1.2.1 Overview
Before constructing the model, we performed Exploratory Data Analysis (EDA) on the data to observe the relationships and potential focus points between the data. During this stage, the potential relationships between individual independent variables and the response variable were closely examined.
Prior to formal commencement, here are the descriptive statistics of each variable.
## age gest sex smokes weight
## Min. :15.00 Min. :224 Female:153 Heavy: 13 Min. : 40.90
## 1st Qu.:24.00 1st Qu.:273 Male :174 Light: 14 1st Qu.: 52.30
## Median :27.00 Median :280 No :300 Median : 56.80
## Mean :27.68 Mean :276 Mean : 58.29
## 3rd Qu.:31.00 3rd Qu.:280 3rd Qu.: 63.00
## Max. :42.00 Max. :301 Max. :100.00
## rate bwt
## Min. :-0.040020 Min. :1640
## 1st Qu.: 0.004185 1st Qu.:3170
## Median : 0.026110 Median :3440
## Mean : 0.028550 Mean :3461
## 3rd Qu.: 0.048285 3rd Qu.:3760
## Max. : 0.124320 Max. :5680The normality of the response variable (Bwt) was also checked using Quantile-Quantile Plot here.
Figure 1
Based on the figure, it can be seen that the Bwt largely
conforms to a normal distribution.
Variables are categorized here as numerical or categorical for the purpose of discussion.
1.2.2 numerical variables
Figure 2
From the Lower Triangle of Correlation Matrix, we can see that
bwt has a strong positive correlation with gest and weight, while the
relationships between other variables are not significant.
Due to the substantial impact of multicollinearity on multiple regression analysis, here we also utilize the Variance Inflation Factor (VIF) for detection.
## age gest weight rate
## 1.073711 1.037887 1.062426 1.027513
The results show no significant collinearity issues.
Figure 3
Further investigation was conducted using a scatterplot matrix
to study the potential relationships between the independent variables
and the response variable. From the graphs, it can be observed that
there is a strong positive correlation between Bwt and gest as well as
weight, a weak positive correlation between Bwt and age, rate, which
requires further processing and exploration. This finding is consistent
with the conclusion obtained from the Lower Triangle of Correlation
Matrix.
1.2.3 Category variables
We use box plots and pirate plots to display the relationship between categorical variables and response variables.
Figure 4
For the categorical variable ‘sex’, it is observed that the
difference in means between the male group and the female group is not
significant, with males being slightly higher overall than females.
There are some outliers in both male and female groups, removing these
outliers before building the model may improve the predictive
performance of the model.
Figure 5
For the categorical variable “smokes”, the mean differences
among the three groups are also not significant. The differences in
sample sizes between the “Heavy” and “Light” groups compared to the “No”
group contribute to the varying shapes of the histograms. There are some
outliers in both the “Light” and “No” groups, and considering removing
these outliers before building the model may improve the predictive
performance of the model.
To further validate, we conducted an analysis using ANOVA. Prior to this, the normality of the response variable Bwt had been confirmed. Here, the homoscedasticity of variance among data groups was examined first, followed by the ANOVA analysis.
## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 586635 586635 2.297 0.131
## Residuals 325 83010399 255417## Df Sum Sq Mean Sq F value Pr(>F)
## smokes 2 1025781 512891 2.013 0.135
## Residuals 324 82571253 254850
The results of the two homogeneity of variance tests show that
the p-values are both greater than 0.05. Therefore, there is not enough
evidence to reject the null hypothesis, indicating that there is no
significant difference in variance between these data groups
statistically.
Moreover, ANOVA analysis further supports our conclusion: the p-values are 0.131 and 0.135 respectively, both greater than 0.05. This also indicates that we do not have enough evidence to reject the null hypothesis, namely,
- Gender does not have a significant impact on birth weight.
- Smoking status also does not have a significant impact on birth weight.
1.3 Model Building and Diagnostic
We only used the training set to build the model, but tested several potential models using the test set to find the best-performing model.
1.3.1 Data Cleaning
Before the formal model building, we processed and cleaned the data. Since the response variable bwt in the training set roughly follows a normal distribution, we calculated Z-scores to detect outliers and removed them; Additionally, we also eliminated outliers in the two categorical variables mentioned in 1.2.2. The cleaned dataset Births_cleaned3 was obtained after these steps.
1.3.2 Models Building
From EDA, we learned that gest and weight have a significant correlation with bwt. Therefore, we start from here to build a linear regression model lm_reduced, which serves as our candidate model 1.
##
## Call:
## lm(formula = bwt ~ gest + weight, data = Births_cleaned3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -942.57 -259.93 -1.72 222.12 1200.50
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1726.138 592.031 -2.916 0.00381 **
## gest 17.213 2.139 8.048 1.81e-14 ***
## weight 7.348 2.219 3.311 0.00104 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 373.6 on 310 degrees of freedom
## Multiple R-squared: 0.2101, Adjusted R-squared: 0.205
## F-statistic: 41.23 on 2 and 310 DF, p-value: < 2.2e-16
We use the Stepwise method for variable selection based on the
potential Model 1, and the model selected by Stepwise serves as our
potential Model 2.
## Start: AIC=3710.98
## bwt ~ gest + weight
##
## Df Sum of Sq RSS AIC
## + rate 1 830413 42448993 3706.9
## + age 1 492316 42787090 3709.4
## <none> 43279406 3711.0
## + sex 1 168840 43110565 3711.8
## + smokes 2 261354 43018052 3713.1
## - weight 1 1530415 44809821 3719.9
## - gest 1 9042569 52321974 3768.4
##
## Step: AIC=3706.91
## bwt ~ gest + weight + rate
##
## Df Sum of Sq RSS AIC
## + age 1 320671 42128321 3706.5
## <none> 42448993 3706.9
## + sex 1 167634 42281359 3707.7
## + smokes 2 316816 42132177 3708.6
## - rate 1 830413 43279406 3711.0
## - weight 1 1405744 43854737 3715.1
## - gest 1 9250715 51699708 3766.6
##
## Step: AIC=3706.54
## bwt ~ gest + weight + rate + age
##
## Df Sum of Sq RSS AIC
## <none> 42128321 3706.5
## - age 1 320671 42448993 3706.9
## + sex 1 155092 41973229 3707.4
## + smokes 2 333998 41794324 3708.0
## - rate 1 658768 42787090 3709.4
## - weight 1 1107024 43235345 3712.7
## - gest 1 9569507 51697829 3768.6##
## Call:
## lm(formula = bwt ~ gest + weight + rate + age, data = Births_cleaned3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1005.7 -252.2 -8.2 231.0 1183.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2095.077 615.143 -3.406 0.000747 ***
## gest 17.982 2.150 8.364 2.12e-15 ***
## weight 6.382 2.243 2.845 0.004740 **
## rate 1460.948 665.702 2.195 0.028939 *
## age 6.161 4.024 1.531 0.126758
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 369.8 on 308 degrees of freedom
## Multiple R-squared: 0.2311, Adjusted R-squared: 0.2211
## F-statistic: 23.15 on 4 and 308 DF, p-value: < 2.2e-16
In addition, due to the small training data set and our desire
for higher prediction stability, one possible technique is to use ridge
regression for modeling. We utilized two ridge regression models. In the
first ridge regression model, the variables were selected by Stepwise to
form variable combinations, and this model represents our potential
Model 3.
## [1] "cv_ridge1$lambda.min = 37.5904970855684"
Figure 6
The second type of ridge regression model, which is our latent
model 4, uses all variables as variable combinations.
## [1] "cv_ridge2$lambda.min = 17.8585315552614"
Figure 7
1.3.3 Model Selection
We trained our models using a training set and tested them using a test set, obtaining the Mean Squared Error (MSE), \(R^2\) (coefficient of determination), and Mean Absolute Error (MAE) for each model.
Below are the MSE, \(R^2\), and MAE for these four models.
## RMSE Rsquared MAE
## 388.8764106 0.1645811 316.5204779## RMSE Rsquared MAE
## 383.2145472 0.1852196 303.7155070## RMSE Rsquared MAE
## 381.7145796 0.1852965 303.2859112## RMSE Rsquared MAE
## 376.8733851 0.1962969 296.9827591Based on the results, we chose potential Model 4 as the final best model. This is because its RMSE and MAE are the smallest, and its \(R^2\) is the highest among the four models, indicating that this model has the smallest absolute difference between predicted values and actual values, as well as the strongest ability to explain data variability.
1.3.4 Model Diagnostic
Now we are conducting model diagnosis and inspection for the best model.
Firstly, residual plots are generated, which helps to visually inspect the distribution of residuals. Ideally, residuals should be randomly distributed around 0 without any apparent pattern.
Figure 8
Figure 8 supported our hypothesis that the residuals are
randomly distributed without any apparent pattern.
We also plotted a Q-Q plot to check the normality of the residuals.
Figure 9
From Figure 9, it can be observed that the residual points
roughly follow a straight line, indicating that the residuals are close
to a normal distribution.
We also conducted additional checks for heteroscedasticity.
Figure 10
Heteroskedasticity checks are also passed.
1.4 Summary of Findings and Conclusions
In our report, we studied the impact of factors such as weight and smoking habits on birthweight of children and established a predictive model. By analyzing multiple variables, we found that maternal weight and gestational age significantly influence birthweight of children. Additionally, we cannot ignore the impact of maternal age, infant growth rate three months after birth, infant gender, and maternal smoking status on birthweight of children.
These findings underscore the importance of providing targeted health guidance during pregnancy, especially regarding weight management recommendations; at the same time, the study reminds us to pay attention to the factor of gestational period. Based on these findings, we recommend medical professionals to pay more attention to these factors, which may be related to infant health.
Future research should further explore possible causal relationships and interactions among these factors, and evaluate the effects of certain interventions.
2.Part(b)
2.1 Model Prediction and Evaluation Using Test Data
We used the best model selected in section 1.3.3 for prediction and output its Mean Squared Error (MSE), \(R^2\) (coefficient of determination), and Mean Absolute Error (MAE).
## RMSE Rsquared MAE
## 376.8733851 0.1962969 296.9827591In addition, scatter plot and bar plot of predicted values against actual values are drawn to visually demonstrate the predictive accuracy of the model.
Figure 11
Figure 12
2.2 Comparison with Full Model
Compared with the Full Model, the best model has a smaller Mean Squared Error (MSE), Mean Absolute Error (MAE), and a larger \(R^2\) (coefficient of determination), which indicates that the best model performs better in prediction than the full model.
## RMSE Rsquared MAE
## 376.8733851 0.1962969 296.9827591## RMSE Rsquared MAE
## 397.5261441 0.2020091 317.8163177
Figure 13
A scatter plot comparing the predicted values and actual values
of the two models is plotted for visual illustration.
Appendix
1.1
## Input Data
Births <- read.table("BirthTrain.txt", header = TRUE)
Births$sex <- as.factor(Births$sex)
Births$smokes <- as.factor(Births$smokes)
str(Births)
Births[1:10,]1.2
1.2.1
## Summary statistics for numerical variables
summary(Births)qqnorm(Births$bwt, main="Q-Q Plot of Response Variable")
qqline(Births$bwt, col="red")1.2.2
## Heatmap Matrix
library(ggplot2)
library(reshape2)
num_vars <- data.frame(
age = Births$age,
gest = Births$gest,
weight = Births$weight,
rate = Births$rate,
bwt = Births$bwt
)
cor_matrix <- cor(num_vars)
cor_matrix[upper.tri(cor_matrix)] <- NA
cor_melted <- melt(cor_matrix, na.rm = TRUE)
ggplot(cor_melted, aes(Var2, Var1, fill = value)) +
geom_tile(color = "white") +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1, 1), space = "Lab",
name = "Pearson\nCorrelation") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
axis.text.y = element_text(size = 12)) +
labs(x = "", y = "", title = "Lower Triangle of Correlation Matrix")## Check Multicollinearity
library(car)
model <- lm(bwt ~ age + gest + weight + rate, data = Births)
vif(model)## Scatterplot Matrix
library(ggplot2)
library(GGally)
ggpairs(Births,
columns = c("bwt", "age", "gest", "weight", "rate"),
title = "Enhanced Pairs Plot with GGally",
upper = list(continuous = wrap("smooth", method = "lm", colour = "lightblue", size = 1, se = TRUE)),
lower = list(continuous = wrap("cor", size = 3)),
diag = list(continuous = "barDiag") ) 1.2.3
###Box plot and pirate plot - sex
library(ggplot2)
library(gridExtra)
library(yarrr)
par(mfrow = c(1, 2))
plot1 <- boxplot(bwt ~ sex, data = Births, main = "Birthweight by Sex")
plot2 <- pirateplot(formula = bwt ~ sex, data = Births,
inf.method = "ci", inf.disp = "line",
main = "Birthweight by Sex",
point.o = 0.4)
par(mfrow = c(1, 1))###Box plot and pirate plot - smokes
library(ggplot2)
library(gridExtra)
library(yarrr)
par(mfrow = c(1, 2))
boxplot(bwt ~ smokes, data = Births, main = "Birthweight by Smoking Status")
plot2 <- pirateplot(formula = bwt ~ smokes, data = Births,
inf.method = "ci", inf.disp = "line",
main = "Birthweight by Smoking Status",
point.o = 0.4)
par(mfrow = c(1, 1))library(car)
# Levene's Test
car::leveneTest(bwt ~ sex, data=Births)
car::leveneTest(bwt ~ smokes, data=Births)
# ANOVA
anova_model1 <- aov(bwt ~ sex, data=Births)
summary(anova_model1)
anova_model2 <- aov(bwt ~ smokes, data=Births)
summary(anova_model2)1.3
1.3.1
##Data Cleaning
Births <- read.table("BirthTrain.txt", header = TRUE)
Births$sex <- as.factor(Births$sex)
Births$smokes <- as.factor(Births$smokes)
z_scores <- scale(Births$bwt)
outliers <- which(abs(z_scores) > 3)
Births[outliers, ]
Births_cleaned1 <- Births[-outliers, ]
library(dplyr)
library(ggplot2)
# SEX
outlier_bounds <- Births_cleaned1 %>%
group_by(sex) %>%
summarise(
Q1 = quantile(bwt, 0.25),
Q3 = quantile(bwt, 0.75),
IQR = IQR(bwt),
Lower_Bound = Q1 - 1.5 * IQR,
Upper_Bound = Q3 + 1.5 * IQR
)
Births_cleaned2 <- Births_cleaned1 %>%
left_join(outlier_bounds, by = "sex") %>%
filter(bwt >= Lower_Bound & bwt <= Upper_Bound)
# Smokes
outlier_bounds1 <- Births_cleaned2 %>%
group_by(smokes) %>%
summarise(
Q1 = quantile(bwt, 0.25),
Q3 = quantile(bwt, 0.75),
IQR = IQR(bwt),
Lower_Bound1 = Q1 - 1.5 * IQR,
Upper_Bound1 = Q3 + 1.5 * IQR
)
Births_cleaned3 <- Births_cleaned2 %>%
left_join(outlier_bounds1, by = "smokes") %>%
filter(bwt >= Lower_Bound1 & bwt <= Upper_Bound1) 1.3.2
## EDA MODEL
library(stats)
lm_reduced <- lm(bwt ~ gest + weight, data = Births_cleaned3)
summary(lm_reduced)## Step-wise MODEL
step_model <- step(lm_reduced, scope = bwt ~ gest + weight + age + rate + sex + smokes)
summary(step_model)## Ridge Model - Based on Step-wise Model
library(glmnet)
x_matrix1 <- model.matrix(bwt ~ gest + weight + rate + age, data = Births_cleaned3)[,-1]
y_vector1 <- Births_cleaned3$bwt
ridge_model1 <- glmnet(x_matrix1, y_vector1, alpha = 0)
cv_ridge1 <- cv.glmnet(x_matrix1, y_vector1, alpha = 0)
print(paste("cv_ridge1$lambda.min =", cv_ridge1$lambda.min))
best_ridge_model1 <- glmnet(x_matrix1, y_vector1, alpha = 0, lambda = cv_ridge1$lambda.min)
plot(cv_ridge1)
abline(v = log(cv_ridge1$lambda.min), col = "red")## Ridge Model - Based on FUll Model
x_matrix2 <- model.matrix(bwt ~ gest + weight + age + rate + sex + smokes, data = Births_cleaned3)[,-1]
y_vector2 <- Births_cleaned3$bwt
ridge_model2 <- glmnet(x_matrix2, y_vector2, alpha = 0)
cv_ridge2 <- cv.glmnet(x_matrix2, y_vector2, alpha = 0)
print(cv_ridge2$lambda.min)
best_ridge_model2 <- glmnet(x_matrix2, y_vector2, alpha = 0, lambda = cv_ridge2$lambda.min)
plot(cv_ridge2)
abline(v = log(cv_ridge2$lambda.min), col = "lightgreen")1.3.3
library(caret)
# Model Selection
BirthsTest <- read.table("BirthTest.txt", header = TRUE)
BirthsTest$sex <- as.factor(BirthsTest$sex)
BirthsTest$smokes <- as.factor(BirthsTest$smokes)
predictions1 <- predict(lm_reduced, newdata = BirthsTest)
predictions2 <- predict(step_model, newdata = BirthsTest)
x_test1 <- model.matrix(bwt ~ gest + weight + rate + age, data = BirthsTest)[,-1]
predictions3 <- predict(best_ridge_model1, newx = x_test1, s = cv_ridge1$lambda.min)
x_test2 <- model.matrix(bwt ~ gest + weight + age + rate + sex + smokes, data = BirthsTest)[,-1]
predictions4 <- predict(best_ridge_model2, newx = x_test2, s = cv_ridge2$lambda.min)
#Result
postResample(predictions1, BirthsTest$bwt)
postResample(predictions2, BirthsTest$bwt)
postResample(predictions3, BirthsTest$bwt)
postResample(predictions4, BirthsTest$bwt)1.3.4
# Calculate the residuals using the predicted values from the ridge regression model.
predictions_ridge <- predict(best_ridge_model2, newx = x_matrix2, s = cv_ridge2$lambda.min)
residuals_ridge <- y_vector2 - predictions_ridge# Residual plot
plot(residuals_ridge, main = "Residuals of Ridge Regression", ylab = "Residuals", xlab = "Index")
abline(h = 0, col = "red")# Q-Q diagram to test normality
qqnorm(residuals_ridge)
qqline(residuals_ridge, col = "red")# Heteroscedasticity
plot(predictions_ridge, residuals_ridge, main = "Residuals vs Fitted", xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, col = "red")2.1
postResample(predictions4, BirthsTest$bwt)plot(BirthsTest$bwt, predictions4 , main = "Predicted vs. Actual Birth Weights",
xlab = "Actual Birth Weight", ylab = "Predicted Birth Weight")
abline(0, 1, col = "red")
par(mfrow = c(1, 1))
barplot(rbind(predictions4, BirthsTest$bwt), beside = TRUE, col = c("blue", "gray"),
main = "Comparison of Predicted and Actual Birth Weights")
legend(x = 0.85, y = 1.15, c("Predicted", "Actual"), fill = c("blue", "gray"), xpd = TRUE, horiz = TRUE)2.2
full_model <- lm(bwt ~ gest + weight + age + rate + sex + smokes, data = Births)
predictions5 <- predict(full_model, newdata = BirthsTest)
postResample(predictions4, BirthsTest$bwt)
postResample(predictions5, BirthsTest$bwt)plot(BirthsTest$bwt, predictions4,
main = "Comparison of Predictions vs. Actual Birth Weights",
xlab = "Actual Birth Weights", ylab = "Predicted Birth Weights",
pch = 19, col = "blue", xlim = range(BirthsTest$bwt), ylim = range(c(predictions4, predictions5)))
points(BirthsTest$bwt, predictions5, pch = 17, col = "red")
legend("topright", legend = c("Best Model", "Full Model"),
col = c("blue", "red"), pch = c(19, 17), cex = 0.8)
abline(0, 1, col = "black")