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     F
labels_stats = c("Min", "Max", "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) {
  
  column_data <- data[, i]
  column_data <- na.omit(column_data)
  
  unique_vals <- unique(column_data)
  if (length(unique_vals) == 2 || length(unique_vals) == 3) {
    next
  }

  if (!is.numeric(column_data)) {
    next
  }

  descriptive_statistics[[i]] <- c(
    min(column_data),
    max(column_data),
    median(column_data),
    quantile(column_data, 0.25),
    mean(column_data),
    quantile(column_data, 0.75),
    max(column_data) - min(column_data),
    IQR(column_data),
    var(column_data),
    sd(column_data)
  )
}

sapply(descriptive_statistics, length)
## [1] 10 10  0 10 10 10 10
descriptive_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)
Statistics
Min Max Median Quantile 25% Mean Quantile 75% Interval IQR Variance Standard deviation
Anni madre 0 46 28 25 28.16 32 46 7 27.81 5.27
N gravidanze 0 12 1 0 0.98 1 12 1 1.64 1.28
Gestazione 25 43 39 38 38.98 40 18 2 3.49 1.87
Peso 830 4930 3300 2990 3284.08 3620 4100 630 275665.68 525.04
Lunghezza neonato 310 565 500 480 494.69 510 255 30 692.67 26.32
Diametro cranio 235 390 340 330 340.03 350 155 20 269.79 16.43

Come possiamo notare, per “Anni madre” è presente un minimo pari a 0, il che è chiaramente dissonante rispetto a quanto atteso. Esploriamo allora i dati relativi a tale parametro:

ggplot(data)+
  geom_bar(aes(x = Anni.madre)) +
  labs(title = "Anni madre", y = "Frequenza assoluta")

Osserviamo che sono presenti valori per 0 e 1. Decidiamo di rimuovere le osservazioni in cui Anni madre è pari a 0 o 1.

data = data[data$Anni.madre > 1, ]

Dopo aver filtrato i dati, osserviamo di nuovo la distribuzione tramite l’istogramma:

ggplot(data)+
  geom_bar(aes(x = Anni.madre)) +
  labs(title = "Anni madre", y = "Frequenza assoluta")

I dati sono stati correttamente rimossi.

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.5778

Con un tale p-value, ben più alto della soglia di significatività di 0.05, siamo fuori dalla zona di rifiuto dell’ipotesi di nulla (cioè che le due variabili sono indipendenti). In altre parole, non possiamo rifiutare l’ipotesi che siano indipendenti e quindi non c’è evidenza che in un dato ospedale ci siano più parti cesarei.

La 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

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.081
t.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.692

Poiché i p-value dei test t sono entrambi inferiori alla soglia di significatività di 0.05, rifiutiamo l’ipotesi nulla e concludiamo che la media del peso e della lunghezza del campione di neonati è significativamente diversa da quella della popolazione di riferimento. In dettaglio, per il peso il campione ha una media di 3284.08 g (intervallo di confidenza al 95%: 3263.49 – 3304.67 g), rispetto alla media della popolazione di 3250 g. Per la lunghezza, la media del campione è di 494.69 mm (intervallo di confidenza al 95%: 493.66 – 495.72 mm), mentre la media della popolazione è di 500 mm (50 cm).

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.215
t.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.6672
t.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.4486

Poiché i p-value di tutti i test t per la differenza delle misure antropometriche tra i sessi sono molto inferiori alla soglia di significatività di 0.05, rifiutiamo l’ipotesi nulla di uguaglianza delle medie e concludiamo che peso, lunghezza e circonferenza cranica sono significativamente diversi tra maschi e femmine.

In particolare:

  1. Il peso medio nei maschi (3408.22 g) è significativamente più alto di quello nelle femmine (3161.13 g), con una differenza media stimata tra -287.11 g e -207.06 g (intervallo di confidenza al 95%).

  2. La lunghezza media nei maschi (499.67 mm) è significativamente maggiore di quella nelle femmine (489.76 mm), con una differenza compresa tra -11.93 mm e -7.88 mm.

  3. La circonferenza cranica media nei maschi (342.45 mm) è significativamente superiore a quella delle femmine (337.63 mm), con una differenza stimata tra -6.09 mm e -3.54 mm.

Ciò riflette l’esperienza comune secondo cui biologicamente il genere maschile presenta dimensioni corporee maggiori di quello femminile.

Regressione lineare multipla

Creazione del modello

shapiro.test(Peso)
## 
##  Shapiro-Wilk normality test
## 
## data:  Peso
## W = 0.97066, p-value < 2.2e-16
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r)
}

numeric_data <- data[sapply(data, is.numeric)]
pairs(numeric_data, upper.panel = panel.smooth, lower.panel = panel.cor)

Nel grafico si osserva una correlazione positiva tra Anni madre e N. gravidanze, un risultato atteso e comprensibile: con l’aumentare dell’età, generalmente aumenta anche la stabilità socio-economica e la maturità necessarie per pianificare e sostenere una gravidanza. Questo porta a una maggiore probabilità di avere più gravidanze in età avanzata.

La correlazione tra la durata della Gestazione e variabili come Peso, Lunghezza e Cranio è altrettanto intuitiva e significativa: al progredire della gestazione aumentano le dimensioni del feto, con incrementi proporzionali in peso e misure corporee. Inoltre, le dimensioni del feto sono fortemente correlate tra loro (es. peso e lunghezza, lunghezza e cranio), il che riflette una crescita armonica e proporzionata.

modello1 = lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio  + 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 + Sesso + Lunghezza:Cranio + I(Cranio^2) + 
##     I(Lunghezza^2) + I(Gestazione^2) + Gestazione:Cranio + Gestazione:Lunghezza)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1188.72  -181.73   -12.22   166.25  1321.12 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -1.693e+02  1.192e+03  -0.142  0.88704    
## Anni.madre            5.876e-01  1.103e+00   0.533  0.59428    
## N.gravidanze          1.409e+01  4.538e+00   3.106  0.00192 ** 
## Fumatrici            -2.274e+01  2.679e+01  -0.849  0.39599    
## Gestazione            1.635e+02  7.516e+01   2.175  0.02969 *  
## Lunghezza            -4.685e+00  5.919e+00  -0.792  0.42869    
## Cranio               -2.270e+01  1.019e+01  -2.228  0.02596 *  
## SessoM                7.335e+01  1.093e+01   6.712 2.37e-11 ***
## I(Cranio^2)           4.934e-02  1.874e-02   2.633  0.00851 ** 
## I(Lunghezza^2)        5.616e-02  6.512e-03   8.623  < 2e-16 ***
## I(Gestazione^2)      -4.527e+00  1.614e+00  -2.805  0.00506 ** 
## Lunghezza:Cranio     -8.479e-02  1.672e-02  -5.071 4.25e-07 ***
## Gestazione:Cranio     1.064e+00  2.083e-01   5.106 3.55e-07 ***
## Gestazione:Lunghezza -2.787e-01  1.969e-01  -1.415  0.15714    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 266.5 on 2486 degrees of freedom
## Multiple R-squared:  0.7438, Adjusted R-squared:  0.7424 
## F-statistic: 555.1 on 13 and 2486 DF,  p-value: < 2.2e-16

Escludiamo a priori il Tipo di parto e l’Ospedale poichè per logica sono variabile considerate non influenti per le previsioni del peso di un neonato. Consideriamo il resto delle variabili, qualche possibile interazioni tra variabili ritenute probabilmente significative e similmente i loro termini non lineari.

modelloSTEP = stepAIC(modello1, 
                   direction = "backward",
                   k = log(nrow(data)))
## Start:  AIC=28021.6
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + 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     20149 176528300 28014
## - Fumatrici             1     51172 176559323 28015
## - Gestazione:Lunghezza  1    142194 176650344 28016
## - I(Cranio^2)           1    492383 177000533 28021
## <none>                              176508151 28022
## - I(Gestazione^2)       1    558826 177066977 28022
## - N.gravidanze          1    685024 177193175 28024
## - Lunghezza:Cranio      1   1825740 178333890 28040
## - Gestazione:Cranio     1   1850869 178359019 28040
## - Sesso                 1   3198844 179706994 28059
## - I(Lunghezza^2)        1   5279929 181788079 28088
## 
## Step:  AIC=28014.07
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Cranio^2) + I(Lunghezza^2) + I(Gestazione^2) + 
##     Lunghezza:Cranio + Gestazione:Cranio + Gestazione:Lunghezza
## 
##                        Df Sum of Sq       RSS   AIC
## - Fumatrici             1     51614 176579914 28007
## - Gestazione:Lunghezza  1    143435 176671735 28008
## - I(Cranio^2)           1    496581 177024881 28013
## <none>                              176528300 28014
## - I(Gestazione^2)       1    561775 177090075 28014
## - N.gravidanze          1    894029 177422329 28019
## - Lunghezza:Cranio      1   1822930 178351230 28032
## - Gestazione:Cranio     1   1844486 178372786 28032
## - Sesso                 1   3204986 179733286 28051
## - I(Lunghezza^2)        1   5304012 181832312 28080
## 
## Step:  AIC=28006.97
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 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    143326 176723240 28001
## - I(Cranio^2)           1    502035 177081948 28006
## <none>                              176579914 28007
## - I(Gestazione^2)       1    560992 177140906 28007
## - N.gravidanze          1    874849 177454763 28012
## - Lunghezza:Cranio      1   1830196 178410109 28025
## - Gestazione:Cranio     1   1843593 178423507 28025
## - Sesso                 1   3192614 179772528 28044
## - I(Lunghezza^2)        1   5317932 181897845 28073
## 
## Step:  AIC=28001.18
## 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 28092

Progressivamente le variabili meno significative sono state eliminate secondo il criterio BIC. Come atteso, variabili come la gestazione e le dimensioni del feto risultano significative nel modello di predizione del peso. Infatti, al progredire della gravidanza, il feto cresce e dunque anche le sue dimensioni aumentano. Coerentemente, un feto più grande mostra anche un peso maggiore. Sono state mantenute anche interazioni tra tali variabili significative. In sintesi, abbiamo sfoltito il modello mantenendo solo le variabili, le interazioni e i termini non lineari signficativi.

summary(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-16

Il modello finale, selezionato tramite procedura backward basata su BIC, spiega il 74% della variabilità del Peso, visto che R² aggiustato = 0.742, ed è globalmente significativo, poichè p < 2.2e-16. Le variabili significative includono: N.gravidanze, Gestazione, Cranio, Sesso, termini quadratici e interazioni tra Lunghezza:Cranio e Gestazione:Cranio, evidenziando relazioni lineari e non lineari complesse. Il termine lineare di Lunghezza non risulta significativo, ma viene mantenuto per la presenza di suoi effetti quadratici e interattivi.

vif(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.162089

I valori elevati dei VIF sono dovuti alla presenza di termini quadratici e interazioni, che per definizione sono correlati alle variabili originali e tra loro.

Bontà del modello

Bisogna valutare R² (coefficiente di determinazione) e RMSE (Root Mean Squared Error):

summary(modelloSTEP)$r.squared
## [1] 0.7434661
residui = residuals(modelloSTEP)
rmse = sqrt(mean(residui^2))
rmse
## [1] 265.8746

Il 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-11

Il 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 e bisogna prestare attenzioni ad essi in fase di analisi dati.

Outliers

plot(rstudent(modelloSTEP))
abline(h = c(-2, 2), col = "blue", lty = 2, lwd = 3)
abline(h = c(-3, 3), col = "red", lty = 2, lwd = 3)

Come si può osservare, sono presenti outliers sia oltre la soglia (-2,2) sia oltre (-3,3). Questi ultimi sono outliers forti e si ipotizza siano errori di misura oppure condizioni molto particolari di gravidanza, associabili a patologie o condizioni cliniche pericolose e insolite.

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.0283780

Previsioni

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. I coefficienti positivi per N.gravidanze, Gestazione e SessoM indicano che un maggior numero di gravidanze, una gestazione più lunga e il sesso maschile sono associati a un aumento del peso neonatale. Il coefficiente negativo per Cranio e Lunghezza è più che compensato dal coefficiente positivo per il loro termine quadratico, variabili per cui quindi vale quando detto per le tre variabili precedenti N. gravidanze, Gestazione e SessoM.

# 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()

Il grafico mostra l’effetto combinato delle settimane di gestazione e del diametro cranico sul peso previsto del neonato. Si osserva che un aumento sia della durata della gestazione sia del diametro cranico è associato a un incremento del peso alla nascita. L’interazione tra le due variabili è evidente: i pesi maggiori si concentrano nelle zone con valori elevati per entrambi i predittori, mentre i pesi più bassi si trovano nei casi di gestazione precoce e cranio più piccolo.