# Carica le librerie necessarie
library(tidyverse) # ggplot2 è incluso in tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(caret) # Per machine learning e analisi statistica
## Caricamento del pacchetto richiesto: lattice
##
## Caricamento pacchetto: 'caret'
##
## Il seguente oggetto è mascherato da 'package:purrr':
##
## lift
library(ggplot2) # Per la creazione di grafici
library(car) # Per l'analisi di regressione
## Caricamento del pacchetto richiesto: carData
##
## Caricamento pacchetto: 'car'
##
## Il seguente oggetto è mascherato da 'package:dplyr':
##
## recode
##
## Il seguente oggetto è mascherato da 'package:purrr':
##
## some
library(lmtest) # Per il test di ipotesi nei modelli
## Caricamento del pacchetto richiesto: zoo
##
## Caricamento pacchetto: 'zoo'
##
## I seguenti oggetti sono mascherati da 'package:base':
##
## as.Date, as.Date.numeric
library(knitr) # Per l'integrazione di risultati nel documento
library(conflicted)
conflict_prefer("filter", "dplyr")
## [conflicted] Will prefer dplyr::filter over any other package.
conflict_prefer("lag", "dplyr")
## [conflicted] Will prefer dplyr::lag over any other package.
conflict_prefer("recode", "dplyr")
## [conflicted] Will prefer dplyr::recode over any other package.
# Carica il dataset in formato CSV
dataset <- read.csv("C:/Users/carbaren/Desktop/R/neonati.csv")
# Verifica i nomi delle colonne del dataset
colnames(dataset)
## [1] "Anni.madre" "N.gravidanze" "Fumatrici" "Gestazione" "Peso"
## [6] "Lunghezza" "Cranio" "Tipo.parto" "Ospedale" "Sesso"
# Mostra le prime righe del dataset
head(dataset)
## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## 1 26 0 0 42 3380 490 325 Nat
## 2 21 2 0 39 3150 490 345 Nat
## 3 34 3 0 38 3640 500 375 Nat
## 4 28 1 0 41 3690 515 365 Nat
## 5 20 0 0 38 3700 480 335 Nat
## 6 32 0 0 40 3200 495 340 Nat
## Ospedale Sesso
## 1 osp3 M
## 2 osp1 F
## 3 osp2 M
## 4 osp2 M
## 5 osp3 F
## 6 osp2 F
# Statistiche descrittive più complete per avere un'idea generale del dataset
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
## Min. : 830 Min. :310.0 Min. :235 Length:2500
## 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Class :character
## Median :3300 Median :500.0 Median :340 Mode :character
## Mean :3284 Mean :494.7 Mean :340
## 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
## Max. :4930 Max. :565.0 Max. :390
## Ospedale Sesso
## Length:2500 Length:2500
## Class :character Class :character
## Mode :character Mode :character
##
##
##
sapply(dataset, function(x) c(mean = mean(x, na.rm = TRUE),
sd = sd(x, na.rm = TRUE),
median = median(x, na.rm = TRUE)))
## Warning in mean.default(x, na.rm = TRUE): l'argomento non è numerico o logico:
## si restituisce NA
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm =
## na.rm): NA introdotti per coercizione
## Warning in mean.default(sort(x, partial = half + 0L:1L)[half + 0L:1L]):
## l'argomento non è numerico o logico: si restituisce NA
## Warning in mean.default(x, na.rm = TRUE): l'argomento non è numerico o logico:
## si restituisce NA
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm =
## na.rm): NA introdotti per coercizione
## Warning in mean.default(sort(x, partial = half + 0L:1L)[half + 0L:1L]):
## l'argomento non è numerico o logico: si restituisce NA
## Warning in mean.default(x, na.rm = TRUE): l'argomento non è numerico o logico:
## si restituisce NA
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm =
## na.rm): NA introdotti per coercizione
## Warning in mean.default(sort(x, partial = half + 0L:1L)[half + 0L:1L]):
## l'argomento non è numerico o logico: si restituisce NA
## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza
## mean 28.164000 0.981200 0.0416000 38.980400 3284.0808 494.69200
## sd 5.273578 1.280587 0.1997133 1.868639 525.0387 26.31864
## median 28.000000 1.000000 0.0000000 39.000000 3300.0000 500.00000
## Cranio Tipo.parto Ospedale Sesso
## mean 340.02920 NA NA NA
## sd 16.42533 NA NA NA
## median 340.00000 NA NA NA
# Verifica la presenza di valori mancanti nel dataset
sapply(dataset, function(x) sum(is.na(x)))
## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza
## 0 0 0 0 0 0
## Cranio Tipo.parto Ospedale Sesso
## 0 0 0 0
# Controllo per età madre improbabili
table(dataset$Anni.madre)
##
## 0 1 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
## 1 1 1 2 6 13 18 24 45 66 74 100 115 131 180 184 197 172 174 200
## 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
## 147 159 110 96 66 64 41 38 27 19 13 8 2 4 1 1
dataset$Anni.madre[dataset$Anni.madre < 15 | dataset$Anni.madre > 60] <- NA
# Imputazione dei valori mancanti con la media per le variabili numeriche
dataset_imputed <- dataset %>%
mutate(across(where(is.numeric), ~ifelse(is.na(.), mean(., na.rm = TRUE), .)))
# Creazione della distribuzione del peso del neonato
ggplot(dataset_imputed, aes(x = Peso)) +
geom_histogram(binwidth = 50, fill = "blue", color = "black") +
theme_minimal() +
labs(title = "Distribuzione del peso del neonato")
# Distribuzione della durata della gravidanza
ggplot(dataset_imputed, aes(x = Gestazione)) +
geom_histogram(binwidth = 1, fill = "green", color = "black") +
theme_minimal() +
labs(title = "Distribuzione della durata della gravidanza")
# Creiamo una tabella delle frequenze tra ospedale e tipo di parto
hospital_parto_cesareo <- table(dataset_imputed$Ospedale, dataset_imputed$Tipo.parto)
hospital_parto_cesareo
##
## Ces Nat
## osp1 242 574
## osp2 254 595
## osp3 232 603
# Eseguiamo un test chi-quadrato per verificare se la distribuzione del tipo di parto dipende dall'ospedale
chisq.test(hospital_parto_cesareo)
##
## Pearson's Chi-squared test
##
## data: hospital_parto_cesareo
## X-squared = 1.0972, df = 2, p-value = 0.5778
# Test t per verificare se il peso è significativamente diverso dalla media della popolazione (3200g)
t.test(dataset_imputed$Peso, mu = 3200)
##
## One Sample t-test
##
## data: dataset_imputed$Peso
## t = 8.0071, df = 2499, p-value = 1.782e-15
## alternative hypothesis: true mean is not equal to 3200
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081
# Test t per verificare se la lunghezza è significativamente diversa dalla media della popolazione (50 cm)
t.test(dataset_imputed$Lunghezza, mu = 50)
##
## One Sample t-test
##
## data: dataset_imputed$Lunghezza
## t = 844.82, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 50
## 95 percent confidence interval:
## 493.6598 495.7242
## sample estimates:
## mean of x
## 494.692
# Matrice di correlazione per le variabili predittive
cor(dataset_imputed %>% select(Peso, Anni.madre, N.gravidanze, Gestazione, Lunghezza, Cranio))
## Peso Anni.madre N.gravidanze Gestazione Lunghezza
## Peso 1.00000000 -0.02403702 0.00240730 0.5917687 0.79603676
## Anni.madre -0.02403702 1.00000000 0.38335003 -0.1361348 -0.06634028
## N.gravidanze 0.00240730 0.38335003 1.00000000 -0.1014919 -0.06040371
## Gestazione 0.59176872 -0.13613476 -0.10149194 1.0000000 0.61892045
## Lunghezza 0.79603676 -0.06634028 -0.06040371 0.6189204 1.00000000
## Cranio 0.70480151 0.01786895 0.03900707 0.4608289 0.60334097
## Cranio
## Peso 0.70480151
## Anni.madre 0.01786895
## N.gravidanze 0.03900707
## Gestazione 0.46082894
## Lunghezza 0.60334097
## Cranio 1.00000000
# Creazione di un modello di regressione lineare con diverse variabili predittive
modello <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = dataset_imputed)
# Riepilogo del modello per ottenere informazioni sui coefficienti e significatività
summary(modello)
##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = dataset_imputed)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1123.41 -181.25 -14.88 160.77 2612.63
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6735.4798 141.4534 -47.616 < 2e-16 ***
## Anni.madre 0.8074 1.1521 0.701 0.4835
## N.gravidanze 11.4050 4.6661 2.444 0.0146 *
## Fumatrici -30.1408 27.5397 -1.094 0.2739
## Gestazione 32.5308 3.8182 8.520 < 2e-16 ***
## Lunghezza 10.2956 0.3007 34.239 < 2e-16 ***
## Cranio 10.4713 0.4261 24.573 < 2e-16 ***
## Tipo.partoNat 29.5056 12.0847 2.442 0.0147 *
## Ospedaleosp2 -11.2377 13.4391 -0.836 0.4031
## Ospedaleosp3 28.1091 13.4967 2.083 0.0374 *
## SessoM 77.5355 11.1781 6.936 5.11e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 669.1 on 10 and 2489 DF, p-value: < 2.2e-16
# Selezione del modello passo per passo (stepwise) escludendo variabili irrilevanti come Tipo.parto e Ospedale
modello_selezionato <- step(modello, direction = "both")
## Start: AIC=28075.39
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 36851 186809099 28074
## - Fumatrici 1 89883 186862131 28075
## <none> 186772248 28075
## - Tipo.parto 1 447324 187219571 28079
## - N.gravidanze 1 448298 187220546 28079
## - Ospedale 2 687324 187459572 28081
## - Sesso 1 3610391 190382639 28121
## - Gestazione 1 5446911 192219158 28145
## - Cranio 1 45311010 232083258 28616
## - Lunghezza 1 87967488 274739735 29038
##
## Step: AIC=28073.88
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 90897 186899996 28073
## <none> 186809099 28074
## + Anni.madre 1 36851 186772248 28075
## - Tipo.parto 1 448222 187257321 28078
## - Ospedale 2 692738 187501837 28079
## - N.gravidanze 1 633756 187442855 28080
## - Sesso 1 3618736 190427835 28120
## - Gestazione 1 5412879 192221978 28143
## - Cranio 1 45588236 232397335 28618
## - Lunghezza 1 87950050 274759149 29036
##
## Step: AIC=28073.1
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 186899996 28073
## + Fumatrici 1 90897 186809099 28074
## + Anni.madre 1 37866 186862131 28075
## - Tipo.parto 1 440684 187340680 28077
## - Ospedale 2 701680 187601677 28079
## - N.gravidanze 1 610840 187510837 28079
## - Sesso 1 3602797 190502794 28119
## - Gestazione 1 5346781 192246777 28142
## - Cranio 1 45632149 232532146 28617
## - Lunghezza 1 88355030 275255027 29039
# Riepilogo del modello ottimale
summary(modello_selezionato)
##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = dataset_imputed)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1113.18 -181.16 -16.58 161.01 2620.19
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.4293 135.9438 -49.340 < 2e-16 ***
## N.gravidanze 12.3619 4.3325 2.853 0.00436 **
## Gestazione 31.9909 3.7896 8.442 < 2e-16 ***
## Lunghezza 10.3086 0.3004 34.316 < 2e-16 ***
## Cranio 10.4922 0.4254 24.661 < 2e-16 ***
## Tipo.partoNat 29.2803 12.0817 2.424 0.01544 *
## Ospedaleosp2 -11.0227 13.4363 -0.820 0.41209
## Ospedaleosp3 28.6408 13.4886 2.123 0.03382 *
## SessoM 77.4412 11.1756 6.930 5.36e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared: 0.7287, Adjusted R-squared: 0.7278
## F-statistic: 836.3 on 8 and 2491 DF, p-value: < 2.2e-16
# Calcoliamo i criteri AIC e BIC per il modello selezionato
AIC(modello_selezionato)
## [1] 35169.79
BIC(modello_selezionato)
## [1] 35228.03
# Calcoliamo il valore di R² per valutare la bontà di adattamento del modello
rsq <- summary(modello_selezionato)$r.squared
cat("R²:", rsq)
## R²: 0.7286934
# Calcoliamo il RMSE (Root Mean Squared Error) per valutare la qualità del modello
pred <- predict(modello_selezionato)
# Calcoliamo i residui
residui <- dataset_imputed$Peso - pred
# Calcoliamo il RMSE
rmse <- sqrt(mean(residui^2))
# Visualizza il valore del RMSE
cat("RMSE:", rmse)
## RMSE: 273.4227
# Diagnosi dei residui per verificare eventuali pattern non catturati dal modello
plot(residui, main = "Residui del modello", xlab = "Indice dei dati", ylab = "Residui")
abline(h = 0, col = "red")
# Distribuzione dei residui
hist(residui, main = "Distribuzione dei residui", xlab = "Residui", col = "gray", border = "black")
# Calcolo del VIF (Variance Inflation Factor) per verificare la multicollinearità
vif(modello_selezionato)
## GVIF Df GVIF^(1/(2*Df))
## N.gravidanze 1.025244 1 1.012544
## Gestazione 1.670237 1 1.292376
## Lunghezza 2.081921 1 1.442886
## Cranio 1.626507 1 1.275346
## Tipo.parto 1.003872 1 1.001934
## Ospedale 1.003038 2 1.000759
## Sesso 1.040337 1 1.019969
# Creazione di un nuovo dato per la previsione
nuovo_dato <- data.frame(
Anni.madre = 28,
N.gravidanze = 3,
Fumatrici = 0,
Gestazione = 39,
Lunghezza = 50,
Cranio = 34,
Tipo.parto = "Nat", # Esempio: "Nat" per parto naturale
Ospedale = "osp1", # Esempio: "osp1" per ospedale 1
Sesso = "M" # Sesso: "M" per maschio
)
# Assicuriamoci che le variabili categoriche siano trattate correttamente come fattori
nuovo_dato$Tipo.parto <- factor(nuovo_dato$Tipo.parto, levels = levels(dataset_imputed$Tipo.parto))
nuovo_dato$Ospedale <- factor(nuovo_dato$Ospedale, levels = levels(dataset_imputed$Ospedale))
nuovo_dato$Sesso <- factor(nuovo_dato$Sesso, levels = levels(dataset_imputed$Sesso))
# Previsione del peso del neonato per il nuovo caso
previsione <- predict(modello_selezionato, nuovo_dato)
# Stampa il risultato
cat("Il peso previsto del neonato è:", previsione, "grammi")
## Il peso previsto del neonato è: NA grammi
# Creiamo un grafico per visualizzare l'effetto della durata della gravidanza sul peso
ggplot(dataset_imputed, aes(x = Gestazione, y = Peso)) +
geom_point() +
geom_smooth(method = "lm", color = "blue") +
labs(title = "Impatto della durata della gravidanza sul peso del neonato",
x = "Durata della gravidanza (settimane)",
y = "Peso del neonato (grammi)")
## `geom_smooth()` using formula = 'y ~ x'
# Creiamo un boxplot per visualizzare l'effetto del fumo sul peso
ggplot(dataset_imputed, aes(x = Fumatrici, y = Peso, fill = factor(Fumatrici))) +
geom_boxplot() +
labs(title = "Differenze nel peso tra fumatrice e non fumatrice")
