Project 3
Francesco Filoseta
2025-07-10
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = "C:/Users/Francesco/Desktop/ProfessionAI")
library(moments)
library(dplyr)##
## Caricamento pacchetto: 'dplyr'## I seguenti oggetti sono mascherati da 'package:stats':
##
## filter, lag## I seguenti oggetti sono mascherati da 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(patchwork)
library(MASS)##
## Caricamento pacchetto: 'MASS'## Il seguente oggetto è mascherato da 'package:patchwork':
##
## area## Il seguente oggetto è mascherato da 'package:dplyr':
##
## selectlibrary(car)## Caricamento del pacchetto richiesto: carData##
## Caricamento pacchetto: 'car'## Il seguente oggetto è mascherato da 'package:dplyr':
##
## recode## Importing DataFrame
getwd()## [1] "C:/Users/Francesco/Desktop/ProfessionAI"data = read.csv("neonatiCSV.csv", sep = ",") # Importing
attach(data)Descriptive statistics
head(data)## 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 Flabels_stats = c("Median", "Quantile 25%", "Mean", "Quantile 75%", "Interval", "IQR", "Variance", "Standard deviation")number_of_variables = length(names(data))
descriptive_statistics = list()
for (i in 1:number_of_variables) {
unique_vals = unique(na.omit(data[,i]))
if (length(unique_vals) == 2 || length(unique_vals) == 3) {
next # skip this column
}
descriptive_statistics[[i]] = c(
median(data[,i]),
quantile(data[,i], .25),
mean(data[,i]),
quantile(data[,i], .75),
max(data[,i]) - min(data[,i]),
IQR(data[,i]),
var(data[,i]),
sd(data[,i]))
}
sapply(descriptive_statistics, length)## [1] 8 8 0 8 8 8 8descriptive_statistics = descriptive_statistics[sapply(descriptive_statistics, length) == length(labels_stats)]df_all_stats = t(data.frame(row.names = labels_stats, descriptive_statistics))
rownames(df_all_stats) = c("Anni madre", "N gravidanze", "Gestazione", "Peso", "Lunghezza neonato", "Diametro cranio") knitr::kable(df_all_stats, caption = "Statistics", digits = 2)| Median | Quantile 25% | Mean | Quantile 75% | Interval | IQR | Variance | Standard deviation | |
|---|---|---|---|---|---|---|---|---|
| Anni madre | 28 | 25 | 28.16 | 32 | 46 | 7 | 27.81 | 5.27 |
| N gravidanze | 1 | 0 | 0.98 | 1 | 12 | 1 | 1.64 | 1.28 |
| Gestazione | 39 | 38 | 38.98 | 40 | 18 | 2 | 3.49 | 1.87 |
| Peso | 3300 | 2990 | 3284.08 | 3620 | 4100 | 630 | 275665.68 | 525.04 |
| Lunghezza neonato | 500 | 480 | 494.69 | 510 | 255 | 30 | 692.67 | 26.32 |
| Diametro cranio | 340 | 330 | 340.03 | 350 | 155 | 20 | 269.79 | 16.43 |
data$Fumatrici <- factor(data$Fumatrici, levels = c(1, 0), labels = c("Sì", "No"))
data$Tipo.parto <- factor(data$Tipo.parto, levels = c("Nat", "Ces"), labels = c("Naturale", "Cesareo"))
peso_classi <- cut(data$Peso,
breaks = c(500,1000,1500,2000,2500,3000,3500,4000,4500,5000),
labels = c("500-1000", "1000-1500", "1500-2000", "2000-2500",
"2500-3000", "3000-3500", "3500-4000", "4000-4500", "4500-5000"))
p1 <- ggplot(data)+
geom_bar(aes(x = Fumatrici), stat = "count") +
labs(title = "Fumatrici vs Non fumatrici", y = "Counts")
p2 <- ggplot(data)+
geom_bar(aes(x = Fumatrici, fill = Ospedale), position = "dodge", stat = "count") +
labs(title = "Fumatrici vs Non fumatrici per Ospedale", y = "Counts")
p3 <- ggplot(data)+
geom_bar(aes(x = Tipo.parto), stat = "count") +
labs(title = "Tipo parto", y = "Counts")
p4 <- ggplot(data)+
geom_bar(aes(x = Tipo.parto, fill = Ospedale), position = "dodge", stat = "count") +
labs(title = "Tipo parto per Ospedale", y = "Counts")
p5 <- ggplot(data)+
geom_bar(aes(x = peso_classi)) +
labs(title = "Peso neonati", y = "Frequenza assoluta")
p6 <- ggplot(data)+
geom_bar(aes(x = peso_classi, fill = Fumatrici), position = "dodge") +
labs(title = "Peso neonati per Fumatrici", y = "Frequenza assoluta")
p7 <- ggplot(data)+
geom_bar(aes(x = Gestazione)) +
labs(title = "Gestazione", y = "Frequenza assoluta")
p8 <- ggplot(data)+
geom_bar(aes(x = Anni.madre)) +
labs(title = "Anni madre", y = "Frequenza assoluta")
p9 <- ggplot(data)+
geom_bar(aes(x = peso_classi, fill = Ospedale), position = "dodge") +
labs(title = "Peso per Ospedale", y = "Frequenza assoluta")
p10 <- ggplot(data)+
geom_bar(aes(x = Ospedale)) +
labs(title = "Nascite negli ospedali", y = "Frequenza assoluta")
# Mostro i grafici 2 per riga
(p1 | p2) /
(p3 | p4) /
(p5 | p6) /
(p9 | p10) /
(p7 | p8) 
Come si può notare, la quantità di fumatrici è molto minore di quella delle non fumatrici, e ciò è vero in tutti e tre gli ospedali. In tutti e tre gli ospedali è altrettanto vero che il parto naturale è predominante rispetto a quello cesareo. L’andamento della distribuzione del peso dei neonati nei 3 ospedali è molto simile, con picco collocata nel range 3000-3500 grammi. Le nascite considerate come dati sono equamente distribuite nei 3 ospedali. Sia la gestazione che gli anni della madre seguono un andamento simile alla normale.
Test
Più parti cesarei in un ospedale
chisq.test(table(Ospedale, Tipo.parto))##
## Pearson's Chi-squared test
##
## data: table(Ospedale, Tipo.parto)
## X-squared = 1.0972, df = 2, p-value = 0.5778Più parti cesarei in un ospedale
Poichè il p-value è molto maggiore di 0.05, siamo fuori dalla zona di rifiuto, cioè non rifiutiamo l’ipotesi nulla, ovvero che le due variabili siano indipendenti. La risposta è dunque No.
t.test(Peso, mu = 3250)##
## One Sample t-test
##
## data: Peso
## t = 3.2456, df = 2499, p-value = 0.001188
## alternative hypothesis: true mean is not equal to 3250
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081t.test(Lunghezza, mu = 50)##
## One Sample t-test
##
## data: 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.692La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione
Fonte per dati della popolazione: https://doi.org/10.1111/j.1651-2227.2006.tb02378.x
Poichè il p-value è minore di 0.05, rifiutiamo l’ipotesi che la media del peso e della lunghezza di questo campione di neonati siano significativamente uguali a quelle della popolazione.
Le misure antropometriche sono significativamente diverse tra i due sessi
t.test(Peso ~ Sesso)##
## Welch Two Sample t-test
##
## data: Peso by Sesso
## t = -12.106, df = 2490.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.1051 -207.0615
## sample estimates:
## mean in group F mean in group M
## 3161.132 3408.215t.test(Lunghezza ~ Sesso)##
## Welch Two Sample t-test
##
## data: Lunghezza by Sesso
## t = -9.582, df = 2459.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.929470 -7.876273
## sample estimates:
## mean in group F mean in group M
## 489.7643 499.6672t.test(Cranio ~ Sesso)##
## Welch Two Sample t-test
##
## data: Cranio by Sesso
## t = -7.4102, df = 2491.4, p-value = 1.718e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -6.089912 -3.541270
## sample estimates:
## mean in group F mean in group M
## 337.6330 342.4486Poihcè il p-value è minore di 0.05, allora rifiutiamo l’ipotesi che le misure antropometriche siano significativamente uguali tra i due sessi. Dunque sono diverse.
Regressione lineare multipla
Creazione del modello
shapiro.test(Peso)##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97066, p-value < 2.2e-16panel.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)
}
numeric_data <- data[sapply(data, is.numeric)]
pairs(numeric_data, upper.panel = panel.smooth, lower.panel = panel.cor)
modello1 = lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + Lunghezza:Cranio + I(Cranio^2) + I(Lunghezza^2) + I(Gestazione^2) + Gestazione:Cranio + Gestazione:Lunghezza)
summary(modello1)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + Lunghezza:Cranio +
## I(Cranio^2) + I(Lunghezza^2) + I(Gestazione^2) + Gestazione:Cranio +
## Gestazione:Lunghezza)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1152.50 -180.02 -10.19 161.57 1329.74
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.594e+02 1.189e+03 -0.218 0.82737
## Anni.madre 5.055e-01 1.100e+00 0.459 0.64602
## N.gravidanze 1.409e+01 4.528e+00 3.112 0.00188 **
## Fumatrici -2.244e+01 2.673e+01 -0.840 0.40124
## Gestazione 1.738e+02 7.504e+01 2.316 0.02063 *
## Lunghezza -4.802e+00 5.905e+00 -0.813 0.41615
## Cranio -2.328e+01 1.017e+01 -2.290 0.02211 *
## Tipo.partoNat 2.781e+01 1.173e+01 2.371 0.01781 *
## Ospedaleosp2 -8.046e+00 1.306e+01 -0.616 0.53777
## Ospedaleosp3 3.255e+01 1.311e+01 2.483 0.01311 *
## SessoM 7.288e+01 1.090e+01 6.687 2.80e-11 ***
## I(Cranio^2) 5.026e-02 1.869e-02 2.689 0.00722 **
## I(Lunghezza^2) 5.582e-02 6.495e-03 8.594 < 2e-16 ***
## I(Gestazione^2) -4.815e+00 1.611e+00 -2.988 0.00283 **
## Lunghezza:Cranio -8.548e-02 1.670e-02 -5.120 3.29e-07 ***
## Gestazione:Cranio 1.070e+00 2.079e-01 5.147 2.86e-07 ***
## Gestazione:Lunghezza -2.598e-01 1.965e-01 -1.322 0.18618
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 265.7 on 2483 degrees of freedom
## Multiple R-squared: 0.7455, Adjusted R-squared: 0.7439
## F-statistic: 454.6 on 16 and 2483 DF, p-value: < 2.2e-16modelloSTEP = stepAIC(modello1,
direction = "backward",
k = log(nrow(data)))## Start: AIC=28028.17
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso + Lunghezza:Cranio +
## I(Cranio^2) + I(Lunghezza^2) + I(Gestazione^2) + Gestazione:Cranio +
## Gestazione:Lunghezza
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 14899 175332965 28021
## - Fumatrici 1 49768 175367834 28021
## - Gestazione:Lunghezza 1 123464 175441530 28022
## - Ospedale 2 769516 176087582 28024
## - Tipo.parto 1 396961 175715027 28026
## - I(Cranio^2) 1 510531 175828598 28028
## <none> 175318066 28028
## - I(Gestazione^2) 1 630467 175948533 28029
## - N.gravidanze 1 683925 176001991 28030
## - Lunghezza:Cranio 1 1850718 177168784 28047
## - Gestazione:Cranio 1 1870395 177188461 28047
## - Sesso 1 3157308 178475374 28065
## - I(Lunghezza^2) 1 5215296 180533362 28094
##
## Step: AIC=28020.56
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + I(Cranio^2) + I(Lunghezza^2) +
## I(Gestazione^2) + Lunghezza:Cranio + Gestazione:Cranio +
## Gestazione:Lunghezza
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 50115 175383080 28014
## - Gestazione:Lunghezza 1 124391 175457356 28015
## - Ospedale 2 774608 176107574 28016
## - Tipo.parto 1 397018 175729983 28018
## - I(Cranio^2) 1 514387 175847352 28020
## <none> 175332965 28021
## - I(Gestazione^2) 1 633449 175966414 28022
## - N.gravidanze 1 877460 176210425 28025
## - Lunghezza:Cranio 1 1848517 177181482 28039
## - Gestazione:Cranio 1 1865112 177198077 28039
## - Sesso 1 3162513 178495478 28057
## - I(Lunghezza^2) 1 5236199 180569164 28086
##
## Step: AIC=28013.45
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso + I(Cranio^2) + I(Lunghezza^2) + I(Gestazione^2) +
## Lunghezza:Cranio + Gestazione:Cranio + Gestazione:Lunghezza
##
## Df Sum of Sq RSS AIC
## - Gestazione:Lunghezza 1 124265 175507344 28007
## - Ospedale 2 781949 176165029 28009
## - Tipo.parto 1 391622 175774702 28011
## - I(Cranio^2) 1 520065 175903145 28013
## <none> 175383080 28014
## - I(Gestazione^2) 1 632849 176015929 28015
## - N.gravidanze 1 858472 176241552 28018
## - Lunghezza:Cranio 1 1855751 177238831 28032
## - Gestazione:Cranio 1 1864327 177247407 28032
## - Sesso 1 3150297 178533377 28050
## - I(Lunghezza^2) 1 5249843 180632923 28079
##
## Step: AIC=28007.4
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso + I(Cranio^2) + I(Lunghezza^2) + I(Gestazione^2) +
## Lunghezza:Cranio + Gestazione:Cranio
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 792290 176299635 28003
## - Tipo.parto 1 399585 175906930 28005
## <none> 175507344 28007
## - I(Cranio^2) 1 643985 176151329 28009
## - N.gravidanze 1 861613 176368958 28012
## - Gestazione:Cranio 1 1758839 177266184 28025
## - Lunghezza:Cranio 1 2451994 177959339 28034
## - I(Gestazione^2) 1 2781562 178288906 28039
## - Sesso 1 3112565 178619910 28044
## - I(Lunghezza^2) 1 7131190 182638535 28099
##
## Step: AIC=28003.01
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso + I(Cranio^2) + I(Lunghezza^2) + I(Gestazione^2) +
## Lunghezza:Cranio + Gestazione:Cranio
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 423605 176723240 28001
## <none> 176299635 28003
## - I(Cranio^2) 1 622878 176922512 28004
## - N.gravidanze 1 909006 177208641 28008
## - Gestazione:Cranio 1 1745863 178045498 28020
## - Lunghezza:Cranio 1 2471223 178770858 28030
## - I(Gestazione^2) 1 2653463 178953097 28033
## - Sesso 1 3151205 179450840 28040
## - I(Lunghezza^2) 1 7113473 183413107 28094
##
## Step: AIC=28001.19
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## I(Cranio^2) + I(Lunghezza^2) + I(Gestazione^2) + Lunghezza:Cranio +
## Gestazione:Cranio
##
## Df Sum of Sq RSS AIC
## <none> 176723240 28001
## - I(Cranio^2) 1 630928 177354167 28002
## - N.gravidanze 1 877927 177601166 28006
## - Gestazione:Cranio 1 1726823 178450063 28018
## - Lunghezza:Cranio 1 2446318 179169557 28028
## - I(Gestazione^2) 1 2643704 179366943 28031
## - Sesso 1 3151614 179874853 28038
## - I(Lunghezza^2) 1 7117769 183841008 28092summary(modelloSTEP)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Cranio^2) + I(Lunghezza^2) + I(Gestazione^2) +
## Lunghezza:Cranio + Gestazione:Cranio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1167.37 -182.19 -11.71 168.33 1318.88
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.032e+02 1.190e+03 -0.171 0.864479
## N.gravidanze 1.483e+01 4.216e+00 3.516 0.000445 ***
## Gestazione 1.802e+02 7.431e+01 2.425 0.015395 *
## Lunghezza -7.234e+00 5.657e+00 -1.279 0.201045
## Cranio -2.058e+01 1.006e+01 -2.046 0.040877 *
## SessoM 7.275e+01 1.092e+01 6.662 3.30e-11 ***
## I(Cranio^2) 5.483e-02 1.839e-02 2.981 0.002901 **
## I(Lunghezza^2) 5.048e-02 5.042e-03 10.012 < 2e-16 ***
## I(Gestazione^2) -6.298e+00 1.032e+00 -6.102 1.21e-09 ***
## Lunghezza:Cranio -9.269e-02 1.579e-02 -5.870 4.94e-09 ***
## Gestazione:Cranio 1.014e+00 2.056e-01 4.932 8.69e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.5 on 2489 degrees of freedom
## Multiple R-squared: 0.7435, Adjusted R-squared: 0.7424
## F-statistic: 721.3 on 10 and 2489 DF, p-value: < 2.2e-16vif(modelloSTEP)## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif## N.gravidanze Gestazione Lunghezza Cranio
## 1.026122 678.640428 780.105639 960.672823
## Sesso I(Cranio^2) I(Lunghezza^2) I(Gestazione^2)
## 1.049694 1433.526168 561.399700 721.544338
## Lunghezza:Cranio Gestazione:Cranio
## 1932.504954 1638.162089In tal modo abbiamo sfoltito il modello mantenendo solo le variabili, le interazioni e i termini non lineari signficativi.
I coefficienti di multi-collinearità alti sono giustificati dalla presenza di variabili che sono trasformate di altre variabili nel modello.
Bontà del modello
Bisogna valutare R² (coefficiente di determinazione) e RMSE (Root Mean Squared Error):
summary(modelloSTEP)$r.squared## [1] 0.7434661residui = residuals(modelloSTEP)
rmse = sqrt(mean(residui^2))
rmse## [1] 265.8746Il modello interpreta circa il 74% dei dati e in media, le previsioni si discostano dai valori osservati di circa 266 grammi per quanto riguarda il peso dei neonati.
Analisi dei residui
Omoschedasticità
predizioni = predict(modelloSTEP)
plot(predizioni,residui, main = "Residui vs predetti")
abline(h = 0, col = "red")
I residui risultano essere omoschedastici.
QQ PLOT
qqnorm(residui)
qqline(residui, col = "red")
shapiro.test(residui)##
## Shapiro-Wilk normality test
##
## data: residui
## W = 0.99103, p-value = 2.306e-11Il p-value è cosi basso che si rifiuta l’ipotesi nulla secondo cui i residui siano distribuoiti normalmente. Le code nel qqplot, inoltre, si discostano dalla linea rappresentante un andamento normale, a causa della presenza di outlier. ## Leverages
hatvalues <- hatvalues(modelloSTEP)
plot(hatvalues, main = "Leverage", ylab = "Hat values")
abline(h = 2*mean(hatvalues), col = "red")
La maggior parte delle osservazioni ha leverage molto basso, vicino a 0,
che è tipico e atteso. Pochi punti sono chiaramente al di sopra della
linea rossa.
Outliers
plot(rstudent(modelloSTEP))
abline(h = c(-2,2), col = 2)
outlierTest(modelloSTEP)## rstudent unadjusted p-value Bonferroni p
## 1306 4.979022 6.8282e-07 0.0017071
## 155 4.603727 4.3586e-06 0.0108960
## 1399 -4.403743 1.1090e-05 0.0277240
## 1694 4.398650 1.1351e-05 0.0283780Previsioni
Analizziamo cosa succede per le donne alla terza gravidanza quando partoriscono alla 39esima settimana:
coeffs <- summary(modelloSTEP)$coefficients
coeffs_df <- data.frame(
Variabile = rownames(coeffs), # Nome variabili
Coefficiente = coeffs[,1] # Valore numerico coefficienti
)
ggplot(coeffs_df, aes(x = reorder(Variabile, Coefficiente), y = Coefficiente)) +
geom_col(fill = "skyblue") +
# coord_flip() +
labs(title = "Effetto stimato delle variabili sul peso",
x = "", y = "Coefficiente") +
theme_minimal()
Come si nota, la gestazione e il sesso (maschio) hanno un impatto preponderente sul peso del neonato.
# Costruzione del dataset di predizione
nuovi_dati <- data.frame(
Gestazione = seq(30, 42, by = 1),
Anni.madre = mean(data$Anni.madre, na.rm = TRUE), # Non serve ma non fa danni
N.gravidanze = median(data$N.gravidanze, na.rm = TRUE),
Lunghezza = mean(data$Lunghezza, na.rm = TRUE),
Cranio = mean(data$Cranio, na.rm = TRUE),
Sesso = factor("F", levels = c("F", "M")) # Deve avere i livelli giusti
)
# Predizione
nuovi_dati$Predizione <- predict(modelloSTEP, newdata = nuovi_dati)
# Grafico
ggplot(nuovi_dati, aes(x = Gestazione, y = Predizione)) +
geom_line(linewidth = 1.2) +
labs(
title = "Peso previsto vs Gestazione",
x = "Settimane di gestazione",
y = "Peso previsto (grammi)"
) +
theme_minimal()
Si osserva un andamento curvilineo, dato dalla presenza di un termine quadratico nel modello. Ciò ha senso poichè il peso cresce molto nelle prime settimane, per poi stabilizzarsi.
nuovi_casi_cranio = data.frame(Gestazione = 39,
Anni.madre = mean(data$Anni.madre, na.rm = TRUE), # Non serve ma non fa danni
N.gravidanze = median(data$N.gravidanze, na.rm = TRUE),
Lunghezza = mean(data$Lunghezza, na.rm = TRUE),
Cranio = seq(300, 400, by = 1),
Sesso = factor("F", levels = c("F", "M"))) # Deve avere i livelli giusti)
nuovi_casi_cranio$Peso_previsto <- predict(modelloSTEP, newdata = nuovi_casi_cranio)
# Grafico
ggplot(nuovi_casi_cranio, aes(x = Cranio, y = Peso_previsto)) +
geom_line(linewidth = 1.2) +
labs(
title = "Peso previsto vs Cranio",
x = "Cranio",
y = "Peso previsto (grammi)"
) +
theme_minimal()
Il modello suggerisce che il peso del neonato cresce in modo accelerato all’aumentare del cranio, il che è biologicamente plausibile, dato che un cranio più grande tende ad associarsi a un feto più sviluppato e quindi più pesante.
# Grid esteso per cranio e gestazione
nuovi_casi_cranio_gestazione <- expand.grid(
Gestazione = seq(30, 42, by = 1),
Cranio = seq(min(data$Cranio, na.rm = TRUE), max(data$Cranio, na.rm = TRUE), length.out = 20),
Anni.madre = mean(data$Anni.madre, na.rm = TRUE),
N.gravidanze = median(data$N.gravidanze, na.rm = TRUE),
Lunghezza = mean(data$Lunghezza, na.rm = TRUE),
Sesso = factor("F", levels = c("F", "M"))
)
# Predizione
nuovi_casi_cranio_gestazione$Peso_previsto <- predict(modelloSTEP, newdata = nuovi_casi_cranio_gestazione )
# Heatmap
ggplot(nuovi_casi_cranio_gestazione, aes(x = Gestazione,
y = Cranio,
fill = Peso_previsto)) +
geom_tile() +
scale_fill_viridis_c(option = "C") +
labs(title = "Effetto combinato di Gestazione e Cranio sul Peso previsto",
x = "Settimane di Gestazione", y = "Diametro Cranio (cm)", fill = "Peso (g)") +
theme_minimal()