<jupyter_start><jupyter_text>STAT 628 Module 2 BodyFat Group 3 Introduction: Our motivation is to obtain a precise, concise and reasonable regression model to predict body fat percentage of a male based on his body informations. Dataset:Our dataset ```BodyFat``` is a real data set of 252 men with measurements of their percentage of body fat and various body circumference measurements. Background: We assume that percentage of body fat for an individual can be estimated once body density has been determined:$$BodyFat(\%) = \frac{495}{BodyDensity(gm/cm^3)} - 450 \tag{1}$$And the BMI can be calculated by:$$BMI = \frac{703\times Weight(lbs)}{Height(inches)} \tag{2}$$ 1. Exploratory data analysis<jupyter_code># install / import packages
usePackage = function(p) {
    if (!is.element(p,installed.packages()[,1]))
        install.packages(p, dep=TRUE)
    require(p, character.only=TRUE)
}

usePackage("ggplot2")
usePackage("Matrix")
usePackage("foreach")
usePackage("glmnet")
usePackage("MASS")
# Load the data
data = read.csv('BodyFat.csv')[,-1]
head(data)<jupyter_output><empty_output><jupyter_text>We can first intuitively analyze the data by looking at its summary tables and histograms. There are some obvious "abnormal" obsevations. Some results are as following:<jupyter_code># summary and hist
summary(data[, c(1, 4, 5, 6)]) 
options(repr.plot.width=6, repr.plot.height=3)
par(mfrow = c(1,2))
hist(data$BODYFAT, breaks = 30, main = "Hist of BODYFAT", xlab = "BODYFAT")
hist(data$WEIGHT, breaks = 30, main = "Hist of BODYFAT", xlab = "WEIGHT")
par(mfrow = c(1,1))<jupyter_output><empty_output><jupyter_text>From the results above, we can find some "abnormal" obsevations: For example, the minimal ```BODYFAT``` is 0 %, which is unlikely to be true; besides, the maximal ```BODYFAT``` is 36 %, which is too large for a living person. Similarly, the minimal ```HEIGHT```, maximal ```WEIGHT``` and maxinal ```ADIPOSITY``` are all unlikely to be true for a nomal man. We can calculate the predicted values of these "abnormal" data using some formula or just delete these observations. Anyway, data cleaning is needed for this dataset.   2. Data cleaningTo find out outliers, we can first look at the **boxplot:**<jupyter_code>col_name = colnames(data)
n_row_0 = nrow(data)
n_col_0 = ncol(data)
# reshape the data into one-dim
data_all_0 = as.matrix(data)
dim(data_all_0) = c(n_row_0*n_col_0,1)
facs_0 = gl(n_col_0,n_row_0,label=col_name)
data_fac_0 = data.frame(data_all_0,facs_0)
# boxplot
box_0 = ggplot(data_fac_0,aes(x=factor(data_fac_0$facs_0),y=data_fac_0$data_all_0,fill = factor(data_fac_0$facs_0)))+
  geom_boxplot(outlier.color = "red",outlier.shape=21,outlier.size = 1.5)+
  labs(title = "boxplot of original data",xlab="factor",ylab="data value")+
  guides(fill=FALSE)+theme(text = element_text(size=10),axis.text.x = element_text(angle=90, hjust=1))
options(repr.plot.width=3, repr.plot.height=3)
box_0
# ggsave("boxplot_original.png",box_0)<jupyter_output><empty_output><jupyter_text>Calculate the rages of ```WEIGHT```,```HEIGHT```,```DENSITY```,```BODYFAT``` between there 1% quantiles and 99% quantiles respectively, any value out of such ranges can be considered as "abnormal" observations:<jupyter_code>weight_q = c(quantile(data$WEIGHT,0.01),quantile(data$WEIGHT,0.99))
#(125.505,245.720)
height_q = c(quantile(data$HEIGHT,0.01),quantile(data$HEIGHT,0.99))
#(64.3825 ,76.0000 )
density_q = c(quantile(data$DENSITY,0.01),quantile(data$DENSITY,0.99))
#(1.016040 ,1.095393)
bodyfat_q = c(quantile(data$BODYFAT,0.01),quantile(data$BODYFAT,0.99))
#(4.355 ,35.582)
print("Quantile ranges:")
quantiles = rbind(weight_q, height_q, density_q, bodyfat_q)
quantiles<jupyter_output>[1] "Quantile ranges:"<jupyter_text>We can then compare the body fat percentages calculated from formula (1) (```CalBodyFat```) and given by the dataset (```RealBodyFat```) of each observation to see if there is any strong difference. Similarly, we can compare the BMI values calculated from formula (2) (```CalBMI```) and given by the dataset (```RealBMI```) to find some outliers:<jupyter_code>#Bodyfat & BMI from calculation
bodyfat_siri = (495/(data$DENSITY))-450
BMI_cal = (703*data$WEIGHT)/((data$HEIGHT)^2)
# compare the calculated results and the given values
compare_bodyfat = data.frame(RealBodyfat=data$BODYFAT,CalBodyfat=bodyfat_siri)
compare_BMI = data.frame(RealBMI=data$ADIPOSITY,CalBMI=BMI_cal)<jupyter_output><empty_output><jupyter_text>2.1 For BMI:<jupyter_code># plot the results
bmi_plot_1<- ggplot(compare_BMI,aes(x=RealBMI,y=CalBMI))+geom_point()+
  geom_abline(intercept = 0,slope = 1)+
  annotate("text",x=29.9,y=150.6,label="42",vjust=-1)+labs(title = "")
options(repr.plot.width=6, repr.plot.height=3)
bmi_plot_1
# ggsave("bmi_compare_1.png",bmi_plot_1)<jupyter_output><empty_output><jupyter_text>The 42nd observation seems to be an outlier. Check his ```WEIGHT``` and ```HEIGHT``` to find out where is the problem:<jupyter_code>paste("WEIGHT of 42nd observation: ", data[42,]$WEIGHT)
paste("Quantiles of WEIGHT: ")
weight_q
paste("HEIGHT of 42nd observation: ", data[42,]$HEIGHT)
paste("Quantiles of HEIGHT: ")
height_q<jupyter_output><empty_output><jupyter_text>From the results above, we can see that the *42nd* observation’s height (**29.5** inches) is out of quantile range, so we use the formula (2) to calculated the height based on his weight and BMI, it turns out that his height should be **69.43** inches:<jupyter_code># replace the height of the 42nd obs
data_new = data
data_new$HEIGHT[42] = sqrt(703*(data$WEIGHT[42])/data$ADIPOSITY[42])
#after changing the height of point 42
compare_BMI$CalBMI[42]<- (703*data_new$WEIGHT[42])/((data_new$HEIGHT[42])^2)
bmi_plot_2<- ggplot(compare_BMI,aes(x=RealBMI,y=CalBMI))+geom_point()+
  geom_abline(intercept = 0,slope = 1)+
  annotate("text",x=23.4,y=25,label="163",vjust=-1)+
  annotate("text",x=24.5,y=18,label="221",vjust=-1)
options(repr.plot.width=6, repr.plot.height=3)
bmi_plot_2
# ggsave("bmi_compare_2.png",bmi_plot_2)<jupyter_output><empty_output><jupyter_text>After changing the height of the *42nd* observation, we can find that the 163rd and 221st observations still don't fit the line perfectly, but after similar analysis, we find out that their heights and weights is within the quantile range above, so we won’t change them.  2.2 For BodyFat:<jupyter_code>#compare BodyFat, plot the results
bodyfat_plot_1<- ggplot(compare_bodyfat,aes(x=RealBodyfat,y=CalBodyfat))+geom_point()+
  geom_abline(intercept = 0,slope = 1)+
  annotate("text",x=6.4,y=9,label="48",vjust=-1)+
  annotate("text",x=17.3,y=-4,label="96",vjust=-1)+
  annotate("text",x=0,y=-3.7,label="182",vjust=-1)
options(repr.plot.width=6, repr.plot.height=3)
bodyfat_plot_1
# ggsave("bodyfat_compare.png",bodyfat_plot_1)
# the data of abnormal observations
compare_bodyfat[c(182, 48, 96),]
myoutput = rbind(density = data$DENSITY[c(48,96)], bodyfat = data$BODYFAT[c(48,96)])
colnames(myoutput) = c("48th", "96th")
myoutput<jupyter_output><empty_output><jupyter_text>From the graph above, we can see that the *182nd* observation have a **ZERO** calculated body fat percentage and a **NEGATIVE** given body fat percentage, which is unlikely to be true for a normal male, so, we will just remove this observation. The *48th* observation has normal body density and body fat percentage, so we will keep it in spite of its unfitness of the line. However, the *96th* obsrvation has a density far away from the density quantile range, so we should modify the body density of this observation or just drop this observation. Since people can not easily obtain body density in real life, we will just remove this observation (We also drop the column ```DENSITY``` since such data is not feasible in real world application):<jupyter_code># we remove point 182 and 96; drop column DENSITY
data_new = data_new[-c(182,96),-2]<jupyter_output><empty_output><jupyter_text>2.3 Other outliersFinally, we use the 0.01 and 0.99 quantiles to replace the values out of 0.01 ~ 0.99 ranges and save the cleaned data into a ```.csv``` file:<jupyter_code>#deal with other outlier
del_outliers = function(data,lower_b,upper_b){
  
  col_n = ncol(data)
  for(i in 1:col_n){
    temp = data[,i]
    q_lower = quantile(temp,lower_b)
    q_upper = quantile(temp,upper_b)
    data[temp < q_lower,i] = q_lower
    data[temp > q_upper,i] = q_upper
  }
  return(data)
}

data_new = del_outliers(data_new,0.01,0.99)
# write.csv(data_new,"cleaned.csv")<jupyter_output><empty_output><jupyter_text>The **boxplot** of the cleaned data:<jupyter_code># boxplot of the cleaned data
data_all<- as.matrix(data_new)
dim(data_all)<- c(nrow(data_new)*ncol(data_new),1)
col_name_1<- colnames(data_new)
facs<- gl(ncol(data_new),nrow(data_new),label=col_name_1)
data_fac<- data.frame(data_all,facs)

box = ggplot(data_fac,aes(x=factor(data_fac$facs),y=data_fac$data_all,fill = factor(data_fac$facs)))+
  geom_boxplot(outlier.color = "red",outlier.shape=21,outlier.size = 1.5)+
  labs(title = "boxplot of adjusted data",xlab="factor",ylab="data value")+guides(fill=FALSE)+
  theme(text = element_text(size=10),axis.text.x = element_text(angle=90, hjust=1))
options(repr.plot.width=3, repr.plot.height=3)
box
# ggsave("boxplot_adjusted.png",box)<jupyter_output><empty_output><jupyter_text>Obviously, there are less outliers after the data cleaning. Next, we will perform feature selection. 3. Feature Selection 3.1 Cross ValidationWe use a **10-fold cross validation** to compute the average MSE (mean-square error) of each model:$$MSE = E((Y - \hat{Y})^2) \tag{3}$$We shuffle the data before every CV and repeated 1000 times.<jupyter_code>#this is a function to calculate the cross validation
cv = function(data, model){
  num = dim(data)[1]
  var = dim(data)[2]
  l = round(n/10) # length of each fold
  MSE = 1:10
  for(i in 1:10){
    test = (l*(i-1)+1):(l*i)
    if(i == 10){
      test = (l*9+1):num
    }
    model = lm(model, data = data[-test,])
    test_X = data[test, 2:var]
    test_Y = data[test, 1]
    predict_Y = predict(model, test_X)
    MSE[i] = sum((test_Y - predict_Y)^2) / length(test_Y)
  }
  return(mean(MSE))
}

# repeat 10-fold CV 1000 times and then return the mean of MSE
FindCV = function(data,mymodel){
  n=nrow(data)
  cv_mean = 0
  for(i in 1:1000){
    shuffle = sample(1:n, n)
    model_cv = cv(data[shuffle,], mymodel)
    cv_mean = cv_mean + model_cv
  }
  cv_mean = cv_mean/1000
  return(cv_mean)
}<jupyter_output><empty_output><jupyter_text>3.2 Stepwise model selectionPerform stepwise selections: We choose **AIC** and **BIC** methods:Let k denote the number of parameters, let RSS denote the residual sum of squares, let n denote the number of observations, we have:  ** AIC (Akaike Information Criterion): **$$AIC = 2k + n\ln{\frac{RSS}{n}} \tag{4}$$Choose the model with minimal AIC.  ** BIC (Bayesian information criterion): **$$BIC = k\ln{n} + n\ln{\frac{RSS}{n}} \tag{5}$$Choose the model with minimal BIC.<jupyter_code>#aic selection
model_full=lm(BODYFAT~.,data=data_new)
model_null=lm(BODYFAT~1, data=data_new)
aic_selection = step(model_null, scope=list(lower=model_null, upper=model_full), 
                     direction="both", trace=0, k = 2)$call$formula
print(aic_selection)
#bic selection
n = dim(data_new)[1]
bic_selection = step(model_full, direction = "both", trace = 0, k = log(n))$call$formula
print(bic_selection)<jupyter_output>BODYFAT ~ AGE + ABDOMEN + WRIST + HEIGHT<jupyter_text>3.3 LASSOLet $\beta$ denote the parameters, we have:  **LASSO** (least absolute shrinkage and selection operator) regression uses an objective function as:$$f(\beta) = RSS_\beta + \lambda||\beta||_1 \tag{6}$$Find the $\beta$ that minimizes $f(\beta)$. Since LASSO can force some parameters to become zero, we can use LASSO to perform feature selection.<jupyter_code>#lasso
X = model.matrix(model_full)[,-1]
Y = data_new$BODYFAT
set.seed(123)
cv_lasso = cv.glmnet(X, Y, alpha = 1, family = "gaussian", nfolds = 10)
lasso.coef = coef(cv_lasso, s = cv_lasso$lambda.1se)
print(lasso.coef)<jupyter_output>15 x 1 sparse Matrix of class "dgCMatrix"
                       1
(Intercept) -10.81225367
AGE           0.02143009
WEIGHT        .         
HEIGHT       -0.23939696
ADIPOSITY     .         
NECK          .         
CHEST         .         
ABDOMEN       0.60378632
HIP           .         
THIGH         .         
KNEE          .         
ANKLE         .         
BICEPS        .         
FOREARM       .         
WRIST        -0.55902587<jupyter_text>It is not difficult to find that **LASSO** and **BIC** methods give the same feature selection result. Then, we will fit and compare the AIC and BIC models. 3.4 Model Diagnostics AIC model:Summary:<jupyter_code>#Diagnosis
fit_aic = lm(aic_selection, data=data_new)
summary(fit_aic) #r-squared:0.7418<jupyter_output><empty_output><jupyter_text>**Diagnostic plots:**<jupyter_code># save the plot
png(filename="aic.png")
par(mfrow = c(2,2))
plot(fit_aic)
par(mfrow = c(1,1))
dev.off()<jupyter_output><empty_output><jupyter_text>**CV MSE:**<jupyter_code># CV MSE for AIC model
aic_cv_mean = FindCV(data_new,aic_selection)
paste("CV MSE for AIC model: ", aic_cv_mean)<jupyter_output><empty_output><jupyter_text>**Box-Cox:**   Test if we need a transformation or not<jupyter_code>boxcox(fit_aic)<jupyter_output><empty_output><jupyter_text>BIC model:Summary:<jupyter_code>fit_bic = lm(bic_selection,data=data_new)
summary(fit_bic) #r-squared:0.7336<jupyter_output><empty_output><jupyter_text>**Diagnostic plots:**<jupyter_code># save the plot
png(filename="bic.png")
par(mfrow = c(2,2))
plot(fit_bic)
par(mfrow = c(1,1))
dev.off()<jupyter_output><empty_output><jupyter_text>**CV MSE:**<jupyter_code>bic_cv_mean = FindCV(data_new,bic_selection)
paste("CV MSE for BIC model: ", bic_cv_mean)<jupyter_output><empty_output><jupyter_text>**Box-Cox:**  Test if we need a transformation or not<jupyter_code>boxcox(fit_bic)
confint.lm(fit_bic,level = 0.95)<jupyter_output><empty_output><jupyter_text>3.6 Other ModelsWe also try some other models which we think might be useful:  $$1st : BODYFAT \sim WEIGHT + ABDOMEN$$$$2nd : BODYFAT \sim ABDOMEN:WEIGHT$$$$3rd : BODYFAT \sim ABDOMEN$$Why we try these models:Since by using AIC, BIC, LASSO to select variables, the model according to BIC selection is still with 4 variables, so according to the correlation of bodyfat and other variables, abdomen and weight have relatively higher correlation, so we will try to use this two variables and build other models to see what happen.<jupyter_code>#try other models
model_1 = BODYFAT ~ WEIGHT + ABDOMEN
fit_model_1 = lm(model_1,data=data_new)
rsq_1 = summary(fit_model_1)$r.squared
model1_cv = FindCV(data_new, model_1)

model_2 = BODYFAT ~ ABDOMEN:WEIGHT
fit_model_2 = lm(model_2,data=data_new)
rsq_2 = summary(fit_model_2)$r.squared
model2_cv = FindCV(data_new, model_2)

model_3 = BODYFAT ~ ABDOMEN
fit_model_3 = lm(model_3,data=data_new)
rsq_3 = summary(fit_model_3)$r.squared
model3_cv = FindCV(data_new, model_3)

# model informations
rsq = c(rsq_1, rsq_2, rsq_3)
cv_mse = c(model1_cv, model2_cv, model3_cv)
other_models = data.frame(ID = c("1st", "2nd", "3rd"), R_square = rsq, CV_MSE = cv_mse)
other_models<jupyter_output><empty_output><jupyter_text>The 1st model has the minimal CV_MSE, choose this one to perform Diagnostics:<jupyter_code># for the 1st model
# save the plot
png(filename="1st.png")
par(mfrow = c(2,2))
plot(fit_model_1)
par(mfrow = c(1,1))
dev.off()
confint.lm(fit_model_1,level = 0.95)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13260309278350516, 'fraction_numerical': 0.03550257731958763, 'mean_word_length': 5.789588801399825, 'pattern_counts': {'":': 0, '<': 108, '<?xml version=': 0, '>': 99, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '438533', 'license': 'MIT', 'n_tokens_mistral': 6121, 'n_tokens_neox': 5649, 'n_words': 1439, 'path': 'module2_extension.ipynb', 'repo_name': 'LingfengZhu/STAT-628', 'size': '15520'}