Modello statistico per la previsione del peso neonatale
Laura Notari
2024-12-09
L’obiettivo dello studio è la creazione di un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita, basandosi su variabili cliniche raccolte da tre ospedali. Il progetto mira a migliorare la gestione delle gravidanze ad alto rischio, ottimizzare le risorse ospedaliere e garantire migliori risultati e garantire migliori risultati per la salute neonatale.
1. Raccolta dei dati e struttura del dataset
Andiamo ad importare il dataset con i dati oggetto del nostro studio
neonati <- read.csv("neonati.csv",
stringsAsFactors = T,
sep = ",")All’interno della funzione read.csv ho specificato il separatore di colonne del foglio elettronico (virgola) e ho utilizzato questo argomento stringAsFactor = True perché all’interno del dataset ci sono delle variabili qualitative che vanno lette come fattori e non come semplici string
Andiamo a verificare la presenza di dati mancanti (Na e NaN)
colSums(is.na(neonati))## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza
## 0 0 0 0 0 0
## Cranio Tipo.parto Ospedale Sesso
## 0 0 0 0Visto che per ogni colonna abbiamo un valore di 0 sappiamo che non ci sono dati mancanti
Verifichiamo quanti set di dati abbiamo
nrow(neonati)## [1] 2500Visualizziamo le prime 5 righe del nostro dataset
head(neonati)## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## 1 26 0 0 42 3380 490 325 Nat
## 2 21 2 0 39 3150 490 345 Nat
## 3 34 3 0 38 3640 500 375 Nat
## 4 28 1 0 41 3690 515 365 Nat
## 5 20 0 0 38 3700 480 335 Nat
## 6 32 0 0 40 3200 495 340 Nat
## Ospedale Sesso
## 1 osp3 M
## 2 osp1 F
## 3 osp2 M
## 4 osp2 M
## 5 osp3 F
## 6 osp2 FDunque abbiamo un dataset costituito da 2500 osservazioni di 10 variabili.
Ecco le variabili presenti:
• Età della madre: Misura dell’età in anni. (variabile quantitativa continua)
• Numero di gravidanze: Quante gravidanze ha avuto la madre. (variabile quantitativa discreta)
• Fumo materno: Un indicatore binario (0=non fumatrice, 1=fumatrice). (variabile qualitativa dummy)
• Durata della gravidanza: Numero di settimane di gestazione. (variabile quantitativa continua)
• Peso del neonato: Peso alla nascita in grammi. (variabile quantitativa continua)
• Lunghezza e diametro del cranio: Lunghezza del neonato e diametro craniale, misurabili anche durante la gravidanza tramite ecografie. (variabile quantitativa continua) I valori riportati nel dataset sono verosimilmente relativi alla circonferenza del cranio e non al diametro.
• Tipo di parto: Naturale o cesareo. (variabile qualitativa su base nominale)
• Ospedale di nascita: Ospedale 1, 2 o 3. (variabile qualitativa su base nominale)
• Sesso del neonato: Maschio (M) o femmina (F). (variabile qualitativa su base nominale)
Anche la colonna con la variabile dummy Fumatrici (che attualmente è una variabile int) deve essere trasformata in fattori
neonati$Fumatrici <- as.factor(neonati$Fumatrici)2. Analisi e modellizzazione
Analisi preliminare
A questo punto facciamo un summary del nostro dataset per svolgere un’analisi preliminare
summary(neonati)## Anni.madre N.gravidanze Fumatrici Gestazione Peso
## Min. : 0.00 Min. : 0.0000 0:2396 Min. :25.00 Min. : 830
## 1st Qu.:25.00 1st Qu.: 0.0000 1: 104 1st Qu.:38.00 1st Qu.:2990
## Median :28.00 Median : 1.0000 Median :39.00 Median :3300
## Mean :28.16 Mean : 0.9812 Mean :38.98 Mean :3284
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:40.00 3rd Qu.:3620
## Max. :46.00 Max. :12.0000 Max. :43.00 Max. :4930
## Lunghezza Cranio Tipo.parto Ospedale Sesso
## Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1256
## 1st Qu.:480.0 1st Qu.:330 Nat:1772 osp2:849 M:1244
## Median :500.0 Median :340 osp3:835
## Mean :494.7 Mean :340
## 3rd Qu.:510.0 3rd Qu.:350
## Max. :565.0 Max. :390Da questo summary possiamo vedere che c’è certamente un problema con la variabile Anni.madre in quanto ha un valore minimo di 0, cosa ovviamente non possibile. Diamo un’occhiata al suo grafico
plot(neonati$Anni.madre, pch=20, col="royalblue3")
Abbiamo due osservazioni molto basse, una pari a 0 e una certamente minore di 5, che vanno eliminate perchè sono frutto di errore in quanto non è fisiologicamente possibile avere figli a quell’età
neonati <- subset(neonati, Anni.madre > 5)Andiamo ad osservare il summary aggiornato
summary(neonati)## Anni.madre N.gravidanze Fumatrici Gestazione Peso
## Min. :13.00 Min. : 0.0000 0:2394 Min. :25.00 Min. : 830
## 1st Qu.:25.00 1st Qu.: 0.0000 1: 104 1st Qu.:38.00 1st Qu.:2990
## Median :28.00 Median : 1.0000 Median :39.00 Median :3300
## Mean :28.19 Mean : 0.9816 Mean :38.98 Mean :3284
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:40.00 3rd Qu.:3620
## Max. :46.00 Max. :12.0000 Max. :43.00 Max. :4930
## Lunghezza Cranio Tipo.parto Ospedale Sesso
## Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1255
## 1st Qu.:480.0 1st Qu.:330 Nat:1770 osp2:848 M:1243
## Median :500.0 Median :340 osp3:834
## Mean :494.7 Mean :340
## 3rd Qu.:510.0 3rd Qu.:350
## Max. :565.0 Max. :390Dal summary possiamo osservare che le madri:
hanno tra i 13 e i 46 anni con una mediana di 28 e una media di 28.19
hanno avuto tra 0 e 12 gravidanze precedenti con una mediana di 1 e una media di 0.98
non fumatrici sono 2394 e fumatrici sono 104.
Le settimane di gestazione a cui sono nati i neonati variano tra 25 e 43 con una mediana di 39 e una media di 38.98
I neonati hanno: - un peso che varia tra 830 e 4930 grammi con una mediana di 3300 e una media di 3284
una lunghezza che varia tra 310 e 565 millimetri con una mediana di 500 e una media di 494.7
un diametro del cranio che varia tra 235 e 390 millimetri con una mediana e una media di 340
728 neonati sono nati per parto cesareo e 1770 per parto naturale
Ci sono 816 neonati nati nell’ospedale 1, 848 nell’ospedale 2 e 834 nell’ospedale 3
Sesso: 1255 sono neonati femmine e 1243 sono maschi
Salviamo il numero di righe del nostro dataset nella variabile n
n <- nrow(neonati)Le righe del dataset adesso sono 2498
Estraiamo le colonne dal dataframe neonati in modo da non dover richiamare tutte le volte il dataset di riferimento e usare il simbolo dollaro
attach(neonati)Saggiamo l’ipotesi che i valori della media di peso e lunghezza del nostro dataset siano significativamente uguali a quelli della popolazione
Prendiamo come valori di riferimento quelli presenti in questa pagina internet: https://www.ospedalebambinogesu.it/da-0-a-30-giorni-come-si-presenta-e-come-cresce-80012/ ovvero Peso alla nascita 3300g e lunghezza 500 mm
#H0: peso medio = 3300
#H1: peso medio != 3300
t.test(Peso,
mu = 3300,
conf.level = 0.95,
alternative = "two.sided")##
## One Sample t-test
##
## data: Peso
## t = -1.505, df = 2497, p-value = 0.1324
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
## 3263.577 3304.791
## sample estimates:
## mean of x
## 3284.184Non si rifiuta l’ipotesi nulla, dunque il valore della media del peso di questo campione è significativamente uguale alla media del peso della popolazione.
Procediamo allo stesso modo per la lunghezza.
#H0: lunghezza media = 500
#H1: lunghezza media != 500
t.test(Lunghezza,
mu = 500,
conf.level = 0.95,
alternative = "two.sided")##
## One Sample t-test
##
## data: Lunghezza
## t = -10.069, df = 2497, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
## 493.6628 495.7287
## sample estimates:
## mean of x
## 494.6958Si rifiuta l’ipotesi nulla, dunque secondo il t-test il valore della media della lunghezza di questo campione non è significativamente uguale alla media della lunghezza della popolazione. La differenza tra le due medie è però di soli 6 mm quindi la questione va approfondita al di là del risultato del test.
Proviamo a fare lo scatterplot delle diverse variabili quantitative per capire se ci sono altri outliers o anomalie
plot(Peso, pch=20, col="royalblue3")
plot(Anni.madre, pch=20, col="royalblue3")
plot(N.gravidanze, pch=20, col="royalblue3")
plot(Gestazione, pch=20, col="royalblue3")
plot(Lunghezza, pch=20, col="royalblue3")
plot(Cranio, pch=20, col="royalblue3")
Guardando gli scatterplot delle singole variabili non sembra esserci nessuna anomalia specifica come accadeva anche leggendo i dati del summary.
Analizziamo invece i grafici di 2 variabili numeriche alla volta per visualizzare le loro relazioni.
plot(Anni.madre, Peso, pch=20, col="royalblue3")
nessuna anomalia rilevabile
plot(N.gravidanze, Peso, pch=20, col="royalblue3")
nessuna anomalia rilevabile
plot(Gestazione, Peso, pch=20, col="royalblue3")
Qui abbiamo un primo valore un po’ fuori dalla nuvola ovvero peso intorno ai 4500g per settimana di gestazione pari a 35
Troviamo il numero di questa osservazione
which(Gestazione==35 & Peso > 4000)## [1] 1551è la numero 1551
plot(Gestazione, Lunghezza, pch=20, col="royalblue3")
qui invece abbiamo un valore un po’ fuori dalla nuvola ovvero gestazione pari a 38 e lunghezza tra 300 e 350 mm Troviamo il numero di questa osservazione
which(Gestazione==38 & Lunghezza < 350)## [1] 1549è la numero 1549
plot(Lunghezza, Peso, pch=20, col="royalblue3")
Qui troviamo un valore “fuori” dalla nuvola formata dagli altri valori: peso molto alto (tra i 4000 e i 4500g) rispetto alla lunghezza che è vicina a 300mm
Troviamo il numero dell’osservazione:
which(Lunghezza<350 & Peso > 4000)## [1] 1549è sempre la numero 1549, l’osservazione che aveva un valore di lunghezza piccolo per la settimana di gestazione
plot(Cranio, Peso, pch=20, col="royalblue3")
qui abbiamo un paio di osservazioni con valori di peso piccolo rispetto alla circonferenza del cranio.
Troviamoli:
which(Cranio>350 & Peso <2000)## [1] 310 1427Sono le osservazioni 310 e 1427
proviamo a guardare il grafico del peso
ggplot()+
geom_histogram(aes(x=Peso),
stat="density",
col="royalblue3",
fill="royalblue3")+
labs(title = "Distribuzione della variabile peso",
x = "Peso (g)",
y = "Densità")+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))
Dal grafico vediamo che la coda di sinistra è più lunga della coda di destra
Realizziamo i grafici dei regressori
ggplot()+
geom_histogram(aes(x=Anni.madre),
col="royalblue4",
fill="royalblue3",
bins = 34)+
labs(title = "Distribuzione della variabile anni madre",
x = "Anni madre",
y = "Frequenze assolute")+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))
la distribuzione ha la coda di destra più lunga
ggplot()+
geom_histogram(aes(x=N.gravidanze),
col="royalblue4",
fill="royalblue3",
bins = 13)+
labs(title = "Distribuzione della variabile numero gravidanze",
x = "Numero gravidanze precedenti",
y = "Frequenze assolute")+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))
Il numero di unità statistiche relative a donne che hanno avuto più di 3 gravidanze precedenti è molto basso.
ggplot()+
geom_histogram(aes(x=Gestazione),
col="royalblue4",
fill="royalblue3",
bins = 19)+
labs(title = "Distribuzione della variabile gestazione",
x = "Gestazione (settimane)",
y = "Frequenze assolute")+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))
la coda di sinistra è molto più lunga di quella di destra
ggplot()+
geom_histogram(aes(x=Lunghezza),
stat="density",
col="royalblue3",
fill="royalblue3")+
labs(title = "Distribuzione della variabile lunghezza",
x = "Lunghezza (mm)",
y = "Densità")+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))
La coda di sinistra è molto più lunga di quella di destra.
ggplot()+
geom_histogram(aes(x=Cranio),
stat="density",
col="royalblue3",
fill="royalblue3")+
labs(title = "Distribuzione della variabile cranio",
x = "Circonferenza cranio (mm)",
y = "Densità")+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))
Anche qui la coda di sinistra è molto più lunga di quella di destra.
ggplot(neonati)+
geom_bar(aes(x=Fumatrici),
fill="royalblue3",
position = "dodge",
stat = "count",
col = "royalblue4")+
geom_text(aes(x = Fumatrici,
label = after_stat(count)),
stat = "count",
vjust = -0.5,
size = 4,
color = "royalblue4") +
labs(title = "Distribuzione della variabile fumatrici",
x = "Fumatrici",
y = "Frequenze assolute")+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))+
guides(fill=guide_legend(title=""))+
scale_x_discrete(labels = c("0" = "No", "1" = "Sì"))
Le mamme sono prevalentemente non fumatrici. Solo il 4.16% sono fumatrici.
ggplot(neonati)+
geom_bar(aes(x=Tipo.parto),
fill="royalblue3",
position = "dodge",
stat = "count",
col = "royalblue4")+
geom_text(aes(x = Tipo.parto,
label = ..count..),
stat = "count",
vjust = -0.5,
size = 4,
color = "royalblue4") +
labs(title = "Distribuzione della variabile tipo parto",
x = "Tipo parto",
y = "Frequenze assolute")+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))+
guides(fill=guide_legend(title=""))+
scale_x_discrete(labels = c("Ces" = "Cesareo", "Nat" = "Naturale"))
Quasi un terzo (29.14%) dei parti è cesareo.
ggplot(neonati)+
geom_bar(aes(x=Ospedale),
fill="royalblue3",
position = "dodge",
stat = "count",
col = "royalblue4")+
labs(title = "Distribuzione della variabile ospedale",
x = "",
y = "Frequenze assolute")+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))+
guides(fill=guide_legend(title=""))+
scale_x_discrete(labels = c("osp1" = "Ospedale 1", "osp2" = "Ospedale 2", "osp3" = "Ospedale 3"))
Possiamo dire che circa un terzo viene dall’ospedale 1, un terzo dall’ospedale 2 e un terzo dall’ospedale 3
ggplot(neonati)+
geom_bar(aes(x=Sesso),
fill="royalblue3",
position = "dodge",
stat = "count",
col = "royalblue4")+
labs(title = "Distribuzione della variabile sesso",
x = "Sesso",
y = "Frequenze assolute")+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))+
guides(fill=guide_legend(title=""))+
scale_x_discrete(labels = c("F" = "Femmina", "M" = "Maschio"))
Circa metà del campione è di sesso femminile e metà di sesso maschile.
ggplot(data = neonati)+
geom_half_boxplot(aes(x=Tipo.parto,y=Peso),
side="l",fill="violet")+
geom_half_violin(aes(x=Tipo.parto,y=Peso),
side="r",fill="pink")+
labs(x="Tipo parto",
y="Peso (g)",
title = "Distribuzione del peso \nin relazione al tipo di parto")+
scale_y_continuous(breaks = seq(750,5000,by=250))+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))+
scale_x_discrete(labels = c("Ces" = "Cesareo", "Nat" = "Naturale"))
Non si notano differenze sostanziali nel peso in base al tipo di parto
ggplot(data = neonati)+
geom_half_boxplot(aes(x=Tipo.parto,y=Peso,fill=Ospedale),
side="l", alpha=0.5)+
geom_half_violin(aes(x=Tipo.parto,y=Peso,fill=Ospedale),
side="r", alpha=0.5)+
labs(x="Tipo parto",
y="Peso (g)",
title = "Distribuzione del peso in relazione \nal tipo di parto e all'ospedale")+
scale_y_continuous(breaks = seq(750,5000,by=250))+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))+
scale_fill_discrete(labels = c("osp1" = "Ospedale 1", "osp2" = "Ospedale 2", "osp3" = "Ospedale 3"))+
guides(fill=guide_legend(title="Ospedale"))+
scale_x_discrete(labels = c("Ces" = "Cesareo", "Nat" = "Naturale"))
Ci sono lievissime differenze nel peso in base al tipo di parto e all’ospedale in cui si è partorito
Proviamo a fare una tabella di contingenza tra le variabili tipo parto e ospedale per vedere se ci sono ospedali che fanno più cesarei di altri
Visto che nel nostro campione le unità statistiche prese dall’ospedale 1 sono in numero diverso rispetto a quelle dell’ospedale 2 e 3 andiamo a valutare la proporzione del tipo di parto e non il numero in sè.
table_proporzioni_ospedale_tipo.parto <-round(prop.table(table(Ospedale, Tipo.parto), margin=1)*100)table_proporzioni_ospedale_tipo.parto## Tipo.parto
## Ospedale Ces Nat
## osp1 30 70
## osp2 30 70
## osp3 28 72C’è una lieve differenza percentuale nell’ospedale 3
Andiamo ad effettuare il test del chi quadro per capire se il numero di parti cesarei è condizionato dall’ospedale in cui avviene il parto
#H0: Il numero di parti cesarei non è condizionato dall'ospedale in cui avviene il parto (ipotesi di indipendenza)
#H1: Il numero di parti cesarei è condizionato dall'ospedale in cui avviene il parto (ipotesi di dipendenza)
chisq.test(table_proporzioni_ospedale_tipo.parto)##
## Pearson's Chi-squared test
##
## data: table_proporzioni_ospedale_tipo.parto
## X-squared = 0.12864, df = 2, p-value = 0.9377Non si rifiuta l’ipotesi di indipendenza, quindi non si può dire che il tipo parto dipenda dall’ospedale in cui avviene.
ggplot(data = neonati)+
geom_half_boxplot(aes(fill=Ospedale,y=Peso),
side="l", alpha=0.5)+
geom_half_violin(aes(fill=Ospedale,y=Peso),
side="r", alpha=0.5)+
labs(x="",
y="Peso (g)",
title = "Distribuzione del peso \nin relazione all'ospedale")+
scale_y_continuous(breaks = seq(750,5000,by=250))+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face = "bold"),
axis.text.x = element_blank(),
axis.ticks.x = element_blank()) +
scale_fill_discrete(
labels = c("osp1" = "Ospedale 1", "osp2" = "Ospedale 2", "osp3" = "Ospedale 3"),
name = "Ospedale")
Come ci si poteva aspettare non ci sono differenze evidenti tra la distribuzione del peso in relazione ai tre ospedali
Analizziamo il grafico relativo alla distribuzione del peso in relazione all’essere fumatrice della madre
ggplot(data = neonati)+
geom_half_boxplot(aes(fill=Fumatrici,y=Peso),
side="l")+
geom_half_violin(aes(fill=Fumatrici,y=Peso),
side="r")+
labs(x="",
y="Peso (g)",
title = "Distribuzione del peso \nin relazione all'essere fumatrice della madre")+
scale_y_continuous(breaks = seq(750,5000,by=250))+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face = "bold"),
axis.text.x = element_blank(),
axis.ticks.x = element_blank()) +
scale_fill_manual(
values = c("0" = "royalblue3", "1" = "gray"),
labels = c("0" = "No", "1" = "Sì"),
name = "Fumatrici")
Dal grafico sembra esserci una differenza tra le due distribuzioni e sembra che l’essere fumatrice porti a un peso neonatale minore.
Leggiamo i valori di peso medio per le due categorie di madri.
neonati_no_fumatrici <- neonati %>%
filter(Fumatrici==0)
neonati_fumatrici <- neonati %>%
filter(Fumatrici==1)
mean(neonati_no_fumatrici$Peso)## [1] 3286.262mean(neonati_fumatrici$Peso)## [1] 3236.346La media per madri non fumatrici è 3286.26 e per madri fumatrici 3236.35
mean(neonati_no_fumatrici$Peso) - mean(neonati_fumatrici$Peso)## [1] 49.91617La differenza è di quasi 50g
Proviamo a saggiare l’ipotesi che questa differenza sia significativa
#H0: peso medio fumatrici = peso medio non fumatrici
#H1: peso medio fumatrici != peso medio non fumatrici
t.test(neonati_fumatrici$Peso,
mu = mean(neonati_no_fumatrici$Peso),
conf.level = 0.95,
alternative = "two.sided")##
## One Sample t-test
##
## data: neonati_fumatrici$Peso
## t = -1.0632, df = 103, p-value = 0.2902
## alternative hypothesis: true mean is not equal to 3286.262
## 95 percent confidence interval:
## 3143.232 3329.460
## sample estimates:
## mean of x
## 3236.346Non si rifiuta l’ipotesi nulla, dunque il valore della media del peso di neonati nati da madri fumatrici è significativamente uguale alla media del peso di neonati nati da madri non fumatrici.
Proviamo però ad approfondire un attimo la questione.
Vi è una classificazione utilizzata per stabilire il grado di prematurità e quindi il rischio alla nascita di un neonato ed è la seguente. Nati: - Estremamente pretermine, prima della 28^ settimana di gestazione. - Molto pretermine, tra le 28 e 31+6 settimane. - Moderatamente pretermine, tra le 32 e 33+6 settimane. - Lievemente pretermine, tra le 34 e le 36+6 settimane.
Proviamo ad analizzare il peso in relazione alla settimana di gestazione e all’essere fumatrice della madre. Dividiamo in 5 classi Gestazione
neonati$Gestazione_CL <- cut(Gestazione,
breaks = c(25,28,32,34,37,43),
include.lowest = TRUE,
right = FALSE,
labels = c("Estremamente pret.","Molto pret.","Moderatamente pret.","Lievemente pret.","A termine"))Visualizziamo le frequenze delle madri fumatrici per ogni classe
table(neonati$Gestazione_CL,Fumatrici)## Fumatrici
## 0 1
## Estremamente pret. 4 0
## Molto pret. 20 0
## Moderatamente pret. 26 1
## Lievemente pret. 110 1
## A termine 2234 102ggplot(data = neonati, aes(y=Peso, x=Gestazione_CL, fill=Fumatrici))+
geom_boxplot()+
theme_classic()+
scale_y_continuous(breaks = seq(750,5000,by=250))+
scale_fill_manual(
values = c("0" = "royalblue3", "1" = "gray"),
labels = c("0" = "No", "1" = "Sì"),
name = "Fumatrici")+
labs(x="Grado di prematurità",
y="Peso (g)",
title = "Distribuzione del peso nei nati a termine\nin relazione all'essere fumatrice della madre")
Per le classi relative ai parti pretermine c’è una scarsissima
(praticamente nulla) rappresentanza di madri fumatrici. Ciò potrebbe
essere dovuto ad esempio alle dimensioni ridotte del campione. Visto ciò
andiamo a rifare l’analisi precedente valutando esclusivamente i dati
relativi a Gestazione maggiore o uguale a 37 ovvero per i nati a
termine.
Analizziamo il grafico relativo alla distribuzione del peso in relazione all’essere fumatrice della madre
ggplot(data = filter(neonati, Gestazione_CL=="A termine"))+
geom_half_boxplot(aes(fill=Fumatrici,y=Peso),
side="l")+
geom_half_violin(aes(fill=Fumatrici,y=Peso),
side="r")+
labs(x="",
y="Peso (g)",
title = "Distribuzione del peso nei nati a termine\nin relazione all'essere fumatrice della madre")+
scale_y_continuous(breaks = seq(750,5000,by=250))+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face = "bold"),
axis.text.x = element_blank(),
axis.ticks.x = element_blank()) +
scale_fill_manual(
values = c("0" = "royalblue3", "1" = "gray"),
labels = c("0" = "No", "1" = "Sì"),
name = "Fumatrici")
Dal grafico sembra esserci una differenza tra le due distribuzioni e sembra che l’essere fumatrice porti a un peso neonatale mediamente minore.
Leggiamo i valori di peso medio per le due categorie di madri.
neonati_no_fumatrici_a_termine <- neonati %>%
filter(Gestazione_CL=="A termine", Fumatrici==0)
neonati_fumatrici_a_termine <- neonati %>%
filter(Gestazione_CL=="A termine", Fumatrici==1)
mean(neonati_no_fumatrici_a_termine$Peso)## [1] 3353.81mean(neonati_fumatrici_a_termine$Peso)## [1] 3256.373La media per madri non fumatrici è 3353.81 e per madri fumatrici 3256.37
mean(neonati_no_fumatrici_a_termine$Peso) - mean(neonati_fumatrici_a_termine$Peso)## [1] 97.43766La differenza è di 97.4g
Proviamo a saggiare l’ipotesi che questa differenza sia significativa
#H0: peso medio fumatrici = peso medio non fumatrici
#H1: peso medio fumatrici != peso medio non fumatrici
t.test(neonati_fumatrici_a_termine$Peso,
mu = mean(neonati_no_fumatrici_a_termine$Peso),
conf.level = 0.95,
alternative = "two.sided")##
## One Sample t-test
##
## data: neonati_fumatrici_a_termine$Peso
## t = -2.1527, df = 101, p-value = 0.03373
## alternative hypothesis: true mean is not equal to 3353.81
## 95 percent confidence interval:
## 3166.582 3346.164
## sample estimates:
## mean of x
## 3256.373Si rifiuta l’ipotesi nulla, dunque secondo il t-test il valore della media del peso di neonati nati da madri fumatrici non è significativamente uguale alla media del peso di neonati nati da madri non fumatrici, la media del peso dei nati da madri non fumatrici si trova fuori dall’intervallo di confidenza.
D’ora in avanti come richiesto si considererà di nuovo il dataset per intero e non un subset dei nati a termine, ritenevo però importante evidenziare questo aspetto.
Ora andiamo ad analizzare la distribuzione delle variabili Peso, Lunghezza e Cranio condizionatamente alla variabile Sesso.
ggplot(data = neonati)+
geom_half_boxplot(aes(y=Peso,fill=Sesso),
side="l",alpha=0.5)+
geom_half_violin(aes(y=Peso,fill=Sesso),
side="r",alpha=0.5)+
labs(x="",
y="Peso (g)",
title = "Distribuzione del peso \nin relazione al sesso del neonato")+
scale_y_continuous(breaks = seq(750,5000,by=250))+
scale_x_continuous(breaks=seq(5,20,by=5))+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))+
scale_fill_discrete(
labels = c("F" = "Femmina", "M" = "Maschio"),
name = "Sesso")
Il peso dei neonati rispetto a quello delle neonate risulta essere maggiore
ggplot(data = neonati)+
geom_half_boxplot(aes(y=Lunghezza,fill=Sesso),
side="l", alpha=0.5)+
geom_half_violin(aes(y=Lunghezza,fill=Sesso),
side="r", alpha=0.5)+
labs(x="",
y="Lunghezza (mm)",
title = "Distribuzione della lunghezza \nin relazione al sesso del neonato")+
scale_y_continuous(breaks = seq(310,570,by=10))+
scale_x_continuous(breaks=seq(5,20,by=5))+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))+
scale_fill_discrete(
labels = c("F" = "Femmina", "M" = "Maschio"),
name = "Sesso")
La lunghezza dei neonati rispetto a quello delle neonate risulta essere maggiore
ggplot(data = neonati)+
geom_half_boxplot(aes(y=Cranio,fill=Sesso),
side="l", alpha=0.5)+
geom_half_violin(aes(y=Cranio,fill=Sesso),
side="r", alpha=0.5)+
labs(x="",
y="Cranio (mm)",
title = "Distribuzione della circonferenza del cranio \nin relazione al sesso del neonato")+
scale_y_continuous(breaks = seq(230,390,by=10))+
scale_x_continuous(breaks=seq(5,20,by=5))+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))+
scale_fill_discrete(
labels = c("F" = "Femmina", "M" = "Maschio"),
name = "Sesso")
La circonferenza del cranio dei neonati rispetto a quello delle neonate risulta essere maggiore
Prima di procedere alla realizzazione della matrice di correlazione e realizzare il primo modello è necessario fare una valutazione. Potremmo pensare di escludere la variabile Ospedale dal modello perché il primo pensiero che viene è “come potrebbe un ospedale diverso cambiare il peso del neonato?” e anche i grafici sopra suggeriscono la stessa cosa. Ritengo però che non avendo informazioni sulla tipologia di ospedali presi in esame, il livello di assistenza fornito, il tipo di terapia intensiva neonatale presente e così via, sia importante prendere in considerazione anche la variabile ospedale. Inoltre al variare delle caratteristiche elencate sopra un certo ospedale potrebbe essere preferito in caso di gravidanza a rischio che porta a un parto prematuro. Quindi potrebbe variare la distribuzione delle settimane di gestazione al variare dell’ospedale. Proviamo a vederne il grafico:
ggplot(data = neonati)+
geom_half_boxplot(aes(y=Gestazione, fill = Ospedale),
side="l",alpha=0.5)+
geom_half_violin(aes(y=Gestazione, fill = Ospedale),
side="r",,alpha=0.5, adjust=2)+
labs(x="",
y="Gestazione (settimane)",
title = "Distribuzione delle settimane di gestazione \nin relazione all'ospedale")+
scale_y_continuous(breaks = seq(26,44,by=2))+
scale_x_continuous(breaks=seq(5,20,by=5))+
scale_fill_discrete(labels = c("osp1" = "Ospedale 1", "osp2" = "Ospedale 2", "osp3" = "Ospedale 3"))+
theme_classic()+
theme(plot.title = element_text(size=14, hjust = 0.5, face = "bold"))+
theme(legend.title = element_text(face="bold"))
Graficamente la mediana è identica e anche l’andamento simile quindi sembra non esserci correlazione tra la variabile Gestazione e la variabile Ospedale, variano un po’ gli estremi per i quali in realtà c’è poca statistica come possiamo vedere qui:
table(Gestazione)## Gestazione
## 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
## 1 1 2 4 3 5 8 9 18 16 33 62 192 437 580 741 328 56 2Possiamo fare table perché nonostante Gestazione sia una variabile quantitativa continua nel nostro dataset assume solo un numero finito di valori interi.
Da quanto visto sinora non sembra essere utile mantenere la variabile Ospedale ma lo decideremo direttamente con i test successivi.
Creazione del modello di regressione
Procediamo alla realizzazione della matrice di correlazione sul dataset privo della colonna Gestazione_CL.
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor=0.8, ...)
{
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(r>=0.1) cex.cor <- 1.4
text(0.5, 0.5, txt, cex = cex.cor)
}
#correlazioni
pairs(x=neonati[,1:10],lower.panel=panel.cor, upper.panel=panel.smooth)
Dalla matrice di correlazione possiamo vedere che i regressori che hanno una correlazione importante con la variabile Peso sono: Lunghezza, Cranio, Gestazione e in piccola parte Sesso. Si rileva anche una possibile correlazione tra le variabili Lunghezza e Cranio, Lunghezza e Gestazione, Cranio e Gestazione. Non si rileva correlazione tra Ospedale e Gestazione. Negli scatterplot si rileva una non linearità nella relazione tra Peso e Gestazione, Lunghezza e Cranio. Partiamo da un modello con tutte le variabili ma senza termini quadratici.
summary(lm(Peso~Anni.madre+N.gravidanze+Fumatrici+Gestazione+Lunghezza+Cranio+Tipo.parto+Ospedale+Sesso, data=neonati))##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1123.26 -181.53 -14.45 161.05 2611.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6735.7960 141.4790 -47.610 < 2e-16 ***
## Anni.madre 0.8018 1.1467 0.699 0.4845
## N.gravidanze 11.3812 4.6686 2.438 0.0148 *
## Fumatrici1 -30.2741 27.5492 -1.099 0.2719
## Gestazione 32.5773 3.8208 8.526 < 2e-16 ***
## Lunghezza 10.2922 0.3009 34.207 < 2e-16 ***
## Cranio 10.4722 0.4263 24.567 < 2e-16 ***
## Tipo.partoNat 29.6335 12.0905 2.451 0.0143 *
## Ospedaleosp2 -11.0912 13.4471 -0.825 0.4096
## Ospedaleosp3 28.2495 13.5054 2.092 0.0366 *
## SessoM 77.5723 11.1865 6.934 5.18e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274 on 2487 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 668.7 on 10 and 2487 DF, p-value: < 2.2e-16Analizziamo la tabella ottenuta con questo summary Nella tabella dei coefficienti vengono mostrate le stime dei coefficienti beta di tutte le variabili che rappresentano gli effetti marginali di ogni singola variabile sulla variabile risposta Per ogni stima bisogna guardare segno, valore e significatività
Anni.madre: coefficiente positivo quindi effetto positivo anche se abbastanza basso ma poco significativo perchè p-value di 0.48 N.gravidanze: il numero di gravidanze precedenti della madre ha effetto positivo e elevato, con p-value sotto la soglia del 5% quindi effetto significativo. Dunque all’aumentare delle gravidanze aumenta il peso del neonato Fumatrici: l’essere fumatrici (valore della variabile = 1) ha un effetto negativo sul peso con significatività al di sopra del 5% classico quindi tale influenza è poco significativa Gestazione: coefficiente alto e positivo con p-value prossimo allo zero quindi significatività altissima dunque per ogni settimana di gestazione si ha un aumento di peso di 32.6g Lunghezza: coefficiente alto e positivo con p-value prossimo allo zero quindi significatività altissima dunque per ogni centimetro si ha un aumento di peso di circa 10.3g Cranio: coefficiente alto e positivo con p-value prossimo allo zero quindi significatività altissima dunque per ogni centimetro di circonferenza del cranio del neonato si ha un aumento di peso di 10.5g Tipo.parto: il tipo di parto naturale avviene per neonati con peso maggiore (coefficiente positivo alto e p-value sotto la soglia del 5% quindi differenza significativa) Ospedale: l’ospedale 2 ha un coefficiente negativo e alto in valore assoluto ma significatività bassa. Invece l’ospedale 3 ha un coefficiente positivo alto e p-value sotto la soglia del 5% quindi effetto significativo. Quindi la variabile nel suo complesso potrebbe aggiungere informazioni. Sesso: il valore Maschio della variabile sesso ha un coefficiente altissimo e p-value prossimo allo zero quindi effetto molto significativo
Si ha un R2 pari a 0.7289 quindi buono ma ancora più importante è il valore dell’R2 aggiustato che è anch’esso pari a 0.7278
Visto che il tipo parto è qualcosa che avviene alla fine e non è quindi utile in questo tipo di studio in tutti i modelli che realizzeremo questa variabile non sarà presente.
Creiamo il mod1 che avrà tutte le variabili tranne tipo parto e vediamo che differenze ci sono rispetto a quello completo.
mod1 <- lm(Peso~Anni.madre+N.gravidanze+Fumatrici+Gestazione+Lunghezza+Cranio+Ospedale+Sesso, data=neonati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Ospedale + Sesso, data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1142.99 -181.56 -14.41 164.29 2612.53
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6710.3696 141.2400 -47.510 < 2e-16 ***
## Anni.madre 0.8124 1.1479 0.708 0.4792
## N.gravidanze 11.0736 4.6716 2.370 0.0178 *
## Fumatrici1 -28.9707 27.5717 -1.051 0.2935
## Gestazione 32.6215 3.8246 8.530 < 2e-16 ***
## Lunghezza 10.2560 0.3008 34.094 < 2e-16 ***
## Cranio 10.5060 0.4265 24.635 < 2e-16 ***
## Ospedaleosp2 -11.1525 13.4606 -0.829 0.4074
## Ospedaleosp3 28.8129 13.5171 2.132 0.0331 *
## SessoM 77.6195 11.1978 6.932 5.28e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2488 degrees of freedom
## Multiple R-squared: 0.7282, Adjusted R-squared: 0.7273
## F-statistic: 740.8 on 9 and 2488 DF, p-value: < 2.2e-16Ci sono modifiche solamente marginali rispetto al modello precedente. Proviamo a semplificare il modello togliendo le variabili non significative.
mod2<-update(mod1,~.-Anni.madre-Fumatrici)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Ospedale + Sesso, data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1132.57 -183.77 -16.22 163.78 2620.13
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6682.3557 135.7530 -49.224 < 2e-16 ***
## N.gravidanze 12.0569 4.3372 2.780 0.00548 **
## Gestazione 32.0899 3.7962 8.453 < 2e-16 ***
## Lunghezza 10.2693 0.3005 34.174 < 2e-16 ***
## Cranio 10.5256 0.4259 24.716 < 2e-16 ***
## Ospedaleosp2 -10.9554 13.4579 -0.814 0.41570
## Ospedaleosp3 29.3428 13.5084 2.172 0.02993 *
## SessoM 77.5211 11.1952 6.924 5.55e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2490 degrees of freedom
## Multiple R-squared: 0.7281, Adjusted R-squared: 0.7273
## F-statistic: 952.4 on 7 and 2490 DF, p-value: < 2.2e-16Le stime sono cambiate soltanto lievemente e la significatività è rimasta pressapoco la stessa. Anche r2 aggiustato è rimasto fisso al 72.73% quindi possiamo dire che le variabili Anni.madre e Fumatrici erano inutili alla spiegazione della variabile risposta
Facciamo una versione del mod2 senza la variabile ospedale
mod2_1<-update(mod2,~.-Ospedale)
summary(mod2_1)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = neonati)
##
## 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-16L’R2 adj è praticamente identico quindi teniamo questa versione come buona
Sulla base di quanto visto nello scatterplot della matrice di correlazione proviamo ad aggiungere il termine quadratico di Gestazione
summary(update(mod2_1,~.+I(Gestazione^2)))##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2), data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1144.0 -181.5 -12.9 165.8 2661.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4646.7158 898.6322 -5.171 2.52e-07 ***
## N.gravidanze 12.5489 4.3381 2.893 0.00385 **
## Gestazione -81.2309 49.7402 -1.633 0.10257
## Lunghezza 10.3502 0.3040 34.045 < 2e-16 ***
## Cranio 10.6376 0.4282 24.843 < 2e-16 ***
## SessoM 75.7563 11.2435 6.738 1.99e-11 ***
## I(Gestazione^2) 1.5168 0.6621 2.291 0.02206 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.5 on 2491 degrees of freedom
## Multiple R-squared: 0.7276, Adjusted R-squared: 0.7269
## F-statistic: 1109 on 6 and 2491 DF, p-value: < 2.2e-16L’aggiunta di I(Gestazione^2) ha “rotto” la variabile Gestazione che adesso non è più significativa e R2 adj non ha avuto un aumento significativo. I(Gestazione^2) non va quindi aggiunto al modello
Anche l’utilizzo delle variabili Gestazione * Lunghezza o Lunghezza * Cranio porta agli stessi risultati
summary(update(mod2_1,~.+Gestazione*Lunghezza))##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + Gestazione:Lunghezza, data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1133.41 -179.98 -11.52 168.93 2652.65
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.991e+03 9.206e+02 -2.163 0.030631 *
## N.gravidanze 1.303e+01 4.321e+00 3.015 0.002594 **
## Gestazione -9.391e+01 2.481e+01 -3.785 0.000157 ***
## Lunghezza -8.476e-02 2.028e+00 -0.042 0.966661
## Cranio 1.076e+01 4.264e-01 25.234 < 2e-16 ***
## SessoM 7.225e+01 1.121e+01 6.445 1.38e-10 ***
## Gestazione:Lunghezza 2.729e-01 5.298e-02 5.151 2.79e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.3 on 2491 degrees of freedom
## Multiple R-squared: 0.7299, Adjusted R-squared: 0.7292
## F-statistic: 1122 on 6 and 2491 DF, p-value: < 2.2e-16Qui a “rompersi” è la variabile Lunghezza. R2 adj non ha avuto un aumento significativo.
summary(update(mod2_1,~.+Lunghezza*Cranio))##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + Lunghezza:Cranio, data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1150.65 -180.93 -13.48 165.99 2865.46
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.803e+03 1.018e+03 -1.771 0.0767 .
## N.gravidanze 1.293e+01 4.323e+00 2.991 0.0028 **
## Gestazione 3.815e+01 3.967e+00 9.616 < 2e-16 ***
## Lunghezza -3.060e-01 2.203e+00 -0.139 0.8895
## Cranio -4.755e+00 3.192e+00 -1.490 0.1365
## SessoM 7.324e+01 1.120e+01 6.537 7.59e-11 ***
## Lunghezza:Cranio 3.157e-02 6.531e-03 4.835 1.41e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.5 on 2491 degrees of freedom
## Multiple R-squared: 0.7296, Adjusted R-squared: 0.7289
## F-statistic: 1120 on 6 and 2491 DF, p-value: < 2.2e-16Qui a “rompersi” sono le variabili Lunghezza e Cranio. R2 adj non ha avuto un aumento significativo.
mod3 <- update(mod2_1,~.+I(Lunghezza^2))
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Lunghezza^2), data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1169.62 -181.77 -12.79 163.77 1786.03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 212.288548 723.852095 0.293 0.769336
## N.gravidanze 14.085464 4.266175 3.302 0.000975 ***
## Gestazione 42.551398 3.876629 10.976 < 2e-16 ***
## Lunghezza -20.267001 3.162718 -6.408 1.76e-10 ***
## Cranio 10.651783 0.418894 25.428 < 2e-16 ***
## SessoM 69.968733 11.038797 6.338 2.75e-10 ***
## I(Lunghezza^2) 0.031655 0.003267 9.690 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.7 on 2491 degrees of freedom
## Multiple R-squared: 0.7369, Adjusted R-squared: 0.7363
## F-statistic: 1163 on 6 and 2491 DF, p-value: < 2.2e-16L’aggiunta del termine quadratico di Lunghezza non aumenta in maniera significativa l’R2 adj. Quindi per il rasoio di Occam preferiamo la soluzione più semplice
summary(update(mod2_1,~.+I(Cranio^2)))##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Cranio^2), data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1138.6 -179.4 -14.8 163.4 2622.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.10118 1151.77280 0.073 0.94180
## N.gravidanze 12.76356 4.31259 2.960 0.00311 **
## Gestazione 38.90540 3.93291 9.892 < 2e-16 ***
## Lunghezza 10.48745 0.30157 34.776 < 2e-16 ***
## Cranio -31.79371 7.16973 -4.434 9.63e-06 ***
## SessoM 73.10236 11.16590 6.547 7.11e-11 ***
## I(Cranio^2) 0.06262 0.01059 5.915 3.77e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 272.8 on 2491 degrees of freedom
## Multiple R-squared: 0.7308, Adjusted R-squared: 0.7301
## F-statistic: 1127 on 6 and 2491 DF, p-value: < 2.2e-16L’aggiunta del termine quadratico di Cranio non aumenta in maniera significativa l’R2 adj. Quindi per il rasoio di Occam preferiamo la soluzione più semplice
Sinora il modello migliore è il mod2_1
Proviamo a semplificare ulteriormente il modello 2_1:
mod4<-lm(Peso~Gestazione+Lunghezza+Cranio)
summary(mod4)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1105.68 -183.52 -12.69 166.41 2622.95
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6778.1583 135.7089 -49.946 <2e-16 ***
## Gestazione 31.7554 3.8247 8.303 <2e-16 ***
## Lunghezza 10.4210 0.3022 34.486 <2e-16 ***
## Cranio 10.7911 0.4285 25.186 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 277.8 on 2494 degrees of freedom
## Multiple R-squared: 0.7207, Adjusted R-squared: 0.7203
## F-statistic: 2145 on 3 and 2494 DF, p-value: < 2.2e-16mod4 rispetto a mod2 ha un r2 aggiustato quasi identico e un’altissima significatività dei regressori quindi per il rasoio di Occam sinora è preferibile il mod4
Riaggiungere singolarmente fumatrici, sesso, anni madre o n gravidanze migliora di pochi millesimi il r2 adj quindi sempre per il rasoio di occam non li riaggiungiamo e scegliamo il modello più semplice ovvero il mod 4
Fumatrici
summary(update(mod4,~.+Fumatrici))##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Fumatrici)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1106.12 -184.29 -13.02 167.10 2620.16
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6778.3648 135.7173 -49.945 <2e-16 ***
## Gestazione 31.9396 3.8313 8.337 <2e-16 ***
## Lunghezza 10.4097 0.3025 34.412 <2e-16 ***
## Cranio 10.7900 0.4285 25.181 <2e-16 ***
## Fumatrici1 -23.3289 27.8768 -0.837 0.403
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 277.8 on 2493 degrees of freedom
## Multiple R-squared: 0.7207, Adjusted R-squared: 0.7203
## F-statistic: 1608 on 4 and 2493 DF, p-value: < 2.2e-16R2 adj identica e non significatività della variabile Fumatrici
Sesso
summary(update(mod4,~.+Sesso))##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1138.1 -184.4 -17.4 163.3 2626.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6651.6732 135.5952 -49.055 < 2e-16 ***
## Gestazione 31.3262 3.7884 8.269 < 2e-16 ***
## Lunghezza 10.2024 0.3009 33.909 < 2e-16 ***
## Cranio 10.6706 0.4247 25.126 < 2e-16 ***
## SessoM 79.1027 11.2205 7.050 2.31e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 275.1 on 2493 degrees of freedom
## Multiple R-squared: 0.7261, Adjusted R-squared: 0.7257
## F-statistic: 1652 on 4 and 2493 DF, p-value: < 2.2e-16R2 adj aumentata solo di mezzo punto
Anni madre
summary(update(mod4,~.+Anni.madre))##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Anni.madre)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1132.88 -185.78 -14.93 167.11 2614.28
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6855.0490 141.4871 -48.450 <2e-16 ***
## Gestazione 32.7288 3.8565 8.487 <2e-16 ***
## Lunghezza 10.4319 0.3021 34.534 <2e-16 ***
## Cranio 10.7189 0.4299 24.934 <2e-16 ***
## Anni.madre 2.0609 1.0792 1.910 0.0563 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 277.6 on 2493 degrees of freedom
## Multiple R-squared: 0.7211, Adjusted R-squared: 0.7206
## F-statistic: 1611 on 4 and 2493 DF, p-value: < 2.2e-16R2 adj quasi identica e Anni.madre al limite della significatività
Numero gravidanze
summary(update(mod4,~.+N.gravidanze))##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1118.36 -181.20 -12.15 167.93 2636.44
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6808.8062 135.8418 -50.123 < 2e-16 ***
## Gestazione 32.8947 3.8360 8.575 < 2e-16 ***
## Lunghezza 10.4643 0.3020 34.651 < 2e-16 ***
## Cranio 10.6486 0.4302 24.752 < 2e-16 ***
## N.gravidanze 13.5089 4.3800 3.084 0.00206 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 277.3 on 2493 degrees of freedom
## Multiple R-squared: 0.7217, Adjusted R-squared: 0.7213
## F-statistic: 1616 on 4 and 2493 DF, p-value: < 2.2e-16R2 adj quasi identica
Come vediamo sotto, anche aggiungere dei termini che tengono conto delle interazioni al modello più semplice non aiuta
summary(lm(Peso~Gestazione*Lunghezza+Cranio))##
## Call:
## lm(formula = Peso ~ Gestazione * Lunghezza + Cranio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1137.31 -180.01 -7.48 164.76 2637.67
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.559e+03 9.269e+02 -1.682 0.0928 .
## Gestazione -1.086e+02 2.494e+01 -4.352 1.4e-05 ***
## Lunghezza -1.069e+00 2.041e+00 -0.524 0.6006
## Cranio 1.103e+01 4.279e-01 25.783 < 2e-16 ***
## Gestazione:Lunghezza 3.030e-01 5.323e-02 5.692 1.4e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 276 on 2493 degrees of freedom
## Multiple R-squared: 0.7242, Adjusted R-squared: 0.7238
## F-statistic: 1637 on 4 and 2493 DF, p-value: < 2.2e-16Qui si rompe la variabile Lunghezza e l’R2 adj è quasi identica
summary(lm(Peso~Gestazione+Lunghezza*Cranio))##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza * Cranio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1108.76 -181.95 -7.75 165.45 2872.86
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.377e+03 1.026e+03 -1.342 0.180
## Gestazione 3.806e+01 3.985e+00 9.550 < 2e-16 ***
## Lunghezza -1.257e+00 2.219e+00 -0.566 0.571
## Cranio -6.116e+00 3.212e+00 -1.904 0.057 .
## Lunghezza:Cranio 3.490e-02 6.571e-03 5.311 1.19e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 276.3 on 2493 degrees of freedom
## Multiple R-squared: 0.7238, Adjusted R-squared: 0.7233
## F-statistic: 1633 on 4 and 2493 DF, p-value: < 2.2e-16Qui si rompe la variabile Lunghezza e Cranio e l’R2 adj è quasi identica
summary(lm(Peso~Gestazione*Lunghezza*Cranio))##
## Call:
## lm(formula = Peso ~ Gestazione * Lunghezza * Cranio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1109.52 -180.55 -10.26 164.14 2583.63
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.596e+04 8.389e+03 4.287 1.88e-05 ***
## Gestazione -1.179e+03 2.380e+02 -4.953 7.81e-07 ***
## Lunghezza -7.186e+01 1.975e+01 -3.638 0.00028 ***
## Cranio -1.216e+02 2.764e+01 -4.400 1.13e-05 ***
## Gestazione:Lunghezza 2.330e+00 5.322e-01 4.377 1.25e-05 ***
## Gestazione:Cranio 3.720e+00 7.606e-01 4.891 1.07e-06 ***
## Lunghezza:Cranio 2.537e-01 6.253e-02 4.058 5.11e-05 ***
## Gestazione:Lunghezza:Cranio -7.129e-03 1.656e-03 -4.304 1.74e-05 ***
## ---
## 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.5 on 7 and 2490 DF, p-value: < 2.2e-16R2 adj quasi identica
Proviamo ad aggiungere il termine quadratico di Lunghezza
mod3_1 <- update(mod4,~.+I(Lunghezza^2))
summary(mod3_1)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + I(Lunghezza^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1128.48 -178.32 -8.38 164.12 1738.23
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 378.389165 730.115696 0.518 0.604
## Gestazione 42.098677 3.892279 10.816 < 2e-16 ***
## Lunghezza -21.206080 3.186356 -6.655 3.46e-11 ***
## Cranio 10.911086 0.420425 25.952 < 2e-16 ***
## I(Lunghezza^2) 0.032782 0.003288 9.969 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 272.4 on 2493 degrees of freedom
## Multiple R-squared: 0.7314, Adjusted R-squared: 0.7309
## F-statistic: 1697 on 4 and 2493 DF, p-value: < 2.2e-16R2 adj aumentata di un punto.
Sinora in base ai dati dell’R2 adj e il rasoio di Occam il miglior modello è il mod4.
Visto però che l’azienda richiede di fare delle stime per, ad esempio, il peso di una neonata considerando una madre alla terza gravidanza, non fumatrice, che partorirà alla 39esima settimana creiamo anche un modello che includa esattamente queste variabili oltre a quelle fortemente significative (Lunghezza e Cranio):
mod_var_richieste <- lm(Peso~N.gravidanze+Fumatrici+Gestazione+Lunghezza+Cranio+Sesso)
summary(mod_var_richieste)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Sesso)
##
## 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 **
## Fumatrici1 -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-16Se non includessimo Lunghezza e Cranio l’R2 adj crollerebbe a meno di 0.38!
summary(lm(Peso~N.gravidanze+Fumatrici+Gestazione+Sesso))##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Sesso)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1498.59 -274.99 -14.02 265.74 1887.17
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3152.312 175.346 -17.978 < 2e-16 ***
## N.gravidanze 24.214 6.511 3.719 0.000204 ***
## Fumatrici1 -111.697 41.521 -2.690 0.007189 **
## Gestazione 162.523 4.499 36.127 < 2e-16 ***
## SessoM 165.433 16.720 9.894 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 413.8 on 2493 degrees of freedom
## Multiple R-squared: 0.3804, Adjusted R-squared: 0.3794
## F-statistic: 382.6 on 4 and 2493 DF, p-value: < 2.2e-16Proviamo a vedere quale modello ci suggerisce R con la libreria MASS
stepwise.mod <- stepAIC(mod1,
direction = "both",
k=log(n)) #criterio BIC## Start: AIC=28116.81
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 37687 187231953 28110
## - Fumatrici 1 83067 187277334 28110
## - Ospedale 2 710947 187905214 28111
## - N.gravidanze 1 422758 187617025 28115
## <none> 187194266 28117
## - Sesso 1 3615113 190809379 28157
## - Gestazione 1 5473810 192668077 28181
## - Cranio 1 45659609 232853876 28654
## - Lunghezza 1 87455591 274649858 29067
##
## Step: AIC=28109.49
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 83954 187315907 28103
## - Ospedale 2 717552 187949505 28103
## <none> 187231953 28110
## - N.gravidanze 1 602588 187834541 28110
## + Anni.madre 1 37687 187194266 28117
## - Sesso 1 3622317 190854270 28150
## - Gestazione 1 5438554 192670508 28173
## - Cranio 1 45914823 233146777 28650
## - Lunghezza 1 87445208 274677161 29059
##
## Step: AIC=28102.79
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Ospedale +
## Sesso
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 726146 188042054 28097
## - N.gravidanze 1 581341 187897249 28103
## <none> 187315907 28103
## + Fumatrici 1 83954 187231953 28110
## + Anni.madre 1 38573 187277334 28110
## - Sesso 1 3607033 190922941 28143
## - Gestazione 1 5375559 192691466 28166
## - Cranio 1 45954185 233270092 28643
## - Lunghezza 1 87854020 275169927 29056
##
## Step: AIC=28096.81
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 188042054 28097
## - N.gravidanze 1 621053 188663107 28097
## + Ospedale 2 726146 187315907 28103
## + Fumatrici 1 92548 187949505 28103
## + Anni.madre 1 45366 187996688 28104
## - Sesso 1 3650790 191692844 28137
## - Gestazione 1 5477493 193519547 28161
## - Cranio 1 46098547 234140601 28637
## - Lunghezza 1 87532691 275574744 29044summary(stepwise.mod)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = neonati)
##
## 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-16Il modello che si ottiene dal metodo stepwise è identico al mod2_1, ha 5 variabili ma ha un R2 adj di poco più di mezzo punto percentuale maggiore del nostro mod4 quindi scegliamo il mod4
Altro modo per valutare se aggiungere o meno un regressore è utilizzare anova() ovvero l’analisi della varianza
anova(mod4,stepwise.mod,mod2,mod1)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 3: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Ospedale +
## Sesso
## Model 4: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Ospedale + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2494 192424285
## 2 2492 188042054 2 4382231 29.1221 3.149e-13 ***
## 3 2490 187315907 2 726146 4.8256 0.008097 **
## 4 2488 187194266 2 121641 0.8084 0.445703
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Secondo il test anova() tra questi modelli è preferibile il mod2
Confrontiamo ora altri modelli
anova(mod3_1,mod3)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + I(Lunghezza^2)
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## I(Lunghezza^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2493 185047535
## 2 2491 181211596 2 3835940 26.365 4.67e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1il mod3 è migliore del mod3_1
anova(stepwise.mod, mod_var_richieste)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 188042054
## 2 2491 187949505 1 92548 1.2266 0.2682lo stepwise.mod è migliore del mod_var_richieste
anova(stepwise.mod,mod3)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## I(Lunghezza^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 188042054
## 2 2491 181211596 1 6830458 93.894 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1il mod3 è migliore dello stepwise.mod
Quindi secondo il test anova i migliore modelli sono il mod3 e il mod2. Non possiamo confrontare questi due modelli in quanto contengono variabili differenti.
Possiamo creare però un modello ad hoc come il mod2 con aggiunta del termine quadratico Lunghezza per confrontare entrambi con esso e vedere se si stabilisce in maniera indiretta un ordine.
mod2_plus <- update(mod2,~.+I(Lunghezza^2))
anova(mod2,mod2_plus)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Ospedale +
## Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Ospedale +
## Sesso + I(Lunghezza^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2490 187315907
## 2 2489 180500632 1 6815275 93.979 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mod3,mod2_plus)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## I(Lunghezza^2)
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Ospedale +
## Sesso + I(Lunghezza^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2491 181211596
## 2 2489 180500632 2 710964 4.9019 0.007504 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Neanche questo confronto è risolutivo, il test anova ci dice che il modello migliore in assoluto è proprio quest’ultimo mod2_plus
A questo punto facciamo anche un summary per questo modello
summary(mod2_plus)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Ospedale + Sesso + I(Lunghezza^2), data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1153.52 -179.76 -12.29 164.18 1767.87
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 203.769407 722.726866 0.282 0.77801
## N.gravidanze 13.695246 4.261771 3.214 0.00133 **
## Gestazione 42.251614 3.871801 10.913 < 2e-16 ***
## Lunghezza -20.211869 3.158064 -6.400 1.85e-10 ***
## Cranio 10.636865 0.418287 25.430 < 2e-16 ***
## Ospedaleosp2 -9.543044 13.214305 -0.722 0.47025
## Ospedaleosp3 29.979854 13.263177 2.260 0.02388 *
## SessoM 69.531232 11.022732 6.308 3.33e-10 ***
## I(Lunghezza^2) 0.031622 0.003262 9.694 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.3 on 2489 degrees of freedom
## Multiple R-squared: 0.738, Adjusted R-squared: 0.7371
## F-statistic: 876.2 on 8 and 2489 DF, p-value: < 2.2e-16Il modello preferito dal test anova ha un R2 adj circa uguale a quello del mod3, di 1.5 punti maggiore dell’R2 adj del mod4. Parliamo però di 6 regressori e un termine quadratico contro 3 regressori, una bella differenza in termini di complessità a fronte di un aumento relativamente piccolo di R2 adj.
Selezione del modello ottimale
Proviamo ad utilizzare altri criteri di valutazione, AIC e BIC. Il modello con i criteri di informazione più bassi è di solito il modello migliore
AIC(mod1,mod2_plus,mod2,mod3,mod3_1,mod4,mod_var_richieste,stepwise.mod) #mod2_plus## df AIC
## mod1 11 35149.60
## mod2_plus 10 35056.64
## mod2 9 35147.22
## mod3 8 35062.46
## mod3_1 6 35110.79
## mod4 5 35206.43
## mod_var_richieste 8 35153.66
## stepwise.mod 7 35152.89BIC(mod1,mod2_plus,mod2,mod3,mod3_1,mod4,mod_var_richieste,stepwise.mod) #mod3## df BIC
## mod1 11 35213.65
## mod2_plus 10 35114.87
## mod2 9 35199.63
## mod3 8 35109.04
## mod3_1 6 35145.72
## mod4 5 35235.55
## mod_var_richieste 8 35200.24
## stepwise.mod 7 35193.65Per AIC il modello migliore è il mod2_plus Per BIC invece è il mod3.
A questo punto, considerando che AIC tende a preferire modelli sovraparametrizzati prendiamo in considerazione il risultato di BIC.
Controlliamo la multicollinearità
car::vif(mod1)## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.190224 1 1.090974
## N.gravidanze 1.188419 1 1.090146
## Fumatrici 1.007051 1 1.003519
## Gestazione 1.695637 1 1.302166
## Lunghezza 2.081857 1 1.442864
## Cranio 1.629339 1 1.276456
## Ospedale 1.003827 2 1.000955
## Sesso 1.040740 1 1.020167car::vif(mod2_plus)## GVIF Df GVIF^(1/(2*Df))
## N.gravidanze 1.026144 1 1.012988
## Gestazione 1.802961 1 1.342744
## Lunghezza 238.051106 1 15.428905
## Cranio 1.626150 1 1.275206
## Ospedale 1.002774 2 1.000693
## Sesso 1.046278 1 1.022877
## I(Lunghezza^2) 230.061562 1 15.167780car::vif(mod2)## GVIF Df GVIF^(1/(2*Df))
## N.gravidanze 1.024531 1 1.012191
## Gestazione 1.670812 1 1.292599
## Lunghezza 2.077805 1 1.441459
## Cranio 1.624925 1 1.274726
## Ospedale 1.002651 2 1.000662
## Sesso 1.040428 1 1.020014car::vif(mod2_1)## N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1.023462 1.669779 2.075747 1.624568 1.040184car::vif(mod3)## N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1.025055 1.801815 238.007682 1.625780 1.046053
## I(Lunghezza^2)
## 230.033474car::vif(mod3_1)## Gestazione Lunghezza Cranio I(Lunghezza^2)
## 1.780167 236.760730 1.605028 228.434188car::vif(mod4)## Gestazione Lunghezza Cranio
## 1.653674 2.048591 1.603713car::vif(mod_var_richieste)## N.gravidanze Fumatrici Gestazione Lunghezza Cranio Sesso
## 1.026101 1.006621 1.676199 2.079728 1.624705 1.040400car::vif(stepwise.mod)## N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1.023462 1.669779 2.075747 1.624568 1.040184Non ci sono problemi di multicollinearità (a parte quella prevedibile tra Lunghezza e il suo termine quadratico I(Lunghezza^2))
Tra il mod3 (6 regressori e un termine quadratico) e il mod4 (3 regressori) abbiamo sempre un R2 adj che differisce di circa 1.5 punti e non vale la pena aggiungere 3 regressori e un termine quadratico per avere solo un punto percentuale di varianza in più spiegata, nonostante alcuni test sopra riportati valutino come migliore il mod3 o il mod2_plus (7 regressori e un termine quadratico)
Tenendo in considerazione il mod4 che risulta essere il migliore in termini di semplicità e rasoio di Occam per i valori di R2 adj, il mod3 segnalato come migliore dal test BIC e il mod_var_richieste andiamo a valutare l’MSE, l’RMSE e il MAE
MSE è la media dei quadrati degli errori tra i valori osservati e quelli predetti da un modello. RMSE è la radice quadrata della media dei quadrati degli errori tra i valori osservati e quelli predetti da un modello quindi è nella stessa unità di misura della variabile peso. Sono entrambi sensibili agli outlier MAE è la media degli errori assoluti tra i valori osservati e quelli predetti. Indica, in media, di quanto il modello si discosta dai valori osservati, nella stessa unità della variabile dipendente. Il MAE è meno sensibile agli outlier
MSE_mod4 <- mean(mod4$residuals^2)
RMSE_mod4 <- sqrt(mean(mod4$residuals^2))
MAE_mod4 <- mean(abs(mod4$residuals))
MSE_mod3 <- mean(mod3$residuals^2)
RMSE_mod3 <- sqrt(mean(mod3$residuals^2))
MAE_mod3 <- mean(abs(mod3$residuals))
MSE_mod_var_richieste <- mean(mod_var_richieste$residuals^2)
RMSE_mod_var_richieste <- sqrt(mean(mod_var_richieste$residuals^2))
MAE_mod_var_richieste <- mean(abs(mod_var_richieste$residuals))
tabella_risultati <- data.frame(
Modello = c("mod4", "mod3", "mod_var_richieste"),
MSE = c(MSE_mod4, MSE_mod3, MSE_mod_var_richieste),
RMSE = c(RMSE_mod4, RMSE_mod3, RMSE_mod_var_richieste),
MAE = c(MAE_mod4, MAE_mod3, MAE_mod_var_richieste)
)
kable(tabella_risultati, caption = "Metriche di Valutazione dei Modelli")| Modello | MSE | RMSE | MAE |
|---|---|---|---|
| mod4 | 77031.34 | 277.5452 | 214.1401 |
| mod3 | 72542.67 | 269.3375 | 208.8608 |
| mod_var_richieste | 75239.99 | 274.2991 | 210.9321 |
I valori di MSE e RMSE del mod4 risultano essere maggiori e questo significa un errore maggiore
Analisi della qualità del modello
mod4
par(mfrow=c(2,2))
plot(mod4, pch=20, col="royalblue3")
Residuals vs Fitted
I punti devono presentarsi sparsi casualmente attorno alla media di zero.
La media di zero sembra essere rispettata per valori più alti mentre i valori più bassi tendono ad essere maggiori di 0. Le osservazioni 1549, 155 e 1305 sono fuori dalla nuvola
Q-Q Residuals
Vengono messi in relazione i residui con i quantili di una distribuzione normale.
Se i punti giacciono pressappoco sulla bisettrice del grafico allora i residui seguono la distribuzione normale.
In questo caso sembra che i punti siano correttamente disposti sopra la bisettrice, a parte qualche problemino nella coda inferiore e soprattutto quella superiore. Ci sono di nuovo le osservazioni 1549, 155 e 1305 che sono più lontane dalle altre, principalmente la 1549.
Scale-Location
Non si devono visualizzare dei patterns ma soltanto una nuvola casuale grossomodo orizzontale attorno a un valore di y che indicherà appunto una varianza costante.
Anche qui si nota una curvatura per i valori più piccoli.
Residuals vs Leverage
Si visualizzano i potenziali valori influenti, ovvero quei valori molto leverage, molto outliers o una combinazione di entrambe le cose.
In generale quando alcuni residui mostrano una distanza di Cook che supera queste linee tratteggiate, vengono considerati residui di osservazioni potenzialmente influenti sulle stime di regressione.
La soglia a 0,5 è quella di avvertimento mentre la soglia fissata a 1 è quella di allarme.
In questo caso l’osservazione 1549 è addirittura oltre la soglia di allarme.
L’osservazione 1549 si era distinta già nell’analisi descrittiva
Calcoliamo ora i valori di leva per il mod4, identifichiamo il valore soglia che classificherà i valori come valori di leva o meno e andiamo infine a visualizzare la soglia anche nel grafico con una linea orizzontale azzurra
lev_mod4<-hatvalues(mod4)
plot(lev_mod4, pch=20, col="royalblue3")
p_mod4 = sum(lev_mod4)
soglia_mod4 = 2*p_mod4/n
abline(h=soglia_mod4,col="turquoise3",lwd=1.5)
Ci sono molti valori leverage ovvero oltre il valore soglia definito.
Possiamo identificare questi valori, ma essendo così tanti guardiamone solo il numero.
length(lev_mod4[lev_mod4>soglia_mod4])## [1] 173Per capire meglio l’impatto di questi leverage andiamo avanti con l’analisi dei residui
plot(rstudent(mod4), pch=20, col="royalblue3")
abline(h=c(-2,2),col="turquoise3",lwd=1.5)
Ci sono molti punti che superano i valori soglia di meno 2 e 2
Andiamo ad effettuare anche il test statistico per gli outlier
car::outlierTest(mod4)## rstudent unadjusted p-value Bonferroni p
## 1549 9.866111 1.5024e-22 3.7531e-19
## 155 5.120997 3.2715e-07 8.1722e-04
## 1305 4.567113 5.1853e-06 1.2953e-02Le tre osservazioni che vedevamo dai grafici mostrano avere qualche problema anche da questo test.
Andiamo anche a visualizzare la distanza di Cook
cook_mod4<-cooks.distance(mod4)
plot(cook_mod4,ylim = c(0,1.2), pch=20, col="royalblue3")
Come avevamo visto il valore massimo della distanza di Cook è maggiore di 1, circa 1.2. Le altre osservazioni sono tutte sotto la soglia di attenzione di 0.5
Andiamo a verificare esattamente a quanto ammonta il massimo della distanza di Cook
max(cook_mod4)## [1] 1.195621Il valore massimo raggiunto delle distanze di cook è di circa 1.2 quindi c’è un valore particolarmente influente sulle stime di regressione.
par(mfrow=c(1,2))
plot(residuals(mod4), pch=20, col="royalblue3")
abline(h=mean(residuals(mod4)), col="turquoise3",lwd=1.5)
plot(density(residuals(mod4)), col="royalblue3")
Possiamo vedere che la distribuzione dei residui assomiglia abbastanza alla normale ma ha una coda a sinistra più schiacciata e una piccola doppia “punta” centrale
Andiamo a fare un test di Shapiro Wilk, un test di normalità in cui l’ipotesi nulla è proprio la distribuzione normale della variabile di riferimento
shapiro.test(mod4$residuals)##
## Shapiro-Wilk normality test
##
## data: mod4$residuals
## W = 0.97559, p-value < 2.2e-16p value molto piccolo quindi si rifiuta l’ipotesi di normalità
Andiamo a saggiare queste altre due ipotesi
RESIDUI A VARIANZA COSTANTE (omoschedasticità) –> TEST DI BREUSCH PAGAN
RESIDUI NON CORRELATI TRA DI LORO –> TEST DI DURBIN WATSON
bptest(mod4)##
## studentized Breusch-Pagan test
##
## data: mod4
## BP = 91.818, df = 3, p-value < 2.2e-16il p value è molto piccolo quindi si rifiuta l’ipotesi di omoschedasticità (residui a varianza non costante)
dwtest(mod4)##
## Durbin-Watson test
##
## data: mod4
## DW = 1.9617, p-value = 0.1687
## alternative hypothesis: true autocorrelation is greater than 0il p value è di 0.1687 (valore più grande della soglia dello 0.05) quindi non si rifiuta l’ipotesi di indipendenza
Visualizziamo questa osservazione 1549 che è quella con la distanza di Cook maggiore di 1
oss_1549 <- neonati[1549, ]
oss_1549## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1551 35 1 0 38 4370 315 374
## Tipo.parto Ospedale Sesso Gestazione_CL
## 1551 Nat osp3 F A termineEscludendo l’osservazione 1549 verifichiamo alcuni valori tenendo presente che la nostra osservazione riporta: Gestazione: 38 Peso: 4370 Lunghezza: 315 Cranio:374
Minima lunghezza del neonato per una gestazione maggiore o uguale a 38
neonati[-1549,] %>%
filter(Gestazione>=38)%>%
summarize(min(Lunghezza))## min(Lunghezza)
## 1 410Massima settimana di gestazione per una lunghezza del neonato compresa tra 310 e 320 mm
neonati[-1549, ] %>%
filter(Lunghezza>=310&Lunghezza<=320)%>%
summarize(max(Gestazione))## max(Gestazione)
## 1 28Massimo peso per una lunghezza del neonato compresa tra 310 e 320 mm
neonati[-1549, ] %>%
filter(Lunghezza>=310&Lunghezza<=320)%>%
summarize(max(Peso))## max(Peso)
## 1 980Massima dimensione del cranio del neonato per una lunghezza compresa tra 310 e 320 mm
neonati[-1549,] %>%
filter(Lunghezza>=310&Lunghezza<=320)%>%
summarize(max(Cranio))## max(Cranio)
## 1 265Lunghezza media per un peso tra i 4300g e i 4400g
neonati[-1549,] %>%
filter(Peso >= 4300 & Peso <= 4400) %>%
summarise(mean(Lunghezza))## mean(Lunghezza)
## 1 534Da questi 4 valori possiamo vedere quanto si discosti questa osservazione dalle altre del dataset.
Proviamo a creare un dataset in cui escludiamo l’osservazione 1549 per vedere come varia il modello, potrebbe anche essere un errore di battitura.
neonati_no_1549 <- neonati[-1549,]Proviamo a creare un modello identico a mod4 applicandolo a questo nuovo dataset
mod4_no_1549 <- lm(Peso~Gestazione+Lunghezza+Cranio, data = neonati_no_1549)
summary(mod4_no_1549)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio, data = neonati_no_1549)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1121.19 -181.67 -12.06 166.83 1452.03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6778.9476 133.1613 -50.908 <2e-16 ***
## Gestazione 28.9674 3.7635 7.697 2e-14 ***
## Lunghezza 11.0595 0.3035 36.441 <2e-16 ***
## Cranio 10.1808 0.4249 23.958 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 272.6 on 2493 degrees of freedom
## Multiple R-squared: 0.7307, Adjusted R-squared: 0.7304
## F-statistic: 2255 on 3 and 2493 DF, p-value: < 2.2e-16summary(mod4)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1105.68 -183.52 -12.69 166.41 2622.95
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6778.1583 135.7089 -49.946 <2e-16 ***
## Gestazione 31.7554 3.8247 8.303 <2e-16 ***
## Lunghezza 10.4210 0.3022 34.486 <2e-16 ***
## Cranio 10.7911 0.4285 25.186 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 277.8 on 2494 degrees of freedom
## Multiple R-squared: 0.7207, Adjusted R-squared: 0.7203
## F-statistic: 2145 on 3 and 2494 DF, p-value: < 2.2e-16Il summary del nuovo modello mostra valori buoni, rispetto al mod4 c’è un miglioramento nei residui e nel r2 adj (1%). I coefficienti tra i due modelli variano ma non di tantissimo.
mod4_no_1549
par(mfrow=c(2,2))
plot(mod4_no_1549, pch=20, col="royalblue3")
La situazione è rimasta analoga alla precedente senza ovviamente l’osservazione 1549 Restano in evidenza principalmente le osservazioni 1305 e 155. Compare anche l’osservazione 1399.
Calcoliamo ora i valori di leva per il mod4_no_1549, identifichiamo il valore soglia che classificherà i valori come valori di leva o meno e andiamo infine a visualizzare la soglia anche nel grafico con una linea orizzontale azzurra
lev_mod4_no_1549<-hatvalues(mod4_no_1549)
plot(lev_mod4_no_1549, pch=20, col="royalblue3")
p_mod4_no_1549 = sum(lev_mod4_no_1549)
n_mod4_no_1549 <- nrow(neonati_no_1549)
soglia_mod4_no_1549 = 2*p_mod4_no_1549/n_mod4_no_1549
abline(h=soglia_mod4_no_1549,col="turquoise3",lwd=1.5)
Ci sono molti valori leverage ovvero oltre il valore soglia definito.
Possiamo identificare questi valori, ma essendo così tanti guardiamone solo il numero.
length(lev_mod4_no_1549[lev_mod4_no_1549>soglia_mod4_no_1549])## [1] 176Per capire meglio l’impatto di questi leverage andiamo avanti con l’analisi dei residui
plot(rstudent(mod4_no_1549), pch=20, col="royalblue3")
abline(h=c(-2,2),col="turquoise3",lwd=1.5)
Ci sono molti punti che superano i valori soglia di meno 2 e 2
Andiamo ad effettuare anche il test statistico per gli outlier
car::outlierTest(mod4_no_1549)## rstudent unadjusted p-value Bonferroni p
## 155 5.373143 8.4549e-08 0.00021112
## 1306 4.671130 3.1546e-06 0.00787710Le due osservazioni che vedevamo dai grafici mostrano avere qualche problema anche da questo test.
Andiamo anche a visualizzare la distanza di Cook
cook_mod4_no_1549<-cooks.distance(mod4_no_1549)
plot(cook_mod4_no_1549,ylim = c(0,1), pch=20, col="royalblue3")
Come avevamo visto dai grafici non abbiamo valori al di sopra della soglia di attenzione della distanza di Cook pari a 0.5.
Andiamo a verificare esattamente a quanto ammonta il massimo della distanza di Cook
max(cook_mod4_no_1549)## [1] 0.09804879il valore massimo raggiunto delle distanze di cook è di circa 0.1 quindi non ci sono valori particolarmente influenti sulle stime di regressione.
par(mfrow=c(1,2))
plot(residuals(mod4_no_1549), pch=20, col="royalblue3")
abline(h=mean(residuals(mod4_no_1549)), col="turquoise3",lwd=1.5)
plot(density(residuals(mod4_no_1549)), col="royalblue3")
Possiamo vedere che la distribuzione dei residui assomiglia abbastanza alla normale ma ha una coda a sinistra più schiacciata e più corta rispetto a quella di destra e una piccola doppia “punta” centrale
Andiamo a fare un test di Shapiro Wilk, un test di normalità in cui l’ipotesi nulla è proprio la distribuzione normale della variabile di riferimento
shapiro.test(mod4_no_1549$residuals)##
## Shapiro-Wilk normality test
##
## data: mod4_no_1549$residuals
## W = 0.98955, p-value = 1.559e-12p value molto piccolo quindi si rifiuta l’ipotesi di normalità
Andiamo a saggiare queste altre due ipotesi
RESIDUI A VARIANZA COSTANTE (omoschedasticità) –> TEST DI BREUSCH PAGAN
RESIDUI NON CORRELATI TRA DI LORO –> TEST DI DURBIN WATSON
bptest(mod4_no_1549)##
## studentized Breusch-Pagan test
##
## data: mod4_no_1549
## BP = 10.634, df = 3, p-value = 0.01388il p value è minore della soglia di 0.05 quindi si rifiuta l’ipotesi di omoschedasticità (residui a varianza non costante)
dwtest(mod4_no_1549)##
## Durbin-Watson test
##
## data: mod4_no_1549
## DW = 1.9623, p-value = 0.173
## alternative hypothesis: true autocorrelation is greater than 0il p value è di 0.173 (valore più grande della soglia dello 0.05) quindi non si rifiuta l’ipotesi di indipendenza
mod3
par(mfrow=c(2,2))
plot(mod3, pch=20, col="royalblue3")
Residuals vs Fitted
La media di zero sembra essere rispettata per valori più alti mentre i valori più bassi tendono ad essere maggiori di 0. Le osservazioni 1551, 155 e 1306 sono fuori dalla nuvola
Q-Q Residuals
In questo caso sembra che i punti siano correttamente disposti sopra la bisettrice, a parte qualche problemino nella coda inferiore e soprattutto quella superiore. Ci sono di nuovo le osservazioni 1551, 155 e 1306 che sono più lontane dalle altre, principalmente la 1551.
Scale-Location
Anche qui si nota una curvatura per i valori più piccoli e le solite 3 osservazioni più lontane.
Residuals vs Leverage
Qui notiamo l’osservazione 1551 oltre la soglia di avvertimento: è un’osservazione che si era già distinta durante l’analisi preliminare dei dati, per un valore di peso alto rispetto alla settimana di gestazione
Calcoliamo ora i valori di leva per lo mod3, identifichiamo il valore soglia che classificherà i valori come valori di leva o meno e andiamo infine a visualizzare la soglia anche nel grafico con una linea orizzontale azzurra
lev_mod3<-hatvalues(mod3)
plot(lev_mod3, pch=20, col="royalblue3")
p_mod3 = sum(lev_mod3)
soglia_mod3 = 2*p_mod3/n
abline(h=soglia_mod3,col="turquoise3",lwd=1.5)
Ci sono molti valori leverage ovvero oltre il valore soglia definito.
Possiamo identificare questi valori, ma essendo così tanti guardiamone solo il numero.
length(lev_mod3[lev_mod3>soglia_mod3])## [1] 142Per capire meglio l’impatto di questi leverage andiamo avanti con l’analisi dei residui
plot(rstudent(mod3), pch=20, col="royalblue3")
abline(h=c(-2,2),col="turquoise3",lwd=1.5)
Ci sono molti punti che superano i valori soglia di meno 2 e 2
Andiamo ad effettuare anche il test statistico per gli outlier
car::outlierTest(mod3)## rstudent unadjusted p-value Bonferroni p
## 1551 7.279956 4.4561e-13 1.1131e-09
## 1306 4.831997 1.4338e-06 3.5817e-03
## 155 4.746031 2.1919e-06 5.4754e-03
## 1399 -4.357506 1.3689e-05 3.4194e-02
## 1694 4.315725 1.6529e-05 4.1290e-02Le tre osservazioni che vedevamo dai grafici mostrano avere qualche problema anche da questo test. Se ne presentano però anche altre 2.
Andiamo anche a visualizzare la distanza di Cook
cook_mod3<-cooks.distance(mod3)
plot(cook_mod3,ylim = c(0,1.4), pch=20, col="royalblue3")
Come avevamo visto il valore massimo della distanza di Cook è maggiore di 1, quasi 1.4. Le altre osservazioni sono tutte sotto la soglia di attenzione di 0.5
Andiamo a verificare esattamente a quanto ammonta il massimo della distanza di Cook
max(cook_mod3)## [1] 1.36394il valore massimo raggiunto delle distanze di cook è di circa 1.36 quindi c’è un valore oltre la soglia di attenzione.
par(mfrow=c(1,2))
plot(residuals(mod3), pch=20, col="royalblue3")
abline(h=mean(residuals(mod3)), col="turquoise3",lwd=1.5)
plot(density(residuals(mod3)), col="royalblue3")
Possiamo vedere che la distribuzione dei residui assomiglia abbastanza alla normale ma ha una coda a destra meno schiacciata e più lunga
Andiamo a fare un test di Shapiro Wilk, un test di normalità in cui l’ipotesi nulla è proprio la distribuzione normale della variabile di riferimento
shapiro.test(mod3$residuals)##
## Shapiro-Wilk normality test
##
## data: mod3$residuals
## W = 0.98572, p-value = 3.811e-15p value molto piccolo quindi si rifiuta l’ipotesi di normalità
Andiamo a saggiare queste altre due ipotesi
RESIDUI A VARIANZA COSTANTE (omoschedasticità) –> TEST DI BREUSCH PAGAN
RESIDUI NON CORRELATI TRA DI LORO –> TEST DI DURBIN WATSON
bptest(mod3)##
## studentized Breusch-Pagan test
##
## data: mod3
## BP = 127.61, df = 6, p-value < 2.2e-16il p value è molto piccolo quindi si rifiuta l’ipotesi di omoschedasticità (residui a varianza non costante)
dwtest(mod3)##
## Durbin-Watson test
##
## data: mod3
## DW = 1.9478, p-value = 0.09601
## alternative hypothesis: true autocorrelation is greater than 0il p value è di 0.09601 (valore più grande della soglia dello 0.05) quindi non si rifiuta l’ipotesi di indipendenza
Visualizziamo l’osservazione 1551, ovvero quella oltre la soglia di attenzione di Cook
oss_1551 <- neonati[1551, ]
oss_1551## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1553 30 4 0 35 4520 520 360
## Tipo.parto Ospedale Sesso Gestazione_CL
## 1553 Nat osp2 F Lievemente pret.Escludendo l’osservazione 1551 verifichiamo alcuni valori tenendo presente che la nostra osservazione riporta:
Gestazione: 35
Peso: 4520
Lunghezza: 520
Cranio:360
Massimo peso del neonato per una gestazione minore o uguale a 35
neonati[-1551,] %>%
filter(Gestazione<=35)%>%
summarize(max(Peso))## max(Peso)
## 1 3140Massima dimensione cranio per una gestazione minore o uguale a 35
neonati[-1551,] %>%
filter(Gestazione<=35)%>%
summarize(max(Cranio))## max(Cranio)
## 1 379Massima lunghezza del neonato per una gestazione minore o uguale a 35
neonati[-1551,] %>%
filter(Gestazione<=35)%>%
summarize(max(Lunghezza))## max(Lunghezza)
## 1 500Peso medio dei neonati alti tra i 515mm e i 525mm
neonati[-1551,] %>%
filter(Lunghezza >= 515 & Lunghezza <= 525) %>%
summarise(mean(Peso))## mean(Peso)
## 1 3670.432Il peso ha un valore totalmente distante dai valori del resto del campione. Potrebbe essere un errore di battitura.
Proviamo a creare un dataset in cui escludiamo l’osservazione 1551 per vedere come varia il modello, potrebbe anche essere un errore di battitura.
neonati_no_1551 <- neonati[-1551,]Proviamo a creare un modello identico a mod3 applicandolo a questo nuovo dataset e leggiamone i due summary a confronto.
mod3_no_1551 <- lm(Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
Sesso + I(Lunghezza^2), data = neonati_no_1551)
summary(mod3_no_1551)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Lunghezza^2), data = neonati_no_1551)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1166.7 -180.5 -13.3 164.1 1779.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 221.05907 722.34369 0.306 0.75961
## N.gravidanze 13.49263 4.26086 3.167 0.00156 **
## Gestazione 43.42167 3.87707 11.200 < 2e-16 ***
## Lunghezza -20.35444 3.15621 -6.449 1.35e-10 ***
## Cranio 10.62032 0.41812 25.400 < 2e-16 ***
## SessoM 70.88786 11.01907 6.433 1.49e-10 ***
## I(Lunghezza^2) 0.03170 0.00326 9.724 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.2 on 2490 degrees of freedom
## Multiple R-squared: 0.7376, Adjusted R-squared: 0.7369
## F-statistic: 1166 on 6 and 2490 DF, p-value: < 2.2e-16summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Lunghezza^2), data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1169.62 -181.77 -12.79 163.77 1786.03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 212.288548 723.852095 0.293 0.769336
## N.gravidanze 14.085464 4.266175 3.302 0.000975 ***
## Gestazione 42.551398 3.876629 10.976 < 2e-16 ***
## Lunghezza -20.267001 3.162718 -6.408 1.76e-10 ***
## Cranio 10.651783 0.418894 25.428 < 2e-16 ***
## SessoM 69.968733 11.038797 6.338 2.75e-10 ***
## I(Lunghezza^2) 0.031655 0.003267 9.690 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.7 on 2491 degrees of freedom
## Multiple R-squared: 0.7369, Adjusted R-squared: 0.7363
## F-statistic: 1163 on 6 and 2491 DF, p-value: < 2.2e-16Il summary del nuovo modello rispetto al mod3 mostra un R2 adj praticamente identico. I coefficienti tra i due modelli variano ma non di tantissimo.
Passiamo all’analisi dei residui di questo nuovo modello
mod3_no_1551
par(mfrow=c(2,2))
plot(mod3_no_1551, pch=20, col="royalblue3")
La situazione è rimasta analoga alla precedente ma si è aggiunta l’osservazione 1551 del nuovo dataset
Restano in evidenza principalmente le osservazioni 1305 e 155.
Calcoliamo ora i valori di leva per il mod3_no_1551, identifichiamo il valore soglia che classificherà i valori come valori di leva o meno e andiamo infine a visualizzare la soglia anche nel grafico con una linea orizzontale azzurra
lev_mod3_no_1551<-hatvalues(mod3_no_1551)
plot(lev_mod3_no_1551, pch=20, col="royalblue3")
p_mod3_no_1551 = sum(lev_mod3_no_1551)
n_mod3_no_1551 <- nrow(neonati_no_1551)
soglia_mod3_no_1551 = 2*p_mod3_no_1551/n_mod3_no_1551
abline(h=soglia_mod3_no_1551,col="turquoise3",lwd=1.5)
Ci sono molti valori leverage ovvero oltre il valore soglia definito.
Possiamo identificare questi valori, ma essendo così tanti guardiamone solo il numero.
length(lev_mod3_no_1551[lev_mod3_no_1551>soglia_mod3_no_1551])## [1] 141Per capire meglio l’impatto di questi leverage andiamo avanti con l’analisi dei residui
plot(rstudent(mod3_no_1551), pch=20, col="royalblue3")
abline(h=c(-2,2),col="turquoise3",lwd=1.5)
Ci sono molti punti che superano i valori soglia di meno 2 e 2
Andiamo ad effettuare anche il test statistico per gli outlier
car::outlierTest(mod3_no_1551)## rstudent unadjusted p-value Bonferroni p
## 1551 7.269113 4.8212e-13 1.2039e-09
## 1306 4.840484 1.3745e-06 3.4321e-03
## 155 4.747261 2.1788e-06 5.4404e-03
## 1399 -4.355687 1.3802e-05 3.4464e-02
## 1694 4.331325 1.5409e-05 3.8476e-02Le tre osservazioni che vedevamo dai grafici mostrano avere qualche problema anche da questo test. Ci sono anche altre 2 osservazioni.
Andiamo anche a visualizzare la distanza di Cook
cook_mod3_no_1551<-cooks.distance(mod3_no_1551)
plot(cook_mod3_no_1551,ylim = c(0,1.4), pch=20, col="royalblue3")
Come avevamo visto dai grafici abbiamo valori al di sopra della soglia di allarme della distanza di Cook. Un’osservazione ha una distanza di Cook pari a circa 1.4 (l’osservazione 1551 del nuovo dataset)
Andiamo a verificare esattamente a quanto ammonta il massimo della distanza di Cook
max(cook_mod3_no_1551)## [1] 1.360464Il valore massimo raggiunto delle distanze di Cook è di circa 1.36 quindi ci sono valori particolarmente influenti sulle stime di regressione.
par(mfrow=c(1,2))
plot(residuals(mod3_no_1551), pch=20, col="royalblue3")
abline(h=mean(residuals(mod3_no_1551)), col="turquoise3",lwd=1.5)
plot(density(residuals(mod3_no_1551)), col="royalblue3")
Possiamo vedere che la distribuzione dei residui assomiglia abbastanza alla normale ma ha una coda a sinistra più schiacciata e più corta rispetto a quella di destra.
Andiamo a fare un test di Shapiro Wilk, un test di normalità in cui l’ipotesi nulla è proprio la distribuzione normale della variabile di riferimento
shapiro.test(mod3_no_1551$residuals)##
## Shapiro-Wilk normality test
##
## data: mod3_no_1551$residuals
## W = 0.98616, p-value = 7.296e-15p value molto piccolo quindi si rifiuta l’ipotesi di normalità
Andiamo a saggiare queste altre due ipotesi
RESIDUI A VARIANZA COSTANTE (omoschedasticità) –> TEST DI BREUSCH PAGAN
RESIDUI NON CORRELATI TRA DI LORO –> TEST DI DURBIN WATSON
bptest(mod3_no_1551)##
## studentized Breusch-Pagan test
##
## data: mod3_no_1551
## BP = 128.99, df = 6, p-value < 2.2e-16il p value è minore della soglia di 0.05 quindi si rifiuta l’ipotesi di omoschedasticità (residui a varianza non costante)
dwtest(mod3_no_1551)##
## Durbin-Watson test
##
## data: mod3_no_1551
## DW = 1.946, p-value = 0.08843
## alternative hypothesis: true autocorrelation is greater than 0il p value è di 0.08843 (valore più grande della soglia dello 0.05) quindi non si rifiuta l’ipotesi di indipendenza
Il mod3_no_1551 ha un outlier molto influente (osservazione 1551 del subset neonati_no_1551) che andrebbe indagato maggiormente prima di utilizzare il modello.
mod_var_richieste
par(mfrow=c(2,2))
plot(mod_var_richieste, pch=20, col="royalblue3")
Residuals vs Fitted
La media di zero sembra essere rispettata per valori più alti mentre i valori più bassi tendono ad essere maggiori di 0. Le osservazioni 1549, 155 e 1305 sono fuori dalla nuvola
Q-Q Residuals
In questo caso sembra che i punti siano correttamente disposti sopra la bisettrice, a parte qualche problemino nella coda inferiore e soprattutto quella superiore. Ci sono di nuovo le osservazioni 1549, 155 e 1305 che sono più lontane dalle altre, principalmente la 1549.
Scale-Location
Anche qui si nota una curvatura per i valori più piccoli e le solite tre osservazioni più distanti dalla nuvola.
Residuals vs Leverage
In questo caso l’osservazione 1549 risulta essere oltre la linea di attenzione.
Calcoliamo ora i valori di leva per il mod_var_richieste, identifichiamo il valore soglia che classificherà i valori come valori di leva o meno e andiamo infine a visualizzare la soglia anche nel grafico con una linea orizzontale azzurra
lev_mod_var_richieste<-hatvalues(mod_var_richieste)
plot(lev_mod_var_richieste, pch=20, col="royalblue3")
p_mod_var_richieste = sum(lev_mod_var_richieste)
soglia_mod_var_richieste = 2*p_mod_var_richieste/n
abline(h=soglia_mod_var_richieste,col="turquoise3",lwd=1.5)
Ci sono moltissimi valori leverage ovvero oltre il valore soglia definito.
Possiamo identificare questi valori.
length(lev_mod_var_richieste[lev_mod_var_richieste>soglia_mod_var_richieste])## [1] 211Per capire meglio l’impatto di questi leverage andiamo avanti con l’analisi dei residui
plot(rstudent(mod_var_richieste), pch=20, col="royalblue3")
abline(h=c(-2,2),col="turquoise3",lwd=1.5)
Ci sono molti punti che superano i valori soglia di meno 2 e 2
Andiamo ad effettuare anche il test statistico per gli outlier
car::outlierTest(mod_var_richieste)## rstudent unadjusted p-value Bonferroni p
## 1549 10.033980 2.9696e-23 7.4182e-20
## 155 5.019634 5.5428e-07 1.3846e-03
## 1305 4.820812 1.5158e-06 3.7865e-03Le tre osservazioni che vedevamo dai grafici mostrano avere qualche problema anche da questo test.
Andiamo anche a visualizzare la distanza di Cook
cook_mod_var_richieste<-cooks.distance(mod_var_richieste)
plot(cook_mod_var_richieste,ylim = c(0,1), pch=20, col="royalblue3")
Come avevamo visto il valore massimo della distanza di Cook è maggiore di 0.5, circa 0.7. Le altre osservazioni sono tutte sotto la soglia di attenzione di 0.5
Andiamo a verificare esattamente a quanto ammonta il massimo della distanza di Cook
max(cook_mod_var_richieste)## [1] 0.7114668il valore massimo raggiunto delle distanze di cook è di circa 0.71.
par(mfrow=c(1,2))
plot(residuals(mod_var_richieste), pch=20, col="royalblue3")
abline(h=mean(residuals(mod_var_richieste)), col="turquoise3",lwd=1.5)
plot(density(residuals(mod_var_richieste)), col="royalblue3")
Possiamo vedere che la distribuzione dei residui assomiglia abbastanza alla normale ma ha una coda a coda a destra meno schiacciata e più lunga
Andiamo a fare un test di Shapiro Wilk, un test di normalità in cui l’ipotesi nulla è proprio la distribuzione normale della variabile di riferimento
shapiro.test(mod_var_richieste$residuals)##
## Shapiro-Wilk normality test
##
## data: mod_var_richieste$residuals
## W = 0.97416, p-value < 2.2e-16p value molto piccolo quindi si rifiuta l’ipotesi di normalità
Andiamo a saggiare queste altre due ipotesi
RESIDUI A VARIANZA COSTANTE (omoschedasticità) –> TEST DI BREUSCH PAGAN
RESIDUI NON CORRELATI TRA DI LORO –> TEST DI DURBIN WATSON
bptest(mod_var_richieste)##
## studentized Breusch-Pagan test
##
## data: mod_var_richieste
## BP = 89.841, df = 6, p-value < 2.2e-16il p value è molto piccolo quindi si rifiuta l’ipotesi di omoschedasticità (residui a varianza non costante)
dwtest(mod_var_richieste)##
## Durbin-Watson test
##
## data: mod_var_richieste
## DW = 1.9539, p-value = 0.1243
## alternative hypothesis: true autocorrelation is greater than 0il p value è di 0.1243 (valore più grande della soglia dello 0.05) quindi non si rifiuta l’ipotesi di indipendenza
A seguito di questi risultati teniamo in considerazione solo il mod_var_richieste e il mod4_no_1549. Come analisi ci fermiamo qui e consideriamo accettabile il mod4_no_1549 ma andrebbero approfondite alcune situazioni problematiche prima di scegliere realmente il modello o si dovrebbero provare dei GLM. Portiamo avanti anche il mod_var_richieste solo perché come detto sopra nella previsione ci è stato chiesto di prevedere il peso di un neonato di cui ci venivano fornite alcune specifiche informazioni.
Previsioni e risultati
Andiamo ora a stimare il peso di una neonata considerando una madre alla terza gravidanza (quindi N.gravidanze=2), non fumatrice, che partorirà alla 39esima settimana.
dati.test <- data.frame(N.gravidanze = 2,
Gestazione = 39,
Lunghezza = mean(Lunghezza), # lunghezza media del neonato
Cranio = mean(Cranio), # circonferenza cranica media
Fumatrici = factor("0", levels = levels(Fumatrici)),
Sesso = factor("F", levels = levels(Sesso)))predict(mod_var_richieste, newdata = dati.test)## 1
## 3260.165Con il mod_var_richieste il peso previsto è 3260.165g
dati.test_no_1549 <- data.frame(N.gravidanze = 2,
Gestazione = 39,
Lunghezza = mean(neonati_no_1549$Lunghezza), # lunghezza media del neonato
Cranio = mean(neonati_no_1549$Cranio)) # circonferenza cranica mediapredict(mod4_no_1549, dati.test_no_1549) # circonferenza cranica media## 1
## 3284.329Con il mod4_no_1549 è 3284.33g
Visualizzazioni di modelli
Andiamo a visualizzare uno scatterplot 3D di peso in relazione a lunghezza e cranio differenziando per sesso. Il grafico può essere esplorato con il mouse.
rgl::setupKnitr(autoprint = TRUE)
car::scatter3d(Peso~Lunghezza+Cranio, point.col = Sesso, groups = Sesso, surface = TRUE, axis.col = c("black", "black", "black"))## Caricamento dei namespace richiesti: mgcvrgl::rglwidget()
Andiamo a visualizzare uno scatterplot 3D di peso in relazione a lunghezza e gestazione differenziando per sesso. Il grafico può essere esplorato con il mouse.
rgl::setupKnitr(autoprint = TRUE)
car::scatter3d(Peso~Lunghezza+Gestazione, point.col = Sesso, groups=Sesso, surface = TRUE, axis.col = c("black", "black", "black"))
rgl::rglwidget()
Adesso visualizziamo la relazione tra lunghezza e peso del neonato differenziata per sesso
ggplot(data=neonati)+
geom_point(aes(x=Lunghezza,
y=Peso,
col=Sesso),position = "jitter", alpha=0.5)+
geom_smooth(aes(x=Lunghezza,
y=Peso,
col=Sesso),
se=F,method = "lm",alpha=0.5)+
geom_smooth(aes(x=Lunghezza,
y=Peso,
linetype = "Totale"),
method="lm", color="royalblue3", fill="royalblue1", se=TRUE)+
labs(x="Lunghezza (mm)",
y="Peso (g)",
title = "Relazione tra Lunghezza e Peso \ncon Regressioni Lineari per Sesso e Totale",
linetype = "Tipo di Regressione")+
theme_classic()+
theme(legend.position = "right",
legend.box = "vertical",
plot.title = element_text(size=14, hjust = 0.5, face = "bold"),
legend.title = element_text(face="bold"),
text = element_text(size=10))+
scale_x_continuous(breaks = seq(300,600,20))+
scale_y_continuous(breaks = seq(500,5000,250))+
scale_linetype_manual(values = c("Totale" = "dashed"))+
scale_color_manual(
labels = c("F" = "Femmina", "M" = "Maschio"),
values = c("F" = "#F8766D", "M" = "#00BFC4"))+
guides(
linetype = guide_legend(order = 1),
col = guide_legend(title = NULL, order = 2)
)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Adesso visualizziamo la relazione tra gestazione e peso del neonato differenziata per sesso
ggplot(data=neonati)+
geom_point(aes(x=Gestazione,
y=Peso,
col=Sesso),position = "jitter", alpha=0.5)+
geom_smooth(aes(x=Gestazione,
y=Peso,
col=Sesso),
se=F,method = "lm",alpha=0.5)+
geom_smooth(aes(x=Gestazione,
y=Peso,
linetype = "Totale"),
method="lm", color="royalblue3", fill="royalblue1", se=TRUE)+
labs(x="Gestazione (settimane)",
y="Peso (g)",
title = "Relazione tra Gestazione e Peso \ncon Regressioni Lineari per Sesso e Totale",
linetype = "Tipo di Regressione")+
theme_classic()+
theme(legend.position = "right",
legend.box = "vertical",
plot.title = element_text(size=14, hjust = 0.5, face = "bold"),
legend.title = element_text(face="bold"),
text = element_text(size=10))+
scale_x_continuous(breaks = seq(25,43,1))+
scale_y_continuous(breaks = seq(500,5000,250))+
scale_linetype_manual(values = c("Totale" = "dashed"))+
scale_color_manual(
labels = c("F" = "Femmina", "M" = "Maschio"),
values = c("F" = "#F8766D", "M" = "#00BFC4"))+
guides(
linetype = guide_legend(order = 1),
col = guide_legend(title = NULL, order = 2)
)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Adesso visualizziamo la relazione tra cranio e peso del neonato differenziata per sesso
ggplot(data=neonati)+
geom_point(aes(x=Cranio,
y=Peso,
col=Sesso),position = "jitter", alpha=0.5)+
geom_smooth(aes(x=Cranio,
y=Peso,
col=Sesso),
se=F,method = "lm",alpha=0.5)+
geom_smooth(aes(x=Cranio,
y=Peso,
linetype = "Totale"),
method="lm", color="royalblue3", fill="royalblue1", se=TRUE)+
labs(x="Cranio (mm)",
y="Peso (g)",
title = "Relazione tra Cranio e Peso \ncon Regressioni Lineari per Sesso e Totale",
linetype = "Tipo di Regressione")+
theme_classic()+
theme(legend.position = "right",
legend.box = "vertical",
plot.title = element_text(size=14, hjust = 0.5, face = "bold"),
legend.title = element_text(face="bold"),
text = element_text(size=10))+
scale_x_continuous(breaks = seq(230,390,10))+
scale_y_continuous(breaks = seq(500,5000,250))+
scale_linetype_manual(values = c("Totale" = "dashed"))+
scale_color_manual(
labels = c("F" = "Femmina", "M" = "Maschio"),
values = c("F" = "#F8766D", "M" = "#00BFC4"))+
guides(
linetype = guide_legend(order = 1),
col = guide_legend(title = NULL, order = 2)
)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Infine visualizziamo la relazione tra gestazione e peso del neonato differenziata per sesso
ggplot(data=neonati)+
geom_smooth(aes(x=Gestazione,
y=Peso,
linetype = "Totale"),
method="lm", color="red", size=1.5,fill="plum1", se=TRUE)+
geom_point(aes(x=Gestazione,
y=Peso,
col=Fumatrici),position = "jitter",alpha=0.8)+
geom_smooth(aes(x=Gestazione,
y=Peso,
col=Fumatrici),
se=F,method = "lm",alpha=0.5)+
labs(x="Gestazione (settimane)",
y="Peso (g)",
title = "Relazione tra Gestazione e Peso \ncon Regressioni Lineari per \nMadri fumatrici vs non fumatrici e Totale",
linetype = "Tipo di Regressione")+
theme_classic()+
theme(legend.position = "right",
legend.box = "vertical",
plot.title = element_text(size=14, hjust = 0.5, face = "bold"),
legend.title = element_text(face="bold"),
text = element_text(size=10))+
scale_x_continuous(breaks = seq(25,43,1))+
scale_y_continuous(breaks = seq(500,5000,250))+
scale_linetype_manual(values = c("Totale" = "solid"))+
scale_color_manual(
labels = c("0" = "Non Fumatrici", "1" = "Fumatrici"),
values = c("0" = "royalblue3", "1" = "darkorange"))+
guides(
linetype = guide_legend(order = 1),
col = guide_legend(title = NULL, order = 2)
)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
In questo grafico la linea della regressione delle madri non fumatrici è quasi completamente sovrapposta alla linea della regressione totale.