Setup

library(tidyverse)
library(ggpubr)
library(car)
library(mfp)
library(corrplot)
library(papeR)
library(knitr)
library(kableExtra)
library(xtable)
library(olsrr)

data(bodyfat)
bf <- bodyfat

General summary

Data analysis process begins with understanding the data. First important observation is that height is measured in inches and weight in pounds. Other body measurements are in centimeters.

summary(bodyfat)

##       case            brozek           siri          density     
##  Min.   :  1.00   Min.   : 0.00   Min.   : 0.00   Min.   :0.995  
##  1st Qu.: 63.75   1st Qu.:12.80   1st Qu.:12.47   1st Qu.:1.041  
##  Median :126.50   Median :19.00   Median :19.20   Median :1.055  
##  Mean   :126.50   Mean   :18.94   Mean   :19.15   Mean   :1.056  
##  3rd Qu.:189.25   3rd Qu.:24.60   3rd Qu.:25.30   3rd Qu.:1.070  
##  Max.   :252.00   Max.   :45.10   Max.   :47.50   Max.   :1.109  
##       age            weight          height           neck      
##  Min.   :22.00   Min.   :118.5   Min.   :29.50   Min.   :31.10  
##  1st Qu.:35.75   1st Qu.:159.0   1st Qu.:68.25   1st Qu.:36.40  
##  Median :43.00   Median :176.5   Median :70.00   Median :38.00  
##  Mean   :44.88   Mean   :178.9   Mean   :70.15   Mean   :37.99  
##  3rd Qu.:54.00   3rd Qu.:197.0   3rd Qu.:72.25   3rd Qu.:39.42  
##  Max.   :81.00   Max.   :363.1   Max.   :77.75   Max.   :51.20  
##      chest           abdomen            hip            thigh      
##  Min.   : 79.30   Min.   : 69.40   Min.   : 85.0   Min.   :47.20  
##  1st Qu.: 94.35   1st Qu.: 84.58   1st Qu.: 95.5   1st Qu.:56.00  
##  Median : 99.65   Median : 90.95   Median : 99.3   Median :59.00  
##  Mean   :100.82   Mean   : 92.56   Mean   : 99.9   Mean   :59.41  
##  3rd Qu.:105.38   3rd Qu.: 99.33   3rd Qu.:103.5   3rd Qu.:62.35  
##  Max.   :136.20   Max.   :148.10   Max.   :147.7   Max.   :87.30  
##       knee           ankle          biceps         forearm          wrist      
##  Min.   :33.00   Min.   :19.1   Min.   :24.80   Min.   :21.00   Min.   :15.80  
##  1st Qu.:36.98   1st Qu.:22.0   1st Qu.:30.20   1st Qu.:27.30   1st Qu.:17.60  
##  Median :38.50   Median :22.8   Median :32.05   Median :28.70   Median :18.30  
##  Mean   :38.59   Mean   :23.1   Mean   :32.27   Mean   :28.66   Mean   :18.23  
##  3rd Qu.:39.92   3rd Qu.:24.0   3rd Qu.:34.33   3rd Qu.:30.00   3rd Qu.:18.80  
##  Max.   :49.10   Max.   :33.9   Max.   :45.00   Max.   :34.90   Max.   :21.40

Check if there are any NA values:

any(is.na(bodyfat))

## [1] FALSE

Conversion of height and weight to centimeters and kilograms.

bf$height <- bf$height * 2.54
bf$weight <- bf$weight * 0.4536

Adjustments on data set

I drop ‘brozek’ column which was the body fat percentage using Brozek’s equation 457/Density - 414.2.
Preferred body fat percentage formula will be obtained by the Siri’s equation 495/Density - 450.
Column ‘case’ also can be dropped since it won’t be used as indexing tool.

bf <- subset(bf, select = -c(brozek))
bf <- subset(bf, select = -c(case))
bf<- rename(bf, fat = siri)

There is a summary for the adjusted data set:

summary(bf)

##       fat           density           age            weight      
##  Min.   : 0.00   Min.   :0.995   Min.   :22.00   Min.   : 53.75  
##  1st Qu.:12.47   1st Qu.:1.041   1st Qu.:35.75   1st Qu.: 72.12  
##  Median :19.20   Median :1.055   Median :43.00   Median : 80.06  
##  Mean   :19.15   Mean   :1.056   Mean   :44.88   Mean   : 81.16  
##  3rd Qu.:25.30   3rd Qu.:1.070   3rd Qu.:54.00   3rd Qu.: 89.36  
##  Max.   :47.50   Max.   :1.109   Max.   :81.00   Max.   :164.72  
##      height            neck           chest           abdomen      
##  Min.   : 74.93   Min.   :31.10   Min.   : 79.30   Min.   : 69.40  
##  1st Qu.:173.35   1st Qu.:36.40   1st Qu.: 94.35   1st Qu.: 84.58  
##  Median :177.80   Median :38.00   Median : 99.65   Median : 90.95  
##  Mean   :178.18   Mean   :37.99   Mean   :100.82   Mean   : 92.56  
##  3rd Qu.:183.51   3rd Qu.:39.42   3rd Qu.:105.38   3rd Qu.: 99.33  
##  Max.   :197.49   Max.   :51.20   Max.   :136.20   Max.   :148.10  
##       hip            thigh            knee           ankle          biceps     
##  Min.   : 85.0   Min.   :47.20   Min.   :33.00   Min.   :19.1   Min.   :24.80  
##  1st Qu.: 95.5   1st Qu.:56.00   1st Qu.:36.98   1st Qu.:22.0   1st Qu.:30.20  
##  Median : 99.3   Median :59.00   Median :38.50   Median :22.8   Median :32.05  
##  Mean   : 99.9   Mean   :59.41   Mean   :38.59   Mean   :23.1   Mean   :32.27  
##  3rd Qu.:103.5   3rd Qu.:62.35   3rd Qu.:39.92   3rd Qu.:24.0   3rd Qu.:34.33  
##  Max.   :147.7   Max.   :87.30   Max.   :49.10   Max.   :33.9   Max.   :45.00  
##     forearm          wrist      
##  Min.   :21.00   Min.   :15.80  
##  1st Qu.:27.30   1st Qu.:17.60  
##  Median :28.70   Median :18.30  
##  Mean   :28.66   Mean   :18.23  
##  3rd Qu.:30.00   3rd Qu.:18.80  
##  Max.   :34.90   Max.   :21.40

Certain data corrections are necessary. It is assumed that height of 75 cm is not possible for 92 kg man aged 44. Therefore, height is corrected by adding 100 cm. After the change height 175 is close to the median.

bf[42,6] <- bf[42,6] + 100

New variable BMI is introduced - Body Mass Index is calculated by the formula:
BMI = Weight(kg) / Height(m)

bf<- bf %>% mutate(bmi = weight / ((height)/100)^2)

Correlations evaluation

It can be guessed that certain body parts size informs about person’s silhouette. Initially I thought that there is a strong relationship between abdomen, hip, and chest and fat percentage. On the other hand I suspected not much correlation between age and body fat.
Indeed my guess was correct, but it is not enough to select variables for regression model.

abd <- ggplot(bf, aes(x=abdomen, y = fat)) + geom_point() + geom_smooth(method = "lm", color="green")
chst <- ggplot(bf, aes(x=chest, y = fat)) + geom_point()+ geom_smooth(method = "lm", color="red")
hp <- ggplot(bf, aes(x=hip, y = fat)) + geom_point()+ geom_smooth(method = "lm", color="blue")
ag <- ggplot(bf, aes(x=age, y = fat)) + geom_point()+ geom_smooth(method = "lm", color="gold")

ggarrange(abd, chst, hp, ag)

The best way to evaluate correlations in the data set will be to use correlation matrix.
It shows how variables are related to each other.

corrplot(cor(bf), method = "pie", type = "upper", diag =FALSE, title = "Correlation matrix - body fat data")

It can be concluded that in face chest, abdomen and hip are strongly correlated with the body fat percentage.
However, it is necessary to minimize multicollinearity in the regression model. On the plot it can be observed that for example chest and abdomen measurements do not only influence fat level but also are strongly correlated with each other. Possibly only one of them should be included int the model.
Correlation between abdomen and chest:

attach(bf)
cor(abdomen, chest)

## [1] 0.9158277

Model selection using Akaike Information Criterion (AIC)

Second general model with introduced bmi variable.

bf.mod2 <- lm(fat ~ density + age + weight + +height + neck + chest + abdomen + hip + thigh + knee + ankle + biceps + forearm + wrist + bmi)
summary(bf.mod2)

## 
## Call:
## lm(formula = fat ~ density + age + weight + +height + neck + 
##     chest + abdomen + hip + thigh + knee + ankle + biceps + forearm + 
##     wrist + bmi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.4273 -0.3765 -0.1143  0.2050 15.1315 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.476e+02  1.206e+01  37.122   <2e-16 ***
## density     -4.113e+02  8.319e+00 -49.445   <2e-16 ***
## age          1.251e-02  9.669e-03   1.293    0.197    
## weight       8.483e-03  4.203e-02   0.202    0.840    
## height       7.917e-03  2.871e-02   0.276    0.783    
## neck        -1.230e-02  6.937e-02  -0.177    0.859    
## chest        2.906e-02  3.054e-02   0.952    0.342    
## abdomen      2.051e-02  3.269e-02   0.627    0.531    
## hip          2.134e-02  4.408e-02   0.484    0.629    
## thigh       -1.508e-02  4.378e-02  -0.345    0.731    
## knee        -7.513e-03  7.313e-02  -0.103    0.918    
## ankle       -8.288e-02  6.625e-02  -1.251    0.212    
## biceps      -5.370e-02  5.138e-02  -1.045    0.297    
## forearm      3.143e-02  5.949e-02   0.528    0.598    
## wrist        1.406e-04  1.610e-01   0.001    0.999    
## bmi          1.906e-02  6.389e-02   0.298    0.766    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.277 on 236 degrees of freedom
## Multiple R-squared:  0.9781, Adjusted R-squared:  0.9767 
## F-statistic: 702.9 on 15 and 236 DF,  p-value: < 2.2e-16

From the general model we can extract which variables to choose for the model. Optimized model can be obtained by selecting predictors by backward selection.
Step function evaluates Akaike Information Criterion and selects proper predictors.

backward2 <- step(bf.mod2, scope = list(lower ~ density), trace = 0)
backward2

## 
## Call:
## lm(formula = fat ~ density + age + chest + hip + ankle)
## 
## Coefficients:
## (Intercept)      density          age        chest          hip        ankle  
##   449.47601   -413.56974      0.01517      0.04054      0.03339     -0.08120

Obviously selected model fits pretty well to the general model, but it still uses many predictors which may lead to overfitting.

plot(fitted(bf.mod2) ~ fitted(backward2))
abline(0,1)

cor(fitted(bf.mod2), fitted(backward2))

## [1] 0.999874

Third model

From the second model summary it can be concluded that density variable has a large t-value and standard error. Moreover it is strongly related with chest. For that reason I will not use it as a predictor in a third model. Ankle won’t be taken into consideration for now.

Proposed model has age, chest and hip as predictors.

bf.mod3 <- lm(fat ~ age + chest + hip)
summary(bf.mod3)

## 
## Call:
## lm(formula = fat ~ age + chest + hip)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8710  -4.4784   0.0518   4.1996  15.1019 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -63.42655    5.25875 -12.061  < 2e-16 ***
## age           0.15266    0.03098   4.928 1.52e-06 ***
## chest         0.42596    0.08278   5.145 5.43e-07 ***
## hip           0.32809    0.09601   3.417 0.000739 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.685 on 248 degrees of freedom
## Multiple R-squared:  0.5441, Adjusted R-squared:  0.5386 
## F-statistic: 98.66 on 3 and 248 DF,  p-value: < 2.2e-16

AIC of proposed model:

AIC(bf.mod3)

## [1] 1596.955

Linear model coefficients:

backward3 <- step(bf.mod3, scope = list(lower ~ age), trace = 0)
backward3

## 
## Call:
## lm(formula = fat ~ age + chest + hip)
## 
## Coefficients:
## (Intercept)          age        chest          hip  
##    -63.4265       0.1527       0.4260       0.3281

Fit of the full and subset models is compared by looking at the corresponding fitted values.

plot(fitted(bf.mod3) ~ fitted(backward3))
abline(0,1)

qqPlot(bf.mod3)

## [1]  39 216

Model and its subset provide the same response variable values.

cor(fitted(backward3), fitted(bf.mod3))

## [1] 1

For now the fitted slope is described by the equation:
fat = -63.42655 + 0.15266 * age + 0.42596 * chest + 0.32809 * hip

Residuals normality evaluation: Distribution seems to be roughly symmetrical and normal. On the Q-Q plot residuals are also really close to the normal line where they are expected to be.

hist(bf.mod3$residuals)

qqnorm(bf.mod3$residuals)
qqline(bf.mod3$residuals)

From the output, the p-value > 0.05 implying that the distribution of the data are not significantly different from normal distribution. In other words, we can assume the normality.

shapiro.test(bf.mod3$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  bf.mod3$residuals
## W = 0.9942, p-value = 0.4471

Global Validation of Linear Models Assumptions:

library(gvlma)
gvmodel3 <- gvlma(bf.mod3) 
summary(gvmodel3)

## 
## Call:
## lm(formula = fat ~ age + chest + hip)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8710  -4.4784   0.0518   4.1996  15.1019 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -63.42655    5.25875 -12.061  < 2e-16 ***
## age           0.15266    0.03098   4.928 1.52e-06 ***
## chest         0.42596    0.08278   5.145 5.43e-07 ***
## hip           0.32809    0.09601   3.417 0.000739 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.685 on 248 degrees of freedom
## Multiple R-squared:  0.5441, Adjusted R-squared:  0.5386 
## F-statistic: 98.66 on 3 and 248 DF,  p-value: < 2.2e-16
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = bf.mod3) 
## 
##                     Value p-value                   Decision
## Global Stat        9.2525 0.05509    Assumptions acceptable.
## Skewness           0.3161 0.57394    Assumptions acceptable.
## Kurtosis           1.7417 0.18692    Assumptions acceptable.
## Link Function      6.2044 0.01274 Assumptions NOT satisfied!
## Heteroscedasticity 0.9902 0.31970    Assumptions acceptable.

As we can see the model obtained linearity, nearly normal residuals, and constant variability (heteroscedasticity). One assumption is still to be met.

Outliers, leverage, influential observations

Below I examine which observations have significant influence in the fitted regression line. In order to to that I start with method which evaluates the difference in fits.

ols_plot_dffits(bf.mod3)

On the plot it can be seen that observation number 39 is extremely influential. Further investigation of this observation is performed:

kable(bf[which(dffits(bf.mod3) < -0.5),]) %>% kable_styling(bootstrap_options = c("bordered", "striped"))

outlierTest(bf.mod3)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##     rstudent unadjusted p-value Bonferroni p
## 39 -2.995523          0.0030185      0.76067

Cook’s distance informs that deleting observation #39 will significantly change fitted slope.

ols_plot_cooksd_bar(bf.mod3)

Studentized residuals vs leverage plot shows which obs. are influencial and which are outliers.

ols_plot_resid_lev(bf.mod3)

In the table we can see three selected outliers, but only one has extremely large Cook’s distance.

influencePlot(bf.mod3)

##       StudRes        Hat      CookD
## 39  -2.995523 0.21284208 0.58767636
## 59   1.035629 0.06166681 0.01761635
## 216  2.723931 0.02422004 0.04488034

Further residuals plots confirm that observation #39 should be adjusted or removed.

ols_plot_resid_stud(bf.mod3)

ols_plot_resid_stand(bf.mod3)

Since there is only one extremely high leverage outlier data set won’t be transformed.
Observation #39 will be removed from the data set.