Modello Statistico per la Previsione del Peso Neonatale

#Importiamo le librerie
library(knitr)
library(moments)
library(ggplot2)
library(car)

## Caricamento del pacchetto richiesto: carData

library(lmtest)

## Caricamento del pacchetto richiesto: zoo

## 
## Caricamento pacchetto: 'zoo'

## I seguenti oggetti sono mascherati da 'package:base':
## 
##     as.Date, as.Date.numeric

library(dplyr)

## 
## Caricamento pacchetto: 'dplyr'

## Il seguente oggetto è mascherato da 'package:car':
## 
##     recode

## I seguenti oggetti sono mascherati da 'package:stats':
## 
##     filter, lag

## I seguenti oggetti sono mascherati da 'package:base':
## 
##     intersect, setdiff, setequal, union

1. Raccolta dei Dati e Struttura del Dataset:

Per costruire il modello predittivo, sono stati raccolti dati su 2500 neonati provenienti da tre ospedali.

#Importiamo il dataset
data<-read.csv("neonati.csv",sep=",", stringsAsFactors = T)
kable(data[1:5, ], caption = "Dataset - Tabella (primi 5 elementi)", digits=3, align="r")

Dataset - Tabella (primi 5 elementi)
Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
26	0	42	3380	490	325	Nat	osp3	M
21	2	39	3150	490	345	Nat	osp1	F
34	3	38	3640	500	375	Nat	osp2	M
28	1	41	3690	515	365	Nat	osp2	M
20	0	38	3700	480	335	Nat	osp3	F

E la struttura del dataset si presenta così:

str(data)

## 'data.frame':    2500 obs. of  10 variables:
##  $ Anni.madre  : int  26 21 34 28 20 32 26 25 22 23 ...
##  $ N.gravidanze: int  0 2 3 1 0 0 1 0 1 0 ...
##  $ Fumatrici   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Gestazione  : int  42 39 38 41 38 40 39 40 40 41 ...
##  $ Peso        : int  3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
##  $ Lunghezza   : int  490 490 500 515 480 495 480 510 500 510 ...
##  $ Cranio      : int  325 345 375 365 335 340 345 349 335 362 ...
##  $ Tipo.parto  : Factor w/ 2 levels "Ces","Nat": 2 2 2 2 2 2 2 2 1 1 ...
##  $ Ospedale    : Factor w/ 3 levels "osp1","osp2",..: 3 1 2 2 3 2 3 1 2 2 ...
##  $ Sesso       : Factor w/ 2 levels "F","M": 2 1 2 2 1 1 1 2 1 1 ...

L’obbiettivo principale è identificare quali delle variabili del dataset sono più predittive del peso alla nascita, con un particolare focus sull’impatto del fumo materno e delle settimane di gestazione.

2. Analisi e Modellizzazione:

Analisi Preliminare:

Analisi delle variabili del dataset “neonati”:

Età della madre (anni): [Anni.madre] variabile quantitativa continua
Numero di gravidanze: [N.gravidanze] variabile quantitativa discreta
Fumo materno: [Fumatrici] variabile qualitativa nominale dicotomica
Durata della gravidanza (settimane): [Gestazione] variabile quantitativa continua
Peso del neonato (g): [Peso] variabile quantitativa continua
Lunghezza del neonato (mm): [Lunghezza] variabile quantitativa continua
Diametro cranio (mm): [Cranio] variabile quantitativa continua
Tipo di parto: [Tipo.parto] variabile qualitativa nominale
Ospedale di nascita: [Ospedale] variabile qualitativa nominale
Sesso del neonato: [Sesso] variabile qualitativa nominale dicotomica

La variabile Peso sarà la variabile risposta del nostro modello mentre le altre variabili sono le variabili esplicative.

Analisi Descrittiva:

Trasformo la variabile Fumatrici e rinomino i campi, per renderla più parlante:

data$Fumatrici<-as.factor(data$Fumatrici)
levels(data$Fumatrici)<-c("Non fumatrici","Fumatrici")

Distribuzione di frequenze:

Iniziamo la nostra analisi descrittiva esaminando le distribuzioni di frequenze, relative alle variabili qualitative: Fumatrici, Tipo.parto, Ospedale e Sesso, riportate nelle seguenti tabelle:

Distribuzione di frequenze per Fumatrici:

n<-nrow(data)
attach(data)
distr_freq_fumatrici<-as.data.frame(
  cbind(
    ni=table(Fumatrici),
    fi=table(Fumatrici)/n
  ))

Distribuzione frequenze Fumatrici - Tabella
Paziente	ni	fi
Non fumatrici	2396	0.958
Fumatrici	104	0.042

Considerando il campione di riferimento, circa il 96% delle madri risulta essere “Non fumatrice”, mentre solo il 4% è fumatrice.

Distribuzione di frequenze per Tipo.parto:

distr_freq_parto<-as.data.frame(
  cbind(
    ni=table(Tipo.parto),
    fi=table(Tipo.parto)/n
  ))

Distribuzione frequenze Tipo di parto - Tabella
Paziente	ni	fi
Ces	728	0.291
Nat	1772	0.709

Per quanto riguarda il tipo di parto (cesareo o naturale), risulta che, nel campione considerato, il 70% delle pazienti hanno avuto un parto naturale, mentre il 30% ha avuto un parto cesareo.

Distribuzione di frequenze per Ospedale:

distr_freq_ospedale<-as.data.frame(
  cbind(
    ni=table(Ospedale),
    fi=table(Ospedale)/n
  ))

Distribuzione frequenze Ospedale - Tabella
Paziente	ni	fi
osp1	816	0.326
osp2	849	0.340
osp3	835	0.334

Per quanto riguarda l’ospedale, possiamo vedere che il campione è omogeneo. Le pazienti sono state selezionate sui tre ospedali con più o meno la stessa frequenza.

Distribuzione di frequenze per Sesso:

distr_freq_sesso<-as.data.frame(
  cbind(
    ni=table(Sesso),
    fi=table(Sesso)/n
  ))

Distribuzione frequenze Sesso - Tabella
Paziente	ni	fi
F	1256	0.502
M	1244	0.498

Considerando il campione di riferimento, si ha un numero uguale di neonati per ambo i sessi.

Indici di posizione:

Nelle tabelle seguenti, invece, vengono esaminati gli indici di posizione delle variabili quantitative:

attach(data)

## I seguenti oggetti sono mascherati da data (pos = 3):
## 
##     Anni.madre, Cranio, Fumatrici, Gestazione, Lunghezza, N.gravidanze,
##     Ospedale, Peso, Sesso, Tipo.parto

position_index<-as.data.frame(cbind(summary(Anni.madre),
                                    summary(N.gravidanze),
                                    summary(Gestazione),
                                    summary(Peso),
                                    summary(Lunghezza),
                                    summary(Cranio)))

Indici di posizione - Tabella
Indici	Anni madre (età)	N° gravidanze (unità)	Gestazione (settimane)	Peso (grammi)	Lunghezza (mm)	Cranio (mm)
Min.	0.00	0.00	25.00	830.00	310.00	235.00
1st Qu.	25.00	0.00	38.00	2990.00	480.00	330.00
Median	28.00	1.00	39.00	3300.00	500.00	340.00
Mean	28.16	0.98	38.98	3284.08	494.69	340.03
3rd Qu.	32.00	1.00	40.00	3620.00	510.00	350.00
Max.	46.00	12.00	43.00	4930.00	565.00	390.00

Guardando gli indici della variabile Anni.madre si nota la presenza di valori anomali (Min. 0 anni), quindi controllo graficamente:

plot(Anni.madre,pch=20,col="lightblue")

Anche graficamente risultano due valori anomali, che elimino scegliendo come età minima della madre un valore arbitrario di 13 anni e correggo gli indici di posizione:

data[data$Anni.madre < 13, "Anni.madre"]

## [1] 1 0

data <-data[data$Anni.madre >= 13, ]

attach(data)

## I seguenti oggetti sono mascherati da data (pos = 3):
## 
##     Anni.madre, Cranio, Fumatrici, Gestazione, Lunghezza, N.gravidanze,
##     Ospedale, Peso, Sesso, Tipo.parto

## I seguenti oggetti sono mascherati da data (pos = 4):
## 
##     Anni.madre, Cranio, Fumatrici, Gestazione, Lunghezza, N.gravidanze,
##     Ospedale, Peso, Sesso, Tipo.parto

position_index<-as.data.frame(cbind(summary(Anni.madre),
                                    summary(N.gravidanze),
                                    summary(Gestazione),
                                    summary(Peso),
                                    summary(Lunghezza),
                                    summary(Cranio)))

Indici di posizione corretti - Tabella
Indici	Anni madre (anni)	N° gravidanze (unità)	Gestazione (settimane)	Peso (grammi)	Lunghezza (mm)	Cranio (mm)
Min.	13.00	0.00	25.00	830.00	310.0	235.00
1st Qu.	25.00	0.00	38.00	2990.00	480.0	330.00
Median	28.00	1.00	39.00	3300.00	500.0	340.00
Mean	28.19	0.98	38.98	3284.18	494.7	340.03
3rd Qu.	32.00	1.00	40.00	3620.00	510.0	350.00
Max.	46.00	12.00	43.00	4930.00	565.0	390.00

Indici di variabilità:

#Calcolo indici di variabilità e di forma
cv<-function(x) {
  return(sd(x)/mean(x)*100)
}
var_name<-c("Anni.madre","N.gravidanze","Gestazione","Peso","Lunghezza","Cranio")
var_index <- data.frame(
  Indici= c("R","IQR","sigma2","sigma","CV [%]"))
index_form<-data.frame(
  Indici=c("skewness","kurtosis-3"))

for (varx in var_name) {
  var_index[[varx]] <- c(
    formatC(diff(range(data[[varx]])), format = "f", digits = 2),
    formatC(IQR(data[[varx]]), format = "f", digits = 2),
    formatC(var(data[[varx]]), format = "f", digits = 2),
    formatC(sd(data[[varx]]), format = "f", digits = 2),
    formatC(cv(data[[varx]]), format = "f", digits = 2)
  )
  
  index_form[[varx]]<-c(
    formatC(skewness(data[[varx]]),format = "f",digits=3),
    formatC(kurtosis(data[[varx]])-3,format = "f",digits=3)
  )
}

Indici di variabilità - Tabella
Indici	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio
R	33.00	12.00	18.00	4100.00	255.00	155.00
IQR	7.00	1.00	2.00	630.00	30.00	20.00
sigma2	27.22	1.64	3.49	275865.90	693.21	269.93
sigma	5.22	1.28	1.87	525.23	26.33	16.43
CV [%]	18.51	130.50	4.79	15.99	5.32	4.83

Distribuzione di frequenze di N.gravidanze:

La variabile del numero di gravidanze mostra un valore di CV molto alto, del 130.5%. E’ una variabile che varia da 0 a 12. Indaghiamo graficamente la sua distribuzione di frequenze:

n<-nrow(data)
distr_freq_gravidanze<-as.data.frame(
  cbind(
    ni=table(N.gravidanze),
    fi=table(N.gravidanze)/n
  ))  
barplot(distr_freq_gravidanze$fi,
        main="Distribuzione di frequenze di N.gravidanze",
        xlab="N° gravidanze",
        ylab="Frequenze relative",
        ylim=c(0,1),
        col="yellow",
        names.arg=c("0","1","2","3","4","5","6","7","8","9","10","11","12"))

Indici di forma:

Indici di forma - Tabella
Indici	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio
skewness	0.151	2.513	-2.065	-0.647	-1.515	-0.785
kurtosis-3	-0.106	10.982	8.256	2.029	6.481	2.945

Considerazioni:

Per quanto riguarda le variabili quantitative e gli indici osservati, si può osservare che:

Anni.madre: presenta un CV non trascurabile (18.51%), con il valore che varia da 13 a 46 anni (un range di 33 anni). L’età media delle gestanti si aggira intorno ai 28 anni (mediana 28 e media 28.19).Per quanto riguarda gli indici di forma presenta valori intorno allo 0, che fanno pensare ad una distrubuzione bene o male simmetrica (leggermente asimmetrica positiva) ed un’indice di curtosi<0 che fa pensare ad una distribuzione leggermente platicurtica.
N.gravidanze: dall’elevato CV (130.51%) si può evincere che ha una grandissima variabilità, di fatti è la variabile con il CV più alto. Il numero di gravidanze varia da 0 a 12. Come possiamo notare dalla distribuzione di frequenze, nel campione preso in esame, la maggior parte delle donne è alla prima gravidanza (circa il 44%), circa il 33% è alla seconda ed il 14% è alla terza. Quindi circa il 91% delle donne prese in esame ha avuto massimo due gravidanze. Aumentando il numero delle gravidanze abbiamo percentuali sempre più basse, con valori per 11 e 12 gravidanze di un solo soggetto ciascuno. Per quanto riguarda invece gli indici di forma, è possibile notare che presenta un’asimmetria positiva (con valori bassi più frequenti, come visto) ed un’indice di curtosi molto>0, che ci fa pensare ad una distribuzione leptocurtica.
Gestazione: il coefficiente di variabilità (4.8%) è piuttosto basso. Abbiamo un valore minimo di 25 settimane (metà del sesto mese), che indica un parto prematuro, fino a 43 settimane (gravidanza post-termine o tardiva). La media (39 settimane) e la mediana (38.98 settimane) si aggirano intorno a 40 settimane, che è il valore normale della gestazione umana. Per quanto riguarda invece gli indici di forma, è possibile notare che presenta un’asimmetria negativa (con valori elevati più frequenti) ed un’indice di curtosi>0 che ci fa pensare ad una distribuzione leptocurtica.
Lunghezza: il coefficiente di variabilità (5.32%) risulta piuttosto basso. Il valore varia da 310 cm (decisamente inferiore anche al primo quartile 480 cm) a 565 cm. Per quanto riguarda invece gli indici di forma, è possibile notare che presenta un indice di simmetria<0, che ci indica un’asimmetria negativa ed un’indice di curtosi>0 che ci fa pensare ad una distribuzione leptocurtica.
Cranio: il coefficiente di variabilità (4.83%) risulta piuttosto basso. Il valore varia da 235 cm a 390 cm. Per quanto riguarda invece gli indici di forma, è possibile notare che presenta un indice di simmetria<0, che ci indica un’asimmetria negativa ed un’indice di curtosi>0 che ci fa pensare ad una distribuzione leptocurtica.
Peso: presenta un CV non trascurabile (circa 16%), da cui è possibile evincere che c’è una variabilità. Il peso infatti varia da 830g, che indica un neonato sottopeso, con un peso ben inferiore alla mediana (3300 g) e alla media (3284.18) a 4930 g (quasi 5 kg).
Per quanto riguarda invece gli indici di forma, è possibile notare che presenta una leggere asimmetria negativa ed un’indice di curtosi>0 che ci fa pensare ad una distribuzione leptocurtica. Trattandosi della nostra variabile risposta andremo ad indagare meglio la sua distribuzione, effettuando uno Shapiro Test.

shapiro.test(Peso)

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

Il test Shapiro-Wilk sulla variabile Peso ci restituisce un p-value<2.2e-16, quindi un p-value molto piccolo e inferiore al livello di significatività alfa scelto a priori del 5%. Pertanto rifiutiamo l’ipotesi nulla di distribuzione normale della variabile Peso, come suggerito anche dai valori di Skewness e di Kurtosis. Ne vediamo anche la distribuzione graficamente:

plot(density(Peso),lwd=1.5, col="tomato3")

Viene confermato che la variabile presenta un’asimmetria negativa (con valori elevati più frequenti) e leptocurtica. Questo dato potrebbe essere dovuto alla presenza di outliers e si potrebbere riflettere sul comportamento dei residui una volta stimato il modello di regressione lineare, in quanto essendo un modello lineare esso non riesce bene a filtrare variabili non normali. Indagheremo meglio il comportamento dei residui una volta stimato il modello di regressione lineare.

Outliers:

Tramite l’utilizzo di Boxplot, andiamo intanto a indagare la presenza di outliers nelle variabili quantitative:

Anni.madre:

ggplot(data)+
  geom_boxplot(aes(y=Anni.madre),
                   color="black",
                   fill="lightblue")+
  labs(title="boxplot di Anni.madre",
       y="Età madre [anni]")+
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

Il boxplot della variabile Anni.madre evidenzia alcuni outliers, sia nella parte bassa che nella parte alta del grafico. Essi rappresentando valori inferiori e superiori ai limiti attesi. Abbiamo delle gravidanza in età superiore ai 42-43 anni e delle gravidanze dove l’età della madre è inferiore ai 15 anni.

N.gravidanze:

ggplot(data)+
  geom_boxplot(aes(y=N.gravidanze),
                   color="black",
                   fill="yellow")+
  labs(title="boxplot di N.gravidanze",
       y="Numero di gravidanze [unità]")+
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

Il boxplot della variabile N.gravidanze evidenzia alcuni outliers nella parte superiore del grafico. Essi ci confermano quanto già esaminato, ossia la rarità di un numero di gravidanze precedenti maggiore di due. Sono casi rari, ma potenzialmente influenti sull’analisi statistica.

Gestazione:

ggplot(data)+
  geom_boxplot(aes(y=Gestazione),
                   color="black",
                   fill="lightcoral")+
  labs(title="boxplot di Gestazione",
       y="Tempo di gestazione [settimane]")+
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

Il boxplot della variabile Gestazione evidenzia alcuni outliers nella parte inferiore del grafico. Considerando una gestazione media nell’intorno delle 40 settimane, i dati mettono in luce diversi casi inferiori alle 35 settimane, valori che indicherebbero la presenza di parti prematuri, con probabili conseguenze sul peso del neonato.

Lunghezza:

ggplot(data)+
  geom_boxplot(aes(y=Lunghezza),
                   color="black",
                   fill="mediumaquamarine")+
  labs(title="boxplot di Lunghezza",
       y="Lunghezza del neonato [cm]")+
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

Il boxplot della variabile Lunghezza evidenzia alcuni outliers nella parte superiore, ma soprattutto nella parte inferiore del grafico, anche se la maggioranza dei valori rientrano nella normalità. Questi dati potrebbero essere dati da errori nella registrazione o potrebbero riferirsi a parti prematuri.

Peso:

ggplot(data)+
  geom_boxplot(aes(y=Peso),
                   color="black",
                   fill="tomato3")+
  labs(title="boxplot di Peso",
       y="Peso del neonato [g]")+
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

Il boxplot della variabile Peso evidenzia alcuni outliers nella parte superiore, ma soprattutto nella parte inferiore del grafico. Questi dati potrebbero essere dati da errori nella registrazione o potrebbero riferirsi a parti prematuri. Questi dati risultano essere rilevanti per il nostro studio, poiché un basso peso alla nascita è un indicatore critico per la salute neonatale e richiede cure specifiche e l’uso di attrezzature specifiche, come incubatrici.

Cranio:

ggplot(data)+
  geom_boxplot(aes(y=Cranio),
                   color="black",
                   fill="orange2")+
  labs(title="boxplot di Cranio",
       y="Diametro cranio del neonato [cm]")+
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

Il boxplot della variabile Cranio evidenzia alcuni outliers nella parte superiore, ma soprattutto nella parte inferiore del grafico, anche se la maggioranza dei valori rientrano nella normalità. Questi dati potrebbero essere dati da errori nella registrazione o potrebbero riferirsi a parti prematuri.

Ipotesi:

A questo punto possiamo procedere a saggiare le seguenti ipotesi:

1. in alcuni ospedali si fanno più parti casarei

Per rispondere a questa ipotesi, bisogna indagare l’indipendenza tra le variabili qualitative Tipo.parto e Ospedale.

Avremo quindi il seguente sistema di ipotesi:

H0: le variabili sono indipendenti
H1: le variabili presentano associazione

e per saggiare l’indipendenza useremo il test chi_quadrato.

tabella<-table(Tipo.parto,Ospedale)
tabella_contingenza <- as.data.frame.matrix(tabella)

Tabella di contingenza tra Opedale e Tipo.parto
	osp1	osp2	osp3
Ces	242	254	232
Nat	574	594	602

test_indipendenza<-chisq.test(tabella_contingenza)
test_indipendenza

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

Eseguendo il test chi_quadro otteniamo un valore di p-value=0.5819, quindi un valore di p-value maggiore del valore di significatività α=0.05 per tanto non rifiutiamo l’ipotesi nulla H0 di indipendenza tra le due variabili. Non sembra quindi esserci una relazione tra l’ospedale ed il numero di parti cesarei eseguiti.

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

Per poter rispondere a questa domanda, sono stati cercati i dati inerenti alla media della popolazione sia per la lunghezza che per il peso, essendo in possesso dei soli dati inerenti al campione. Le informazioni sono state prese dal seguente link e si riferiscono ai dati dell’Ospedale Pediatrico Bambin Gesù. Dal sito ci viene detto che la lunghezza media per la popolazione italiana è di 50 cm mentre il peso medio è di 3300 g con un’oscillazione tra maschi e femmine di 150 g.

Per scegliere il test migliore per il confronto bisogna vedere la distribuzione per le variabili Peso e Lunghezza. Sappiamo già che il Peso non ha una distribuzione normale, andiamo a vedere la variabile Lunghezza eseguendo un test di Shapiro-Wilk:

shapiro.test(Lunghezza)

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

Anche in questo caso, con un p-value<2.2e-16 rifiutiamo l’ipotesi di normalità, come del resto avevamo intuito dai valori di Skewness e Kurtosis. La vediamo graficamente:

plot(density(Lunghezza),lwd=1.5,col="mediumaquamarine")

Useremo in entrambi i casi quindi il test non parametrico del t test, il test di Wilcoxon. Avremo quindi il seguente sistema di ipotesi:

H0: differenza tra le medie=0
H1: differenza tra le medie!=0

Lunghezza_Popolazione<-500
Peso_Popolazione<-3300
wilcox.test(Lunghezza,mu=Lunghezza_Popolazione)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  Lunghezza
## V = 875939, p-value < 2.2e-16
## alternative hypothesis: true location is not equal to 500

wilcox.test(Peso,mu=Peso_Popolazione)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  Peso
## V = 1493714, p-value = 0.9482
## alternative hypothesis: true location is not equal to 3300

Dal test di Wilcoxon possiamo evincere che per la Lunghezza, che ha un valore di p-value significativamente minore rispetto ad α=0.05, rifiutiamo l’ipotesi nulla H0 e possiamo ritenere il campione significativamente diverso rispetto alla popolazione. Nel caso, invece, del Peso, che ha un valore di p-value > α non rifiutiamo l’ipotesi nulla H0; quindi il campione non è significativamente diverso dalla popolazione.

3. Le misure antropometriche sono significativamente tra i due sessi

Le misure antropometriche sono: Lunghezza, Peso e Cranio. Anche nel caso della variabile Cranio vediamo se è distribuita come una normale:

shapiro.test(Cranio)

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

Con un p-value<2.2e-16, rifiutiamo l’ipotesi di normalità; essendo tutte e tre le variabili distribuite in modo non normale, anche in questo caso useremo il test di Wilcoxon.

Differenze in base al sesso del diametro del cranio:

Avremo quindi il seguente sistema di ipotesi:

H0: differenza tra le medie=0
H1: differenza tra le medie!=0

Valutiamo la differenza sia graficamente, tramite boxplot:

ggplot(data)+
  geom_boxplot(aes(x=Sesso,
                   y=Cranio,
                   fill=Sesso),col="black")+
  scale_fill_manual(values=c("violetred","royalblue"))+
  labs(title="Boxplot diametro cranio e sesso",
       x="Sesso",
       y="Diametro Cranio [mm]")+
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

E tramite il test di Wilcoxon:

wilcox.test(Cranio~Sesso)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Cranio by Sesso
## W = 639872, p-value = 7.141e-15
## alternative hypothesis: true location shift is not equal to 0

Il valore molto piccolo di p-value=7.141e-15, ci porta ad escludere l’ipotesi nulla. Tra maschi e femmine c’è una significativa differenza nelle dimensioni del cranio.

Differenze in base al sesso della Lunghezza:

Avremo quindi il seguente sistema di ipotesi:

H0: differenza tra le medie=0
H1: differenza tra le medie!=0

Valutiamo la differenza sia graficamente, tramite boxplot:

ggplot(data)+
  geom_boxplot(aes(x=Sesso,
                   y=Lunghezza,
                   fill=Sesso),col="black")+
  scale_fill_manual(values=c("violetred","royalblue"))+
  labs(title="Boxplot lunghezza e sesso",
       x="Sesso",
       y="Lunghezza [mm]")+
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

E tramite il test di Wilcoxon:

wilcox.test(Lunghezza~Sesso)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Lunghezza by Sesso
## W = 593412, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

Il valore molto piccolo di p-value<2.2e-16, ci porta ad escludere l’ipotesi nulla. Tra maschi e femmine c’è una significativa differenza nella lunghezza.

Differenze in base al sesso del Peso:

Avremo quindi il seguente sistema di ipotesi:

H0: differenza tra le medie=0
H1: differenza tra le medie!=0

Valutiamo la differenza sia graficamente, tramite boxplot:

ggplot(data)+
  geom_boxplot(aes(x=Sesso,
                   y=Peso,
                   fill=Sesso),col="black")+
  scale_fill_manual(values=c("violetred","royalblue"))+
  labs(title="Boxplot peso e sesso",
       x="Sesso",
       y="Peso [g]")+
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

E tramite il test di Wilcoxon:

wilcox.test(Peso~Sesso)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Peso by Sesso
## W = 537495, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

Il valore molto piccolo di p-value<2.2e-16, ci porta ad escludere l’ipotesi nulla. Tra maschi e femmine c’è una significativa differenza nel peso.

Creazione del Modello di Regressione:

Abbiamo già visto come la variabile risposta Peso non è distribuita normalmente. Andiamo a vedere la correlazione esistente tra la variabile risposta e i regressori e tra i regressori stessi, creando la Matrice di Correlazione:

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = 1.5)
}
pairs(data,upper.panel = panel.smooth, lower.panel = panel.cor)

Dalla matrice di correlazione possiamo vedere che la variabile risposta Peso presenta una correlazione con le variabili quantitative Gestazione, Lunghezza e Cranio. Ed abbiamo, nello specifico:

Peso~Gestazione: 0.59, valore alto e >0, confermato anche dallo scatterplot dove la nuvola sembra allinearsi e andare verso l’alto, pur presentando una non linearità per valori bassi di gestazione.

plot(Peso,Gestazione,pch=20, col="lightcoral")

Peso~Lunghezza: 0.80, valore alto e >0, confermato anche dallo scatterplot, che mostra una chiara tendenza, seppur sembra esserci non linearità.

plot(Peso,Lunghezza, pch=20, col="mediumaquamarine")

Peso~Cranio: 0.70, valore alto e >0, confermato anche qui dallo scatterplot. Potrebbe esserci anche qui una leggera non linearità.

plot(Peso, Cranio, pch=20, col="orange2")

Tra le tre variabili potrebbe esserci multicollinearità.

Per quanto riguarda le variabili Anni.madre e N.gravidanze, non sembrerebbe esserci un’evidente correlazione:

Peso~Anni.madre: 0.024, valore basso e >0, confermato da scatterplot dove la nuvola è sparsa con tendenza a linea orizzontale.
Peso~N.gravidanze: 0.0023, valore basso e >0, confermato da scatterplot dove la nuvola è sparsa con tendenza a linea orizzontale.

Variabili Qualitative:

Per quanto riguarda le variabili qualitative, sappiamo che la matrice di correlazione ci dà poche informazioni, quindi andiamo ad indagarle tramite altri metodi statistici.

Peso~Tipo.parto:

Andiamo a vedere graficamente tramite boxplot:

boxplot(Peso~Tipo.parto, main="Relazione tra peso e tipo di parto",col=c("wheat2","lightblue"))

wilcox.test(Peso~Tipo.parto)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Peso by Tipo.parto
## W = 633851, p-value = 0.5244
## alternative hypothesis: true location shift is not equal to 0

Dal boxplot appare evidente che il tipo di parto non influenzi il peso del neonato. Questo dato ci viene confermato anche dal wilcoxon test, dove il p-value=0.5244 è maggiore del livello di significatività α=0.05 quindi non rifiuto H0, ipotesi di uguaglianza. La variabile Tipo.parto non sembra essere significativamente rilevante, almeno per il nostro campione.

Peso~Sesso:

Andiamo a vedere graficamente tramite boxplot:

boxplot(Peso~Sesso, main="Relazione tra peso e sesso",col=c("violetred","royalblue"))

wilcox.test(Peso~Sesso)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Peso by Sesso
## W = 537495, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

Dal boxplot appare evidente che il sesso influenzi il peso del neonato. I maschi pesano di più rispetto alle femmine. Questo dato ci viene confermato anche dal wilcoxon test, dove il p-value<2.2e-16 è molto vicino a 0, quindi rifiuto H0, ipotesi di uguaglianza. La variabile Sesso sembra essere significativamente rilevante, per tanto mi aspetto un beta di regressione significativo.

Peso~Ospedale:

Andiamo a vedere graficamente tramite boxplot:

boxplot(Peso~Ospedale, main="Relazione tra peso e ospedale",col=c("gold","gold2","gold4"))

pairwise.wilcox.test(Peso,Ospedale,paired = FALSE, 
                     pool.sd = T,p.adjust.method = "bonferroni")

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  Peso and Ospedale 
## 
##      osp1 osp2
## osp2 1.00 -   
## osp3 0.71 0.26
## 
## P value adjustment method: bonferroni

Come già visto, i parti sono equamente distribuiti tra i 3 ospedali. Anche la matrice di confronto del test presenta p-value elevati, che ci portano ad accettare l’ipotesi di uguaglianza. Per tanto, la variabile non sembra essere statisticamente rilevante, almeno per il nostro campione.

Peso~Fumatrici:

Andiamo a vedere graficamente tramite boxplot:

boxplot(Peso~Fumatrici, main="Relazione tra il peso ed il fumo", col=c("powderblue","lightslategrey"))

wilcox.test(Peso~Fumatrici)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Peso by Fumatrici
## W = 138069, p-value = 0.05928
## alternative hypothesis: true location shift is not equal to 0

Dal boxplot appare evidente che il fumo influenzi il peso del neonato. I neonati figli di chi non è fumatrice pesano di più rispetto ai neonati di chi fuma. Nel test, però, il p-value è al limite, con un valore di 0.05928, che ci porterebbe a non rifiutare l’uguaglianza. Questo dato appare in disaccordo con la letteratura scientifica e potrebbe essere dovuto al campione, dove le gestanti fumatrici sono del 30%.

Creazione del Modello di Regressione:

MODELLO 1:

Dalle analisi effettuate, abbiamo già visto che le variabili Tipo.parto e Ospedale non influenzano il peso del neonato. Per tanto, facciamo un primo modello di regressione con tutte le variabili tranne queste due variabili.

mod1<-lm(Peso~. -(Tipo.parto+Ospedale), data= data )
summary(mod1)

## 
## Call:
## lm(formula = Peso ~ . - (Tipo.parto + Ospedale), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1160.6  -181.3   -15.7   163.6  2630.7 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -6712.2405   141.3339 -47.492  < 2e-16 ***
## Anni.madre             0.8803     1.1491   0.766    0.444    
## N.gravidanze          11.3789     4.6767   2.433    0.015 *  
## FumatriciFumatrici   -30.3958    27.6080  -1.101    0.271    
## Gestazione            32.9472     3.8288   8.605  < 2e-16 ***
## Lunghezza             10.2316     0.3011  33.979  < 2e-16 ***
## Cranio                10.5198     0.4271  24.633  < 2e-16 ***
## SessoM                78.0787    11.2132   6.963 4.24e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2490 degrees of freedom
## Multiple R-squared:  0.7272, Adjusted R-squared:  0.7264 
## F-statistic: 948.3 on 7 and 2490 DF,  p-value: < 2.2e-16

Come ci aspettavamo, le variabili Gestazione, Lunghezza, Cranio e Sesso hanno tutte valore di p-value vicino allo 0 e risultano essere variabili molto significative. Anche la variabile N.gravidanze sembra essere rilevante. All’aumentare del numero delle gravidanze, il peso del neonato sembra aumentare. Anni.madre con un valore di p-value=0.444, quindi maggiore del nostro α conferma di non essere una variabile significativa. La stessa cosa avviene per la variabile Fumatrice, che però è importante al fine dello studio e che probabilmente non risulta rilevante a causa del campionamento. Il valore di Adjusted R_squared è 0.7264 quindi valido.

MODELLO 2:

Creiamo un secondo modello, rimuovendo la variabile Anni.madre.

mod2<-update(mod1,~. -Anni.madre, data= data)
summary(mod2)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1150.24  -181.32   -15.73   162.95  2635.69 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -6682.2637   135.7983 -49.207  < 2e-16 ***
## N.gravidanze          12.6996     4.3470   2.921  0.00352 ** 
## FumatriciFumatrici   -30.5728    27.6048  -1.108  0.26818    
## Gestazione            32.6437     3.8079   8.573  < 2e-16 ***
## Lunghezza             10.2309     0.3011  33.979  < 2e-16 ***
## Cranio                10.5366     0.4265  24.707  < 2e-16 ***
## SessoM                78.1596    11.2117   6.971 4.01e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2491 degrees of freedom
## Multiple R-squared:  0.7271, Adjusted R-squared:  0.7265 
## F-statistic:  1106 on 6 and 2491 DF,  p-value: < 2.2e-16

Le variabili mantengono la loro significatività e l’Adjusted R_squared è di 0.7265 quindi resta invariato. Migliora la significatività di N.gravidanze. Il modello potrebbe essere valido.

MODELLO 3:

La variabile Fumatrici non è molto significativa. Proviamo a vedere l’effetto sul modello, eliminandola:

mod3<-update(mod2,~. -Fumatrici, data= data)
summary(mod3)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.37  -180.98   -15.57   163.69  2639.09 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.7251   135.8036 -49.201  < 2e-16 ***
## N.gravidanze    12.4554     4.3416   2.869  0.00415 ** 
## Gestazione      32.3827     3.8008   8.520  < 2e-16 ***
## Lunghezza       10.2455     0.3008  34.059  < 2e-16 ***
## Cranio          10.5410     0.4265  24.717  < 2e-16 ***
## SessoM          77.9807    11.2111   6.956 4.47e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1327 on 5 and 2492 DF,  p-value: < 2.2e-16

Le variabili mantengono la loro significatività e non abbiamo effetto sull’Adjusted R_squared. Essendo però importante ai fini dello studio ed essendo probabilmente poco influente a causa del campione, la teniamo.

MODELLO 4:

Abbiamo visto che, per questo campione, la variabile Fumatrici non è significativa. Potrebbe essere interessante indagare l’effetto di interazione tra Fumatrici e la variabile Gestazione:

mod4<-update(mod2,~. +Fumatrici*Gestazione, data= data)
summary(mod4)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso + Fumatrici:Gestazione, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.83  -181.58   -16.95   163.76  2635.68 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -6699.6947   136.7290 -49.000  < 2e-16 ***
## N.gravidanze                     12.7452     4.3470   2.932   0.0034 ** 
## FumatriciFumatrici              795.7016   757.5315   1.050   0.2936    
## Gestazione                       33.2007     3.8418   8.642  < 2e-16 ***
## Lunghezza                        10.2252     0.3011  33.957  < 2e-16 ***
## Cranio                           10.5313     0.4265  24.693  < 2e-16 ***
## SessoM                           78.7436    11.2241   7.016 2.94e-12 ***
## FumatriciFumatrici:Gestazione   -21.0469    19.2830  -1.091   0.2752    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2490 degrees of freedom
## Multiple R-squared:  0.7273, Adjusted R-squared:  0.7265 
## F-statistic: 948.6 on 7 and 2490 DF,  p-value: < 2.2e-16

Le variabili significative restano tali, la variabile Fumatrici resta poco significativa e anche l’interazione Fumatrici:Gestazione non risulta essere significativa. L’Adjusted R_squared non varia.

Dagli scatterplot è emerso l’effetto quadratico relativo alle variabili Gestazione, Cranio e Lunghezza. Andiamo ad indagare il loro effetto sui modelli:

MODELLO 5:

Indaghiamo l’effetto quadratico della variabile Gestazione:

mod5<-update(mod2,~. +I(Gestazione^2), data= data)
summary(mod5)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso + I(Gestazione^2), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1144.90  -181.16   -12.71   165.06  2658.43 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -4665.6268   898.7841  -5.191 2.26e-07 ***
## N.gravidanze          12.7825     4.3435   2.943  0.00328 ** 
## FumatriciFumatrici   -29.3545    27.5870  -1.064  0.28740    
## Gestazione           -79.9534    49.7533  -1.607  0.10818    
## Lunghezza             10.3353     0.3043  33.960  < 2e-16 ***
## Cranio                10.6326     0.4282  24.830  < 2e-16 ***
## SessoM                75.9482    11.2447   6.754 1.78e-11 ***
## I(Gestazione^2)        1.5031     0.6622   2.270  0.02331 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.5 on 2490 degrees of freedom
## Multiple R-squared:  0.7277, Adjusted R-squared:  0.7269 
## F-statistic: 950.7 on 7 and 2490 DF,  p-value: < 2.2e-16

Aggiungendo l’effetto quadratico della variabile Gestazione, vediamo che la variabile Gestazione perde di significatività e l’Adjusted R_squared non cambia in maniera significativa.

MODELLO 6:

Indaghiamo l’effetto quadratico della variabile Cranio:

mod6<-update(mod2,~. +I(Cranio^2), data= data)
summary(mod6)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso + I(Cranio^2), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1139.39  -178.60   -14.03   163.57  2619.67 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          56.88134 1152.12461   0.049  0.96063    
## N.gravidanze         12.97552    4.31813   3.005  0.00268 ** 
## FumatriciFumatrici  -26.69191   27.42789  -0.973  0.33056    
## Gestazione           39.10747    3.93843   9.930  < 2e-16 ***
## Lunghezza            10.47375    0.30190  34.692  < 2e-16 ***
## Cranio              -31.63012    7.17177  -4.410 1.08e-05 ***
## SessoM               73.27785   11.16748   6.562 6.45e-11 ***
## I(Cranio^2)           0.06238    0.01059   5.890 4.39e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 272.8 on 2490 degrees of freedom
## Multiple R-squared:  0.7309, Adjusted R-squared:  0.7301 
## F-statistic: 966.1 on 7 and 2490 DF,  p-value: < 2.2e-16

Aggiungendo l’effetto quadratico della variabile Cranio, vediamo che le variabili restano significative e l’Adjusted R_squared con un valore di 0.7301 migliora leggermente. Bisogna capire se, per il principio del resoio di Occam, è un miglioramento tale da aggiungere un’altra variabile.

MODELLO 7:

Indaghiamo l’effetto quadratico della variabile Lunghezza:

mod7<-update(mod2,~. +I(Lunghezza^2), data= data)
summary(mod7)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso + I(Lunghezza^2), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1170.26  -182.00   -11.21   163.18  1785.31 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        196.075358 724.115115   0.271 0.786583    
## N.gravidanze        14.273458   4.271642   3.341 0.000846 ***
## FumatriciFumatrici -24.005963  27.115144  -0.885 0.376062    
## Gestazione          42.733042   3.882223  11.007  < 2e-16 ***
## Lunghezza          -20.208568   3.163544  -6.388 2.00e-10 ***
## Cranio              10.648125   0.418933  25.417  < 2e-16 ***
## SessoM              70.127565  11.040734   6.352 2.52e-10 ***
## I(Lunghezza^2)       0.031583   0.003268   9.664  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.7 on 2490 degrees of freedom
## Multiple R-squared:  0.737,  Adjusted R-squared:  0.7363 
## F-statistic: 996.9 on 7 and 2490 DF,  p-value: < 2.2e-16

Aggiungendo l’effetto quadratico della variabile Lunghezza, vediamo che le variabili restano significative e l’Adjusted R_squared con un valore di 0.7363 migliora leggermente. Bisogna capire se, per il principio del resoio di Occam, è un miglioramento tale da aggiungere un’altra variabile.

MODELLO 8:

Indaghiamo l’effetto quadratico della variabile Gestazione e della variabile Lunghezza sul modello:

mod8<-update(mod7,~. + I(Gestazione^2), data = data)
summary(mod8)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso + I(Lunghezza^2) + I(Gestazione^2), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1191.66  -180.67   -13.04   165.68  1402.73 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -2.379e+03  9.060e+02  -2.625 0.008712 ** 
## N.gravidanze        1.466e+01  4.255e+00   3.446 0.000579 ***
## FumatriciFumatrici -2.464e+01  2.700e+01  -0.913 0.361527    
## Gestazione          3.366e+02  6.274e+01   5.365 8.82e-08 ***
## Lunghezza          -3.207e+01  4.039e+00  -7.941 3.02e-15 ***
## Cranio              1.044e+01  4.195e-01  24.898  < 2e-16 ***
## SessoM              7.276e+01  1.101e+01   6.610 4.69e-11 ***
## I(Lunghezza^2)      4.362e-02  4.143e-03  10.528  < 2e-16 ***
## I(Gestazione^2)    -3.872e+00  8.250e-01  -4.693 2.83e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 268.6 on 2489 degrees of freedom
## Multiple R-squared:  0.7393, Adjusted R-squared:  0.7385 
## F-statistic: 882.4 on 8 and 2489 DF,  p-value: < 2.2e-16

Le variabili rimangono significative, tranne la variabile Fumatrici. Tutte le variabili riducono di molto il valore del coefficiente di regressione, come vediamo dai valori di “Estimate” e per Lunghezza abbiamo un parametro negativo. In compenso abbiamo un miglioramento ulteriore di Adjusted R_squared che ha valore 0.7385.

Selezione del modello ottimale:

Andiamo a vedere il modello migliore tra gli otto modelli creati, utilizzando i criteri BIC e AIC:

AIC(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8)

##      df      AIC
## mod1  9 35155.07
## mod2  8 35153.66
## mod3  7 35152.89
## mod4  9 35154.46
## mod5  9 35150.49
## mod6  9 35121.09
## mod7  9 35063.67
## mod8 10 35043.66

BIC(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8)

##      df      BIC
## mod1  9 35207.48
## mod2  8 35200.24
## mod3  7 35193.65
## mod4  9 35206.87
## mod5  9 35202.90
## mod6  9 35173.50
## mod7  9 35116.08
## mod8 10 35101.90

Per entrambi i criteri, il modello migliore risulta essere il MODELLO 8, che considera gli effetti quadratici sia della Lunghezza che della Grandezza. Prendendo in considerazione il criterio BIC, che penalizza un numero troppo alto di variabili, possiamo notare che, tolti i modelli relativi agli effetti quadratici (MODELLI 6-7-8) ed escludendo il MODELLO 3, dove non consideravamo la variabile Fumatrici, il modello più semplice (con meno variabili) che risulta il più valido è il MODELLO 2. Usiamo a questo punto il test ANOVA per verificare se è presente un aumento significativo della varianza spiegata tra il MODELLO 2 ed il MODELLO 8.

anova(mod2,mod8)

## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Sesso
## Model 2: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Lunghezza^2) + I(Gestazione^2)
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1   2491 187949505                                  
## 2   2489 179565443  2   8384063 58.107 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Il valore del p-value<2.2e-16 ci dice che abbiamo un significativo aumento della varianza spiegata, quindi scelgo il MODELLO 8.

Controlliamo a questo punto la presenza di multicollinearità tra le variabili:

vif(mod8)

##    N.gravidanze       Fumatrici      Gestazione       Lunghezza          Cranio 
##        1.027979        1.007279      475.939678      391.486566        1.643729 
##           Sesso  I(Lunghezza^2) I(Gestazione^2) 
##        1.049062      373.045477      453.588126

I valori sono <5, quindi non abbiamo multicollinearità, tranne nel caso di Lunghezza e Gestazione, dove gli alti valori sono dovuti alla presenza dei loro effetti quadratici.

Analisi della Qualità del Modello:

Il modello scelto, ossia il MODELLO 8, presenta un valore di coefficiente di determinazione R_squared = 0.7393, mentre di Adjusted R_squared = 0.7385. Abbiamo quindi che quasi il 74% della variabilità osservata sul peso neonatale è spiegata dai regressori inclusi nel modello. Le variabili più significative, come già visto, sono:

N.gravidanze: ad ogni gravidanza precedente è associato un aumento medio di 14.66 g sul peso del neonato.
Gestazione: per ogni settimana avremo un incremento di 336.6 g, con un rallentamento nelle settimane avanzate (effetto quadratico significativo)
Lunghezza: Relazione curvilinea, con un impatto positivo che si riduce per lunghezze maggiori.
SessoM: I neonati maschi pesano mediamente 72.76 g in più rispetto alle femmine.
Cranio: Incremento medio di 10.44 g.

La variabile Fumatrici resta poco significativa, ma come detto questo potrebbe essere dovuto al campione ed essendo oggetto dello studio, la manteniamo.

La presenza di termini quadratici nel modello migliora la capacità di rappresentare relazioni non lineari, garantendo una descrizione più accurata delle dinamiche tra le variabili.

RMSE_mod8<-sqrt(mean(mod8$residuals^2))
RMSE_mod8

## [1] 268.1113

L’ RMSE ha un valore di 268.1113.

Analisi dei residui e presenza si outliers e leverages:

I residui devono essere puliti e rispettare le seguenti assunzioni:

media =0
varianza costante
distribuzione normale

E’ inoltre importante osservare la presenza di eventuali Outliers e Leverages, valori anomali che potrebbero influenzare le previsioni.

Andiamo a effettuare una prima analisi grafica:

par(mfrow=c(2,2))
plot(mod8,pch=20,col="limegreen")

Dai grafici precedenti è possibile effettuare alcune considerazioni:

Il primo grafico ci mostra che effettivamente la media è intorno allo 0.
Il secondo grafico, che riguarda la distribuzione normale dei residui, ci dice che la distribuzione dei residui è tendenzialmente normale nella parte centrale, mentre tende ad allontanarsi dalla distribuzione normale sulle due code.
Il terzo grafico ci mostra una distribuzione della nuvola abbastanza casuale, che ci farebbe pensare ad una varianza costante dei residui.
Il quarto grafico evidenzia la presenza di un valore leverage critico (1551).

Andiamo ad indagare i risultati grafici attraverso diversi test statistici:

bptest(mod8)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod8
## BP = 98.102, df = 8, p-value < 2.2e-16

dwtest(mod8)

## 
##  Durbin-Watson test
## 
## data:  mod8
## DW = 1.9496, p-value = 0.1041
## alternative hypothesis: true autocorrelation is greater than 0

shapiro.test(residuals(mod8))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod8)
## W = 0.98914, p-value = 7.704e-13

Il test Breusch-Pagan, che valuta l’omoschedasticità, ci restituisce un valore di p-value<2.2e-16. Per tanto l’assunzione di varianza costante non sembra essere rispettata. Detto ciò, unendo le considerazioni grafiche a quelle statistiche, possiamo notare che non ci sono pattern forti ed è probabile che l’eteroschedasticità sia minima e non comprometta la capacità predittiva del modello. Il test Durbin-Watson, che valuta la correlazione tra i residui, ha un valore di p-value=0.1041 maggiore di α=0.05, Per tanto possiamo dire che i residui non sono autocorrelati. Il test Shapiro-Wilk, che valuta la distribuzione normale, ci restituisce un valore di p-value=7.704e-13. Anche qui l’assunzione di distribuzione normale non sembra essere rispettata. Detto ciò, unendo le considerazioni grafiche a quelle statistiche, che mostrano delle anomalie sulle code, potremmo ritenerla accettabile. Il grafico della distribuzione dei residui si presenterà così:

plot(density(residuals(mod8)),lwd=1.5,col="limegreen")

Outliers e Leverages:

Indaghiamo i valori leverages:

lev<-hatvalues(mod8)
plot(lev,pch=20,col="limegreen")
p=sum(lev)
soglia=2*p/n
abline(h=soglia,col="red3", lwd=2)

length(lev[lev>soglia])

## [1] 213

Abbiamo 213 valori leverages.

Indaghiamo anche i valori outliers:

plot(rstudent(mod8),pch=20,col="limegreen")
abline(h=c(-2,2),col="red3", lwd=2)

outlierTest(mod8)

##       rstudent unadjusted p-value Bonferroni p
## 1551  6.063273         1.5369e-09   3.8391e-06
## 1306  4.953796         7.7664e-07   1.9400e-03
## 155   4.466320         8.3131e-06   2.0766e-02
## 1399 -4.459600         8.5758e-06   2.1422e-02
## 1694  4.302183         1.7565e-05   4.3879e-02

Abbiamo 5 valori outliers. Andiamo a vedere la distanza di Cook:

cook <- cooks.distance(mod8)
plot(cook,pch=20,col="limegreen")

max(cook)

## [1] 1.324176

Come visto anche nell’analisi grafica, abbiamo un’osservazione (n° 1551,che avendo eliminato due dati per N.madre diventa la n° 1549) la cui distanza di Cook supera la soglia di allarme, con una distanza di Cook di 1.3. Questo valore influenzerà il modello. A questo punto indaghiamo i dati corrispondenti a questo valore, per capire di che tipo di dati si tratta e se potrebbero essere frutto di una registrazione errata o se, effettivamente, potrebbe trattarsi di dati reali.

oss_1549<-data[1549, ]
oss_1549

##      Anni.madre N.gravidanze     Fumatrici Gestazione Peso Lunghezza Cranio
## 1551         35            1 Non fumatrici         38 4370       315    374
##      Tipo.parto Ospedale Sesso
## 1551        Nat     osp3     F

I dati ci dicono che abbiamo un neonato a fine gestazione (38 settimane), con un peso elevato e sopra la media (4370g). Facciamo ulteriori indagini:

data[-1549,]%>%
  filter(Gestazione>=38)%>%
  summarize(min(Lunghezza))

##   min(Lunghezza)
## 1            410

data[-1549,]%>%
  filter(Gestazione>=38)%>%
  summarize(mean(Peso))

##   mean(Peso)
## 1   3379.273

Considerando un valore di gestazione di 38 settimane, la lunghezza minima attesa sarebbe di 410 mm, quindi molto superiore rispetto al nostro dato. Mentre il peso medio è di 3380 g, molto inferiore rispetto a quello rilevato nei dati. Facendo altre indagini:

data[-1549,]%>%
  filter(Lunghezza>=310&Lunghezza<=320)%>%
  summarize(max(Gestazione))

##   max(Gestazione)
## 1              28

data[-1549,]%>%
  filter(Lunghezza>=310&Lunghezza<=320)%>%
  summarize(max(Peso))

##   max(Peso)
## 1       980

data[-1549,]%>%
  filter(Lunghezza>=310&Lunghezza<=320)%>%
  summarize(max(Cranio))

##   max(Cranio)
## 1         265

Possiamo notare che la lunghezza indicata (315 mm), ce l’aspetteremmo per una gestazione di 28 settimane. In relazione sempre alla lunghezza, il peso atteso ha un valore nettamente inferiore (980 g), rispetto a quello rilevato (4370 g), così come il diametro del cranio.

data[-1549,]%>%
  filter(Peso>=4300&Peso<=4400)%>%
  summarize(mean(Lunghezza))

##   mean(Lunghezza)
## 1             534

Anche la lunghezza media attesa (534 mm) per il peso rilevato risulta essere molto superiore, rispetto a quella registrata (315 mm). Per tanto, il dato risulta anomalo e potrebbe essere frutto di rilevazioni scorrette, non solo rispetto al campione, ma rispetto alla popolazione (lunghezza=500 mm e peso=3300) per tanto decido di escluderlo dal dataset.

data_new<-data[-1549,]
mod8_new<-update(mod7,~. + I(Gestazione^2), data = data_new)
summary(mod8_new)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso + I(Lunghezza^2) + I(Gestazione^2), data = data_new)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1188.3  -180.5   -11.3   163.5  1319.1 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -2.829e+03  9.026e+02  -3.134 0.001743 ** 
## N.gravidanze        1.457e+01  4.224e+00   3.448 0.000574 ***
## FumatriciFumatrici -2.379e+01  2.681e+01  -0.887 0.374956    
## Gestazione          2.005e+02  6.622e+01   3.028 0.002484 ** 
## Lunghezza          -1.922e+01  4.537e+00  -4.235 2.37e-05 ***
## Cranio              1.009e+01  4.205e-01  24.000  < 2e-16 ***
## SessoM              7.355e+01  1.093e+01   6.728 2.12e-11 ***
## I(Lunghezza^2)      3.085e-02  4.621e-03   6.676 3.03e-11 ***
## I(Gestazione^2)    -2.135e+00  8.679e-01  -2.460 0.013980 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 266.7 on 2488 degrees of freedom
## Multiple R-squared:  0.7427, Adjusted R-squared:  0.7418 
## F-statistic: 897.6 on 8 and 2488 DF,  p-value: < 2.2e-16

Eliminato il dato e guardando il summary, notiamo che seppur il valore di Adjusted R_squared sia aumentato (0.7418), le variabili Gestazione e il suo valore quadratico hanno perso significatività. Scarto il modello.

Avendo scartato il modello 8, riprendo in considerazione il secondo modello migliore secondo i valori di BIC, ossia il modello 7, che considerava il solo effetto quadratico della variabile Lunghezza. Consideriamo sempre il dataset con l’osservazione rimossa.

mod7_new<-update(mod2,~. +I(Lunghezza^2), data= data_new)
summary(mod7_new)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso + I(Lunghezza^2), data = data_new)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1177.33  -179.14   -10.48   163.31  1326.26 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -1.625e+03  7.591e+02  -2.141 0.032403 *  
## N.gravidanze        1.436e+01  4.228e+00   3.397 0.000691 ***
## FumatriciFumatrici -2.336e+01  2.684e+01  -0.871 0.384052    
## Gestazione          3.795e+01  3.898e+00   9.735  < 2e-16 ***
## Lunghezza          -1.166e+01  3.344e+00  -3.488 0.000495 ***
## Cranio              1.015e+01  4.203e-01  24.136  < 2e-16 ***
## SessoM              7.236e+01  1.093e+01   6.619 4.41e-11 ***
## I(Lunghezza^2)      2.323e-02  3.432e-03   6.768 1.62e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 267 on 2489 degrees of freedom
## Multiple R-squared:  0.7421, Adjusted R-squared:  0.7413 
## F-statistic:  1023 on 7 and 2489 DF,  p-value: < 2.2e-16

Il modello si semplifica, avendo rimosso un regressore, ossia l’effetto quadratico della variabile Gestazione. Tutti i regressori restano significativi, a meno della variabile Fumatrici, come già visto. Il valore di Adjusted R_squared rimane pressoché identico (0.7413).

Esaminiamo l’eventuale multicollinearità:

vif(mod7_new)

##   N.gravidanze      Fumatrici     Gestazione      Lunghezza         Cranio 
##       1.027603       1.007247       1.859536     266.527752       1.668233 
##          Sesso I(Lunghezza^2) 
##       1.046737     255.645668

Anche qui i valori restano <5, tranne per Lunghezza ed il suo effetto quadratico, quindi non abbiamo multicollinearità.

Vediamo l’RMSE del modello:

RMSE_mod7_new<-sqrt(mean(mod7_new$residuals^2))
RMSE_mod7_new

## [1] 266.5289

L’RMSE ha un valore di 266.5289.

Effettuiamo nuovamente l’analisi dei residui per il nuovo modello, partendo dall’analisi grafica:

par(mfrow=c(2,2))
plot(mod7_new, pch=20,col="limegreen")

Dai grafici precedenti è possibile notare che:

Il primo grafico ci mostra che effettivamente la media è intorno allo 0.
Il secondo grafico, che riguarda la distribuzione normale dei residui, ci dice che la distribuzione dei residui è tendenzialmente normale nella parte centrale, mentre tende ad allontanarsi dalla distribuzione normale sulle due code, come visto nel modello precedente.
Il terzo grafico ci mostra una distribuzione della nuvola abbastanza casuale, che ci farebbe pensare ad una varianza costante dei residui.
Il quarto grafico non presenta più la criticità dell’osservazione 1551, che abbiamo eliminato.

Andiamo ad indagare i risultati grafici attraverso diversi i test statistici:

bptest(mod7_new)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod7_new
## BP = 15.81, df = 7, p-value = 0.02691

dwtest(mod7_new)

## 
##  Durbin-Watson test
## 
## data:  mod7_new
## DW = 1.9499, p-value = 0.1053
## alternative hypothesis: true autocorrelation is greater than 0

shapiro.test(residuals(mod7_new))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod7_new)
## W = 0.98995, p-value = 3.172e-12

In questo caso, il test Breusch-Pagan, che valuta l’omoschedasticità, ci restituisce un valore di p-value=0.02691. L’assunzione di varianza costante, quindi, non sembra comunque essere rispettata. Detto ciò, unendo le considerazioni grafiche a quelle statistiche, possiamo notare che anche in questo caso non ci sono pattern forti ed è probabile che l’eteroschedasticità sia minima e non comprometta la capacità predittiva del modello. Il test Durbin-Watson, che valuta la correlazione tra i residui, ha un valore di p-value=0.1053 maggiore di α=0.05. Per tanto possiamo dire che i residui non sono autocorrelati. Il test Shapiro-Wilk, che valuta la distribuzione normale, ci restituisce un valore di p-value=3.172e-12. Anche qui l’assunzione di distribuzione normale non sembra essere rispettata. Detto ciò, unendo le considerazioni grafiche a quelle statistiche, che mostrano delle anomalie sulle code, potremmo ritenerla accettabile. Il grafico della distribuzione dei residui si presenterà così:

plot(density(residuals(mod7_new)), lwd=1.5, col="limegreen")

Andiamo ad indagare i valori leverages:

lev<-hatvalues(mod7_new)
plot(lev,pch=20,col="limegreen")
p=sum(lev)
soglia=2*p/n
abline(h=soglia,col="red3",lwd=2)

length(lev[lev>soglia])

## [1] 214

Abbiamo 214 valori leverages. Andiamo a indagare i valori outliers:

plot(rstudent(mod7_new),pch=20,col="limegreen")
abline(h=c(-2,2),col="red3",lwd=2)

outlierTest(mod7_new)

##       rstudent unadjusted p-value Bonferroni p
## 155   5.016215         5.6416e-07    0.0014087
## 1306  4.915885         9.4135e-07    0.0023506
## 1399 -4.432176         9.7331e-06    0.0243030
## 1694  4.352331         1.4014e-05    0.0349920

Abbiamo 4 valori outliers. Andiamo a verificare la distanza di Cook:

cook <- cooks.distance(mod7_new)
plot(cook,pch=20,col="limegreen")

max(cook)

## [1] 0.05027286

In questo caso la distanza di Cook è di 0.05, quindi minore anche della soglia di allarme.

A seguito di tutte le analisi effettuate, decidiamo di usare il mod7_new, seppur alcune delle assunzioni non siano rispettate (omoschedasticità e distribuzione normale). Probabilmente con i GLM potrebbe essere possibile trovare un modello migliore. Scegliamo anche, per motivi di studio, di tenere la variabile Fumatrici seppur poco significativa, come già esplicato, per motivi di studio.

3. Previsione e Risultati:

Utilizziamo il modello scelto per effettuare una stima del Peso per una neonata considerando una madre non fumatrice, alla terza gravidanza, che partorirà alla 39esima settimana. Non avendo il valore di Lunghezza e Cranio prendiamo come valore la media dei valori del nostro dataset.

lunghezza_media=mean(data_new$Lunghezza)
cranio_media=mean(data_new$Cranio)

neonata<- data.frame(
  N.gravidanze = 2,
  Fumatrici="Non fumatrici",
  Gestazione = 39,
  Lunghezza = lunghezza_media,
  Cranio = cranio_media,
  Sesso = "F",
  I.Lunghezza.2 = lunghezza_media^2
)

Dati neonata - Tabella
N.gravidanze	Fumatrici	Gestazione	Lunghezza	Cranio	Sesso	I.Lunghezza.2
2	Non fumatrici	39	494.768	340.016	F	244795.1

peso_previsto <-predict(mod7_new, newdata = neonata)
peso_previsto

##        1 
## 3248.289

La stima della variabile Peso sarà di 3248.289g. Il valore è valido.

4. Visualizzazioni:

Visualizziamo graficamente i risultati del modello per mostrare le relazioni tra la nostra variabile risposta Peso e le variabili più significative.

Peso e N.gravidanze:

ggplot(data_new,aes(x=N.gravidanze,
                    y=Peso,
                    group=Sesso))+
  geom_point(aes(color=Sesso), lwd=3, position="Jitter")+
  geom_smooth(aes(color=Sesso),lwd=2,se=F,method="lm")+
  scale_color_manual(name="",
                     breaks = c("F","M"),
                     values=c("violetred2","royalblue2"),
                     labels=c("F","M"))+
  labs(
    title = "Relazione tra il numero di gravidanze e il peso del neonato",
    x = "Settimane di gestazione",
    y = "Peso del neonato [g]",
    color = "Sesso"
  ) +
  scale_x_continuous(breaks=seq(0,12,1))+
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

## Warning in geom_point(aes(color = Sesso), lwd = 3, position = "Jitter"):
## Ignoring unknown parameters: `linewidth`

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

Nel grafico possiamo notare una differenza tra maschi e femmine. All’aumentare del numero delle gravidanze, la retta mostra che le femmine tendono a diminuire leggermente di peso, mentre i maschi ad aumentarlo leggermente. Quindi tra i due sessi sembra esserci una tendenza diversa.

Peso e Gestazione:

ggplot(data_new,aes(x=Gestazione,
                    y=Peso,
                    group=Sesso))+
  geom_point(aes(color=Sesso), lwd=3, position="Jitter")+
  geom_smooth(aes(color=Sesso),lwd=2,se=F,method="lm")+
  scale_color_manual(name="",
                     breaks = c("F","M"),
                     values=c("violetred2","royalblue2"),
                     labels=c("F","M"))+
  labs(
    title = "Relazione tra le settimane di gestazione e il peso del neonato",
    x = "Settimane di gestazione",
    y = "Peso del neonato [g]",
    color = "Sesso"
  ) +
  scale_x_continuous(breaks=seq(25,43,1))+
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

## Warning in geom_point(aes(color = Sesso), lwd = 3, position = "Jitter"):
## Ignoring unknown parameters: `linewidth`

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

Il grafico mostra come all’aumentare delle settimane di gestazione, per entrambi i sessi, il peso aumenta. Si noti anche la differenza tra maschi e femmine, i primi, come già visto con i boxplot, pesano più delle seconde.

Peso e Lunghezza:

ggplot(data_new,aes(x=Lunghezza,
                    y=Peso,
                    group=Sesso))+
  geom_point(aes(color=Sesso), lwd=3, position="Jitter")+
  geom_smooth(aes(color=Sesso),lwd=2,se=F,method="lm")+
  scale_color_manual(name="",
                     breaks = c("F","M"),
                     values=c("violetred2","royalblue2"),
                     labels=c("F","M"))+
  labs(
    title = "Relazione tra la lunghezza [mm] ed il peso del neonato",
    x = "Lunghezza [mm]",
    y = "Peso del neonato [g]",
    color = "Sesso"
  ) +
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

## Warning in geom_point(aes(color = Sesso), lwd = 3, position = "Jitter"):
## Ignoring unknown parameters: `linewidth`

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

Il grafico mostra come all’aumentare della lunghezza, aumenta anche il peso. Anche qui abbiamo una leggera differenza tra i due sessi.

Peso e Cranio:

ggplot(data_new,aes(x=Cranio,
                    y=Peso,
                    group=Sesso))+
  geom_point(aes(color=Sesso), lwd=3, position="Jitter")+
  geom_smooth(aes(color=Sesso),lwd=2,se=F,method="lm")+
  scale_color_manual(name="",
                     breaks = c("F","M"),
                     values=c("violetred2","royalblue2"),
                     labels=c("F","M"))+
  labs(
    title = "Relazione tra il diametro del Cranio [mm] ed il peso del neonato",
    x = "Diametro Cranio [mm]",
    y = "Peso del neonato [g]",
    color = "Sesso"
  ) +
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

## Warning in geom_point(aes(color = Sesso), lwd = 3, position = "Jitter"):
## Ignoring unknown parameters: `linewidth`

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

Il grafico mostra come all’aumentare della dimensione del cranio, aumenta anche il peso. Anche qui abbiamo una leggera differenza tra i due sessi.

Peso-Gestazione-Fumo:

Nell’ultimo grafico vediamo la relazione tra le settimane di gestazione, il peso ed il fumo.

ggplot(data_new,aes(x=Gestazione,
                    y=Peso,
                    group=Fumatrici))+
  geom_point(aes(color=Fumatrici), lwd=3, position="Jitter")+
  geom_smooth(aes(color=Fumatrici),lwd=2,se=F,method="lm")+
  scale_color_manual(name="",
                     breaks = c("Non fumatrici","Fumatrici"),
                     values=c("powderblue","lightslategrey"),
                     labels=c("Non fumatrici","Fumatrici"))+
  labs(
    title = "Relazione tra gestazione-peso-fumo",
    x = "Settimane di gestazione",
    y = "Peso del neonato [g]",
    color = "Fumatrici"
  )+
  scale_x_continuous(breaks=seq(25,43,1))+
  theme_minimal() + 
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), 
    axis.title.x = element_text(size = 12, face = "bold"),  
    axis.title.y = element_text(size = 12, face = "bold"), 
    legend.title = element_text(size = 12, face = "bold"),  
    legend.position = "right"  
  )

## Warning in geom_point(aes(color = Fumatrici), lwd = 3, position = "Jitter"):
## Ignoring unknown parameters: `linewidth`

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

Osservando il grafico, possiamo vedere come in effetti nel nostro campione, essendo le gestanti non fumatrici decisamente più numerose rispetto alle fumatrici, è difficile trovare un’evidenza netta. Nonostante ciò, possiamo notare una tendenza del peso, intorno alla fine della gravidanza, ad essere inferiore nel caso di madre fumatrice.

Modello Statistico per la Previsione del Peso Neonatale

Scionti_Maria

2025-02-13

1. Raccolta dei Dati e Struttura del Dataset:

2. Analisi e Modellizzazione:

Analisi Preliminare:

Analisi delle variabili del dataset “neonati”:

Analisi Descrittiva:

Distribuzione di frequenze:

Distribuzione di frequenze per Fumatrici:

Distribuzione di frequenze per Tipo.parto:

Distribuzione di frequenze per Ospedale:

Distribuzione di frequenze per Sesso:

Indici di posizione:

Indici di variabilità:

Distribuzione di frequenze di N.gravidanze:

Indici di forma:

Outliers:

Anni.madre:

N.gravidanze:

Gestazione:

Lunghezza:

Peso:

Cranio:

Ipotesi:

Creazione del Modello di Regressione:

Variabili Qualitative:

Peso~Tipo.parto:

Peso~Sesso:

Peso~Ospedale:

Peso~Fumatrici:

Creazione del Modello di Regressione:

MODELLO 1:

MODELLO 2:

MODELLO 3:

MODELLO 4:

MODELLO 5:

MODELLO 6:

MODELLO 7:

MODELLO 8:

Selezione del modello ottimale:

Analisi della Qualità del Modello:

Analisi dei residui e presenza si outliers e leverages:

Outliers e Leverages:

3. Previsione e Risultati:

4. Visualizzazioni:

Peso e N.gravidanze:

Peso e Gestazione:

Peso e Lunghezza:

Peso e Cranio:

Peso-Gestazione-Fumo: