List of required packages
Package Version
knitr 1.51
ggplot2 4.0.3
moments 0.14.1
dplyr 1.2.1
kableExtra 1.4.0
tidyr 1.3.2
broom 1.0.13
car 3.1.5

Summary

Objective: Create a statistical model capable of accurately predicting newborn birth weight, based on clinical variables collected from three hospitals. The project aims to improve the management of high-risk pregnancies, optimize hospital resources, and ensure better neonatal health outcomes.

The project is part of a growing focus on preventing neonatal complications. Predicting newborn birth weight represents a key opportunity to improve clinical planning and reduce the risks associated with problematic births, such as premature births or low birth weight infants. The main benefits this project will bring to the company and the healthcare sector are:

  1. Improved clinical predictions

    • Newborn weight is a key indicator of their health. Having an accurate predictive model allows medical staff to intervene promptly in case of abnormalities, reducing perinatal complications such as breathing difficulties or hypoglycemia.

  2. Optimizing hospital resources

    • Knowing in advance which newborns may require intensive care helps hospitals organize their human and technological resources efficiently. This translates into reduced operating costs and better planning for neonatal intensive care unit (NICU) utilization.

  3. Prevention and identification of risk factors

    • The model will be able to identify the factors that most negatively impact newborn weight (such as maternal smoking, multiple pregnancies, or advanced maternal age). This information is valuable for prevention and personalized pregnancy management, allowing for proactive interventions in cases of high risk.

  4. Evaluation of hospital practices

    • Through a comparative analysis between the three hospitals involved, the company will be able to identify any differences in clinical outcomes, such as a higher incidence of cesarean sections at a given facility. This allows for monitoring the quality of practices and harmonizing protocols across different hospitals, improving the consistency of care.

  5. Supporting Strategic Planning

    • Data analysis and forecasting can be used to make informed decisions not only at the clinical level but also at the strategic level. The company will be able to leverage this information to implement new public health policies, ensuring a positive impact on neonatal mortality and morbidity rates.

1. Data collection and dataset structure

The table below shows a subset of the data in the dataset:

Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
199 24 0 0 40 3550 500 348 Nat osp2 M
215 29 1 0 38 3300 500 340 Nat osp3 F
802 22 0 0 40 3150 505 322 Ces osp3 F
810 32 0 0 40 3450 515 344 Nat osp3 M
2118 27 0 0 40 3370 505 332 Nat osp1 M

There are no NA values within the dataset:

##   Anni.madre N.gravidanze    Fumatrici   Gestazione         Peso    Lunghezza 
##            0            0            0            0            0            0 
##       Cranio   Tipo.parto     Ospedale        Sesso 
##            0            0            0            0

The dataset consists of 2,500 observations on 10 variables, that can be described as below:

Description of the study variables
Variable Description Type
Mother’s age Measure of age in years Quantitative discrete
Number of pregnancies How many pregnancies the mother has had Quantitative discrete
Maternal smoking A binary indicator (0 = non-smoker, 1 = smoker) Qualitative nominal
Length of pregnancy Number of weeks of gestation Quantitative discrete
Newborn weight Newborn’s weight in grams Quantitative continuous
Newborn length Newborn’s length in centimeters Quantitative continuous
Skull diameter Newborn’s cranial diameter, which can also be measured during pregnancy via ultrasound Quantitative continuous
Type of birth Natural or C-section Qualitative nominal
Hospital of birth Hospital 1, 2, or 3 Qualitative nominal
Newborn’s gender Male (M) or Female (F) Qualitative nominal

2. Analysis and modeling

2.1 Preliminary Analysis

2.1.1 Analysis for quantitative variables
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
Min. : 0.00 Min. : 0.0000 Min. :0.0000 Min. :25.00 Min. : 830 Min. :310.0 Min. :235 Length:2500 Length:2500 Length:2500
1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:38.00 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Class :character Class :character Class :character
Median :28.00 Median : 1.0000 Median :0.0000 Median :39.00 Median :3300 Median :500.0 Median :340 Mode :character Mode :character Mode :character
Mean :28.16 Mean : 0.9812 Mean :0.0416 Mean :38.98 Mean :3284 Mean :494.7 Mean :340 NA NA NA
3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:40.00 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350 NA NA NA
Max. :46.00 Max. :12.0000 Max. :1.0000 Max. :43.00 Max. :4930 Max. :565.0 Max. :390 NA NA NA

Looking at the table above, you can see an anomaly for the variable Anni.madre, which has a minimum value of 0. We then proceed with checking for outliers to identify those due to measurement or transcription errors within the dataset.

Outliers for variable Anni.madre
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
138 13 0 0 38 2760 470 325 Nat osp2 F
205 45 2 0 38 3850 505 384 Nat osp3 M
230 43 1 0 35 2050 455 305 Nat osp2 F
260 44 1 0 40 3500 480 346 Nat osp1 F
335 44 0 1 38 3150 465 335 Nat osp3 F
855 43 0 0 38 3600 510 336 Nat osp3 M
1075 14 1 0 39 3510 490 365 Nat osp2 M
1106 46 5 0 36 2710 470 347 Nat osp2 M
1152 1 1 0 41 3250 490 350 Nat osp2 F
1380 0 0 0 39 3060 490 330 Nat osp3 M
1532 14 0 0 39 3550 500 355 Ces osp1 M
2026 44 0 0 40 3050 505 345 Nat osp1 M
2098 44 5 0 38 3280 475 346 Ces osp1 F

The values 1 and 0 for the mothers’ ages (ID 1152 and ID 1380) can be removed as they are biologically impossible and represent obvious data entry errors.

On the contrary, all other identified outliers (adolescent mothers and those older than 43 years) should be considered real values as they are biologically plausible and reflect the variability of the collected data.

# remove wrong values for variable Anni.madre
df <- df %>% filter(Anni.madre >= 13)
summary(df$Anni.madre)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13.00   25.00   28.00   28.19   32.00   46.00

The table below summarizes each quantitative variable using measures of position (mean, median, quartiles, min and max), dispersion (variance, standard deviation, IQR and CV), and shape (skewness and excess kurtosis).

Descriptive statistics for quantitative variables
Variable Mean Median Q1 Q3 Min Max Variance SD IQR CV Skewness Excess Kurtosis
Anni.madre 28.19 28 25 32 13 46 27.22 5.22 7 0.19 0.15 -0.11
N.gravidanze 0.98 1 0 1 0 12 1.64 1.28 1 1.30 2.51 10.98
Gestazione 38.98 39 38 40 25 43 3.49 1.87 2 0.05 -2.07 8.26
Peso 3284.18 3300 2990 3620 830 4930 275865.90 525.23 630 0.16 -0.65 2.03
Lunghezza 494.70 500 480 510 310 565 693.21 26.33 30 0.05 -1.51 6.48
Cranio 340.03 340 330 350 235 390 269.93 16.43 20 0.05 -0.79 2.94

Based on density graph below and the above table, we can assume Anni.madre is clearly close to normal distribution with mean ≈ median, skewness ≈ 0.15 and excess kurtosis ≈ −0.11. Variables Peso and Cranio approach a bell shape curve but still skewness.

2.1.2 Analysis for qualitative variables

In order to analyze the remaining qualitative variables related the data frame, it is appropriate to compute frequency distributions.

As shown in the frequency tables below, the distribution of qualitative variables is balanced across the main categories.

Regarding the newborn’s sex, a balanced distribution is observed between the two categories, with percentages of ≈ 49.76% and ≈ 50.24% respectively. The distribution across hospitals is also homogeneous, with each facility contributing approximately one-third of the observations.

Frequency distribution - Sesso
ni fi Ni Fi
F 1255 0.5024 1255 0.5024
M 1243 0.4976 2498 1.0000
Frequency distribution - Ospedale
ni fi Ni Fi
osp1 816 0.3267 816 0.3267
osp2 848 0.3395 1664 0.6661
osp3 834 0.3339 2498 1.0000

Regarding the type of birth, a prevalence of vaginal births with ≈ 71% is observed compared to cesarean sections (≈ 29%), consistent with what is generally observed in clinical practice.

Frequency distribution - Tipo.parto
ni fi Ni Fi
Ces 728 0.2914 728 0.2914
Nat 1770 0.7086 2498 1.0000

Finally, the variable relating to maternal smoking highlights a clear predominance of non-smokers with ≈ 96%, while the share of smokers is ≈ 4%.

Frequency distribution - Fumatrici
ni fi Ni Fi
0 2394 0.9584 2394 0.9584
1 104 0.0416 2498 1.0000

2.2 Hypotheses tests

In this section, we will test the following hypotheses with appropriate tests.

2.2.1 Some hospitals perform more cesarean sections

The following hypotheses are to be tested:

\[ H_0: \text{The type of delivery is independent of the hospital where it occurs} \] \[ H_1: \text{The type of delivery depends on the hospital where it occurs} \]

res_chi2 <- chisq.test(df$Tipo.parto, df$Ospedale)
print(res_chi2)
## 
##  Pearson's Chi-squared test
## 
## data:  df$Tipo.parto and df$Ospedale
## X-squared = 1.083, df = 2, p-value = 0.5819

p-value = 0.5819 is above the 0.05 significance level, by suggesting that the type of delivery is independent of the hospital where it occurs.

2.2.2 The average weight and length of this sample of newborns are significantly similar to those of the general population

Based on data from the World Health Organization Child Growth Standards (https://www.cdc.gov/growthcharts/who-growth-charts.htm) and the Bambin Gesù Hospital in Rome (https://www.ospedalebambinogesu.it/da-0-a-30-giorni-come-si-presenta-e-come-cresce-80012), we can define reference values for average neonatal weight and length to use as a comparison with our data. For this purpose, we used an average weight of 3300 grams and an average length of 500 millimeters.

Let’s proceed with the below hypotheses:

\[ H_0: \mu = 3300 \quad \] \[ H_1: \mu \neq 3300 \quad \]

mu_newborn_weight <- 3300
# Student's test for weight
res_weight <- t.test(df$Peso, mu = mu_newborn_weight, alternative = "two.sided")
print(res_weight)
## 
##  One Sample t-test
## 
## data:  df$Peso
## t = -1.505, df = 2497, p-value = 0.1324
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.577 3304.791
## sample estimates:
## mean of x 
##  3284.184

With a p-value of 0.13, it is possible to state that the mean newborn weight does not differ significantly from the theoretical population mean, therefore the null hypothesis is not rejected.

\[ H_0: \mu = 500 \quad \] \[ H_1: \mu \neq 500 \quad \]

mu_newborn_length <- 500
# Student's test for length
res_length <- t.test(df$Lunghezza, mu = mu_newborn_length, alternative = "two.sided")
print(res_length)
## 
##  One Sample t-test
## 
## data:  df$Lunghezza
## t = -10.069, df = 2497, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.6628 495.7287
## sample estimates:
## mean of x 
##  494.6958

The mean newborn length is significantly shorter than the hypothesized mean of 500 mm (p < 0.001). Therefore, the null hypothesis is rejected.

2.2.3 Anthropometric measurements are significantly different between the two sexes

Let us now analyze the differences between the three anthropometric measurements examined in the dataset (newborn weight, length and diameter of the skull) between the two sexes using the Student test for independent groups:

\[ H_0: \mu_M = \mu_F \quad \] \[ H_1: \mu_M \neq \mu_F \quad \]

## 
##  Welch Two Sample t-test
## 
## data:  df$Peso by Sesso
## t = -12.115, df = 2488.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -287.4841 -207.3844
## sample estimates:
## mean in group F mean in group M 
##        3161.061        3408.496
## 
##  Welch Two Sample t-test
## 
## data:  df$Lunghezza by Sesso
## t = -9.5823, df = 2457.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -11.939001  -7.882672
## sample estimates:
## mean in group F mean in group M 
##        489.7641        499.6750
## 
##  Welch Two Sample t-test
## 
## data:  df$Cranio by Sesso
## t = -7.4366, df = 2489.4, p-value = 1.414e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.110504 -3.560417
## sample estimates:
## mean in group F mean in group M 
##        337.6231        342.4586

The results of the tests reported above lead us to the conclusion that the anthropometric measurements considered are significantly different between the two sexes.

Specifically:

  • the mean weight is ≈ 3161 g for females and 3408 g for males, with a p-value < 2.2e-16.

  • the mean length is ≈ 500 mm for males and ≈ 490 mm for females, with a p-value < 2.2e-16

  • the mean skull diameter is ≈ 338 mm for males and ≈ 342 mm for females, with a p-value = 1.414e-13

Therefore, the null hypothesis is rejected.

2.3 Creating the regression model

A multiple linear regression model will be developed that includes all relevant variables. This will allow us to quantify the impact of each independent variable on newborn weight and any interactions. For example, we expect that a longer gestation period would increase newborn weight on average.

2.3.1 Normality test for weight variable

Let’s check that the model’s response variable (newborn weight) is approximately normal, considering the shape indices, as per section 2.1.1, and performing the Shapiro-Wilk normality test:

Skewness and Excess Kurtosis for newborn’s weight
Variable Skewness Excess Kurtosis
Peso -0.65 2.03 
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Peso
## W = 0.97068, p-value < 2.2e-16

With a p-value < 2.2e-16, the test leads us to reject the null hypothesis, so the weight variable does not appear to follow a Normal distribution.

However, given the large sample size, the test becomes sensitive even to small deviations from the Normal distribution.

Therefore, by combining a visual analysis with a histogram graph, it is possible to assert that the weight variable tends to follow a Normal distribution that is moderately left-skewed and leptokurtic.

2.3.2 Correlation matrix for numeric variables

We use the correlation matrix tool to assess the linear relationship between the numerical variables in the dataset (thus excluding categorical variables) according to the Pearson correlation coefficient.

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r)
}

# consider only quantitative variables to be used for the correlation analysis
df_numeric <- df %>% select(all_of(vars_quantitative))
pairs(df_numeric, upper.panel = panel.smooth, lower.panel = panel.cor)

The analysis shows that newborn weight is correlated with both length and skull diameter, suggesting that longer newborns with broader skulls tend to have higher weight. Furthermore, gestational length also shows a positive correlation with weight, indicating that a greater number of weeks of gestation is associated with a slight increase in birth weight.

2.3.3 Model definition

We define the first model by excluding the variables “Type of birth” and “Hospital” which have no direct predictive value with respect to the newborn’s weight:

mod1 <- lm(Peso ~ .-Tipo.parto -Ospedale, data = df)
summary (mod1)
## 
## Call:
## lm(formula = Peso ~ . - Tipo.parto - Ospedale, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1160.6  -181.3   -15.7   163.6  2630.7 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6712.2405   141.3339 -47.492  < 2e-16 ***
## Anni.madre       0.8803     1.1491   0.766    0.444    
## N.gravidanze    11.3789     4.6767   2.433    0.015 *  
## Fumatrici      -30.3958    27.6080  -1.101    0.271    
## Gestazione      32.9472     3.8288   8.605  < 2e-16 ***
## Lunghezza       10.2316     0.3011  33.979  < 2e-16 ***
## Cranio          10.5198     0.4271  24.633  < 2e-16 ***
## SessoM          78.0787    11.2132   6.963 4.24e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2490 degrees of freedom
## Multiple R-squared:  0.7272, Adjusted R-squared:  0.7264 
## F-statistic: 948.3 on 7 and 2490 DF,  p-value: < 2.2e-16

The F-statistic = 948.30, with a p-value < 2.2e-16, confirms that the overall model is highly significant. The value of Adjusted R² = 0.7264 indicates that a large part of the variance of the dependent variable “newborn weight” is explained by the independent variables of the model.

From the model summary, the following variables are statistically relevant, showing a significant linear relationship with newborn’s weight:

  • number of pregnancies (p-value = 0.015)
  • length of pregnancy, length of newborn, skull diameter (p-value < 2.2e-16)
  • male newborn tend to have a higher average weight than female (p-value < 2.2e-16)

The number of previous pregnancies also appears to have some influence, but much less so than the other variables listed above.

# check for multicollinearity
vif(mod1)
##   Anni.madre N.gravidanze    Fumatrici   Gestazione    Lunghezza       Cranio 
##     1.189264     1.187447     1.006692     1.694331     2.079749     1.628987 
##        Sesso 
##     1.040493

The VIF values are below the threshold of 5, indicating there are no multicollinearity effects.

2.3.4 Optimal model selection

We define a second model by removing the variables “mother’s age” and “maternal smoking” which in the previous analysis showed very high p-values, therefore indicating that they do not have a significant impact on the newborn’s weight.

mod2 <- update(mod1, ~ .-Anni.madre -Fumatrici)
summary (mod2)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.37  -180.98   -15.57   163.69  2639.09 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.7251   135.8036 -49.201  < 2e-16 ***
## N.gravidanze    12.4554     4.3416   2.869  0.00415 ** 
## Gestazione      32.3827     3.8008   8.520  < 2e-16 ***
## Lunghezza       10.2455     0.3008  34.059  < 2e-16 ***
## Cranio          10.5410     0.4265  24.717  < 2e-16 ***
## SessoM          77.9807    11.2111   6.956 4.47e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1327 on 5 and 2492 DF,  p-value: < 2.2e-16

The F-statistic = 1327, with a p-value < 2.2e-16, confirms that also this 2nd model is highly significant. The value of Adjusted R² = 0.7265 is substantially unchanged compared to the first model.

# check for multicollinearity
vif(mod2)
## N.gravidanze   Gestazione    Lunghezza       Cranio        Sesso 
##     1.023462     1.669779     2.075747     1.624568     1.040184

Also for this 2nd model, the VIF values are below the threshold of 5, indicating there are no multicollinearity effects.

We proceed with further evaluations by performing the ANOVA (analysis of variance), AIC (akaike information criterion) and BIC (bayesian information criterion) tests:

anova(mod2, mod1)
## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ (Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso) - Tipo.parto - 
##     Ospedale
##   Res.Df       RSS Df Sum of Sq      F Pr(>F)
## 1   2492 188042054                           
## 2   2490 187905214  2    136840 0.9067  0.404
AIC(mod2, mod1)
##      df      AIC
## mod2  7 35152.89
## mod1  9 35155.07
BIC(mod2, mod1)
##      df      BIC
## mod2  7 35193.65
## mod1  9 35207.48

The ANOVA test reports a p-value = 0.404, confirming what was deduced from the comparison between the Adjusted R²: the difference between the two models is not statistically significant and the simplification does not involve a substantial loss of information.

Even in the comparison between the AIC and BIC tests, the lower values for mod2 suggest a preference for the simpler one.

2.3.5 Conclusions

As a result of above tests, we can conclude that mod2 does not imply a loss of prediction compared to mod1, despite being a simplified model. Therefore, mod2 appears to be the preferred model.

2.4 Residue analysis

In this section, we will perform a residual analysis of the chosen regression model (mod2) by studying the differences between the observed and predicted values to verify whether the basic statistical assumptions are met.

We will use graphs to understand the behavior of the residuals:

par(mfrow=c(2,2))
plot(mod2)

  • Residuals vs. Fitted and Scale-Location: a tendency for the residuals and standardized residuals, to become more dispersed with increasing weight values is observed, with a slight curvature, which could indicate potential heteroskedasticity. Additionally, some extreme points are also present.

  • Q-Q residuals: the residuals are arranged along the bisector of the graph, with some deviations in the tails. Combined with the result of the Shapiro-Wilk test below and the density graph, we can assume a slight deviation from perfect normality of the residuals due to the large dataset size.

# check the Normality assumption for the model residuals
shapiro.test(residuals(mod2))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod2)
## W = 0.97414, p-value < 2.2e-16
plot(density(residuals(mod2)))

- Residuals vs. Leverage: there are some observations that show high leverage values based on their Cook’s distance values. In particular the value for the observation #1549 is between the warning and alarm thresholds.These could represent residuals of influential observations.

Let’s proceed by using numerical methods to further investigate leverage values and outliers.

2.4.1 Leverage values

We can check which observations tend to be distant from the others in the regressor space, so those observations with leverage values higher than the threshold value:

\[ \text{threshold value} = \frac{2p}{n} \] where:

  • n: number of observations
  • p: number of estimated parameters including the intercept
# function definition for leverage
leverage.for.model <- function(mod, print_values = TRUE){
  lev <- hatvalues(mod)
  plot(lev)
  # definition for threshold value used to classify values as leverage or not
  p <- sum(lev)
  n <- length(lev)
  threshold.value = 2*p/n
  
  abline(h=threshold.value, col=2)
  
  if(print_values){
    # print leverage values above the threshold
    lev[lev>threshold.value]
  }
}

leverage.for.model(mod2)

##          13          15          34          67          89          96 
## 0.005632888 0.007056483 0.006748974 0.005896189 0.012817563 0.005351869 
##         101         106         131         134         151         155 
## 0.007528094 0.014487904 0.007237755 0.007553514 0.010889886 0.007209736 
##         161         189         190         204         205         206 
## 0.020341133 0.004894598 0.005366831 0.014494112 0.005351828 0.009482503 
##         220         294         305         310         312         315 
## 0.007403130 0.005914307 0.005445189 0.028815300 0.013173723 0.005386836 
##         378         440         442         445         486         492 
## 0.015942080 0.005405733 0.007725708 0.007511227 0.005165714 0.008274392 
##         497         516         582         587         592         614 
## 0.005167522 0.013080707 0.011667393 0.008415872 0.006385013 0.005300091 
##         638         656         657         684         697         702 
## 0.006693153 0.005934494 0.005323517 0.008825948 0.005864920 0.005203163 
##         729         748         750         757         765         805 
## 0.005024448 0.008567254 0.006943755 0.008146487 0.006076412 0.014358658 
##         828         893         895         913         928         946 
## 0.007180179 0.005076266 0.005297643 0.005574016 0.022745493 0.006910965 
##         947         956         985        1008        1014        1049 
## 0.008409518 0.007791530 0.007040035 0.005343994 0.008474086 0.004956561 
##        1067        1091        1106        1130        1165        1180 
## 0.008467035 0.008940030 0.005967511 0.031739177 0.005513581 0.005678043 
##        1187        1199        1218        1237        1247        1272 
## 0.006481533 0.005493009 0.030697778 0.005912496 0.014631359 0.007089553 
##        1290        1292        1310        1320        1324        1355 
## 0.006118589 0.006074423 0.009626391 0.009295311 0.004857340 0.005306940 
##        1356        1383        1393        1398        1400        1409 
## 0.006967713 0.012641828 0.005129491 0.005932524 0.004816704 0.008049539 
##        1418        1426        1427        1448        1503        1549 
## 0.005156670 0.008195421 0.021758961 0.015106684 0.013334256 0.048802841 
##        1551        1554        1571        1591        1604        1608 
## 0.008507417 0.005923458 0.005049368 0.005624961 0.005009312 0.008725821 
##        1615        1617        1626        1684        1691        1699 
## 0.004869997 0.015069038 0.005070717 0.009356578 0.005079185 0.010846383 
##        1710        1716        1725        1733        1778        1779 
## 0.006993461 0.006961226 0.013303625 0.004886226 0.025544997 0.016833696 
##        1807        1825        1866        1890        1960        1965 
## 0.008711082 0.006068655 0.005206137 0.005333985 0.005541287 0.005339716 
##        1975        2035        2038        2044        2084        2087 
## 0.006928794 0.004890003 0.011504028 0.005471894 0.013194769 0.006293791 
##        2096        2112        2113        2118        2138        2144 
## 0.005097113 0.013318873 0.011779730 0.018667270 0.006244737 0.005804705 
##        2146        2147        2155        2173        2198        2213 
## 0.007930323 0.013589469 0.005910248 0.032531736 0.011679031 0.004894044 
##        2214        2218        2219        2222        2223        2242 
## 0.008120291 0.005415586 0.021633907 0.005841119 0.005593576 0.006929530 
##        2255        2305        2315        2316        2335        2357 
## 0.006171101 0.013972901 0.007677230 0.004834368 0.005230839 0.010068259 
##        2406        2420        2434        2435        2450        2456 
## 0.009702475 0.021536269 0.004986609 0.023951342 0.023848222 0.008509142 
##        2469        2476 
## 0.020905930 0.005777111
2.4.2 Outliers

In a similar way, we can proceed to identify the outliers, that is, those extreme values of the response variable.

This procedure involves using the rstudent function on the regression model to build another scatter-plot graph, with horizontal thresholds representing the number of standard deviations away from the regression line:

# outliers for lm function definition

outliers.for.model <- function(model, print_values = FALSE, columns = NULL){
  resz.model <- rstudent(model)
  plot(resz.model)
  # 95% of data between ±2 sd
  abline(h=c(-2,2), col="red")
  # 99.7% of data between ±3 sd
  abline(h=c(-3,3), col="blue")

  if(print_values == TRUE && is.null(columns) == FALSE)
    print(df[which(abs(resz.model) > 3), columns])
}

outliers.for.model(mod2, TRUE, all.vars.mod2)

##      Peso N.gravidanze Gestazione Lunghezza Cranio Sesso
## 119  3410            0         40       550    372     M
## 130  4240            2         39       485    352     M
## 155  3610            0         36       410    330     M
## 310  1560            3         28       420    379     F
## 791  4440            1         41       510    335     M
## 1036 4330            0         40       500    355     F
## 1267 3790            1         40       460    332     F
## 1292 4600            3         38       485    380     M
## 1305 4900            0         41       510    352     F
## 1397 2560            2         38       525    349     M
## 1539 4540            0         38       530    343     M
## 1549 4370            1         38       315    374     F
## 1551 4520            4         35       520    360     F
## 1633 3430            2         39       445    322     F
## 1692 3850            1         36       460    334     F
## 1918 4930            0         39       550    350     F
## 2021 4650            1         39       510    354     M
## 2113 1890            1         32       500    309     F
## 2193 3980            1         40       480    335     F
## 2217 2500            1         39       490    352     M
## 2223 3140            0         35       465    290     F
## 2313 2800            0         42       520    340     M

The point cloud appears to be distributed homogeneously around zero, suggesting homoscedasticity. However, several points exceed the |3| threshold, confirming the presence of isolated outliers.

Given the sample size considered, we perform the formal Bonferroni test to correct for multiple testing bias, that identify #3 statistically significant outliers with Bonferroni p-value << 0.05:

res.outlierTest.mod2 <- outlierTest(mod2)
res.outlierTest.mod2
##       rstudent unadjusted p-value Bonferroni p
## 1549 10.046230         2.6345e-23   6.5810e-20
## 155   5.025345         5.3818e-07   1.3444e-03
## 1305  4.824963         1.4848e-06   3.7092e-03
df[as.numeric(names(res.outlierTest.mod2$rstudent)), all.vars.mod2]
##      Peso N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1549 4370            1         38       315    374     F
## 155  3610            0         36       410    330     M
## 1305 4900            0         41       510    352     F

Let’s continue analyzing the model using Cook’s distance, which indicates how much the entire regression model would change if a single point were removed. A high Cook’s distance implies that the specific observation has a significant and potentially critical impact on the coefficients of the linear regression model.

# cook's distance with graph function definition
cooks.distance.graph <- function(lm){
  cook <- cooks.distance(lm)
  plot(cook, ylim = c(0,1))
  # print index value for top 3 observations
  # based on cook's distances
  top3.cook.indices <- order(cook, decreasing = TRUE)[1:3]
                                                      
  text(x = top3.cook.indices, y = cook[top3.cook.indices], 
       labels = top3.cook.indices, col = "red", pos = 3)
  text(x = top3.cook.indices, y = cook[top3.cook.indices], 
       labels = round(cook[top3.cook.indices],4), col = "blue", pos = 4, cex = 0.6)
}

cooks.distance.graph(mod2)

Based on the tests performed, particular attention is paid to observation 1549:

  • Cook’s distance ≈ 0.83: This value defines the newborn as an outlier with an extremely significant impact on the regression coefficients.
  • rstudent = 10.046230 and Bonferroni p-value = 6.5810e-20: This indicates that the observation is approximately 10 standard deviations away from the regression line estimated by the model, with a probability ≈ 0 that the observation belongs to the distribution of the rest of the sample.

These tests reveal a data entry error, and this is confirmed by observing the data for #1549:

##      Peso N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1549 4370            1         38       315    374     F

The weight = 4370g is typical of newborns of advanced gestational age but biologically possible at 38 weeks of gestation. However, it is not compatible with a length = 315 mm, which is even smaller than the skull diameter (374 mm).

Since it is not possible to correct this error, we proceed with removing the data, redefining the regression model and analyzing the residuals.

2.5 Redefinition of the model by outlier removal

We define a second version for the mod2 by removing the outlier identified in the previous paragraph from data:

df_2.1 <- df[rownames(df) != "1549",]
mod2.1 <- lm(Peso ~ .-Tipo.parto -Ospedale -Anni.madre -Fumatrici, data = df_2.1)
summary (mod2.1)
## 
## Call:
## lm(formula = Peso ~ . - Tipo.parto - Ospedale - Anni.madre - 
##     Fumatrici, data = df_2.1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1165.68  -179.74   -12.42   162.92  1410.68 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6683.8326   133.1602 -50.194  < 2e-16 ***
## N.gravidanze    13.1448     4.2576   3.087  0.00204 ** 
## Gestazione      29.6341     3.7369   7.930 3.27e-15 ***
## Lunghezza       10.8899     0.3019  36.077  < 2e-16 ***
## Cranio           9.9192     0.4227  23.465  < 2e-16 ***
## SessoM          78.1376    10.9929   7.108 1.53e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.3 on 2491 degrees of freedom
## Multiple R-squared:  0.7372, Adjusted R-squared:  0.7367 
## F-statistic:  1398 on 5 and 2491 DF,  p-value: < 2.2e-16

The F-statistic = 1398, with a p-value < 2.2e-16, confirms that also the mod2.1 is highly significant. With the removal of observation 1549, there was an improvement with the adjusted R² increased from 0.7264 to 0.7367

We can now proceed with the residual analysis for the regression model mod2.1, as in section #2.4:

par(mfrow=c(2,2))
plot(mod2.1)

# check the Normality assumption for the model residuals
shapiro.test(residuals(mod2.1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod2.1)
## W = 0.98891, p-value = 5.23e-13
leverage.for.model(mod2.1, FALSE)

outliers.for.model(mod2.1, FALSE)

res.outlierTest.mod2.1 <- outlierTest(mod2.1)
res.outlierTest.mod2.1
##       rstudent unadjusted p-value Bonferroni p
## 155   5.285313         1.3639e-07   0.00034056
## 1305  4.940535         8.3077e-07   0.00207440
## 1397 -4.348609         1.4251e-05   0.03558600
df[as.numeric(names(res.outlierTest.mod2.1$rstudent)), all.vars.mod2]
##      Peso N.gravidanze Gestazione Lunghezza Cranio Sesso
## 155  3610            0         36       410    330     M
## 1305 4900            0         41       510    352     F
## 1397 2560            2         38       525    349     M
cooks.distance.graph(mod2.1)

Removing outlier #1549 appears to have slightly improved the fit of the regression model. However, as highlighted by the Residuals Q-Q plot, as for mod2, the residuals show some deviations in the tails. This, combined with the Shapiro-Wilk test result (W = 0.98891, p-value = 5.23e-13), leads us to confirm a slight deviation from perfect normality of the residuals due to the large sample size.

Then, it’s possible to consider a non-linear relationship between the features.

2.6 Updating the model with non-linear relationship

We then update the model, by adding a quadratic (non-linear) relationship with the variables for length of pregnancy and newborn’s length:

mod2.2 <- update(mod2.1, ~. + I(Gestazione^2) + I(Lunghezza^2))
summary(mod2.2)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Gestazione^2) + I(Lunghezza^2), data = df_2.1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1187.6  -180.9   -12.4   164.2  1319.9 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -2.811e+03  9.024e+02  -3.115 0.001860 ** 
## N.gravidanze     1.438e+01  4.219e+00   3.408 0.000664 ***
## Gestazione       2.000e+02  6.621e+01   3.020 0.002552 ** 
## Lunghezza       -1.925e+01  4.537e+00  -4.244 2.28e-05 ***
## Cranio           1.010e+01  4.204e-01  24.011  < 2e-16 ***
## SessoM           7.339e+01  1.093e+01   6.715 2.32e-11 ***
## I(Gestazione^2) -2.130e+00  8.678e-01  -2.454 0.014198 *  
## I(Lunghezza^2)   3.090e-02  4.621e-03   6.687 2.80e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 266.7 on 2489 degrees of freedom
## Multiple R-squared:  0.7426, Adjusted R-squared:  0.7419 
## F-statistic:  1026 on 7 and 2489 DF,  p-value: < 2.2e-16

The coefficients I(Gestazione^2) (p-value ≈ 0.0142) and I(Lunghezza^2) (p-value ≈ 0) pass the significance tests. The Adjusted R² increase from 0.7367 to 0.7419. This means that the quadratic model explains approximately 0.52% more of the total variability for newborn’s weight (≈ 74.19%).

anova(mod2.2, mod2.1)
## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Gestazione^2) + I(Lunghezza^2)
## Model 2: Peso ~ (Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso) - Tipo.parto - 
##     Ospedale - Anni.madre - Fumatrici
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1   2489 177006781                                  
## 2   2491 180719906 -2  -3713125 26.106 6.018e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(mod2.2, mod2.1)
##        df      AIC
## mod2.2  9 34992.81
## mod2.1  7 35040.64
BIC(mod2.2, mod2.1)
##        df      BIC
## mod2.2  9 35045.21
## mod2.1  7 35081.40

The BIC value for mod2.2 (35045.21) is slightly lower than BIC for mod2.1 (35081.40).

par(mfrow=c(2,2))
plot(mod2.2)

# check the Normality assumption for the model residuals
shapiro.test(residuals(mod2.2))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod2.2)
## W = 0.99048, p-value = 8.393e-12
plot(density(residuals(mod2.2)))

As for the results discussed in section #2.4, also for mod2.2, the density graph combined with the result of the Shapiro-Wilk (W = 0.99048, p-value = 8.393e-12), show some deviations in the tails and a not perfect normality of the residuals.

# studentized Breusch-Pagan test
lmtest::bptest(mod2.2)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod2.2
## BP = 17.125, df = 7, p-value = 0.01661
# Durbin-Watson test
lmtest::dwtest(mod2.2)
## 
##  Durbin-Watson test
## 
## data:  mod2.2
## DW = 1.9496, p-value = 0.104
## alternative hypothesis: true autocorrelation is greater than 0

Also, the graphs Residuals vs fitted, Scale-location and the Breusch-Pagan test (BP = 17.125, p-value = 0.01661), confirm a potential heteroskedasticity.

The result of the Durbin-Watson test on model mod2.2 (DW = 1.9496, p-value = 0.104) indicates that the residuals are not autocorrelated.

leverage.for.model(mod2.2, FALSE)

outliers.for.model(mod2.2, FALSE)

res.outlierTest.mod2.2 <- outlierTest(mod2.2)
res.outlierTest.mod2.2
##       rstudent unadjusted p-value Bonferroni p
## 1305  4.978664         6.8407e-07    0.0017081
## 155   4.845078         1.3434e-06    0.0033545
## 1397 -4.476653         7.9239e-06    0.0197860
## 1692  4.344872         1.4495e-05    0.0361940
df[as.numeric(names(res.outlierTest.mod2.2$rstudent)), all.vars.mod2]
##      Peso N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1305 4900            0         41       510    352     F
## 155  3610            0         36       410    330     M
## 1397 2560            2         38       525    349     M
## 1692 3850            1         36       460    334     F
cooks.distance.graph(mod2.2)

While the previous conclusions on the outliers and leverage charts remain valid, it was shown that outlier #1549 acted as a real distorting factor.

Mod2.2 is therefore better than mod2.1.

3 Predictions and results

Once the model has been validated, we will use it to make practical predictions. For example, we could estimate the weight of a newborn girl, considering a mother in her third pregnancy who is due to give birth at 39 weeks.

row <- data.frame(
  N.gravidanze = 3,
  Gestazione = 39,
  Lunghezza = quantitative_vars_table_results[quantitative_vars_table_results$Variable == 'Lunghezza', 'Mean'],
  Cranio = quantitative_vars_table_results[quantitative_vars_table_results$Variable == 'Cranio', 'Mean'],
  Sesso = "F"
)

prediction.mod2.2 <- predict(mod2.2, newdata = row, interval = "confidence", level = 0.95)

# print as table with columns for prediction renamed
table.prediction.mod2.2 <- cbind(row, prediction.mod2.2)
names(table.prediction.mod2.2)[(ncol(table.prediction.mod2.2) - 2):ncol(table.prediction.mod2.2)] <- c("Fit", "CI.lower", "CI.upper")

kable(table.prediction.mod2.2,
      caption = "Prediction with mod2.2 at 95% CI",
      align = "lllllrrr",
      digits = 2,
      booktabs = TRUE)
Prediction with mod2.2 at 95% CI
N.gravidanze Gestazione Lunghezza Cranio Sesso Fit CI.lower CI.upper
3 39 494.7 340.03 F 3262.66 3239.74 3285.58

4 Visualizations

In this section, graphs will be used to show the impact of the number of weeks of gestation and smoking on the predicted weight.

# df_2.1 is the original dataframe with the
# exclusion of outlier #1549
ggplot(data = df_2.1) +
  geom_point(aes(x = Gestazione,
                 y = Peso,
                 col = factor(Fumatrici)), position = "jitter") +
  geom_smooth(aes(x = Gestazione,
                  y = Peso,
                  col = factor(Fumatrici)), se = FALSE, method = "lm", formula = y ~ x) +
  geom_smooth(aes(x = Gestazione,
                  y = Peso), col="black", se = FALSE, method = "lm", formula = y ~ x) +
  labs(title = "Impact of the number of weeks of gestation and smoking on the predicted weight",
       x = "Pregnancy length (weeks)",
       y = "Newborn's weight (g)",
       color = "Maternal smoking") +
  scale_color_manual(values = c("skyblue", "orange"), 
                     labels = c("False", "True")) +
  theme_minimal()

As with the summary for mod1 in Section 2.3.3, the graph above indicates that the maternal smoking feature has little influence on the newborn’s weight. Indeed, it can be observed that the black line, that represents the relationship between the newborn’s weight with the length of gestation, appears to coincide (save for a very slight difference in slope) with the sky blue line (Maternal smoking = False).

5 Conclusions

Scope of this study was to implement a linear regression model capable of predicting a newborn’s weight based on anthropometric characteristics.

Model mod2.2 explains 74.19% of the total variability in newborn weight by using only five predictors (two of which involve a quadratic relationship), demonstrating excellent predictive capability.

The Shapiro-Wilk test rejected perfect normality of the residuals, however on large datasets this test is overcritical. According to the Central Limit Theorem, as sample size increases, regression coefficient estimates tend to follow a normal distribution, even if the residuals are not normally distributed.

However, regarding the presence of heteroscedasticity, the linear regression model could be replaced by a Generalized Linear Model (GLM), capable to better fitting the data and handling distributions other than the Normal distribution.