progetto statistica inferenziale p

Dettagli del Progetto 1.

Raccolta dei Dati e Struttura del Dataset Per costruire il modello predittivo, abbiamo raccolto dati su 2500 neonati provenienti da tre ospedali. Le variabili raccolte includono:

Età della madre: Misura dell’età in anni.
Numero di gravidanze: Quante gravidanze ha avuto la madre.
Fumo materno: Un indicatore binario (0=non fumatrice, 1=fumatrice).
Durata della gravidanza: Numero di settimane di gestazione.
Peso del neonato: Peso alla nascita in grammi.
Lunghezza e diametro del cranio: Lunghezza del neonato e diametro craniale, misurabili anche durante la gravidanza tramite ecografie.
Tipo di parto: Naturale o cesareo. Ospedale di nascita: Ospedale 1, 2 o 3.
Sesso del neonato: Maschio (M) o femmina (F).

L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature.

suppressPackageStartupMessages(library("ggplot2"))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(car))
suppressPackageStartupMessages(library(lmtest))

data <- read.csv("neonati.csv", stringsAsFactors = T)
data$Fumatrici <- as.factor(data$Fumatrici)
suppressWarnings(attach(data))
summary(data)

##    Anni.madre     N.gravidanze     Fumatrici   Gestazione         Peso     
##  Min.   : 0.00   Min.   : 0.0000   0:2396    Min.   :25.00   Min.   : 830  
##  1st Qu.:25.00   1st Qu.: 0.0000   1: 104    1st Qu.:38.00   1st Qu.:2990  
##  Median :28.00   Median : 1.0000             Median :39.00   Median :3300  
##  Mean   :28.16   Mean   : 0.9812             Mean   :38.98   Mean   :3284  
##  3rd Qu.:32.00   3rd Qu.: 1.0000             3rd Qu.:40.00   3rd Qu.:3620  
##  Max.   :46.00   Max.   :12.0000             Max.   :43.00   Max.   :4930  
##    Lunghezza         Cranio    Tipo.parto Ospedale   Sesso   
##  Min.   :310.0   Min.   :235   Ces: 728   osp1:816   F:1256  
##  1st Qu.:480.0   1st Qu.:330   Nat:1772   osp2:849   M:1244  
##  Median :500.0   Median :340              osp3:835           
##  Mean   :494.7   Mean   :340                                 
##  3rd Qu.:510.0   3rd Qu.:350                                 
##  Max.   :565.0   Max.   :390

colSums(is.na(data))

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

n <- dim(data)[1]
data$Anni.madre <- ifelse(data$Anni.madre <= 12, median(data$Anni.madre), data$Anni.madre)

Analisi e Modellizzazione Analisi Preliminare Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.

statistiche <- function(variabile){
  mean_variabile <- sum(variabile)/length(variabile)
  median_variabile <- median(variabile)
  var_variabile <- var(variabile)
  dev_std_variabile <- sd(variabile)
  stats <- c(Media= mean_variabile,
             Mediana= median_variabile,
             Varianza= var_variabile,
             Deviazione_Standard= dev_std_variabile)
  stats_approx <- round(stats,digits = 2)
  
  quantile_variabile=round(quantile(variabile),digits = 2)
  stats_finali <- c(stats_approx,Quantile=quantile_variabile)
  
  return(stats_finali)
}

variabili_quant <- data[c("Anni.madre","N.gravidanze","Gestazione","Peso","Lunghezza","Cranio")]

risultati_statistiche <- lapply(variabili_quant,statistiche)
df_statistiche <- as.data.frame(do.call(cbind,risultati_statistiche))
kable(df_statistiche)

	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio
Media	28.19	0.98	38.98	3284.08	494.69	340.03
Mediana	28.00	1.00	39.00	3300.00	500.00	340.00
Varianza	27.20	1.64	3.49	275665.68	692.67	269.79
Deviazione_Standard	5.22	1.28	1.87	525.04	26.32	16.43
Quantile.0%	13.00	0.00	25.00	830.00	310.00	235.00
Quantile.25%	25.00	0.00	38.00	2990.00	480.00	330.00
Quantile.50%	28.00	1.00	39.00	3300.00	500.00	340.00
Quantile.75%	32.00	1.00	40.00	3620.00	510.00	350.00
Quantile.100%	46.00	12.00	43.00	4930.00	565.00	390.00

df_num <- data[, sapply(data, is.numeric)]
moments::skewness(df_num)

##   Anni.madre N.gravidanze   Gestazione         Peso    Lunghezza       Cranio 
##    0.1512083    2.5142541   -2.0653133   -0.6470308   -1.5146991   -0.7850527

shapiro.test(Peso)

## 
##  Shapiro-Wilk normality test
## 
## data:  Peso
## W = 0.97066, p-value < 2.2e-16

round(cor(df_num),2)

##              Anni.madre N.gravidanze Gestazione  Peso Lunghezza Cranio
## Anni.madre         1.00         0.38      -0.13 -0.02     -0.06   0.02
## N.gravidanze       0.38         1.00      -0.10  0.00     -0.06   0.04
## Gestazione        -0.13        -0.10       1.00  0.59      0.62   0.46
## Peso              -0.02         0.00       0.59  1.00      0.80   0.70
## Lunghezza         -0.06        -0.06       0.62  0.80      1.00   0.60
## Cranio             0.02         0.04       0.46  0.70      0.60   1.00

L’analisi circa l’asimmetria delle variabili riporta i seguenti risultati: le variabili Peso, Lunghezza, Gestazione e Cranio presentano dei valori negativi, sottolineando un’asimmetria negativa e dunque una coda sinistra pronunciata. N.gravidanze riporta valori altamente positivi, in questo caso la coda sarà verso destra per via dell’asimmetria positiva. L’unica variabile che presenta una distribuzione simmetrica è Anni.Madre.

Il test di normalità di Shapiro-Wilk per la variabile target Peso conferma quanto visto nell’analisi precedente: con un valore del p-value inferiore di 0.05 si rifiuta l’ipotesi nulla di normalità nella distribuzione.

La matrice di correlazione mostra come le variabili Peso, Lunghezza e Cranio siano fortemente correlate. Con un valore più basso da sottolineare è la correlazione tra Gestazione e le variabili sopracitate.

df_only_x <- subset(df_num, select = setdiff(names(df_num), "Peso"))

panel.cor <- function(df_only_x, Peso, digits = 2, prefix = "", cex.cor, ...)
{
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(df_only_x, Peso))
  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)
}
pairs(df_num,upper.panel=panel.smooth, lower.panel = panel.cor)

Lo scatterplot matrix conferma quanto detto precedentemente: le variabili “fisiche” come Peso, Lunghezza e Cranio mostrano una correlazione maggiore, anche graficamente possiamo osservare come i punti nelle relazioni tra queste variabili mostrino un andamento crescente. Anche la variabile Gestazione mostra un andamento positivo se posta in relazione con le variabili fisiche del neonato. discorso differente, sia in termini di correlazione che in termini di relazione con le altre variabili per le variabili N.gravidanze e Anni.madre: i punti appaiono sparsi, mostrando un trend abbastanza piatto ad eccezione della relazione tra le stesse variabili.

Rappresentazioni grafiche

data$fascia_età1 <- cut(
  data$Anni.madre,
  breaks= c(12,20,27,34,41,48),
  labels= c("12-19","20-26","27-33","34-40","41-47"))

ggplot(data=data,aes(x=fascia_età1,y=Peso))+
  geom_boxplot(fill="lightblue",color="darkblue",alpha=0.8)+
  labs(title = "Peso alla nascita per fasce d'età materna",
       y="Peso neonato (grammi)",
       x="Fascia d'età madre (anni)")+
  theme_minimal(base_size = 14)+
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
    axis.title = element_text(face = "bold"))

Il boxplot conferma quanto visto nello scatterplot matrix: il Peso non è influenzato in modo significativo dalla sola età della madre, lo possiamo osservare dalla mediana, questa risulta stabile per tutte le fasce d’età. la presenza di molti outlier sottolinea la necessità di approfondire il fenomeno con ulteriori variabili.

ggplot(data=data,aes(x=Fumatrici,y=Peso))+
  geom_boxplot(fill="lightblue",col="darkblue",alpha=0.8)+
  labs(title="Peso alla nascita rispetto allo stato di fumo delle madri",
       x="Stato di fumo",
       y="Peso (grammi)")+
  theme_minimal(base_size = 14)+
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
    axis.title = element_text(face = "bold"))

il boxplot mostra come sia associato un minor peso per i neonati che hanno una madre fumatrice. però il numero di osservazioni risulta abbastanza ridotto per confermare tale ipotesi.

ggplot(data=data,aes(x=Gestazione,y=Peso))+
  geom_point(fill="lightblue",col="darkblue",alpha=0.8)+
  labs(title="Peso alla nascita rispetto alle settimane di gestazione",
       x="Periodo di gestazione (settimane)",
       y="Peso neonato (in grammi)")+
  theme_minimal(base_size = 14)+
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
    axis.title = element_text(face = "bold"))

Il grafico conferma quanto visto nelle analisi precedenti, ossia la correlazione positvia tra il periodo di gestazione e il peso del neonato.

ggplot(data=data,aes(x=Lunghezza,y=Peso))+
  geom_point(fill="lightblue",col="darkblue",alpha=0.8)+
  labs(title="Peso alla nascita rispetto alla lunghezza del neonato",
       x="Lunghezza neonato (millimetri)",
       y="Peso neonato (grammi)")+
  theme_minimal(base_size = 14)+
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
    axis.title = element_text(face = "bold"))

ggplot(data=data,aes(x=Cranio,y=Peso))+
  geom_point(fill="lightblue",col="darkblue",alpha=0.8)+
  labs(title="Peso alla nascita rispetto al diametro del cranio",
       x="Diametro cranio (millimetri)",
       y="Peso neonato (grammi)")+
  theme_minimal(base_size = 14)+
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
    axis.title = element_text(face = "bold"))

I due grafici trovano conferma con quanto detto precedentemente, il Peso è correlato in modo positivo con la lunghezza e la misura del cranio del neonato.

ggplot(data=data,aes(x=Sesso,y=Peso))+
  geom_boxplot(fill="lightblue",color="darkblue",alpha=0.8)+
  labs(title = "Peso alla nascita rispetto al sesso",
       x="Sesso neonato",
       y="Peso neonato (grammi)")+
  theme_minimal(base_size = 14)+
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
    axis.title = element_text(face = "bold"))

Il boxplot mostra come i neonati maschi tendenzialmente nascano con un peso maggiore rispetto alle femmine.

Inoltre si saggeranno le seguenti ipotesi con i test adatti:

in alcuni ospedali si fanno più parti cesarei

tabella_contingenza <- table(data$Ospedale,data$Tipo.parto)
chiquad <- chisq.test(tabella_contingenza)
chiquad

## 
##  Pearson's Chi-squared test
## 
## data:  tabella_contingenza
## X-squared = 1.0972, df = 2, p-value = 0.5778

con un p-value di gran lunga superiore a 0.05, possiamo accettare l’ipotesi nulla, confermando l’indipendenza tra le due variabili.

La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione

mu_peso=3300 #dati OMS
t.test(data$Peso, mu= mu_peso)

## 
##  One Sample t-test
## 
## data:  data$Peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.490 3304.672
## sample estimates:
## mean of x 
##  3284.081

mu_lunghezza=500 #dati OMS
t.test(data$Lunghezza,mu=mu_lunghezza)

## 
##  One Sample t-test
## 
## data:  data$Lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.6598 495.7242
## sample estimates:
## mean of x 
##   494.692

Con un p-value di 0.1296 (maggiore di 0.05), non rifiutiamo l’ipotesi nulla. Pertanto, la media del peso della popolazione da cui è stato estratto il campione non differisce dalla media dell’OMS (3300g). Con un p-value di 2.2e-16 (inferiore di 0.05), rifiutiamo l’ipotesi nulla. Pertanto, la media della lunghezza della popolazione da cui è stato estratto il campione, differisce dalla media dell’OMS (500mm).

Le misure antropometriche sono significativamente diverse tra i due sessi

t.test(Peso ~ Sesso, data=data)

## 
##  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, data=data)

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

con un p-value di 2.2e-16 (minore di 0.05), rifiutiamo l’ipotesi nulla, confermando quanto visto precedentemente: le medie del peso dei neonati differiscono in base al sesso.

con un p-value di 2.2e-16 (minore di 0.05), rifiutiamo l’ipotesi nulla, confermando quanto visto precedentemente: le medie della lunghezza dei neonati differiscono in base al sesso.

Creazione del Modello di Regressione

Verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti. In questo modo, potremo quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato ed eventuali interazioni.
Ad esempio, ci aspettiamo che una maggiore durata della gestazione aumenterebbe in media il peso del neonato.

modello_completo <- lm(Peso~Anni.madre + N.gravidanze + Sesso + Ospedale + Tipo.parto + Cranio + Lunghezza + Fumatrici + Gestazione, data = data)
summary(modello_completo) #commento sul coefficiente della gestazione

## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Sesso + Ospedale + 
##     Tipo.parto + Cranio + Lunghezza + Fumatrici + Gestazione, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1123.3  -181.2   -14.6   160.7  2612.6 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6735.1677   141.3977 -47.633  < 2e-16 ***
## Anni.madre        0.7983     1.1463   0.696   0.4862    
## N.gravidanze     11.4118     4.6665   2.445   0.0145 *  
## SessoM           77.5473    11.1779   6.938 5.07e-12 ***
## Ospedaleosp2    -11.2216    13.4388  -0.835   0.4038    
## Ospedaleosp3     28.0984    13.4972   2.082   0.0375 *  
## Tipo.partoNat    29.5027    12.0848   2.441   0.0147 *  
## Cranio           10.4725     0.4261  24.580  < 2e-16 ***
## Lunghezza        10.2951     0.3007  34.237  < 2e-16 ***
## Fumatrici1      -30.1567    27.5396  -1.095   0.2736    
## Gestazione       32.5265     3.8179   8.520  < 2e-16 ***
## ---
## 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

ggplot(data, aes(x = Gestazione, y = Peso)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", color = "blue")

## `geom_smooth()` using formula = 'y ~ x'

il modello nel complesso spiega circa il 72.7% della variabilità del Peso dei neonati (osservando Rquadro e Rquadro aggiustato), è altamente significativo con un valore del p_value pari a < 2.2e-16. Nello specifico, le variabili più significative, ossia quelle che hanno ottenuto un p-value inferiore a 0.01 oppure 0.001, sono: sesso maschile con un coefficiente positivo di 77.5, questo singifica che tendenzialmente i neonati maschi nascono con 77.5g in più rispetto alle femmine; gestazione con un coefficiente positivo di 32.5g per ogni settimana in più di gestazione. cranio e lunghezza con un impatto positivo sul peso del bambino. meno significative (al di sotto della soglia di 0.05) sono: tipo di parto, n.gravidanze e ospedale.

il grafico conferma la forte relazione lineare positiva spiegata nelle precedenti analisi tra il peso e le settimane di gestazione, possiamo notare come ci sia un aumento del peso all’aumentare delle settimane di gestazione, sottolineando la singificatività di questa variabile per prevedere il peso dei neonati.

Selezione del Modello Ottimale Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC), selezioneremo il modello più parsimonioso, eliminando le variabili non significative. Verranno considerati anche modelli con interazioni tra le variabili e possibili effetti non lineari.

modello_BIC_step <- MASS::stepAIC(modello_completo,
                              direction = "both",
                              k=log(n))

## Start:  AIC=28139.46
## Peso ~ Anni.madre + N.gravidanze + Sesso + Ospedale + Tipo.parto + 
##     Cranio + Lunghezza + Fumatrici + Gestazione
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     36392 186809099 28132
## - Fumatrici     1     89979 186862686 28133
## - Ospedale      2    686397 187459103 28133
## - Tipo.parto    1    447233 187219939 28138
## - N.gravidanze  1    448762 187221469 28138
## <none>                      186772707 28140
## - Sesso         1   3611588 190384294 28180
## - Gestazione    1   5446558 192219264 28204
## - Cranio        1  45338051 232110758 28675
## - Lunghezza     1  87959790 274732497 29096
## 
## Step:  AIC=28132.12
## Peso ~ N.gravidanze + Sesso + Ospedale + Tipo.parto + Cranio + 
##     Lunghezza + Fumatrici + Gestazione
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     90897 186899996 28126
## - Ospedale      2    692738 187501837 28126
## - Tipo.parto    1    448222 187257321 28130
## <none>                      186809099 28132
## - N.gravidanze  1    633756 187442855 28133
## + Anni.madre    1     36392 186772707 28140
## - Sesso         1   3618736 190427835 28172
## - Gestazione    1   5412879 192221978 28196
## - Cranio        1  45588236 232397335 28670
## - Lunghezza     1  87950050 274759149 29089
## 
## Step:  AIC=28125.51
## Peso ~ N.gravidanze + Sesso + Ospedale + Tipo.parto + Cranio + 
##     Lunghezza + Gestazione
## 
##                Df Sum of Sq       RSS   AIC
## - Ospedale      2    701680 187601677 28119
## - Tipo.parto    1    440684 187340680 28124
## <none>                      186899996 28126
## - N.gravidanze  1    610840 187510837 28126
## + Fumatrici     1     90897 186809099 28132
## + Anni.madre    1     37311 186862686 28133
## - Sesso         1   3602797 190502794 28165
## - Gestazione    1   5346781 192246777 28188
## - Cranio        1  45632149 232532146 28664
## - Lunghezza     1  88355030 275255027 29086
## 
## Step:  AIC=28119.23
## Peso ~ N.gravidanze + Sesso + Tipo.parto + Cranio + Lunghezza + 
##     Gestazione
## 
##                Df Sum of Sq       RSS   AIC
## - Tipo.parto    1    463870 188065546 28118
## <none>                      187601677 28119
## - N.gravidanze  1    651066 188252743 28120
## + Ospedale      2    701680 186899996 28126
## + Fumatrici     1     99840 187501837 28126
## + Anni.madre    1     43845 187557831 28126
## - Sesso         1   3649259 191250936 28160
## - Gestazione    1   5444109 193045786 28183
## - Cranio        1  45758101 233359778 28657
## - Lunghezza     1  88054432 275656108 29074
## 
## Step:  AIC=28117.58
## Peso ~ N.gravidanze + Sesso + Cranio + Lunghezza + Gestazione
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      188065546 28118
## - N.gravidanze  1    623141 188688687 28118
## + Tipo.parto    1    463870 187601677 28119
## + Ospedale      2    724866 187340680 28124
## + Fumatrici     1     91892 187973654 28124
## + Anni.madre    1     45044 188020502 28125
## - Sesso         1   3655292 191720838 28158
## - Gestazione    1   5464853 193530399 28181
## - Cranio        1  46108583 234174130 28658
## - Lunghezza     1  87632762 275698308 29066

summary(modello_BIC_step)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Sesso + Cranio + Lunghezza + 
##     Gestazione, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.44  -180.81   -15.58   163.64  2639.72 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.1445   135.7229 -49.226  < 2e-16 ***
## N.gravidanze    12.4750     4.3396   2.875  0.00408 ** 
## SessoM          77.9927    11.2021   6.962 4.26e-12 ***
## Cranio          10.5402     0.4262  24.728  < 2e-16 ***
## Lunghezza       10.2486     0.3006  34.090  < 2e-16 ***
## Gestazione      32.3321     3.7980   8.513  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1328 on 5 and 2494 DF,  p-value: < 2.2e-16

modello_AIC_step <- MASS::stepAIC(modello_completo,
                                  direction = "both",
                                  k=2)

## Start:  AIC=28075.39
## Peso ~ Anni.madre + N.gravidanze + Sesso + Ospedale + Tipo.parto + 
##     Cranio + Lunghezza + Fumatrici + Gestazione
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     36392 186809099 28074
## - Fumatrici     1     89979 186862686 28075
## <none>                      186772707 28075
## - Tipo.parto    1    447233 187219939 28079
## - N.gravidanze  1    448762 187221469 28079
## - Ospedale      2    686397 187459103 28081
## - Sesso         1   3611588 190384294 28121
## - Gestazione    1   5446558 192219264 28145
## - Cranio        1  45338051 232110758 28617
## - Lunghezza     1  87959790 274732497 29038
## 
## Step:  AIC=28073.88
## Peso ~ N.gravidanze + Sesso + Ospedale + Tipo.parto + Cranio + 
##     Lunghezza + Fumatrici + Gestazione
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     90897 186899996 28073
## <none>                      186809099 28074
## + Anni.madre    1     36392 186772707 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 + Sesso + Ospedale + Tipo.parto + Cranio + 
##     Lunghezza + Gestazione
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      186899996 28073
## + Fumatrici     1     90897 186809099 28074
## + Anni.madre    1     37311 186862686 28075
## - Tipo.parto    1    440684 187340680 28077
## - Ospedale      2    701680 187601677 28078
## - 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

summary(modello_AIC_step)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Sesso + Ospedale + Tipo.parto + 
##     Cranio + Lunghezza + Gestazione, data = data)
## 
## 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 ** 
## SessoM           77.4412    11.1756   6.930 5.36e-12 ***
## Ospedaleosp2    -11.0227    13.4363  -0.820  0.41209    
## Ospedaleosp3     28.6408    13.4886   2.123  0.03382 *  
## Tipo.partoNat    29.2803    12.0817   2.424  0.01544 *  
## Cranio           10.4922     0.4254  24.661  < 2e-16 ***
## Lunghezza        10.3086     0.3004  34.316  < 2e-16 ***
## Gestazione       31.9909     3.7896   8.442  < 2e-16 ***
## ---
## 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

utilizzando sia il BIC che l’AIC otteniamo due modelli differenti, il modello BIC è un modello più piccolo, questo perchè viene penalizzato maggiormente l’aggiunta di variabili nel modello: dunque l’aumento della complessità a un certo punto non giustifica l’aumento minimo dell’adattamento. l’AIC invece impone una penalità migliore alle variabili che vengono aggiunte, di conseguenza avremo un modello con più variabili. Confrontando i risultati dei modelli possiamo osservare come la variabilità spiegata da entrambi sia abbastanza simile: per il modello BIC abbiamo un Rquadro pari a 0.7270 e per il modello AIC 0.7287. inoltre, il modello ottenuto con il BIC ha tutte le variabili con un p-value significativo (tutte le variabili ad eccezione di N.gravidanze hanno un p-value inferiore a 0.001, N.gravidanze ha un p-value inferiore a 0.01).

in questo blocco di codice ho cercato, di trovare il modello più accurato, aggiungendo anche dei termini di interazioni tra le variabili e non linearità.

mod3 <- update(modello_BIC_step,~.-N.gravidanze)
summary(mod3)

## 
## Call:
## lm(formula = Peso ~ Sesso + Cranio + Lunghezza + Gestazione, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.2  -184.3   -17.6   163.3  2627.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6651.1188   135.5172 -49.080  < 2e-16 ***
## SessoM         79.1049    11.2117   7.056 2.22e-12 ***
## Cranio         10.6704     0.4245  25.139  < 2e-16 ***
## Lunghezza      10.2054     0.3007  33.939  < 2e-16 ***
## Gestazione     31.2737     3.7856   8.261 2.31e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2495 degrees of freedom
## Multiple R-squared:  0.7261, Adjusted R-squared:  0.7257 
## F-statistic:  1654 on 4 and 2495 DF,  p-value: < 2.2e-16

mod4 <- lm(Peso ~ Gestazione*Fumatrici+Lunghezza+Cranio+Sesso)
summary(mod4)

## 
## Call:
## lm(formula = Peso ~ Gestazione * Fumatrici + Lunghezza + Cranio + 
##     Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1138.31  -182.30   -18.44   164.51  2624.11 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -6667.9088   136.4437 -48.869  < 2e-16 ***
## Gestazione               32.0177     3.8254   8.370  < 2e-16 ***
## Fumatrici1              777.3046   758.4089   1.025    0.306    
## Lunghezza                10.1865     0.3011  33.836  < 2e-16 ***
## Cranio                   10.6641     0.4245  25.123  < 2e-16 ***
## SessoM                   79.8505    11.2261   7.113 1.48e-12 ***
## Gestazione:Fumatrici1   -20.4706    19.3051  -1.060    0.289    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2493 degrees of freedom
## Multiple R-squared:  0.7263, Adjusted R-squared:  0.7257 
## F-statistic:  1103 on 6 and 2493 DF,  p-value: < 2.2e-16

mod5 <- lm(Peso ~ Gestazione*N.gravidanze+Lunghezza+Cranio)
summary(mod5)

## 
## Call:
## lm(formula = Peso ~ Gestazione * N.gravidanze + Lunghezza + Cranio)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1117.73  -183.37   -12.96   168.49  2636.09 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             -6711.3721   159.2765 -42.137  < 2e-16 ***
## Gestazione                 30.3151     4.4067   6.879 7.58e-12 ***
## N.gravidanze              -68.2071    70.7718  -0.964    0.335    
## Lunghezza                  10.4681     0.3018  34.688  < 2e-16 ***
## Cranio                     10.6525     0.4300  24.772  < 2e-16 ***
## Gestazione:N.gravidanze     2.1064     1.8206   1.157    0.247    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 277.2 on 2494 degrees of freedom
## Multiple R-squared:  0.7218, Adjusted R-squared:  0.7213 
## F-statistic:  1294 on 5 and 2494 DF,  p-value: < 2.2e-16

mod6 <- lm(Peso ~ Gestazione*Sesso+Lunghezza+Cranio)
summary(mod6)

## 
## Call:
## lm(formula = Peso ~ Gestazione * Sesso + Lunghezza + Cranio)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1133.41  -183.10   -17.12   162.35  2625.80 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -6539.1643   166.8151 -39.200  < 2e-16 ***
## Gestazione           28.3239     4.5716   6.196 6.77e-10 ***
## SessoM             -191.1270   235.0952  -0.813    0.416    
## Lunghezza            10.2095     0.3007  33.953  < 2e-16 ***
## Cranio               10.6713     0.4244  25.143  < 2e-16 ***
## Gestazione:SessoM     6.9246     6.0174   1.151    0.250    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2494 degrees of freedom
## Multiple R-squared:  0.7262, Adjusted R-squared:  0.7257 
## F-statistic:  1323 on 5 and 2494 DF,  p-value: < 2.2e-16

mod7 <- lm(Peso ~ Gestazione+I(Gestazione^2)+Lunghezza+Cranio)
summary(mod7)

## 
## Call:
## lm(formula = Peso ~ Gestazione + I(Gestazione^2) + Lunghezza + 
##     Cranio)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1170.92  -182.76   -13.31   166.15  2652.23 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -4245.2988   905.8184  -4.687 2.93e-06 ***
## Gestazione       -109.4206    50.0637  -2.186  0.02894 *  
## I(Gestazione^2)     1.8835     0.6663   2.827  0.00474 ** 
## Lunghezza          10.5477     0.3047  34.621  < 2e-16 ***
## Cranio             10.9058     0.4296  25.384  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 277.3 on 2495 degrees of freedom
## Multiple R-squared:  0.7215, Adjusted R-squared:  0.7211 
## F-statistic:  1616 on 4 and 2495 DF,  p-value: < 2.2e-16

ho raggruppato il tutto per creare degli insights che potessero aiutarmi nella decisione.

statistiche_modelli <- function(model,name){
  r2 <- summary(model)$r.squared
  rmse <- sqrt(mean(resid(model)^2))
  bic <- BIC(model)
  
  data.frame(Modello=name,
             R2=r2,
             RMSE=rmse,
             BIC=bic)}

risultati <- rbind(statistiche_modelli(modello_completo,"Modello Completo"),
                   statistiche_modelli(modello_AIC_step,"modello_AIC_step"),
                   statistiche_modelli(modello_BIC_step,"modello_BIC_step"),
                   statistiche_modelli(mod3,"mod3"),
                   statistiche_modelli(mod4,"mod4"),
                   statistiche_modelli(mod5,"mod5"),
                   statistiche_modelli(mod6,"mod6"),
                   statistiche_modelli(mod6,"mod7"))
kable(risultati)

Modello	R2	RMSE	BIC
Modello Completo	0.7288782	273.3296	35241.97
modello_AIC_step	0.7286934	273.4227	35228.03
modello_BIC_step	0.7270015	274.2740	35220.10
mod3	0.7260969	274.7280	35220.54
mod4	0.7263205	274.6158	35234.15
mod5	0.7218447	276.8523	35266.88
mod6	0.7262423	274.6551	35227.04
mod7	0.7262423	274.6551	35227.04

possiamo notare come l’Rquadro sia abbastanza simile per tutti i modelli, ho scelto di utilizzare il mod3 in quanto utilizza meno predittori mantenendo performance simili ai modelli più grandi, in modo tale da rispettare l’equilibrio tra la predizione e l’interpetabilità dei dati.

anova(modello_completo,modello_AIC_step,modello_BIC_step,mod3,mod4,mod5,mod6,mod7)

## Analysis of Variance Table
## 
## Model 1: Peso ~ Anni.madre + N.gravidanze + Sesso + Ospedale + Tipo.parto + 
##     Cranio + Lunghezza + Fumatrici + Gestazione
## Model 2: Peso ~ N.gravidanze + Sesso + Ospedale + Tipo.parto + Cranio + 
##     Lunghezza + Gestazione
## Model 3: Peso ~ N.gravidanze + Sesso + Cranio + Lunghezza + Gestazione
## Model 4: Peso ~ Sesso + Cranio + Lunghezza + Gestazione
## Model 5: Peso ~ Gestazione * Fumatrici + Lunghezza + Cranio + Sesso
## Model 6: Peso ~ Gestazione * N.gravidanze + Lunghezza + Cranio
## Model 7: Peso ~ Gestazione * Sesso + Lunghezza + Cranio
## Model 8: Peso ~ Gestazione + I(Gestazione^2) + Lunghezza + Cranio
##   Res.Df       RSS Df Sum of Sq       F    Pr(>F)    
## 1   2489 186772707                                   
## 2   2491 186899996 -2   -127290  0.8482  0.428328    
## 3   2494 188065546 -3  -1165550  5.1775  0.001444 ** 
## 4   2495 188688687 -1   -623141  8.3042  0.003989 ** 
## 5   2493 188534660  2    154027  1.0263  0.358478    
## 6   2494 191617989 -1  -3083329 41.0895 1.734e-10 ***
## 7   2494 188588551  0   3029438                      
## 8   2495 191839039 -1  -3250488 43.3172 5.656e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Il test anova sottolinea l’importanza dei predittori lunghezza, cranio, gestazione e sesso. anche la variabile N.gravidanze risulta significativa, ma produce dei miglioramenti minimi, quindi per ragioni di parsimonia viene scartata.

vif(mod3)

##      Sesso     Cranio  Lunghezza Gestazione 
##   1.038813   1.606131   2.069517   1.653502

il tasso di inflazione della varianza mostra come i preditori del mod3 non siano correlati tra loro (vif <5)

Analisi della Qualità del Modello Una volta ottenuto il modello finale, valuteremo la sua capacità predittiva utilizzando metriche come R² e il Root Mean Squared Error (RMSE). Un’attenzione particolare sarà rivolta all’analisi dei residui e alla presenza di valori influenti, che potrebbero distorcere le previsioni, indagando su di essi.

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

Residuals vs fitted plot: il grafico dei residui mostra un buon addattamento generale per quanto riguarda la linearità, osserviamo una leggera eteroschedasticità. QQplot dei residui: il QQplot dei residui standardizzati mostra una distribuzione approssivatiamente normale, ad eccezione delle code (pesanti) e outlier che non seguono la distribuzione normale. si potrebbero “tagliare” le code o eliminare gli outlier per ridurre la distorsione e avvicinarsi a una distribuzione gaussiana. Scale-location: la radice dei residui è relativamente costante lungo i fitted values ad eccezione di un leggero aumento nella parte destra del plot. questo suggerisce e conferma la presenza di una lieve eteroschedasticità. Residuals vs laverage plot: il grafico mostra la presenza di un valore (1551) altamente influente (sia outlier che laverage point), questo può distorcere le stime della regressione.

Analisi outliers e high laverage points

#laverage point
lev <- hatvalues(mod3)
plot(lev)
p <- sum(lev)
soglia <- 2*p/n
abline(h=soglia,col=2)

plot(rstudent(mod3))
abline(h=c(-2,2))

outlierTest(mod3)

##      rstudent unadjusted p-value Bonferroni p
## 1551 9.986149         4.7193e-23   1.1798e-19
## 155  4.951276         7.8654e-07   1.9663e-03
## 1306 4.781188         1.8440e-06   4.6100e-03

cook <- cooks.distance(mod3)
plot(cook)

max(cook)

## [1] 0.9783883

l’osservazione 1551 è sia un outlier che un high laverage point. proverò a rimuoverla e ristimare il modello mod3 per osservare eventuali differenze.

new_data <- data[-1551,]
suppressWarnings(attach(new_data))

## The following objects are masked from data:
## 
##     Anni.madre, Cranio, Fumatrici, Gestazione, Lunghezza, N.gravidanze,
##     Ospedale, Peso, Sesso, Tipo.parto

new_mod3 <- lm(Peso~ Sesso + Cranio + Lunghezza + Gestazione, data= new_data)
summary(new_mod3)

## 
## Call:
## lm(formula = Peso ~ Sesso + Cranio + Lunghezza + Gestazione, 
##     data = new_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1153.79  -181.99   -14.94   163.85  1391.36 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6651.7241   132.9131 -50.046  < 2e-16 ***
## SessoM         79.3075    10.9963   7.212 7.27e-13 ***
## Cranio         10.0591     0.4208  23.906  < 2e-16 ***
## Lunghezza      10.8440     0.3018  35.935  < 2e-16 ***
## Gestazione     28.4861     3.7233   7.651 2.83e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.7 on 2494 degrees of freedom
## Multiple R-squared:  0.7362, Adjusted R-squared:  0.7358 
## F-statistic:  1740 on 4 and 2494 DF,  p-value: < 2.2e-16

summary(mod3)

## 
## Call:
## lm(formula = Peso ~ Sesso + Cranio + Lunghezza + Gestazione, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.2  -184.3   -17.6   163.3  2627.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6651.1188   135.5172 -49.080  < 2e-16 ***
## SessoM         79.1049    11.2117   7.056 2.22e-12 ***
## Cranio         10.6704     0.4245  25.139  < 2e-16 ***
## Lunghezza      10.2054     0.3007  33.939  < 2e-16 ***
## Gestazione     31.2737     3.7856   8.261 2.31e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2495 degrees of freedom
## Multiple R-squared:  0.7261, Adjusted R-squared:  0.7257 
## F-statistic:  1654 on 4 and 2495 DF,  p-value: < 2.2e-16

Osservando l’Rquadro, possiamo notare un leggero aumento dovuto alla rimozione dell’outlier e laverage point che provocava distorsione nella stima.

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

dai grafici non si riescono a osservare delle differenze sostanziali rimuovendo l’osservazione 1551.

ulteriori test

bptest(new_mod3)

## 
##  studentized Breusch-Pagan test
## 
## data:  new_mod3
## BP = 9.3632, df = 4, p-value = 0.05263

bptest(mod3)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod3
## BP = 89.148, df = 4, p-value < 2.2e-16

osservando il risultato del test Breusch-Pagan, possiamo notare come la rimozione dell’osservazione 1551 ha inciso sulla presenza o meno di eteroschedasticità: nel mod3 abbiamo un p-value al di sotto della soglia critica di 0.05 evidenziando la presenza di eteroschedasticità; il new_mod3 presenta un p-value poco al di sopra della soglia critica (0.05263), un dato completamente diverso rispetto al modello precedente, suggerendo che non vi sia presenza di eteroschedasticità.

dwtest(new_mod3)

## 
##  Durbin-Watson test
## 
## data:  new_mod3
## DW = 1.956, p-value = 0.1354
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(mod3)

## 
##  Durbin-Watson test
## 
## data:  mod3
## DW = 1.9557, p-value = 0.1337
## alternative hypothesis: true autocorrelation is greater than 0

il test Durbin-Watson, mostra valori simili per entrambi i modelli, sottolineando che non ci sia autocorrelazione positiva nei residui.

shapiro.test(residuals(new_mod3))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(new_mod3)
## W = 0.98864, p-value = 3.278e-13

shapiro.test(residuals(mod3))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod3)
## W = 0.9742, p-value < 2.2e-16

osservando i risultati dello shapiro-wilk test, possiamo evidenziare come entrambi i modelli non soddisfano l’assunzione di normalità: il modello new_mod3 però presenta una migliore forma dei residui (valore di W) rispetto al modello mod3.

graficamente troviamo conferma:

plot(density(residuals(new_mod3)))

plot(density(residuals(mod3)))

ulteriore analisi dei modelli osservando il mean square error, bias e la varianza

set.seed(100)
training=sample(1:nrow(new_data),size=nrow(new_data)*(2/3))
train_set=new_data[training,]
test_set=new_data[-training,]

evaluate_model <- function(model, test_data, response_var) {
  predictions <- predict(model, newdata = test_data)
  true_values <- test_data[[response_var]]
  mse <- round(mean((true_values - predictions)^2),2)
  bias <- round(mean(predictions - true_values),2)
  variance <- round(var(predictions),2)
  return(list(MSE = mse, Bias = bias, Variance = variance))
}
modelli <- list(
  modello_completo = modello_completo,
  modello_BIC_step = modello_BIC_step,
  mod3 = mod3,
  new_mod3 = new_mod3,
  mod4= mod4,
  mod5= mod5,
  mod6=mod6,
  mod7=mod7)

risultati <- data.frame(sapply(modelli, evaluate_model, test_data = test_set, response_var = "Peso"))
rownames(risultati) <- c("MSE","Bias","Variance")
kable(risultati)

	modello_completo	modello_BIC_step	mod3	new_mod3	mod4	mod5	mod6	mod7
MSE	67699.15	67962.04	68128.92	67995.47	68235.2	68921.68	68188.5	69364.34
Bias	9.64	10.64	11.1	9.97	11.16	12.12	11.31	11.59
Variance	192559.8	192587.5	191783.3	195396.6	191840.7	188929	191644.8	191389.8

confrontando il MSE dei modelli, si osserva che il new_mod3, pur essendo il modello con il minor numero di predittori, presenta uno dei valori più bassi del MSE, confermando quanto detto precedentemente sul compromesso tra semplicità e accuratezza. osservando il BIAS possiamo notare come tutti siano positivi, indicando una tendenza a sovrastimare leggermente il peso. il confronto tra i valori della varianza mostrano come il mod5 sia il modello con meno dispersione dei dati attorno alla media

Previsioni e Risultati Una volta validato il modello, lo useremo per fare previsioni pratiche. Ad esempio, potremo stimare il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.

test <- data.frame(
  Gestazione= 39,
  N.gravidanze= 3,
  Sesso= "F"
)
modello_predict <- lm(Peso~N.gravidanze+Sesso+Gestazione, data=data)
summary(modello_predict)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Sesso + Gestazione, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1494.15  -272.82   -14.43   267.31  1892.39 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3136.835    175.410 -17.883  < 2e-16 ***
## N.gravidanze    23.395      6.509   3.595 0.000331 ***
## SessoM         165.144     16.729   9.872  < 2e-16 ***
## Gestazione     162.025      4.499  36.014  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 414.2 on 2496 degrees of freedom
## Multiple R-squared:  0.3784, Adjusted R-squared:  0.3777 
## F-statistic: 506.5 on 3 and 2496 DF,  p-value: < 2.2e-16

predict(modello_predict,test)

##        1 
## 3252.311

il modello ha predetto un Peso per la neonata alla 39esima settimana di gestazione pari a: 3252g

Visualizzazioni Infine, utilizzeremo grafici e rappresentazioni visive per comunicare i risultati del modello e mostrare le relazioni più significative tra le variabili. Ad esempio, potremmo visualizzare l’impatto del numero di settimane di gestazione e del fumo sul peso previsto.

data_gestazione <- expand.grid(
  Gestazione = seq(min(new_data$Gestazione), max(new_data$Gestazione), length = 100),
  Sesso = c("F", "M"),
  Lunghezza = mean(new_data$Lunghezza),
  Cranio = mean(new_data$Cranio)
)

data_gestazione$Pred <- predict(new_mod3, newdata = data_gestazione)

ggplot(data_gestazione, aes(x = Gestazione, y = Pred, color = Sesso)) +
  geom_line(size = 1.2) +
  theme_minimal() +
  labs(
    title = "Peso in funzione della Gestazione per Sesso",
    x = "Settimane di Gestazione",
    y = "Peso (g)"
  )

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

possiamo osservare come i neonati abbiano un peso maggiore rispetto alle neonate lungo le settiamen di gestazione analizzate. nello specifico, si evidenzia una crescita lineare e simile durante le settimane per entrambi i sessi.

ggplot(data=new_data)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col=Fumatrici),position="jitter")+
  geom_smooth(aes(x=Gestazione, 
                  y=Peso,
                  col=Fumatrici),se=F,method="lm")+
  geom_smooth(aes(x=Gestazione, 
                  y=Peso),col="black",se=F,method="lm")

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

possiamo osservare come la gestazione sia una delle variabili principali nell’influenza del peso del neonato. ai figli delle madri fumatrici è associato a un minor peso rispetto a quelli delle madri non fumatrici.

progetto statistica inferenziale p_AI

Marco Falchi

2025-11-18