Neonatal Health Solutions

## Visualizza le prime righe del dataset
head(neonati)

## # A tibble: 6 × 10
##   Anni.madre N.gravidanze Fumatrici Gestazione  Peso Lunghezza Cranio Tipo.parto
##        <dbl>        <dbl>     <dbl>      <dbl> <dbl>     <dbl>  <dbl> <chr>     
## 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       
## # ℹ 2 more variables: Ospedale <chr>, Sesso <chr>

# Definisco le colonne numeriche e categoriche
numeric_columns <- c("Anni.madre","Gestazione","Peso","Lunghezza","Cranio")
categorical_columns <- c("N.gravidanze", "Fumatrici", "Tipo.parto", "Ospedale")

“Anni.madre” è variabile quantitativa continua.

“N.gravidanze” è variabile quantitativa discreta.

“Gestazione” è variabile quantitativa continua.

“Peso” è variabile quantitativa continua

“Lunghezza” è variabile quantitativa continua

“Cranio” è variabile quantitativa continua

“Fumatrici” è variabile categoriale seppur numerica

“Tipo.parto” è variabile categoriale

“Ospedale” è variabile categoriale

# Caricare i pacchetti necessari
library(knitr)
library(kableExtra)

# Calcola la frequenza delle categorie per le colonne categoriche
frequenze_categoriche <- lapply(neonati[categorical_columns], table)

# Calcolare le frequenze per ciascuna colonna categorica e combinare i risultati
lista_frequenze <- lapply(names(neonati[categorical_columns]), function(col_name) {
  freq_table <- table(neonati[[col_name]])
  freq_df <- as.data.frame(freq_table)
  colnames(freq_df) <- c("Categoria", "Frequenza")
  freq_df$Colonna <- col_name
  return(freq_df)
})

# Combina tutti i data frame in uno solo
frequenze_combinate <- do.call(rbind, lista_frequenze)

# Visualizza la tabella combinata con knitr::kable
library(knitr)
library(kableExtra)

kable(frequenze_combinate, format = "markdown", caption = "Frequenze per le Colonne Categoriche") %>%
  kable_styling(full_width = FALSE)

Frequenze per le Colonne Categoriche

Categoria	Frequenza	Colonna
0	1096	N.gravidanze
1	818	N.gravidanze
2	340	N.gravidanze
3	150	N.gravidanze
4	48	N.gravidanze
5	21	N.gravidanze
6	11	N.gravidanze
7	1	N.gravidanze
8	8	N.gravidanze
9	2	N.gravidanze
10	3	N.gravidanze
11	1	N.gravidanze
12	1	N.gravidanze
0	2396	Fumatrici
1	104	Fumatrici
Ces	728	Tipo.parto
Nat	1772	Tipo.parto
osp1	816	Ospedale
osp2	849	Ospedale
osp3	835	Ospedale

# Calcola gli indici di posizione per le colonne numeriche
indici_di_posizione_numerici <- lapply(neonati[numeric_columns], function(x) {
  if (all(is.na(x))) return(NA)
  list(
    media = mean(x, na.rm = TRUE),
    mediana = median(x, na.rm = TRUE),
    minimo = min(x, na.rm = TRUE),
    massimo = max(x, na.rm = TRUE),
    varianza = var(x, na.rm = TRUE),
    deviazione_standard = sd(x, na.rm = TRUE)
  )
})

# Converti la lista in un dataframe
indici_numerici_df <- do.call(rbind, lapply(indici_di_posizione_numerici, as.data.frame))
rownames(indici_numerici_df) <- names(indici_di_posizione_numerici)

library(knitr)

# Arrotonda tutte le colonne numeriche a 2 cifre decimali
indici_numerici_df <- indici_numerici_df %>%
  mutate(across(where(is.numeric), ~ round(.x, 2)))


# Visualizza la tabella con knitr::kable
kable(indici_numerici_df, caption = "Indici di Posizione per le Colonne Numeriche", format = "markdown")

Indici di Posizione per le Colonne Numeriche

	media	mediana	minimo	massimo	varianza	deviazione_standard
Anni.madre	28.16	28	0	46	27.81	5.27
Gestazione	38.98	39	25	43	3.49	1.87
Peso	3284.08	3300	830	4930	275665.68	525.04
Lunghezza	494.69	500	310	565	692.67	26.32
Cranio	340.03	340	235	390	269.79	16.43

# Calcola indici di variabilità per le colonne numeriche
calcola_variabilita <- function(x) {
  if (all(is.na(x))) return(NA)
  return (list(
    range_val = range(x,na.rm = TRUE),
    varianza = var(x,na.rm = TRUE),
    deviazione_standard = sd(x,na.rm = TRUE),
    coeff_variazione = (sd(x,na.rm = TRUE) / mean(x,na.rm = TRUE)) * 100,
    iqr = IQR(x,na.rm = TRUE)
  ))
}

# Converti la lista in un dataframe
# Applica la funzione alle colonne filtrate e converte in dataframe
indici_variabilita <- lapply(neonati[numeric_columns], calcola_variabilita)
indici_variabilita_df <- do.call(rbind, lapply(indici_variabilita, as.data.frame))
rownames(indici_variabilita_df) <- names(calcola_variabilita)

# Arrotonda tutte le colonne numeriche a 2 cifre decimali
indici_variabilita_df <- indici_variabilita_df %>%
  mutate(across(where(is.numeric), ~ round(.x, 2)))

# Visualizza la tabella con knitr::kable
kable(indici_variabilita_df, caption = "Indici di variabilita per le Colonne Numeriche", format = "markdown")

Indici di variabilita per le Colonne Numeriche

range_val	varianza	deviazione_standard	coeff_variazione	iqr
0	27.81	5.27	18.72	7
46	27.81	5.27	18.72	7
25	3.49	1.87	4.79	2
43	3.49	1.87	4.79	2
830	275665.68	525.04	15.99	630
4930	275665.68	525.04	15.99	630
310	692.67	26.32	5.32	30
565	692.67	26.32	5.32	30
235	269.79	16.43	4.83	20
390	269.79	16.43	4.83	20

“Anni.madre” presenta valori anomali come si evince dal valore minimo uguale a 0.

“Gestazione” presenta valori plausibili anche se con molti outliers.

“Peso” presenta valori plausibili ma con molti Outliers e in generale con Dev standard e Coeff variabilità elevati.

“Lunghezza” e “Cranio” hanno valori plausibili, anch’essi con molti Outliers.

Di seguito i grafici della variabilità delle colonne numeriche:

# Funzione per creare un grafico boxplot
calcola_grafico <- function(x, x_label, y_label) {
  x <- na.omit(x)
  df <- data.frame(Index = seq_along(x), Value = sort(x))
  quantiles <- quantile(x, na.rm = TRUE)
  quantile_labels <- data.frame(
    x = rep(max(df$Index), length(quantiles)),
    y = quantiles,
    label = round(quantiles, 2)
  )
  ggplot(df, aes(x = factor(1), y = Value)) +
    geom_boxplot(width = 0.5) +
    geom_text(data = quantile_labels, aes(x = x, y = y, label = label),
              color = "red", hjust = -0.5, vjust = -0.5) +
    ggtitle("Grafico di Variabilità") +
    xlab(x_label) +
    ylab(y_label) +
    theme_minimal() +
    theme(axis.text.x = element_blank()) +
    coord_cartesian(xlim = c(0.5, 1.5))
}

# Crea e mostra i grafici per le colonne numeriche
grafici <- lapply(numeric_columns, function(column_name) {
  x <- neonati[[column_name]]
  if (all(is.na(x))) return(NULL)
  calcola_grafico(x, paste("Index of", column_name), column_name)
})
grafici <- grafici[!sapply(grafici, is.null)]
grid.arrange(grobs = grafici, ncol = 2)

plot(neonati$Anni.madre, neonati$Peso)

Tramite il grafico sopra possiamo notare come i valori anomali siano sicuramente minori del valore 10(preso per semplicità, dato che è impossibile una mamma abbia meno di 10 anni).

Adesso si vuole identificare e sostituire i valori anomali in Anni.madre

# Crea una tabella con i valori di Anni.madre
anni_madre_tabella <- data.frame(Anni.madre = neonati$Anni.madre)

I valori Anni.madre vengono identificati di seguito, mostrati e rimpiazzati con la media degli anni delle altre madri del dataset. Le righe con i valori anomali hanno dei valori validi.

# Utilizza subset per filtrare il dataset e ottenere i valori anomali
anomalies <- subset(neonati, Anni.madre < 10)

# Stampa le righe anomale
pander(anomalies)

Table continues below

Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio
1	1	0	41	3250	490	350
0	0	0	39	3060	490	330

Tipo.parto	Ospedale	Sesso
Nat	osp2	F
Nat	osp3	M

# Calcola la media dei valori validi di 'Anni.madre'
mean_age <- mean(neonati$Anni.madre[neonati$Anni.madre >= 10], na.rm = TRUE)

# Sostituisci i valori anomali con la media calcolata
neonati$Anni.madre[neonati$Anni.madre < 10] <- mean_age

# Utilizza subset per filtrare il dataset e controllare se ci sono ancora anomalie
check_anomalies <- subset(neonati, Anni.madre < 10)

# Controlla se 'check_anomalies' non è NULL e ha righe
if (!is.null(check_anomalies) && nrow(check_anomalies) > 0) {
  print(check_anomalies)
} else {
  print("Non ci sono valori anomali!")
}

## [1] "Non ci sono valori anomali!"

Si vuole testare la veridicità di alcune tesi: Ipotesi1:“Non è vero che le madri fumatrici abbiano un neonato con peso nettamente diverso madri fumatrici”

#Traccia Precedente
t_test_1 <- t.test(Peso ~ Fumatrici, data = neonati)
pander(t_test_1)

Welch Two Sample t-test: Peso by Fumatrici (continued below)

Test statistic	df	P value	Alternative hypothesis	mean in group 0
1.034	114.1	0.3033	two.sided	3286

mean in group 1
3236

Il test t conferma l’ipotesi nulla. Non esiste una differenza rilevante tra peso neonati nati da madre fumatrice e non.

Ipotesi2: “Non è vero che le madre con gravidanze precedenti abbiano un neonati con peso diverso da quello delle madri di prima gravidanza”

#Si analizza il legame tra madri con gravidanze precedenti e  peso
t_test_2 <- t.test(Peso ~ N.gravidanze>0, data = neonati)
pander(t_test_2)

Welch Two Sample t-test: Peso by N.gravidanze > 0 (continued below)

Test statistic	df	P value	Alternative hypothesis	mean in group FALSE
-0.4306	2376	0.6668	two.sided	3279

mean in group TRUE
3288

L’ipotesi nulla è confermata. Non esiste una differenza rilevante tra peso neonati nati da madre con precedenti gravidanze e non.

Ipotesi3: Peso e età della madre sono correlati.

# Calcola la correlazione tra l'età della madre e il peso del neonato
correlazione <- cor.test(neonati$Anni.madre, neonati$Peso)

# Visualizza i risultati
pander::pander(correlazione)

Pearson’s product-moment correlation: neonati$Anni.madre and neonati$Peso Il test di correlazione indica un quasi nullo legame tra anni madre e peso neonato. L’ipotesi è falsa.

Test statistic	df	P value	Alternative hypothesis	cor
-1.189	2498	0.2346	two.sided	-0.02378

Ipotesi4 :Esiste una differenza di numero parti cesario in base all’ospedale?

Verifichiamolo tramite il test del chi-quadro

#modifichiamo il dataset codificando le variabili categoriche
library(dplyr)

neonati <- neonati %>%
  mutate(
    Tipo.Parto_Naturale = ifelse(Tipo.parto == "Nat", 1, 0),
    Tipo.Parto_Cesario = ifelse(Tipo.parto == "Ces", 1, 0),
    Sesso_M = ifelse(Sesso == "M", 1, 0),
    Sesso_F = ifelse(Sesso == "F", 1, 0),
    Ospedale = case_when(
      Ospedale == "osp1" ~ 1,
      Ospedale == "osp2" ~ 2,
      Ospedale == "osp3" ~ 3,
      TRUE ~ NA_real_
    )
  ) %>%
  dplyr::select(-Tipo.parto, -Sesso)

# Creare una tabella di contingenza
tabella_2 <- table(neonati$Tipo.Parto_Cesario, neonati$Ospedale)
# Eseguire il test del chi-quadro di indipendenza
test_chi_quadro <- chisq.test(tabella_2)
# Visualizzare i risultati
pander(test_chi_quadro)

Pearson’s Chi-squared test: tabella_2 Il p-value è maggiore di 0.05, ciò significa che la tesi che ammette differenze significative di parti cesari in base all’ospedale è falsa.

Test statistic	df	P value
1.097	2	0.5778

test_df = data.frame("Media_Peso_M"=3400 , "Media_Peso_F"= 3200 ,"Lunghezza_Media_M"= 50.5, "Lunghezza_Media_F"=49.5 )

Ipotesi5: La media peso del campione analizzato(il nostro dataset) è uguale a quella della popolazione(la popolazione generale).

Ho: La media del campione è uguale alla media della popolazione H1: La media del campione è diversa dalla media della popolazione

Il test sarà svolto separatamente tra uomini e donne:

#Eseguire il test t per confrontare il peso del campione con la media della popolazione

t_test_peso_m <- t.test(neonati$Peso[neonati$Sesso_M == 1], mu = test_df$Media_Peso_M)

p-value molto superiore a 0.05, (0.56), non si rifiuta H0, non ci sono differenze significative. La media peso maschile del campione è uguale a quella della popolazione.

t_test_peso_f <- t.test(neonati$Peso[neonati$Sesso_F == 1], mu = test_df$Media_Peso_F)

p-value inferiore a 0.05 (0.008), si rifiuta H0, ci sono differenze significative. La media peso femminile del campione non è uguale a quella della popolazione.

Ipotesi6: La media lunghezza neonati del campione analizzato(il nostro dataset) è uguale a quella della popolazione(la popolazione generale).

# Eseguire il test t per confrontare la lunghezza del campione con la media della popolazione
t_test_lunghezza_m <- t.test(neonati$Lunghezza[neonati$Sesso_M == 1], mu = test_df$Lunghezza_Media_M)

p-value inferiore a 0.05, si rifiuta H0, ci sono differenze significative La media lunghezza maschile del campione non è uguale a quella della popolazione.

t_test_lunghezza_f <- t.test(neonati$Lunghezza[neonati$Sesso_F == 1], mu = test_df$Lunghezza_Media_F)

p-value inferiore a 0.05, si rifiuta H0, ci sono differenze significative La media lunghezza femminile del campione non è uguale a quella della popolazione.

# Stampare i risultati del test t
print(t_test_peso_m)

## 
##  One Sample t-test
## 
## data:  neonati$Peso[neonati$Sesso_M == 1]
## t = 0.58679, df = 1243, p-value = 0.5574
## alternative hypothesis: true mean is not equal to 3400
## 95 percent confidence interval:
##  3380.748 3435.683
## sample estimates:
## mean of x 
##  3408.215

print(t_test_peso_f)

## 
##  One Sample t-test
## 
## data:  neonati$Peso[neonati$Sesso_F == 1]
## t = -2.6172, df = 1255, p-value = 0.008971
## alternative hypothesis: true mean is not equal to 3200
## 95 percent confidence interval:
##  3131.997 3190.267
## sample estimates:
## mean of x 
##  3161.132

print(t_test_lunghezza_m)

## 
##  One Sample t-test
## 
## data:  neonati$Lunghezza[neonati$Sesso_M == 1]
## t = 659.05, df = 1243, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 50.5
## 95 percent confidence interval:
##  498.3301 501.0043
## sample estimates:
## mean of x 
##  499.6672

print(t_test_lunghezza_f)

## 
##  One Sample t-test
## 
## data:  neonati$Lunghezza[neonati$Sesso_F == 1]
## t = 566.68, df = 1255, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 49.5
## 95 percent confidence interval:
##  488.2401 491.2885
## sample estimates:
## mean of x 
##  489.7643

#3.Le misure antropometriche sono significativamente diverse tra i due sessi
t_test_antro_l <- t.test(neonati$Lunghezza[neonati$Sesso_M == 1],neonati$Lunghezza[neonati$Sesso_F == 1])

Ipotesi nulla: La lunghezza è significativamente diversa tra i due sessi. Il p-value è molto minore di 0.05, non ci sono evidenze per confutare l’ip nulla, quindi l’ipotesi è vera.

t_test_antro_c <- t.test(neonati$Cranio[neonati$Sesso_M == 1],neonati$Cranio[neonati$Sesso_F == 1])

Ipotesi nulla: Il diametro del cranio è significativamente diversa tra i due sessi. Il p-value è molto minore di 0.05, non ci sono evidenze per confutare l’ip nulla.

print(t_test_antro_l)

## 
##  Welch Two Sample t-test
## 
## data:  neonati$Lunghezza[neonati$Sesso_M == 1] and neonati$Lunghezza[neonati$Sesso_F == 1]
## t = 9.582, df = 2459.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   7.876273 11.929470
## sample estimates:
## mean of x mean of y 
##  499.6672  489.7643

print(t_test_antro_c)

## 
##  Welch Two Sample t-test
## 
## data:  neonati$Cranio[neonati$Sesso_M == 1] and neonati$Cranio[neonati$Sesso_F == 1]
## t = 7.4102, df = 2491.4, p-value = 1.718e-13
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  3.541270 6.089912
## sample estimates:
## mean of x mean of y 
##  342.4486  337.6330

Quali sono le variabili che più influiscono sul peso di un neonato?

Matrice correlazione

correlation_matrix <- cor(neonati, use = "complete.obs")
ggcorrplot(correlation_matrix, method = "circle")

La matrice indica chiaramente come “Sesso_M”,“Gestazione”,“Lunghezza”, “Cranio” siano le variabili più influenti.

Creazione dei modelli di previsione

Si vuole provare a creare un modello con tutte le features:

#Creazione del Modello di Regressione
mod1 <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + Ospedale + Tipo.Parto_Naturale  + Sesso_M , data = neonati)
summary(mod1)

## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Ospedale + Tipo.Parto_Naturale + Sesso_M, 
##     data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1139.63  -182.97   -15.91   163.33  2617.32 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -6758.9675   141.7633 -47.678  < 2e-16 ***
## Anni.madre              0.7956     1.1472   0.693   0.4881    
## N.gravidanze           11.7043     4.6681   2.507   0.0122 *  
## Fumatrici             -30.5268    27.5598  -1.108   0.2681    
## Gestazione             32.6751     3.8201   8.553  < 2e-16 ***
## Lunghezza              10.2754     0.3008  34.161  < 2e-16 ***
## Cranio                 10.4852     0.4263  24.594  < 2e-16 ***
## Ospedale               14.1345     6.7536   2.093   0.0365 *  
## Tipo.Parto_Naturale    29.8162    12.0931   2.466   0.0137 *  
## Sesso_M                77.9335    11.1850   6.968 4.11e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.1 on 2490 degrees of freedom
## Multiple R-squared:  0.7284, Adjusted R-squared:  0.7274 
## F-statistic: 741.8 on 9 and 2490 DF,  p-value: < 2.2e-16

La bontà del modello è decente,ma ci sono tutte le variabili inseriti e come abbiamo visto soltanto alcune di queste sono correlate realmente al Peso, mentre altre non hanno alcun legame con la sua variazione. Si vuole quindi creare un modello più semplice, senza rinunciare al buon funzionamento del modello.

Le variabili più significative sono quelle inserite nel seguente modello

mod2 <- lm(Peso ~ Gestazione + Cranio + Lunghezza + Sesso_M , data = neonati)
summary(mod2)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Cranio + Lunghezza + Sesso_M, 
##     data = neonati)
## 
## 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 ***
## Gestazione     31.2737     3.7856   8.261 2.31e-16 ***
## Cranio         10.6704     0.4245  25.139  < 2e-16 ***
## Lunghezza      10.2054     0.3007  33.939  < 2e-16 ***
## Sesso_M        79.1049    11.2117   7.056 2.22e-12 ***
## ---
## 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

# Estrai i coefficienti
coeffs2 <- coef(mod2)
pander(coeffs2)

(Intercept)	Gestazione	Cranio	Lunghezza	Sesso_M
-6651	31.27	10.67	10.21	79.1

Variazione di y in funzione di x: Per ogni settimana aggiuntiva di gestazione (Gestazione), il peso previsto del neonato aumenta di 31.27 unità. Per ogni centimetro aggiuntivo di Lunghezza (Lunghezza), il peso previsto diminuisce di 10 unità. Probabilmente influenzato dalla combinazione di Lunghezza con Cranio. Per ogni unità di incremento nella misura del cranio (Cranio), il peso previsto aumenta di 10.67 unità. Se il neonato è maschio (Sesso_M), il peso previsto aumenta di 79 unità.

Multiple R squared 0.72, Adjusted R-squared 0.72, il modello sembra buono.

Si può notare come anche senza alcune variabili le prestazioni del modello non varia. Sicuramente questo modello è da preferire al precedente.

Si vuole provare un ulteriore miglioramento combinando Lunghezza e Cranio.

mod3 <- lm(Peso ~ Gestazione + Cranio + Lunghezza + Cranio * Lunghezza + Sesso_M, data = neonati )
summary(mod3)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Cranio + Lunghezza + Cranio * 
##     Lunghezza + Sesso_M, data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1139.05  -181.44   -15.46   163.32  2850.10 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.839e+03  1.019e+03  -1.804   0.0714 .  
## Gestazione        3.692e+01  3.951e+00   9.344  < 2e-16 ***
## Cranio           -4.409e+00  3.194e+00  -1.380   0.1676    
## Lunghezza        -2.013e-01  2.205e+00  -0.091   0.9273    
## Sesso_M           7.449e+01  1.121e+01   6.648 3.64e-11 ***
## Cranio:Lunghezza  3.113e-02  6.537e-03   4.763 2.01e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.8 on 2494 degrees of freedom
## Multiple R-squared:  0.7286, Adjusted R-squared:  0.728 
## F-statistic:  1339 on 5 and 2494 DF,  p-value: < 2.2e-16

Il miglioramento del modello quasi nullo rispetto al precedente.

Si vuole provare a inserire delle variabili quadratiche e studiarne la correlazione

#Creo una nuova Colonna Cranio^2
neonati$Cranio2 <- neonati$Cranio^2

#I valori di Cranio ^2 migliorano la relazione con Peso?

mod5 <- lm(Peso ~  Cranio2 , data = neonati)
summary(mod5) 

## 
## Call:
## lm(formula = Peso ~ Cranio2, data = neonati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2652.4  -240.8    -5.2   245.9  1424.8 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.925e+02  7.969e+01  -7.435 1.43e-13 ***
## Cranio2      3.345e-02  6.846e-04  48.864  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 375.5 on 2498 degrees of freedom
## Multiple R-squared:  0.4887, Adjusted R-squared:  0.4885 
## F-statistic:  2388 on 1 and 2498 DF,  p-value: < 2.2e-16

Non c’è una differenza significativa rispetto alla semplice relaz. Cranio->Peso. Proviamo un ulteriore combinazione di Cranio2 e il suo prodotto * Lunghezza

mod6 <- lm(Peso ~   Sesso_M + Cranio + Lunghezza + Cranio2 +Cranio2 * Lunghezza , data = neonati)
summary(mod6)

## 
## Call:
## lm(formula = Peso ~ Sesso_M + Cranio + Lunghezza + Cranio2 + 
##     Cranio2 * Lunghezza, data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1398.93  -186.25   -15.23   165.99  2640.14 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -1.220e+03  1.364e+03  -0.895  0.37106    
## Sesso_M            7.890e+01  1.137e+01   6.941 4.96e-12 ***
## Cranio            -2.922e+01  1.247e+01  -2.344  0.01918 *  
## Lunghezza          1.518e+01  2.074e+00   7.318 3.37e-13 ***
## Cranio2            7.471e-02  2.610e-02   2.863  0.00423 ** 
## Lunghezza:Cranio2 -3.072e-05  1.791e-05  -1.715  0.08642 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 278 on 2494 degrees of freedom
## Multiple R-squared:  0.7201, Adjusted R-squared:  0.7196 
## F-statistic:  1283 on 5 and 2494 DF,  p-value: < 2.2e-16

Si voleva provare ad aggiungere le variabili al quadrato per aumentare la bontà del modello, pur rispettando il principio di marginalità. Come si può vedere la bontà arriva addirittura a diminuire, nonostante si abbia un modello più completo.

neonati$Lunghezza2 <- neonati$Lunghezza^2
summary(neonati)

##    Anni.madre     N.gravidanze       Fumatrici        Gestazione   
##  Min.   :13.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.19   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       Ospedale    
##  Min.   : 830   Min.   :310.0   Min.   :235   Min.   :1.000  
##  1st Qu.:2990   1st Qu.:480.0   1st Qu.:330   1st Qu.:1.000  
##  Median :3300   Median :500.0   Median :340   Median :2.000  
##  Mean   :3284   Mean   :494.7   Mean   :340   Mean   :2.008  
##  3rd Qu.:3620   3rd Qu.:510.0   3rd Qu.:350   3rd Qu.:3.000  
##  Max.   :4930   Max.   :565.0   Max.   :390   Max.   :3.000  
##  Tipo.Parto_Naturale Tipo.Parto_Cesario    Sesso_M          Sesso_F      
##  Min.   :0.0000      Min.   :0.0000     Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.0000      1st Qu.:0.0000     1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :1.0000      Median :0.0000     Median :0.0000   Median :1.0000  
##  Mean   :0.7088      Mean   :0.2912     Mean   :0.4976   Mean   :0.5024  
##  3rd Qu.:1.0000      3rd Qu.:1.0000     3rd Qu.:1.0000   3rd Qu.:1.0000  
##  Max.   :1.0000      Max.   :1.0000     Max.   :1.0000   Max.   :1.0000  
##     Cranio2         Lunghezza2    
##  Min.   : 55225   Min.   : 96100  
##  1st Qu.:108900   1st Qu.:230400  
##  Median :115600   Median :250000  
##  Mean   :115890   Mean   :245413  
##  3rd Qu.:122500   3rd Qu.:260100  
##  Max.   :152100   Max.   :319225

Si vuole inserire sia Cranio ^2 che Lunghezza ^2 nel modello

mod7 <- lm(Peso ~ Sesso_M + Cranio + Lunghezza + Cranio2 + Lunghezza2 + Cranio * Lunghezza , data = neonati)
summary(mod7)

## 
## Call:
## lm(formula = Peso ~ Sesso_M + Cranio + Lunghezza + Cranio2 + 
##     Lunghezza2 + Cranio * Lunghezza, data = neonati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1660.0  -182.8   -12.5   164.1  1354.4 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -3.427e+03  1.136e+03  -3.018  0.00257 ** 
## Sesso_M           7.778e+01  1.123e+01   6.927 5.47e-12 ***
## Cranio            1.581e+01  9.466e+00   1.670  0.09509 .  
## Lunghezza        -3.779e+00  4.584e+00  -0.824  0.40982    
## Cranio2           4.374e-02  1.855e-02   2.358  0.01844 *  
## Lunghezza2        4.042e-02  5.014e-03   8.061 1.16e-15 ***
## Cranio:Lunghezza -7.017e-02  1.328e-02  -5.283 1.38e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2493 degrees of freedom
## Multiple R-squared:  0.7271, Adjusted R-squared:  0.7264 
## F-statistic:  1107 on 6 and 2493 DF,  p-value: < 2.2e-16

Il modello sembra non essere migliorato, con una Bontà pressochè identica ai modelli molto più semplici.

Proviamo a vedere se la relazione tra Lunghezza, Cranio rispetto a Peso varia se le features vengono elevate al quadrato.

plot(neonati$Peso, neonati$Lunghezza)

plot(neonati$Peso, neonati$Cranio)

plot(neonati$Peso,neonati$Lunghezza2)

plot(neonati$Peso, neonati$Cranio2)

Graficamente non si vedono cambiamenti di forma nella distribuzione tra variabile e variabile al quadrato nei confronti del Peso.

mod8 <- lm(Peso ~ Gestazione +  Sesso_M + Cranio + Cranio2 + Lunghezza + Lunghezza2 + Cranio*Lunghezza, data = neonati)
summary(mod8)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Sesso_M + Cranio + Cranio2 + 
##     Lunghezza + Lunghezza2 + Cranio * Lunghezza, data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1294.23  -182.55   -12.29   167.34  1298.24 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -6.624e+02  1.145e+03  -0.578  0.56313    
## Gestazione        3.972e+01  3.897e+00  10.192  < 2e-16 ***
## Sesso_M           7.311e+01  1.101e+01   6.638 3.89e-11 ***
## Cranio            3.525e+00  9.354e+00   0.377  0.70633    
## Cranio2           5.383e-02  1.820e-02   2.957  0.00313 ** 
## Lunghezza        -1.133e+01  4.553e+00  -2.489  0.01289 *  
## Lunghezza2        4.336e-02  4.922e-03   8.810  < 2e-16 ***
## Cranio:Lunghezza -6.002e-02  1.306e-02  -4.598 4.48e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.1 on 2492 degrees of freedom
## Multiple R-squared:  0.738,  Adjusted R-squared:  0.7373 
## F-statistic:  1003 on 7 and 2492 DF,  p-value: < 2.2e-16

# Estrai i coefficienti
coeffs8 <- coef(mod8)
pander(coeffs8)

Table continues below

(Intercept)	Gestazione	Sesso_M	Cranio	Cranio2	Lunghezza	Lunghezza2
-662.4	39.72	73.11	3.525	0.05383	-11.33	0.04336

Cranio:Lunghezza
-0.06002

Variazione di y in funzione di x: Per ogni settimana aggiuntiva di gestazione (Gestazione), il peso previsto del neonato aumenta di 39.72 unità. Per ogni centimetro aggiuntivo di lunghezza (Lunghezza), il peso previsto diminuisce di 11 unità. Probabilmente influenzato dalla combinazione di lunghezza con Cranio. Per ogni cm^2 di Cranio2 e lunghezza2 si ha un aumento del peso di 0.05. Per ogni unità di incremento nella misura del cranio (Cranio), il peso previsto aumenta di 3.52 unità. Se il neonato è maschio (Sesso_M), il peso previsto aumenta di 73 unità.

Analisi statistica

Verifichiamo innanzitutto la normalità dei dati.

shapiro.test(neonati$Peso)

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

moments::skewness(neonati$Peso)

## [1] -0.6470308

moments::kurtosis(neonati$Peso - 3)

## [1] 5.031532

Il test di Shapiro-Wilk rifiutà l’ipotesi di normalità. L’asimmetria(Skewness) indica una leggere pendenza verso sinistra. La curtosi indica una netta presenza di dati nelle code, indicando una forte presenza di valori estremi.

Proviamo a testare i due modelli( mod2 e mod8) tramite Il Criterio d’informazione di Akaike(AIC) e il Criterio di Informazione Bayesiano(BIC)

AIC(mod2,mod8)

##      df      AIC
## mod2  6 35185.60
## mod8  9 35080.56

BIC(mod2,mod8)

##      df      BIC
## mod2  6 35220.54
## mod8  9 35132.98

Il modello AIC e BIC suggerscono un modello 8 leggermente migliore.

Si vuole testare la collinearità delle varie features tramite il Variance Inflation Factor(vif), indicatore di collinearità.

install.packages("car")

## pacchetto 'car' aperto con successo con controllo somme MD5
## 
## I pacchetti binari scaricati sono in
##  C:\Users\Dario\AppData\Local\Temp\Rtmpq2W63S\downloaded_packages

library(car)

vif(mod2)

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

vif(mod8)

##       Gestazione          Sesso_M           Cranio          Cranio2 
##         1.829395         1.046824       814.499089      1376.745185 
##        Lunghezza       Lunghezza2 Cranio:Lunghezza 
##       495.342315       524.505004      1294.704633

Il test di collinearità del mod8 sfora chiaramente il range di collinearità, dal momento che le variabili al quadrato o combinate devono necessariamente essere correlate.

Il modello2 supera il test di collinearità.

Il miglior modello da utilizzare sarà quindi il mod2.

Il “mod8” è leggermente migliore , tuttavia la differenza rispetto al “mod2” in termini di R2 è minima, per il Principio di Parsimonia si vuole scegliere il “mod2” in quanto più semplice, nonostante la piccola differenza di Bontà.

Il riferimento d’ora in poi sarà quindi il “mod2”, si prova quindi a testarne i residui di seguito.

install.packages("MASS")

## pacchetto 'MASS' aperto con successo con controllo somme MD5

## 
## I pacchetti binari scaricati sono in
##  C:\Users\Dario\AppData\Local\Temp\Rtmpq2W63S\downloaded_packages

library(MASS)

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

Diagnostica grafica dei residui:

Grafico Residuals vs Fitted: Non ci sono pattern particolari seguiti da tanti. Q-Q Residuals: I residui seguono approssimativamente( a meno della coda finale) la direzione della prima bisettrice. Scale Location: Nuvola che accenna leggermente una curva, ma grosso modo è distribuita in modo casuale rispetto a una linea orizzontale. Residuals vs Leverage: Un solo valore è critico secondo la distanza di Cook, il 1551.

#Test sui Residui
install.packages("lmtest")

## pacchetto 'lmtest' aperto con successo con controllo somme MD5
## 
## I pacchetti binari scaricati sono in
##  C:\Users\Dario\AppData\Local\Temp\Rtmpq2W63S\downloaded_packages

library(lmtest)

Si utilizza il test di Shapiro-Wilk per verificare la normalità dei residui.

shapiro.test(residuals(mod2))#Controllo normalità

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

Shapiro-Wilk Test: Scopo: Questo test verifica l’ipotesi di normalità dei residui del modello. Uno degli assunti fondamentali della regressione lineare classica è che gli errori (residui) seguano una distribuzione normale. Se i residui non sono normalmente distribuiti, le stime dei parametri del modello possono essere distorte e i test di significatività potrebbero non essere validi. La normalità è particolarmente importante quando si desidera fare inferenze sui parametri del modello.

Il test dà esito negativo.

Si procede quindi col Test di Breusch-Pagan per verificare l’omoschedasticità.

bptest(mod2)#Test Breusch-Pagan(Omoschedasticità)

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

Scopo: Questo test verifica l’ipotesi di omoschedasticità, che significa che la varianza degli errori sia costante per tutti i valori delle variabili indipendenti. L’omoschedasticità è un altro assunto fondamentale della regressione lineare. Se i residui mostrano eteroschedasticità (varianza non costante), le stime dei coefficienti di regressione possono essere inefficaci, portando a errori standard errati e, di conseguenza, a test di significatività non affidabili. L’eteroschedasticità può anche indicare che il modello non ha catturato bene la struttura dei dati

Il test dà esito negativo.

Si vuole provare il test di Darwin-Watson per testare l’incorrelazione dei residui, qui di seguito.

dwtest(mod2)#Test Darwin-Watson(Incorrelazione)

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

Scopo: Questo test verifica l’ipotesi di incorrelazione dei residui, ovvero che gli errori non siano autocorrelati. L’autocorrelazione nei residui può indicare che ci sono pattern nel tempo o nelle osservazioni non catturati dal modello. Se si viola l’assunto di indipendenza dei residui, può portare a stime dei parametri distorte e a test di ipotesi non validi.

Si vuole mostrare una rappresentazione grafica della distribuzione dei residui:

#normalità dei residui grafica
plot(density(residuals(mod2)))

Pur non essendo perfettamente simmetrica la curva si avvicina molto a una Gaussiana, indicando una buona normalità.

Si vogliono adesso controllare i leverage values:

#Leverage
lev <- hatvalues(mod2)
plot(lev)

p = sum(lev)

n <- length(fitted(mod2)) 
soglia = 2*p/n
summary(lev[lev>soglia])

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.004006 0.004603 0.005559 0.007743 0.008183 0.048522

Nello spazio dei regressori ci sono diversi valori fuori norma, esattamente 143 valori leverage.

#Outliers
plot(rstudent(mod2))
abline(h = c(-2,2),col = 2)

#alternativa con corr Bonferroni
outlierTest(mod2)

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

Si hanno 3 valori Outliers.

#distanza di Cook
cook <- cooks.distance(mod2)
plot(cook)

max_cook <- max(cook)

Sono valori problematici per il modello?

Di seguito si vuole approfondire il valore 1551, il valore più estremo.

# Controlla che il dataset abbia almeno 1553 righe
if (nrow(neonati) >= 1551) {
  # Estrai la 1551ª riga
  osservazione_1551 <- neonati[1551, ]
  
  # Visualizza le caratteristiche della 1551ª osservazione
  View(osservazione_1551)
} else {
  print("Il dataset non ha 1551 righe.")
}

# Visualizzare i valori della distanza di Cook
summary(cook)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000000 0.0000270 0.0001312 0.0008970 0.0003974 0.9783883

# Trovare l'indice dell'osservazione con la distanza di Cook massima
index_max_cook <- which(cook == max_cook)

# Visualizzare il valore massimo e l'indice dell'osservazione
pander("Il valore massimo della distanza di Cook è:", max_cook, "per l'osservazione:", int(index_max_cook), "\n")

Il valore massimo della distanza di Cook è:

La distanza massima è di circa 1, per un solo valore: il 1551.

Proviamo a calcolare gli indici di posizione per i soli neonati femmina , e proviamo quanto il valore 1551 è simile a questi.

Si filtra il dataset per tenere solo le osservazioni con Sesso_F uguale a 1 e se ne calcolano gli indici di posizione.

neonati_femmine <- neonati[neonati$Sesso_F == 1, ]

# Calcola gli indici di posizione per le colonne numeriche nel dataset filtrato
indici_di_posizione_numerici_femmine <- lapply(neonati_femmine[numeric_columns], function(x) {
  if (all(is.na(x))) return(NA)
  list(
    media = mean(x, na.rm = TRUE),
    mediana = median(x, na.rm = TRUE),
    minimo = min(x, na.rm = TRUE),
    massimo = max(x, na.rm = TRUE),
    varianza = var(x, na.rm = TRUE),
    deviazione_standard = sd(x, na.rm = TRUE)
  )
})

# Converti la lista in un dataframe
indici_numerici_df_femmine <- do.call(rbind, lapply(indici_di_posizione_numerici_femmine, as.data.frame))
rownames(indici_numerici_df) <- names(indici_di_posizione_numerici_femmine)

library(knitr)

# Arrotonda tutte le colonne numeriche a 2 cifre decimali
indici_numerici_df_femmine <- indici_numerici_df_femmine %>%
  mutate(across(where(is.numeric), ~ round(.x, 2)))


# Visualizza la tabella con knitr::kable
kable(indici_numerici_df_femmine, caption = "Indici di Posizione per le Colonne Numeriche", format = "markdown")

Indici di Posizione per le Colonne Numeriche

	media	mediana	minimo	massimo	varianza	deviazione_standard
Anni.madre	28.14	28	13	44	26.28	5.13
Gestazione	38.73	39	25	43	3.96	1.99
Peso	3161.13	3160	830	4930	277001.31	526.31
Lunghezza	489.76	490	310	565	758.13	27.53
Cranio	337.63	340	235	390	280.15	16.74

L’osservazione 1551 ha un valore Peso molto maggiore dei valori media, settimane di Gestazione uguali alla media, Lunghezza molto minore(1551 ha valore lunghezza 315 contro una media di 490), Cranio leggermente superiore alla media. Ne risulta un neonato un po’ particolare: pesante , poco lungo e con la diametro del cranio notevole. Trattandosi di un solo valore anomalo, decido di tenerlo nel dataset.

Creiamo il dataset da utilizzare con le sole variabili del modello2 e verifichiamo un altro indicatore: il RMSE.

neonati_subset <- neonati[, c("Sesso_M", "Cranio", "Lunghezza", "Gestazione")]
summary(neonati_subset)

##     Sesso_M           Cranio      Lunghezza       Gestazione   
##  Min.   :0.0000   Min.   :235   Min.   :310.0   Min.   :25.00  
##  1st Qu.:0.0000   1st Qu.:330   1st Qu.:480.0   1st Qu.:38.00  
##  Median :0.0000   Median :340   Median :500.0   Median :39.00  
##  Mean   :0.4976   Mean   :340   Mean   :494.7   Mean   :38.98  
##  3rd Qu.:1.0000   3rd Qu.:350   3rd Qu.:510.0   3rd Qu.:40.00  
##  Max.   :1.0000   Max.   :390   Max.   :565.0   Max.   :43.00

# Ottieni le predizioni del modello
predicted_values <- predict(mod2, newdata = neonati_subset)
View(predicted_values)

# Calcola l'MSE
observed_values <- neonati$Peso
summary(observed_values)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     830    2990    3300    3284    3620    4930

mse <- mean((observed_values - predicted_values)^2)

# Calcola l'RMSE
rmse <- sqrt(mse)

# Stampa l'RMSE
cat("Il Root Mean Squared Error (RMSE) è:", rmse, "\n")

## Il Root Mean Squared Error (RMSE) è: 274.728

summary(mod2)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Cranio + Lunghezza + Sesso_M, 
##     data = neonati)
## 
## 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 ***
## Gestazione     31.2737     3.7856   8.261 2.31e-16 ***
## Cranio         10.6704     0.4245  25.139  < 2e-16 ***
## Lunghezza      10.2054     0.3007  33.939  < 2e-16 ***
## Sesso_M        79.1049    11.2117   7.056 2.22e-12 ***
## ---
## 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

Il RMSE è accettabile(271) in quanto minore della dev. standard(525).

CONCLUSIONE:

Analizzando i residui in maniera grafica il modello darebbe un esito positivo, anche se nelle code si accennano delle irregolarità dovute ad outliers, in particolare il valore 1551 supera la distanza di Cook. Il RMSE dà un risultato accettabile. Il test di Shapiro-Wilk e il Bp-test non indicano una correttezza del modello. Ciò è probabilmente dovuto alle irregolarità sopra citate. Nonostante questo, darò per buono il modello basandomi sulla rappresentazione grafica dei residui.

Prendiamo come ottimale il mod2.

3. Previsioni e Risultati

Si vuole fare un esempio pratico per prevedere il peso di un neonato, con dati non presenti nel dataset di addestramento del modello.

Si consideri una madre di 39 anni alla terza gravidanza e alla 35 settimana di gestazione. Sappiamo che il numero di gravidanze non ha nessun effetto sul modello, quindi possiamo ignorarne le info. L’unico dato utile è dato dalle settimane di gestazione.

Creiamo quindi un neonato con 35 settimane di gestazione, maschio ,con diametro Cranio e Lunghezza medio(media del dataset).

# Calcola la media per le variabili Cranio e Lunghezza
mean_cranio <- mean(neonati$Cranio, na.rm = TRUE)
mean_lunghezza <- mean(neonati$Lunghezza, na.rm = TRUE)

# Crea una nuova riga con la media dei valori e aggiungila al DataFrame
new_row <- data.frame(Sesso_M = 1,  # Supponiamo che '0' indica 'Femmina'
                      Cranio = mean_cranio,
                      Lunghezza = mean_lunghezza,
                      Gestazione = 39)

# Calcola la previsione del peso per la nuova riga
new_prediction <- predict(mod2, newdata = new_row)

# Converti la previsione in intero
new_prediction_int <- as.integer(new_prediction)

# Visualizza il risultato
pander(new_prediction_int)

3324 Il neonato ipotizzato avrà un peso di 3324 grammi , secondolo le previsioni del modello2.

Rappresenterò adesso i dati più significativi graficamente.

plot(neonati$Gestazione, neonati$Peso)

neonati_subset$Peso <- neonati$Peso

Graficamente si direbbe che all’aumentare del numero di settimane di Gestazione aumenta anche il Peso dei Neonati.

library(ggplot2) 

# install.packages("viridis")
install.packages("viridis")

## pacchetto 'viridis' aperto con successo con controllo somme MD5
## 
## I pacchetti binari scaricati sono in
##  C:\Users\Dario\AppData\Local\Temp\Rtmpq2W63S\downloaded_packages

library(viridis)

Il seguente grafico mostra la relazione tra Lunghezza, Cranio e Peso.

ggplot(neonati_subset, aes(x = Cranio, y = Lunghezza, color = Peso)) +
  geom_point(size = 4) +
  scale_color_viridis(option = "D") +  # Utilizza una scala di colori viridis
  labs(title = "Relazione tra Cranio, Lunghezza e Peso",
       x = "Cranio ",
       y = "Lunghezza ",
       color = "Peso (kg)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom")

All’aumentare delle dimensioni di Cranio e Lunghezza aumenta anche il Peso.

Visualizziamo la relazione tra Cranio, Gestazione E Peso Ottenuto.

ggplot(neonati_subset, aes(x = Cranio, y = Gestazione, color = Peso)) +
  geom_point(size = 4) +
  scale_color_viridis(option = "B") +  # Utilizza una scala di colori viridis
  labs(title = "Relazione tra Cranio, Gestazione e Peso",
       x = "Cranio (migliaia di cm)",
       y = "Settimane di Gestazione",
       color = "Peso (kg)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom")

All’aumentare del diametro del Cranio aumenta il Peso del neonato, in maniera meno evidente aumenta anche il numero di settimane di Gestazione.

Visualizziamo la relazione tra Lunghezza, Gestazione E Peso Ottenuto.

ggplot(neonati_subset, aes(x = Lunghezza, y = Gestazione, color = Peso)) +
  geom_point(size = 4) +
  scale_color_viridis(option = "C") +  # Utilizza una scala di colori viridis
  labs(title = "Relazione tra Lunghezza, Gestazione e Peso",
       x = "Lunghezza (migliaia di cm)",
       y = "Settimane di Gestazione",
       color = "Peso (kg)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom")

All’aumentare della lunghezza è correlato un aumento di peso , vale lo stesso all’aumentare delle settimane di Gestazione.

Il seguente grafico vuole mostrare come all’aumentare delle settimane di Gestazione è plausibile un aumento delle dimensioni in lunghezza e diametro del cranio del bambino.

ggplot(neonati_subset, aes(x = Cranio, y = Lunghezza, color = Gestazione)) +
  geom_point(size = 4) +
  scale_color_viridis(option = "A") +  # Utilizza una scala di colori viridis
  labs(title = "Relazione tra Cranio, Lunghezza e Sett Gestazione",
       x = "Cranio (migliaia di cm)",
       y = "Lunghezza (migliaia di cm)",
       color = "Gestazione (settimane)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom" )

Si può confermare che all’aumentare delle settimane di gestazione aumentano lunghezza e peso. ________________________________________________________________________________

Il seguente grafico vuole mostrare come il sesso maschile sia correlato a un aumento della lunghezza e un aumento del peso previsto.

ggplot(neonati_subset, aes(x = Lunghezza , y = Sesso_M, color = Peso)) +
  geom_point(size = 4) +
  scale_color_viridis(option = "A") +  # Utilizza una scala di colori viridis
  labs(title = "Relazione tra Lunghezza, Sesso_M e Peso_Prev",
       x = "Lunghezza (migliaia di cm)",
       y = "Sesso(0 femmina, 1 maschio)",
       color = "Peso previsto") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom" )

Il grafico non evidenzia nettamente la differenza di Peso tra neonati di Sesso M e Sesso F. Anche se il sesso è una features rilevante per il calcolo del peso, in questo grafico non è esplicito.

Si vuole dare una rappresentazione tridimensionale dei dati.

install.packages("plotly")

## pacchetto 'plotly' aperto con successo con controllo somme MD5
## 
## I pacchetti binari scaricati sono in
##  C:\Users\Dario\AppData\Local\Temp\Rtmpq2W63S\downloaded_packages

library(plotly)
#Si vuole dare una rappresentazione tridimensionale dei dati.

plot_ly(data = neonati, x = ~Gestazione, y = ~Lunghezza, z = ~Cranio, color = ~Peso, colors = c('blue', 'red')) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'Settimane di Gestazione'),
                      yaxis = list(title = 'Lunghezza in cm '),
                      zaxis = list(title = 'Cranio in cm ')))

All’aumentare di Cranio, Lunghezza e Gestazione aumenta anche il Peso.

Il seguente grafico presenta delle linee di smussamento per evidenziare le relazioni tra le variabili, diminuendo il rumore.

# Carica il pacchetto ggplot2
library(ggplot2)

# Crea il grafico con geom_smooth per diverse categorie
ggplot(neonati, aes(x = Lunghezza, y = Cranio, color = Peso)) +
  geom_point(size = 3)+ 
  geom_smooth(method = "lm", se = FALSE) +  # Aggiunge una linea di smussamento per ogni categoria
  labs(title = "Andamento Peso neonato al variare delle sue dimensioni",
       x = "Lunghezza",
       y = "Diametro Cranio") +
  theme_minimal()+
  theme(legend.position = "bottom")  

Si ha una retta di regressione abbastanza netta.

Neonatal Health Solutions adesso ha un modello valido di regressione per il peso dei neonati.