Project Specifications

Company: Neonatal Health Solutions

Objective: to create a statistical model capable of accurately predicting newborn weight at birth, 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 outcomes for newborn health.

The dataset contains data on \(2500\) newborns from three hospitals. The variables collected include:

Mother’s age: age of mother, measured in years.
Number of pregnancies: how many pregnancies the mother experienced.
Maternal smoking: binary indicator (0 = non-smoker, 1 = smoker).
Duration of pregnancy: number of weeks of gestation.
Newborn weight: birth weight in grams.
Length: newborn length in millimeters.
Skull diameter: diameter of newborn’s skull in millimeters.
Type of delivery: natural or cesarean.
Birthing hospital: hospital 1, 2 or 3.
Newborn sex: male (M) or female (F).

1 Data collection and dataset structure

First of all, we define a function to display the tables more efficiently.

# --------------------------- TABLE DISPLAYING ---------------------------

show_table <- function(data, digits = 2, scrollable = FALSE){
  # data: input data frame to display
  # digits: number of decimal places to show in numeric columns
  # scrollable: logical, whether to enable horizontal scrolling
  
  numeric_cols <- sapply(data, is.numeric)
  data[numeric_cols] <- lapply(data[numeric_cols], round, digits = digits)
  
  table <- kable(data) %>%
    kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "bordered"),
      full_width = FALSE,
      position = "center"
    )
  
  if (scrollable == TRUE) {
    table <- table %>%
      scroll_box(width = "100%", height = "auto")
  }
  
  return(table)
}

# ------------------------------------------------------------------------

Now we load the dataset, inspect its dimensions, check for missing values and show the first ten records.

# Load the dataset
dataset_url <- "https://drive.google.com/uc?export=download&id=1ChfwftuOSH-WLIto_1AvV-_sQIksGeTq"
dataset <- read.csv(dataset_url)

# Inspect dataset size
records <- nrow(dataset)
variables <- ncol(dataset)

# Ckeck for missing values
cat("The dataset",
    ifelse(anyNA(dataset), "contains", "has not"),
    "missing values.\n")

# Display the first ten records
show_table(head(dataset, 10))

## The dataset has not missing values.

Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
26	0	42	3380	490	325	Nat	osp3	M
21	2	39	3150	490	345	Nat	osp1	F
34	3	38	3640	500	375	Nat	osp2	M
28	1	41	3690	515	365	Nat	osp2	M
20	0	38	3700	480	335	Nat	osp3	F
32	0	40	3200	495	340	Nat	osp2	F
26	1	39	3100	480	345	Nat	osp3	F
25	0	40	3580	510	349	Nat	osp1	M
22	1	40	3670	500	335	Ces	osp2	F
23	0	41	3700	510	362	Ces	osp2	F

The dataset includes 10 variables and 2500 records. Let’s display the dataset structure.

# Show the dataset structure
str(dataset)

## 'data.frame':    2500 obs. of  10 variables:
##  $ Anni.madre  : int  26 21 34 28 20 32 26 25 22 23 ...
##  $ N.gravidanze: int  0 2 3 1 0 0 1 0 1 0 ...
##  $ Fumatrici   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Gestazione  : int  42 39 38 41 38 40 39 40 40 41 ...
##  $ Peso        : int  3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
##  $ Lunghezza   : int  490 490 500 515 480 495 480 510 500 510 ...
##  $ Cranio      : int  325 345 375 365 335 340 345 349 335 362 ...
##  $ Tipo.parto  : chr  "Nat" "Nat" "Nat" "Nat" ...
##  $ Ospedale    : chr  "osp3" "osp1" "osp2" "osp2" ...
##  $ Sesso       : chr  "M" "F" "M" "M" ...

Let’s analyze the types of statistical variables contained in the dataset.

Discrete quantitative variables: Anni.madre, N.gravidanze, Gestazione
Continuous quantitative variables: Peso, Lunghezza, Cranio
Categorical nominal variables: Fumatrici, Tipo.parto, Ospedale, Sesso

We can visualize some basic descriptive statistics.

summary(dataset)

##    Anni.madre     N.gravidanze       Fumatrici        Gestazione   
##  Min.   : 0.00   Min.   : 0.0000   Min.   :0.0000   Min.   :25.00  
##  1st Qu.:25.00   1st Qu.: 0.0000   1st Qu.:0.0000   1st Qu.:38.00  
##  Median :28.00   Median : 1.0000   Median :0.0000   Median :39.00  
##  Mean   :28.16   Mean   : 0.9812   Mean   :0.0416   Mean   :38.98  
##  3rd Qu.:32.00   3rd Qu.: 1.0000   3rd Qu.:0.0000   3rd Qu.:40.00  
##  Max.   :46.00   Max.   :12.0000   Max.   :1.0000   Max.   :43.00  
##       Peso        Lunghezza         Cranio        Tipo.parto        Ospedale   
##  Min.   : 830   Min.   :310.0   Min.   :235   Length   :2500   Length   :2500  
##  1st Qu.:2990   1st Qu.:480.0   1st Qu.:330   N.unique :   2   N.unique :   3  
##  Median :3300   Median :500.0   Median :340   N.blank  :   0   N.blank  :   0  
##  Mean   :3284   Mean   :494.7   Mean   :340   Min.nchar:   3   Min.nchar:   4  
##  3rd Qu.:3620   3rd Qu.:510.0   3rd Qu.:350   Max.nchar:   3   Max.nchar:   4  
##  Max.   :4930   Max.   :565.0   Max.   :390                                    
##        Sesso     
##  Length   :2500  
##  N.unique :   2  
##  N.blank  :   0  
##  Min.nchar:   1  
##  Max.nchar:   1  
##

We can observe that the variable Anni.madre has a minimum of zero: this value does not make sense in the context of the study, so let’s examine the frequency distribution of mothers’ ages.

# Frequency distribution of Anni.madre
freq_Anni.madre <- data.frame(table(dataset$Anni.madre))
names(freq_Anni.madre) <- c("Anni.madre", "Frequency")

freq_Anni.madre_transposed <- as.data.frame(t(freq_Anni.madre))
colnames(freq_Anni.madre_transposed) <- NULL

show_table(freq_Anni.madre_transposed, scrollable = TRUE)


Anni.madre	0	1	13	14	15	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45	46
Frequency	1	1	1	2	6	13	18	24	45	66	74	100	115	131	180	184	197	172	174	200	147	159	110	96	66	64	41	38	27	19	13	8	2	4	1	1

The dataset contains two anomalous records: one mother is \(0\) years old and another is \(1\) years old. Since these values are physically impossible, they are likely due to sampling errors during dataset construction. To prevent them from affecting the analysis, these two records will be removed from the dataset.

# Remove the anomalous records
dataset <- dataset%>%
  filter((Anni.madre != 0) & (Anni.madre != 1))

To make the dataset more readable, we’ll convert the values of Fumatrici from (0, 1) to (“Non-smoker”, Smoker”) in the dataset.

# Transform Fumatrici labels
dataset$Fumatrici <- factor(
  dataset$Fumatrici, 
  levels = c(0, 1), 
  labels = c("Non-smoker", "Smoker")
  )

show_table(head(dataset, 10))

Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
26	0	Non-smoker	42	3380	490	325	Nat	osp3	M
21	2	Non-smoker	39	3150	490	345	Nat	osp1	F
34	3	Non-smoker	38	3640	500	375	Nat	osp2	M
28	1	Non-smoker	41	3690	515	365	Nat	osp2	M
20	0	Non-smoker	38	3700	480	335	Nat	osp3	F
32	0	Non-smoker	40	3200	495	340	Nat	osp2	F
26	1	Non-smoker	39	3100	480	345	Nat	osp3	F
25	0	Non-smoker	40	3580	510	349	Nat	osp1	M
22	1	Non-smoker	40	3670	500	335	Ces	osp2	F
23	0	Non-smoker	41	3700	510	362	Ces	osp2	F

Before proceeding with our analysis, let’s attach the dataset to make all the variables accessible in the workspace.

attach(dataset)

2 Analysis and Modeling

2.1 Descriptive Statistics

In this chapter we will analyze the variables to understand their distribution and identify any outliers.
For quantitative variables, we’ll compute the measures of central tendency, variability and shape, together with the number of detected outliers.
For qualitative variables, we’ll compute the absolute and relative frequencies of their categories.

# ----------- DESCRIPTIVE STATISTICS FOR QUANTITATIVE VARIABLES -----------

quantitative_statistics <- function(df, variables){
  # df: data frame containing the data to analyze
  # variables: names of the variable to describe
  
  results <- list()
  
  for(v in variables){
    x <- df[[v]]
    
    if(is.numeric(x)){
      mu <- mean(x, na.rm = TRUE)
      sigma <- sd(x, na.rm = TRUE)
      min_value <- min(x, na.rm = TRUE)
      max_value <- max(x, na.rm = TRUE)
      q1 = as.numeric(quantile(x, probs = 0.25, na.rm = TRUE))
      q3 = as.numeric(quantile(x, probs = 0.75, na.rm = TRUE))
      iqr_value = IQR(x, na.rm = TRUE)
      lower_bound = q1 - 1.5 * iqr_value
      upper_bound = q3 + 1.5 * iqr_value
      n_outliers = sum(x < lower_bound | x > upper_bound, na.rm = TRUE)
      
      results[[v]] <- c(
        # Measures of central tendency
        mean = mu,
        min = min_value,
        q1 = q1,
        median = median(x, na.rm = TRUE),
        q3 = q3,
        max = max_value,
        
        # Measures of variability
        variance = var(x, na.rm = TRUE),
        std_dev = sigma,
        range = max_value - min_value,
        interquartile_range = iqr_value,
        cv = ifelse(mu == 0, NA, sigma/mu),
  
        # Measures of shape
        skewness_index = skewness(x, na.rm = TRUE),
        kurtosis_index = kurtosis(x, na.rm = TRUE) - 3,
        
        # Number of outliers
        n_outliers = n_outliers
      )
    }
  }
  
  stats <- as.data.frame(results)
  
  rownames(stats) <- c(
  "Mean value", "Minimum", "First quartile (Q1)", "Median value", 
  "Third quartile (Q3)", "Maximum", "Variance", "Standard deviation", 
  "Range", "Interquartile range", "Coefficient of variation",
  "Fisher's Skewness", "Excess Kurtosis", "Number of outliers"
  )
  
  return(stats)
}

# -------------------------------------------------------------------------

# Descroptive statistics for quantitative variables
quantitative_summary <- quantitative_statistics(
  dataset,
  c("Anni.madre", "N.gravidanze", "Gestazione","Peso", "Lunghezza", "Cranio")
  )

show_table(quantitative_summary)

	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio
Mean value	28.19	0.98	38.98	3284.18	494.70	340.03
Minimum	13.00	0.00	25.00	830.00	310.00	235.00
First quartile (Q1)	25.00	0.00	38.00	2990.00	480.00	330.00
Median value	28.00	1.00	39.00	3300.00	500.00	340.00
Third quartile (Q3)	32.00	1.00	40.00	3620.00	510.00	350.00
Maximum	46.00	12.00	43.00	4930.00	565.00	390.00
Variance	27.22	1.64	3.49	275865.90	693.21	269.93
Standard deviation	5.22	1.28	1.87	525.23	26.33	16.43
Range	33.00	12.00	18.00	4100.00	255.00	155.00
Interquartile range	7.00	1.00	2.00	630.00	30.00	20.00
Coefficient of variation	0.19	1.30	0.05	0.16	0.05	0.05
Fisher’s Skewness	0.15	2.51	-2.07	-0.65	-1.51	-0.79
Excess Kurtosis	-0.11	10.98	8.26	2.03	6.48	2.94
Number of outliers	11.00	246.00	67.00	69.00	59.00	48.00

From the descriptive statistics showed above, we can extract some useful information:

All quantitative variables present a modest number of outliers. In any case, these observations represent reasonable quantities in the context of our analysis; therefore, since the sampling methodology is unknown, we cannot exclude any outliers. We will discuss how to handle them in the future analysis.
Anni.madre appears to have an approximately normal distribution, with very low absolute values of skewness and kurtosis, and the lowest number of outliers among the quantitative variables.
N.gravidanze has a modest positive skewness and a very high positive kurtosis, deviating strongly from a normal distribution; moreover, it has the highest value of the coefficient of variation, suggesting high dispersion, and it shows the highest number of outliers among the quantitative variables.
Gestazione, Peso, Lunghezza and Cranio are all left-skewed and show moderate to high positive kurtosis.

Below, we show the density function of Peso, which will be the target variable in the the multiple linear regression model.

# Density plot of newborns weight
ggplot(dataset, aes(x = Peso)) +
  geom_density(col = "black", fill = "skyblue") +
  labs(
    title = "Density plot of newborns weight",
    x = "Weight (g)",
    y = "Density"
  ) +
  scale_x_continuous(breaks = seq(0, max(Peso), by = 500)) +
  theme_light() +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5, face = "bold", color = "red3"),
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title.x = element_text(size = 13, margin = margin(t = 10)),
    axis.title.y = element_text(size = 13, margin = margin(r = 10))
  )

Now let’s take a look at the categorical variables, showing the absolute and relative frequencies of their categories.

# ----------- DESCRIPTIVE STATISTICS FOR CATEGORICAL VARIABLES ----------- 

categorical_statistics <- function(df, variable){
  # df: data frame containing the data to analyze
  # variable: name of the variable to describe
  
  x <- df[[variable]]
  N <- length(x)
  ni <- table(x)
  fi <- ni/N
  
  freq_summary <- data.frame(
    labels = names(ni),
    abs_freq = as.vector(ni),
    rel_freq = as.vector(fi)
  )
  
  colnames(freq_summary) <- c(variable, "Frequency", "Relative frequency")
    
  return(freq_summary)
}

# ------------------------------------------------------------------------

# Table of frequencies of Fumatrici
show_table(categorical_statistics(dataset, "Fumatrici"))

Fumatrici	Frequency	Relative frequency
Non-smoker	2394	0.96
Smoker	104	0.04

The class representing non-smoking mothers accounts for the vast majority of the data.

# Table of frequencies of Tipo.parto
show_table(categorical_statistics(dataset, "Tipo.parto"))

Tipo.parto	Frequency	Relative frequency
Ces	728	0.29
Nat	1770	0.71

We can observe that the natural delivery is the most frequent category.

# Table of frequencies of Ospedale
show_table(categorical_statistics(dataset, "Ospedale"))

Ospedale	Frequency	Relative frequency
osp1	816	0.33
osp2	848	0.34
osp3	834	0.33

The observations are uniformly distributed across the three hospitals, with each hospital accounting for about one-third of the sample.

# Table of frequencies of Sesso
show_table(categorical_statistics(dataset, "Sesso"), 3)

Sesso	Frequency	Relative frequency
F	1255	0.502
M	1243	0.498

The sample is evenly distributed between male and female newborns. This is helpful because sex will be an important control variable in the multiple linear regression model.

The next step is to test the following hypothesis using appropriate statistical tests:

Some hospitals perform more cesarean sections.
The average weight and length of this sample of newborns are significantly similar to those of the overall population.
Anthropometric measures differ significantly between the two sexes.

2.1.1 Do some hospitals perform more cesarean sections?

Let’s start with the first hypothesis test: it translates into testing for independence between the variables Ospedale and Tipo.parto. The system of assumptions is as follows:

H0: Ospedale and Tipo.parto are independent;
H1: Ospedale and Tipo.parto are not independent.

To test this hypothesis, we will use the Chi-square test with a significance level of \(5\)%. First, let’s look at the contingency table.

contingency_Ospedale_Tipo.parto <- table(Tipo.parto, Ospedale)
kable(contingency_Ospedale_Tipo.parto)

	osp1	osp2	osp3
Ces	242	254	232
Nat	574	594	602

# Chi-square test for indipendence between Ospedale and Tipo.parto
indipendence_test <- chisq.test(contingency_Ospedale_Tipo.parto)

indipendence_test_results <- data.frame(
  Chi_square = indipendence_test$statistic,
  p_value = indipendence_test$p.value,
  Degrees_of_freedom = indipendence_test$parameter,
  row.names = NULL
)

show_table(indipendence_test_results, 3)

Chi_square	p_value	Degrees_of_freedom
1.083	0.582	2

The Chi-square test yielded a p-value \(\simeq 0.58\): this means that we don’t have sufficient evidence to reject the null hypothesis of independence with \(5\)% significance level. Therefore, we do not observe a statistically significant association between the variables Ospedale and Tipo.parto, so there is no evidence of differences in the number of cesarean sections among the three hospitals.

2.1.2 Are weight and length of this sample of newborns similar to those of the population?

Now we want to test if the average weight and length of this sample of newborns are significantly similar to those of the overall population. According to the medical literature, the average weight and length of newborns are \(3.3\) kg and \(50\) cm, respectively.

Let’s check if Peso and Lunghezza follow a normal-like distribution, by using a Shapiro-Wilk test.

cat("SHAPIRO-WILK TESTS RESULTS \n")
shapiro_Peso <- shapiro.test(Peso)
cat("p-value (Peso) =", 
    format(shapiro_Peso$p.value, scientific=TRUE, digits=4),
    "\n")
shapiro_Lunghezza <- shapiro.test(Lunghezza)
cat("p-value (Lunghezza) =", 
    format(shapiro_Lunghezza$p.value, scientific=TRUE, digits=4))

## SHAPIRO-WILK TESTS RESULTS 
## p-value (Peso) = 3.379e-22 
## p-value (Lunghezza) = 3.313e-36

The p-values of Shapiro-Wilk tests on both Peso and Lunghezza are extremely low, showing very strong evidence against normality hypothesis. Although the normality assumption is not satisfied, the t-test may still be reasonably robust in our case because of the large sample size of 2498 observations. In fact, the distribution of the sample mean is approximately normal for large samples, as stated by the Central Limit Theorem, and the t-Student distribution tends toward the standard normal distribution as the degrees of freedom increase.

Thus, we will perform both the t-test and the Wilcoxon test, which is non-parametric so it does not require the normality assumption, and then we will compare them and discuss the results.

Let’s start considering the newborns weight, the system of assumptions is as follows:

H0: Mean weight of sample is equal to \(3300\) g ;
H1: Mean weight of sample is different from \(3300\) g.

To test this hypothesis, we will first use the two-sided t-test with a significance level of \(5\)%. Let’s create a function to efficiently summarize the results of a t-test.

# -------------------- SUMMARIZE METRICS FOR A T-TEST --------------------

summarize_t_test <- function(t_test){
  # t_test: t.test to be summarized
  
  p_value <- ifelse(
    t_test$p.value > 10^(-4),
    round(t_test$p.value, 4),
    "< 0.0001"
  )
  conf_int <- paste0("[", round(t_test$conf.int[1], 2), ", ", round(t_test$conf.int[2], 2), "]")
  dof <- round(t_test$parameter, 2)
    
  # One-sample t-test
  if(length(t_test$estimate) == 1){
    population_mean <- round(t_test$null.value, 2)
    sample_mean <- round(t_test$estimate, 2)
    
    summary_table <- data.frame(
    Metric = c("Population mean", "Sample mean", "p-value", 
                "Confidence interval", "Degrees of freedom"),
    Value = c(
      population_mean,
      sample_mean,
      p_value,
      conf_int,
      dof
    )
  )
}
  
  # Two-sample t-test
  else if(length(t_test$estimate) == 2){
    group_names <- names(t_test$estimate)
    mean_group1 <- round(as.numeric(t_test$estimate[1]), 2)
    mean_group2 <- round(as.numeric(t_test$estimate[2]), 2)
    mean_diff <- round(as.numeric(t_test$estimate[1] - t_test$estimate[2]), 2)
    
    summary_table <- data.frame(
    Metric = c(group_names[1], group_names[2], "Mean difference", 
                "p-value", "Confidence interval", "Degrees of freedom"),
    Value = c(
      mean_group1,
      mean_group2,
      mean_diff,
      p_value,
      conf_int,
      dof
    )
  )
  }

  return(summary_table)
}

# -----------------------------------------------------------------------

Now we can compute the t-test described before.

t_test_Peso <- t.test(
  Peso,
  mu = 3300,
  conf.level = 0.95,
  alternative = "two.sided")

show_table(summarize_t_test(t_test_Peso))

Metric	Value
Population mean	3300
Sample mean	3284.18
p-value	0.1324
Confidence interval	[3263.58, 3304.79]
Degrees of freedom	2497

With a p-value \(\simeq 0.13 > 0.05\), we don’t reject the null hypothesis. Let’s check the results using the Wilcoxon test.

wilcox_test_Peso <- wilcox.test(Peso, mu=3300)
cat("WILCOXON TEST RESULTS \n")
cat("p-value = ", round(wilcox_test_Peso$p.value, 3))

## WILCOXON TEST RESULTS 
## p-value =  0.948

Also the Wilcoxon test don’t provide sufficient evidence to reject the null hypothesis. Therefore, there is insufficient evidence to conclude that the newborns weight differs from the population reference value.

Now, let’s consider the newborns length, the system of assumptions is:

H0: Mean length of sample is equal to \(500\) mm ;
H1: Mean length of sample is different from \(500\) mm.

To test this hypothesis, once again, we will use the two-sided t-test with a significance level of \(5\)%.

t_test_Lunghezza <- t.test(
  Lunghezza,
  mu = 500,
  conf.level = 0.95,
  alternative = "two.sided")

show_table(summarize_t_test(t_test_Lunghezza))

Metric	Value
Population mean	500
Sample mean	494.7
p-value	< 0.0001
Confidence interval	[493.66, 495.73]
Degrees of freedom	2497

With a p-value very close to zero, which is below any conventional significance level, we reject the null hypothesis. Let’s check the results using the Wilcoxon test.

wilcox_test_Lunghezza <- wilcox.test(Lunghezza, mu=500)
cat("WILCOXON TEST RESULTS \n")
cat("p-value = ", format(wilcox_test_Lunghezza$p.value, scientific=TRUE, digits=3))

## WILCOXON TEST RESULTS 
## p-value =  1.42e-16

The Wilcoxon test also exhibits a p-value close to zero, providing strong evidence against the null hypothesis. There is sufficient reason to conclude that newborns length differs from the population reference value.

2.1.3 Do anthropometric measures differ significantly by sex?

We aim to test if body measures like weight, length and skull diameter differ significantly between male and female newborns. This translates into considering the following set of hypotheses for each anthropometric variable of interest:

H0: The anthropometric measure is equal between male and female newborns;
H1: The anthropometric measure is different between male and female newborns.

Previously, we observed that Peso and Lunghezza dont’t respect the normality assumption, let’s check if the variable Cranio follow a normal-like distribution using the Shapiro-Wilk test.

cat("SHAPIRO-WILK TEST RESULTS \n")
shapiro_Cranio <- shapiro.test(Cranio)
cat("p-value (Cranio) =", 
    format(shapiro_Cranio$p.value, scientific=TRUE, digits=4))

## SHAPIRO-WILK TEST RESULTS 
## p-value (Cranio) = 1.296e-24

The p-values of Shapiro-Wilk tests on Cranio is close to zero, suggesting very strong evidence against normality hypothesis. For the same reason explained in the previous section, the t-test may still be sufficiently robust so we’ll perform both the t-test and the Wilcoxon test to check if anthropometric measures differ significantly by sex.

Let’s start computing two-sided t-tests for independent samples (males and females), choosing a \(5\)% significance level, one t-test for each anthropometric variable.

t_test_Peso_Sesso <- t.test(Peso ~ Sesso)
show_table(summarize_t_test(t_test_Peso_Sesso))

Metric	Value
mean in group F	3161.06
mean in group M	3408.5
Mean difference	-247.43
p-value	< 0.0001
Confidence interval	[-287.48, -207.38]
Degrees of freedom	2488.67

wilcox_test_Peso_Sesso <- wilcox.test(Peso ~ Sesso)
cat("WILCOXON TEST RESULTS \n")
cat("p-value = ", format(wilcox_test_Peso_Sesso$p.value, scientific=TRUE, digits=3))

## WILCOXON TEST RESULTS 
## p-value =  2.91e-41

Both the t-test and the Wilcoxon test show a p-value smaller than any reasonable significance level, providing strong evidence against the null hypothesis: newborns weight differs significantly between the two sexes. We can observe that males weigh, on average, \(247\) g more than females at birth.

t_test_Lunghezza_Sesso <- t.test(Lunghezza ~ Sesso)
show_table(summarize_t_test(t_test_Lunghezza_Sesso))

Metric	Value
mean in group F	489.76
mean in group M	499.67
Mean difference	-9.91
p-value	< 0.0001
Confidence interval	[-11.94, -7.88]
Degrees of freedom	2457.3

wilcox_test_Lunghezza_Sesso <- wilcox.test(Lunghezza ~ Sesso)
cat("WILCOXON TEST RESULTS \n")
cat("p-value = ", format(wilcox_test_Lunghezza_Sesso$p.value, scientific=TRUE, digits=3))

## WILCOXON TEST RESULTS 
## p-value =  2.66e-25

Both the t-test and the Wilcoxon test show a p-value closer to zero, suggesting strong evidence against the null hypothesis: newborns length differs significantly between the two sexes. We can notice that male newborns length is, on average, \(9.9\) mm higher than female newborns length.

t_test_Cranio_Sesso <- t.test(Cranio ~ Sesso)
show_table(summarize_t_test(t_test_Cranio_Sesso))

Metric	Value
mean in group F	337.62
mean in group M	342.46
Mean difference	-4.84
p-value	< 0.0001
Confidence interval	[-6.11, -3.56]
Degrees of freedom	2489.39

wilcox_test_Cranio_Sesso <- wilcox.test(Cranio ~ Sesso)
cat("WILCOXON TEST RESULTS \n")
cat("p-value = ", format(wilcox_test_Cranio_Sesso$p.value, scientific=TRUE, digits=3))

## WILCOXON TEST RESULTS 
## p-value =  7.14e-15

Both the t-test and the Wilcoxon test show a p-value smaller than any reasonable significance level, providing strong evidence against the null hypothesis: newborns skull diameter differs significantly between the two sexes. We can observe that males have a skull diameter that is, on average, \(4.8\) mm larger than females.

In conclusion, we can say with a high level of confidence that anthropometric measures differ significantly between the two sexes, as we would expect.

2.2 Building of Multiple Linear Regression Model with all relevant variables

The objective of this section is to build a multiple linear regression model to predict newborns weight. We previously verified that the distribution of the target variable Peso is not normal; now let’s display the correlation matrix between the quantitative variables.

# ----------------- SHOW CORRELATION MATRIX WITH PAIRS() -----------------

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  par(usr = c(0, 1, 0, 1))
  r <- 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)
}

# ------------------------------------------------------------------------

# Quantitative variables
dataset_numeric <- dataset[sapply(dataset, is.numeric)]

# Correlation matrix and scatterplots
pairs(dataset_numeric, upper.panel = panel.smooth, lower.panel = panel.cor)

The correlation matrix, evaluated with Pearson’s correlation index, allow us to understand the linear relationship between the quantitative variables.

In particular, we can notice that the target variable Peso shows a strong positive linear correlation with Lunghezza (\(\rho \simeq 0.8\)) and Cranio (\(\rho \simeq 0.7\)), and a moderate positive correlation with Gestazione (\(\rho \simeq 0.59\)), while it doesn’t exhibit a clear linear correlation with Anni.madre and N.gravidanze.

Regarding the independent variables, we observe that Gestazione is moderately correlated with Lunghezza (\(\rho \simeq 0.62\)) and Cranio (\(\rho \simeq 0.46\)), and also Lunghezza and Cranio are moderately correlated (\(\rho \simeq 0.6\)). In the subsequent analysis we will assess whether these correlations between explanatory variables cause multicollinearity problems.

We also note that Peso appears to have a non-linear relationship with the variables Gestazione and Cranio.

Now we build the first multiple linear regression model by including the relevant available predictors, excluding the variable Ospedale, since it cannot be used to predict the newborns weight and it would not be useful for generalizing the model to different samples.

# Initial model
model_1 <- lm(Peso ~. -Ospedale, data = dataset)
summary(model_1)

## 
## Call:
## lm(formula = Peso ~ . - Ospedale, data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1140.10  -181.96   -14.86   160.30  2629.68 
## 
## Coefficients:
##                   Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)     -6738.3356   141.5660 -47.599 < 0.0000000000000002 ***
## Anni.madre          0.8681     1.1479   0.756               0.4496    
## N.gravidanze       11.6900     4.6733   2.501               0.0124 *  
## FumatriciSmoker   -31.7061    27.5836  -1.149               0.2505    
## Gestazione         32.8963     3.8248   8.601 < 0.0000000000000002 ***
## Lunghezza          10.2691     0.3012  34.098 < 0.0000000000000002 ***
## Cranio             10.4850     0.4268  24.564 < 0.0000000000000002 ***
## Tipo.partoNat      30.3855    12.1052   2.510               0.0121 *  
## SessoM             78.0234    11.2013   6.966     0.00000000000417 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.4 on 2489 degrees of freedom
## Multiple R-squared:  0.7279, Adjusted R-squared:  0.727 
## F-statistic: 832.3 on 8 and 2489 DF,  p-value: < 0.00000000000000022

The first model achieved an \(R^2\) and adjusted \(R^2\) of about \(0.73\), which means it can explain approximately \(73\)% of the variability in newborns weight.

Let’s look at the coefficients and the p-values: choosing a significance level of \(5\)% for each predictor to consider it statistically significant, we can observe that six variables are relevant in predicting weight, which are N.gravidanze, Gestazione, Lunghezza, Cranio, Tipo.parto and Sesso, whereas the other two variables, Anni.madre and Fumatrici, show a p-value \(> 0.05\) so we don’t consider them statistically significant.

We can also notice that, holding all other terms constant, male newborns weigh approximately \(78\) grams more than female newborns; moreover, babies born by natural delivery weigh on average about \(30\) grams more than babies born by cesarean section.

2.3 Selecting the Optimal Model

In order to select the optimal multiple linear regression model, we will compare a set of candidate models by evaluating their performance with BIC (Bayesian information criterion), which penalizes the additional parameters more heavily than AIC (Akaike information criterion). The preferred multiple linear regression model will be the one with the lowest BIC value.

As a first step, since Anni.madre and Fumatrici are not statistically significant in the initial model, we consider models obtained by removing them separately and jointly, and we compare their BIC values.

model_2 <- update(model_1, ~. -Anni.madre)
model_3 <- update(model_1, ~. -Fumatrici)
model_4 <- update(model_1, ~. - Anni.madre -Fumatrici)

bic_values <- BIC(model_1, model_2, model_3, model_4)
bic_values <- bic_values[order(bic_values$BIC), ]

show_table(bic_values)

	df	BIC
model_4	8	35195.25
model_2	9	35201.73
model_3	9	35202.49
model_1	10	35208.98

We can see that the model with the lowest BIC value is the one obtained by discarding both Anni.madre and Fumatrici, let’s check its performance.

summary(model_4)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Sesso, data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1129.14  -181.97   -16.26   160.95  2638.18 
## 
## Coefficients:
##                 Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)   -6708.0171   136.0715 -49.298 < 0.0000000000000002 ***
## N.gravidanze     12.7356     4.3385   2.935              0.00336 ** 
## Gestazione       32.3253     3.7969   8.514 < 0.0000000000000002 ***
## Lunghezza        10.2833     0.3009  34.177 < 0.0000000000000002 ***
## Cranio           10.5063     0.4263  24.648 < 0.0000000000000002 ***
## Tipo.partoNat    30.1601    12.1027   2.492              0.01277 *  
## SessoM           77.9171    11.1994   6.957     0.00000000000442 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.4 on 2491 degrees of freedom
## Multiple R-squared:  0.7277, Adjusted R-squared:  0.727 
## F-statistic:  1109 on 6 and 2491 DF,  p-value: < 0.00000000000000022

This model obtained a \(R^2 \simeq R^2_{adj} \simeq 0.73\): removing two variables, the performance of reduced model has not changed compared to the initial model. We can also observe that all the remaining predictors are statistically significant.

Let’s compare the two models using the ANOVA test, which evaluates whether the first complete model provides a significantly better fit than the reduced one.

anova(model_4, model_1)

## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Sesso
## Model 2: Peso ~ (Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso) - Ospedale
##   Res.Df       RSS Df Sum of Sq     F Pr(>F)
## 1   2491 187574428                          
## 2   2489 187430749  2    143679 0.954 0.3853

The ANOVA test shows a p-value greater than \(5\)%: this means that there is not enough evidence to say that model_1 is better than model_4, so the more complex model doesn’t significantly improve the fit. At this point, we prefer the reduced model to the initial model. Let’s continue searching a better solution.

current_model <- model_4
model_5 <- update(current_model, ~. -N.gravidanze)
model_6 <- update(current_model, ~. -Gestazione)
model_7 <- update(current_model, ~. -Lunghezza)
model_8 <- update(current_model, ~. -Cranio)
model_9 <- update(current_model, ~. -Tipo.parto)
model_10 <- update(current_model, ~. -Sesso)

bic_values <- BIC(current_model, model_5, model_6,
                  model_7, model_8, model_9, model_10)
bic_values <- bic_values[order(bic_values$BIC), ]

show_table(bic_values)

	df	BIC
model_9	7	35193.65
current_model	8	35195.25
model_5	7	35196.05
model_10	7	35235.50
model_6	7	35259.08
model_8	7	35732.60
model_7	7	36147.96

We can notice that the model with the lowest BIC value is the one obtained removing Tipo.parto, let’s check its performance.

summary(model_9)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = dataset)
## 
## 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 < 0.0000000000000002 ***
## N.gravidanze    12.4554     4.3416   2.869              0.00415 ** 
## Gestazione      32.3827     3.8008   8.520 < 0.0000000000000002 ***
## Lunghezza       10.2455     0.3008  34.059 < 0.0000000000000002 ***
## Cranio          10.5410     0.4265  24.717 < 0.0000000000000002 ***
## SessoM          77.9807    11.2111   6.956     0.00000000000447 ***
## ---
## 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: < 0.00000000000000022

The performance is approximately unchanged compared to the previous model, let’s compare them with ANOVA test.

anova(model_9, model_4)

## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Sesso
##   Res.Df       RSS Df Sum of Sq      F  Pr(>F)  
## 1   2492 188042054                              
## 2   2491 187574428  1    467626 6.2101 0.01277 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The ANOVA test yields a p-value smaller than \(5\)%, indicating that adding the predictor Tipo.parto significantly improves the fit; on the other hand, we have seen that model_9 has a lower BIC value and the \(R^2_{adj}\) changes very little, so we prefer to choose the more parsimonious model. Moreover, the variable Tipo.parto has a big limitation in real application contests, since it can only be determined at the time of birth, and therefore it cannot be used to predict newborns weight days or week before birth.

Let’s continue our backward procedure.

current_model <- model_9
model_11 <- update(current_model, ~. -N.gravidanze)
model_12 <- update(current_model, ~. -Gestazione)
model_13 <- update(current_model, ~. -Lunghezza)
model_14 <- update(current_model, ~. -Cranio)
model_15 <- update(current_model, ~. -Sesso)

bic_values <- BIC(current_model, model_11, model_12,
                  model_13, model_14, model_15)
bic_values <- bic_values[order(bic_values$BIC), ]

show_table(bic_values)

	df	BIC
current_model	7	35193.65
model_11	6	35194.06
model_15	6	35233.86
model_12	6	35257.55
model_14	6	35733.53
model_13	6	36140.54

We can observe that the current model is the one with the lowest BIC value, let’s do the ANOVA test to compare it with model_11, obtained removing N.gravidanze.

anova(model_11, model_9)

## Analysis of Variance Table
## 
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##   Res.Df       RSS Df Sum of Sq      F   Pr(>F)   
## 1   2493 188663107                                
## 2   2492 188042054  1    621053 8.2304 0.004154 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

With a p-value smaller than \(5\)%, the model with N.gravidanze significantly improves the fit compared to the model without it, so we choose model_9.

Now, let’s try to improve again this model by adding non-linear terms and interaction effects between predictors. More specifically, we will include quadratic terms for Gestazione and Cranio, since we previously observed a possible non-linear relationship between the target variable Peso and these two variables. Moreover, we will add in the model the following interaction terms:

Lunghezza - Cranio: these are anatomical variables that could interact to affect weight through a combined effect.
Gestazione - Cranio: these variables are both related to fetal development, so it is possible that the effect of gestation on weight depends on skull diameter or vice versa.
Gestazione - Lunghezza: also these variables are both related to fetal growth, so their effect on weight may be non-additive.

Please note that we will not consider taking into account more complex non-linear effects, as well as different interaction effects between variables. However, the goal is to obtain a model that can balance generalizability and simplicity.

Let’s start adding only the quadratic terms to model_9.

current_model <- model_9
model_16 <- update(model_9, ~. +I(Gestazione^2))
model_17 <- update(model_9, ~. +I(Cranio^2))
model_18 <- update(model_9, ~. +I(Gestazione^2) +I(Cranio^2))

bic_values <- BIC(current_model, model_16, model_17, model_18)
bic_values <- bic_values[order(bic_values$BIC), ]

show_table(bic_values)

	df	BIC
model_17	8	35166.63
model_18	9	35171.70
current_model	7	35193.65
model_16	8	35196.21

The model obtained adding only the quadratic term in Cranio achieved the lowest BIC value; instead, adding only the quadratic term in Gestazione seems to reduce the model performance. Let’s check the statistics of model_17.

summary(model_17)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Cranio^2), data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.6  -179.4   -14.8   163.4  2622.6 
## 
## Coefficients:
##                Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)    84.10118 1151.77280   0.073              0.94180    
## N.gravidanze   12.76356    4.31259   2.960              0.00311 ** 
## Gestazione     38.90540    3.93291   9.892 < 0.0000000000000002 ***
## Lunghezza      10.48745    0.30157  34.776 < 0.0000000000000002 ***
## Cranio        -31.79371    7.16973  -4.434      0.0000096316623 ***
## SessoM         73.10236   11.16590   6.547      0.0000000000711 ***
## I(Cranio^2)     0.06262    0.01059   5.915      0.0000000037748 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 272.8 on 2491 degrees of freedom
## Multiple R-squared:  0.7308, Adjusted R-squared:  0.7301 
## F-statistic:  1127 on 6 and 2491 DF,  p-value: < 0.00000000000000022

The \(R^2\) and the adjusted \(R^2\) slightly improved, let’s test if this improvement is statistically significant or not.

anova(model_9, model_17)

## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Cranio^2)
##   Res.Df       RSS Df Sum of Sq      F         Pr(>F)    
## 1   2492 188042054                                       
## 2   2491 185437522  1   2604532 34.987 0.000000003775 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

With a p-value near to zero, we can say that adding the quadratic term \(\text{Cranio} ^2\) significantly improves the fit of the previous model. This suggests us that the relationship between the target variable of weight and the predictor of skull diameter is not perfectly linear.

Now, we will introduce also the interaction terms described before.

current_model <- model_17
model_18 <- update(current_model, ~. +Lunghezza:Cranio)
model_19 <- update(current_model, ~. +Gestazione:Cranio)
model_20 <- update(current_model, ~. +Gestazione:Lunghezza)
model_21 <- update(current_model, ~.
                   +Lunghezza:Cranio +Gestazione:Cranio +Gestazione:Lunghezza)

bic_values <- BIC(current_model, model_18, model_19, model_20, model_21)
bic_values <- bic_values[order(bic_values$BIC), ]

show_table(bic_values)

	df	BIC
current_model	8	35166.63
model_19	9	35169.69
model_20	9	35172.28
model_18	9	35174.45
model_21	11	35181.60

We can observe that the current model, model_17, has the lowest BIC value, followed by model_19, the one obtained adding the interaction term between Gestazione and Cranio. Let’s compare them using the ANOVA test.

anova(model_17, model_19)

## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Cranio^2)
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Cranio^2) + Gestazione:Cranio
##   Res.Df       RSS Df Sum of Sq      F  Pr(>F)  
## 1   2491 185437522                              
## 2   2490 185084578  1    352944 4.7483 0.02942 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Since p-value \(<0.05\), adding the interaction term Gestazione-Cranio significantly improves the fit. However, since the model without this term achieved a lower BIC value, we decide to choose that model. This choice is driven by the goal of reducing complexity and improving the interpretability.

2.4 Model Quality Analysis

We have thus obtained our optimal multiple linear regression model. The coefficients of the optimal model, the model_17, are the following.

# Final multiple linear regression model
final_model <- model_17

# Summary of model's coefficients
coef_summary <- data.frame(
  Term = names(final_model$coefficients),
  Coefficient = round(final_model$coefficients, 1),
  row.names = NULL
)

show_table(coef_summary)

Term	Coefficient
(Intercept)	84.1
N.gravidanze	12.8
Gestazione	38.9
Lunghezza	10.5
Cranio	-31.8
SessoM	73.1
I(Cranio^2)	0.1

Therefore, the linear regression equation for predicting neonatal weight is constructed as follows.

\[ \widehat{Peso} \, [\text{g}] = 84.1 + 12.8 \cdot N.gravidanze + 38.9 \cdot Gestazione + 10.5 \cdot Lunghezza - 31.8 \cdot Cranio + 0.1 \cdot Cranio^2 + 73.1 \cdot I(Sesso = M) \] where \(\widehat{Peso}\) is the predicted value of the target variable, measured in grams, and \(I(sesso = M)\) is the indicator function of male sex.

summary(final_model)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Cranio^2), data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.6  -179.4   -14.8   163.4  2622.6 
## 
## Coefficients:
##                Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)    84.10118 1151.77280   0.073              0.94180    
## N.gravidanze   12.76356    4.31259   2.960              0.00311 ** 
## Gestazione     38.90540    3.93291   9.892 < 0.0000000000000002 ***
## Lunghezza      10.48745    0.30157  34.776 < 0.0000000000000002 ***
## Cranio        -31.79371    7.16973  -4.434      0.0000096316623 ***
## SessoM         73.10236   11.16590   6.547      0.0000000000711 ***
## I(Cranio^2)     0.06262    0.01059   5.915      0.0000000037748 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 272.8 on 2491 degrees of freedom
## Multiple R-squared:  0.7308, Adjusted R-squared:  0.7301 
## F-statistic:  1127 on 6 and 2491 DF,  p-value: < 0.00000000000000022

The final model explains about \(73\)% of the variability in newborns weight, and it has a residual standard error of about \(273\) grams. Let’s evaluate the RMSE, which represents the typical mean error of the model’s predictions.

RMSE <- sqrt(mean(residuals(final_model)^2))
round(RMSE, 2)

## [1] 272.46

Now let’s check if there is the problem of multicollinearity among the predictors by evaluating the VIF of the variables.

# Compute VIF values
vif_summary <- data.frame(
  Predictor = names(vif(final_model)),
  VIF = vif(final_model),
  row.names = NULL
)

show_table(vif_summary)

Predictor	VIF
N.gravidanze	1.02
Gestazione	1.81
Lunghezza	2.11
Cranio	465.42
Sesso	1.05
I(Cranio^2)	452.90

From the VIF summary above we can observe that Cranio and its quadratic term show very high VIF values, indicating a strong multicollinearity between these variables: this was expected, since one term is the square of the other. We can continue to use this model, taking into account that it may struggle to distinguish between these two terms, so the individual coefficients may be unstable. For the other predictors, VIF values are smaller than the usual threshold of \(5\), suggesting that there are no multicollinearity issues for them.

Now let’s analyze the residuals of the model.

# Residuals plots
par(mfrow = c(2, 2))
plot(final_model)

Looking at the graph “Residuals vs Fitted”, which shows the relationship between the model’s predicted values and their residuals, the residuals seem randomly distributed around the mean of zero and don’t exhibit a strong systematic pattern.

The “Q-Q Residuals” plot compares the quantiles of the residuals with the quantiles of a normal distribution; we can observe that the points lie on the diagonal line in the central part of the distribution, while they deviate in the tails, especially for the right tail, where they lie above the diagonal. This suggests a deviation from normality for the highest positive residuals.

Studying the “Scale-Location” graph, which shows the relationship between the fitted values and the square root of the absolute standardized residuals, we can see that the dispersion of the points appears to increase significantly in the right half of the graph, suggesting possible heteroscedasticity.

From the plot “Residuals vs Leverage” we can observe that there is only one observations with a Cook’s distance greater than \(0.5\), which we can consider a warning threshold for outliers or leverage points.

That point in particular, the record 1549, is a very high positive residual and deviates considerably from the other observations; later we will examine whether this is a significantly anomalous value.

Now let’s verify our deduction using appropriate statistical tests. Let’s start by checking if the residuals are normally distributed by using the Shapiro-Wilk test with a significance level of \(5\)%.

# Model's residuals
model_residuals <- residuals(final_model)

# Shapiro-Wilk test for normality of residuals
shapiro.test(model_residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model_residuals
## W = 0.97436, p-value < 0.00000000000000022

The Shapiro-Wilk test yields a p-value \(<5\)%, so we have to reject the null hypothesis of normality of residuals. Let’s visualize the residuals density plot and compute the descriptive statistics of residuals distribution.

# Density plot of residuals
ggplot() +
  geom_density(
    aes(x = model_residuals),
    col = "black", 
    fill = "skyblue") +
  labs(
    title = "Density plot of residuals",
    x = "Residuals (g)",
    y = "Density"
  ) +
  theme_light() +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5, face = "bold", color = "red3"),
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title.x = element_text(size = 13, margin = margin(t = 10)),
    axis.title.y = element_text(size = 13, margin = margin(r = 10))
  )

# Residuals' descriptive statistics
residuals_df <- data.frame(Model_Residuals = model_residuals)
show_table(quantitative_statistics(residuals_df, "Model_Residuals"))

	Model_Residuals
Mean value	0.00
Minimum	-1138.58
First quartile (Q1)	-179.37
Median value	-14.80
Third quartile (Q3)	163.37
Maximum	2622.58
Variance	74264.13
Standard deviation	272.51
Range	3761.16
Interquartile range	342.74
Coefficient of variation	-16365868172109000.00
Fisher’s Skewness	0.67
Excess Kurtosis	4.25
Number of outliers	43.00

We observe that the residuals distribution is slightly positively skewed and has a high excess kurtosis, exhibiting heavier tails. The median is slightly negative, suggesting that the model tends to overestimate the newborns weight.

Let’s examine the heteroscedasticity of residuals by using the Breusch-Pagan test.

# Breusch-Pagan test for heteroscedasticity of residuals
bptest(final_model)

## 
##  studentized Breusch-Pagan test
## 
## data:  final_model
## BP = 90.755, df = 6, p-value < 0.00000000000000022

With a p-value near to zero, the test suggests the presence of heteroscedasticity: the variance of residuals is not constant over the range of values. Now let’s examine the presence of autocorrelation among the residuals using the Durbin-Watson test.

# Durbin-Watson test for autocorrelation of residuals
dwtest(final_model)

## 
##  Durbin-Watson test
## 
## data:  final_model
## DW = 1.9523, p-value = 0.1165
## alternative hypothesis: true autocorrelation is greater than 0

With a p-value \(>5\)%, we don’t reject the null hypothesis of absence of autocorrelation among the residuals. The next step is to verify the presence of leverage and outliers, and to test if they have a significant impact on the regression coefficients.

# Leverage values
lev <- hatvalues(final_model)

p <- sum(lev)   # p = number of model's parameters
n <- nrow(dataset)   # n = sample size

# Threshold for leverages
threshold_lev = 2 * p/n

lev_df <- data.frame(
  records = 1:length(lev),
  leverage = lev
)

ggplot(lev_df, aes(
  x = records, 
  y = leverage)) +
  geom_point() +
  geom_hline(yintercept = threshold_lev, color = "red3") +
  labs(
    title = "Scatterplot of leverages",
    x = "Index",
    y = "Leverage"
  ) +
  theme_light() +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5, face = "bold", color = "red3"),
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title.x = element_text(size = 13, margin = margin(t = 10)),
    axis.title.y = element_text(size = 13, margin = margin(r = 10))
  )

# Number of high-leverage points
count_lev <- sum(lev > threshold_lev)
prc_lev <- round((count_lev/n) * 100, 2)
cat(paste0("Number of high leverage observations: ", count_lev, " (", prc_lev, "% of sample)"))

## Number of high leverage observations: 151 (6.04% of sample)

We can see that there is a large number of high-leverage points, which are observations whose leverage value exceeds the \(2 \cdot \frac{p}{n}\) threshold (where \(p\) is the number of parameters of the model, and \(n\) is the sample size).

Now let’s find out if there are some outliers observations.

# Studentized residuals
student_res <- rstudent(final_model)

out_df <- data.frame(
  records = 1:length(student_res),
  studentized_residuals = student_res
)

ggplot(out_df, aes(
  x = records, 
  y = studentized_residuals)) +
  geom_point() +
  geom_hline(yintercept = c(-2, 2), color = "red3") +
  labs(
    title = "Scatterplot of Studentized residuals",
    x = "Index",
    y = "Studentized Residual"
  ) +
  theme_light() +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5, face = "bold", color = "red3"),
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title.x = element_text(size = 13, margin = margin(t = 10)),
    axis.title.y = element_text(size = 13, margin = margin(r = 10))
  )

# Number of outliers points
count_out <- sum(abs(student_res) > 2)
prc_out <- round((count_out/n) * 100, 2)
cat(paste0("Number of outlier observations: ", count_out, " (", prc_out, "% of sample)"))

## Number of outlier observations: 102 (4.08% of sample)

We can observe that there are many outliers, which are observations with absolute studentized residual greater than \(2\), which we can consider a warning threshold.

Now let’s compute the Cook’s distance for the observations, so we can determine whether there are leverage points or outliers that could have a significant impact on the regression estimates.

# Cook's distances
cook_dist <- cooks.distance(final_model)

cook_df <- data.frame(
  records = 1:length(cook_dist),
  cook_distance = cook_dist
)

ggplot(cook_df, aes(
  x = records, 
  y = cook_distance)) +
  geom_point() +
  geom_hline(yintercept = 0.5, color = "red3") +
  labs(
    title = "Scatterplot of Cook's distances",
    x = "Index",
    y = "Cook's distance"
  ) +
  theme_light() +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5, face = "bold", color = "red3"),
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title.x = element_text(size = 13, margin = margin(t = 10)),
    axis.title.y = element_text(size = 13, margin = margin(r = 10))
  )

# High Cook's distance points
high_cook_points <- data.frame(
  Observation = which(cook_dist > 0.5),
  Cook_distance = cook_dist[cook_dist > 0.5],
  row.names = NULL
)

show_table(high_cook_points)

Observation	Cook_distance
1549	0.71

We can see there is only one point that has a Cook’s distance that exceeds the threshold of \(0.5\)%, the observation 1549; that record is precisely the anomalous observation we had noted while studying the model’s residuals. Let’s display the values of leverage and studentized residual for the anomalous record.

anomalous_data_index <- which(cook_dist > 0.5)

cat(paste0(
  "Leverage value of anomalous observation: ", 
  round(lev[anomalous_data_index], 3), 
  " (Threshold = ",
  round(threshold_lev, 3), 
  ") \n\n")
)

cat(paste0(
  "Absolute Studentized residual of anomalous observation: ", 
  round(abs(student_res[anomalous_data_index]), 3),
  " (Threshold = 0.5",
  ") \n")
)

## Leverage value of anomalous observation: 0.049 (Threshold = 0.006) 
## 
## Absolute Studentized residual of anomalous observation: 10.052 (Threshold = 0.5)

We can notice that the anomalous observation 1549 is both a leverage and an outlier point. In fact, it has a leverage value greater than the warning threshold computed brefore, and it also has an absolute Studentized residual that exceeds the warning threshold of \(0.5\). Let’s analyze the impact of this observation on the performance and coefficients of the multiple linear regression model.

# Model without the anomalous observation
model_without_anomaly <- lm(
  Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
    Sesso + I(Cranio^2), data = dataset[-anomalous_data_index, ]
)

coef <- final_model$coefficients
coef_no_anomaly <- model_without_anomaly$coefficients

coef_comparison_df <- data.frame(
  Model_with_anomaly = coef,
  Model_without_anomaly = coef_no_anomaly
)

# Percentage changes in coefficients
coef_comparison_df$Percentage_change <- 100 * (
  (coef_comparison_df$Model_without_anomaly -
   coef_comparison_df$Model_with_anomaly) /
    coef_comparison_df$Model_with_anomaly)

show_table(coef_comparison_df)

	Model_with_anomaly	Model_without_anomaly	Percentage_change
(Intercept)	84.10	-36.27	-143.12
N.gravidanze	12.76	13.44	5.33
Gestazione	38.91	36.06	-7.31
Lunghezza	10.49	11.12	6.07
Cranio	-31.79	-31.67	-0.38
SessoM	73.10	73.34	0.33
I(Cranio^2)	0.06	0.06	-1.75

cat("R^2 of complete model =", round(summary(final_model)$r.squared, 3), "\n\n")
cat("R^2 of model without anomalous record =", 
    round(summary(model_without_anomaly)$r.squared, 3))

## R^2 of complete model = 0.731 
## 
## R^2 of model without anomalous record = 0.741

We can notice that, excluding the anomalous observation, the multiple linear regression model shows a \(1\)% improvement in the \(R^2\) value. Moreover, looking at the percentage variations in the predictors’ coefficients, we observe that there are modest changes for the variables N.gravidanze, Gestazione and Lunghezza. We can say that the observation 1549 is a significantly influential record for regression estimates.

Let’s observe the values of variables for this particular observation.

show_table(dataset[anomalous_data_index, ])

	Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
1549	35	1	Non-smoker	38	4370	315	374	Nat	osp3	F

This observation seems to assume plausible values for the variables; we can notice that, relative to the very low length of the female newborn (\(31.5\) cm, much lower than the first quartile of Lunghezza, which is \(48\) cm), the observed weight (\(4.37\) kg) is significantly greater than the median weight of newborns (\(3.30\) kg). Despite that, there are not sufficient reasons to remove this record from the dataset, since it represents a real and valid observation.

3 Predictions and Results

Let’s use the trained model to make practical predictions. For example, we want to estimate the weight of a female newborn, assuming a mother in her third pregnancy who will deliver at \(39\) weeks. Our final model requires the predictors N.gravidanze, Gestazione, Sesso - whose values are know in this case - and also Lunghezza and Cranio, whose values are unknown. To infer the weight in the example, we will use the mean values for Lunghezza and Cranio.

# New data to estimate
data_to_predict <- data.frame(
  N.gravidanze = 3,
  Gestazione = 39,
  Sesso = "F",
  Lunghezza = mean(Lunghezza),
  Cranio = mean(Cranio)
)

# Model's prediction
predicted_weight <- predict(final_model, newdata = data_to_predict)

cat(paste0(
  "Predicted weight = ",
  round(predicted_weight, 2),
  " grams"
))

## Predicted weight = 3257.47 grams

Let’s make another example: we aim to predict the weight of a male newborn, assuming a mother in her first pregnancy who will deliver at \(28\) weeks.

# New data to estimate
data_to_predict <- data.frame(
  N.gravidanze = 1,
  Gestazione = 28,
  Sesso = "M",
  Lunghezza = mean(Lunghezza),
  Cranio = mean(Cranio)
)

# Model's prediction
predicted_weight <- predict(final_model, newdata = data_to_predict)

cat(paste0(
  "Predicted weight = ",
  round(predicted_weight, 2),
  " grams"
))

## Predicted weight = 2877.08 grams

4 Data Visualizations

In this chapter, we will create some graphs to show the most significant relationships between the variables.

For example, we could visualize the impact of the number of weeks of gestation on the predicted weight, grouped by smoking class.

# Scatterplot of Gestazione vs Peso by Fumatrici
ggplot(data=dataset) +
  geom_point(aes(
    x = Gestazione,
    y = Peso,
    color = Fumatrici
  )) + 
geom_smooth(aes(
  x = Gestazione,
  y = Peso,
  color = Fumatrici),
  se = FALSE,
  method= "lm") +
  labs(
    title = "Gestation vs Weight by Smoking",
    x = "Gestation (weeks)",
    y = "Weight (g)"
  ) +
  theme_light() +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5, face = "bold", color = "red3"),
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title.x = element_text(size = 13, margin = margin(t = 10)),
    axis.title.y = element_text(size = 13, margin = margin(r = 10))
  )

Considering only Gestazione and Fumatrici as predictors, we can observe that weight increases as the number of weeks of gestation increases. We can also see a small difference in the linear regression’s slopes between smoker and non-smoker mothers. However, taking into consideration all the variables, our multiple linear regression model excluded smoking as a significant predictor.

Now let’s visualize the relationship between the number of pregnancies and the weight, grouping by newborns sex.

# Scatterplot of N.gravidanze vs Peso by Sesso
ggplot(data = dataset) +
  geom_point(aes(
    x = N.gravidanze,
    y = Peso,
    color = Sesso
  )) + 
  geom_smooth(aes(
    x = N.gravidanze,
    y = Peso,
    color = Sesso),
    se = FALSE,
    method = "lm"
  ) + 
  labs(
    x = "Number of pregnancies",
    y = "Weight (g)",
    title = "Pregnancies vs Weight by Sex"
  ) +
  theme_light() +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5, face = "bold", color = "red3"),
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title.x = element_text(size = 13, margin = margin(t = 10)),
    axis.title.y = element_text(size = 13, margin = margin(r = 10))
  )

Considering only N.gravidanze and Sesso, we can observe that the linear fits results in approximately horizontal lines for both male and female newborns, which means that the number of pregnancies and sex, when considered on their own, are not useful to predict the weight.

Let’s display a scatterplot showing the relationship between newborns length and weight, grouping by sex.

# Scatterplot of Lunghezza vs Peso by Sesso
ggplot(data = dataset) +
  geom_point(aes(
    x = Lunghezza,
    y = Peso,
    color = Sesso
  )) + 
  geom_smooth(aes(
    x = Lunghezza,
    y = Peso,
    color = Sesso),
    se = FALSE,
    method = "lm"
  ) + 
  labs(
    x = "Length (mm)",
    y = "Weight (g)",
    title = "Length vs Weight by Sex"
  ) +
  theme_light() +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5, face = "bold", color = "red3"),
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title.x = element_text(size = 13, margin = margin(t = 10)),
    axis.title.y = element_text(size = 13, margin = margin(r = 10))
  )

From the graph above, we can notice an upward trend: weight increases as length increases, and this growth trend is very similar among male and female newborns.

Let’s visualize the relationship between newborns skull diameter and weight, grouping by sex.

# Scatterplot of Cranio vs Peso by Sesso
ggplot(data = dataset) +
  geom_point(aes(
    x = Cranio,
    y = Peso,
    color = Sesso
  )) + 
  geom_smooth(aes(
    x = Cranio,
    y = Peso,
    color = Sesso),
    se = FALSE,
    method = "lm"
  ) + 
  labs(
    x = "Skull diameter (mm)",
    y = "Weight (g)",
    title = "Skull diameter vs Weight by Sex"
  ) +
  theme_light() +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5, face = "bold", color = "red3"),
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title.x = element_text(size = 13, margin = margin(t = 10)),
    axis.title.y = element_text(size = 13, margin = margin(r = 10))
  )

Once again, we can observe an upward trend: weight increases as skull diameter increases, with a similar growth rate and a systematic difference between male and female newborns. This is consistent with the results obtained by the multiple linear regression model we trained: under equal conditions, male newborns statistically weigh more than female newborns.

Finally, let’s visualize the impact of mother’s age on the predicted weight, grouped by sex.

# Scatterplot of Anni.madre vs Peso by Sesso
ggplot(data = dataset) +
  geom_point(aes(
    x = Anni.madre,
    y = Peso,
    color = Sesso
  )) + 
  geom_smooth(aes(
    x = Anni.madre,
    y = Peso,
    color = Sesso),
    se = FALSE,
    method = "lm"
  ) + 
  labs(
    x = "Mother's age (years)",
    y = "Weight (g)",
    title = "Mother's age vs Weight by Sex"
  ) +
  theme_light() +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5, face = "bold", color = "red3"),
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title.x = element_text(size = 13, margin = margin(t = 10)),
    axis.title.y = element_text(size = 13, margin = margin(r = 10))
  )

We notice that linear fits have a slope close to zero, suggesting that there isn’t a clear association between mother’s age and neonatal weight for both sexes.

5 Conclusions

We built a multiple linear regression model to predict neonatal weight based on the variables N.gravidanze, Gestazione, Lunghezza, Cranio and Sesso.

Studying its performance and analyzing the residuals, we observed that the model explains about \(73\)% of variability of weight in the training dataset, but its residuals don’t follow a normal distribution and exhibit heteroscedasticity. Despite that, we can still consider it a fairly good model for predicting newborn weight.

Moreover, we noticed that the variable Fumatrici is not statistically significant in the model, suggesting that there are no appreciable differences in birth weight between babies born to smoking and non-smoking mothers.

Finally, also Anni.madre is not statistically significant in the model, so mother’s age does not seem to be a good predictor of neonatal weight.

Statistical Model for Predicting Neonatal Weight

Leonardo Arduini

2026-08-14