Modello_predittivo_peso
Vittoria
2024-09-09
data<-read.csv("neonati.csv")Il dataset oggetto di studio contiene dati riferiti alle nascite in 3 ospedali; in particolare sono descritte le principali caratteristiche fisiche dei neonati ed alcune variabili, riferite alle madri e al parto, che possano ragionevolmente influenzare le caratteristiche del neonato. L’obiettivo dello studio infatti e’ di prevedere il peso dei futuri neonati sulla base di questi dati. Variabili: età della madre: quantitativa discreta numero di gravidanze sostenute: quantitativa discreta Madre fumatrice: qualitativa su scala nominale (espressa come fattore) N° di settimane di gestazione: quantitativa discreta peso in grammi del neonato: quantitativa continua Lunghezza in mm del neonato: quantitativa continua Diametro in mm del cranio del neonato: quantitativa continua Tipo di parto: qualitativa su scala nominale Ospedale: qualitativa su scala nominale Sesso del neonato: qualitativa su scala nominale
Caratteristiche dei neonati:
summary(data$Peso)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 830 2990 3300 3284 3620 4930summary(data$Lunghezza)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 310.0 480.0 500.0 494.7 510.0 565.0summary(data$Cranio)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 235 330 340 340 350 390table(data$Sesso)##
## F M
## 1256 1244library(ggplot2)
ggplot(data, aes(x=Peso))+
geom_density()+
geom_vline(aes(xintercept=mean(Peso)), color="red")+
labs(title="Distribuzione dei valori di peso (g)", x="", y="densita'" )+
theme_minimal()
ggplot(data, aes(x=Lunghezza))+
geom_density()+
geom_vline(aes(xintercept=mean(Lunghezza)), color="red")+
labs(title="Distribuzione dei valori di lunghezza (mm)", x="", y="densita'" )+
theme_minimal()
ggplot(data, aes(x=Cranio))+
geom_density()+
geom_vline(aes(xintercept=mean(Cranio)), color="red")+
labs(title="Distribuzione dei valori del diametro del cranio (mm)", x="", y="densita'" )+
theme_minimal()
Caratteristiche delle madri:
summary(data$Anni.madre)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 25.00 28.00 28.16 32.00 46.00summary(data$N.gravidanze)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 1.0000 0.9812 1.0000 12.0000table(data$Fumatrici)##
## 0 1
## 2396 104Salta all’occhio la presenza di valori spuri per gli anni della madre che potrebbero inficiare l’analisi (min=0)
table(data$Anni.madre)##
## 0 1 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
## 1 1 1 2 6 13 18 24 45 66 74 100 115 131 180 184 197 172 174 200
## 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
## 147 159 110 96 66 64 41 38 27 19 13 8 2 4 1 1I valori da correggere sono 0 e 1. Sebbene si possa provare a considerare i corrispondenti valori di N.gravidanze per “dedurre” i valori corretti, penso che sostituirli con il valore mediano (che peraltro, essendoci tutti valori interi, praticamente coincide con quello medio) sia una soluzione piu’ pratica e piuttosto robusta:
data$Anni.madre[data$Anni.madre==0] <- 28
data$Anni.madre[data$Anni.madre==1] <- 28
table(data$Anni.madre)##
## 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
## 1 2 6 13 18 24 45 66 74 100 115 131 180 184 197 174 174 200 147 159
## 33 34 35 36 37 38 39 40 41 42 43 44 45 46
## 110 96 66 64 41 38 27 19 13 8 2 4 1 1Altre caratteristiche:
table(data$Tipo.parto)##
## Ces Nat
## 728 1772table(data$Ospedale)##
## osp1 osp2 osp3
## 816 849 835summary(data$Gestazione)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 25.00 38.00 39.00 38.98 40.00 43.00Tra queste variabili, creo una rappresentazione grafica solo per quelle numeriche, per visualizzarle meglio:
plot(table(data$Anni.madre),
main="Distribuzione di frequenza dell'eta' delle madri", ylab="")
plot(table(data$N.gravidanze),
main="Distribuzione di frequenza del numero di gravidanze", ylab="")
plot(table(data$Gestazione),
main="Distribuzione di frequenza delle settimane di gestazione", ylab="")
- Per saggiare l’ipotesi che la media del peso e della lunghezza dei neonati del campione siano uguali a quelle della popolazione, utilizzo un t-test (ho preso i valori medi dal sito di un ospedale di Roma) ed imposto il livello di significativita’ a 0.05. Ma, prima di farlo, eseguo un test di normalita’ dei dati:
shapiro.test(data$Peso)##
## Shapiro-Wilk normality test
##
## data: data$Peso
## W = 0.97066, p-value < 2.2e-16p-value<0.05: i dati non seguono una distribuzione normale, per cui eseguo un test non parametrico per confrontare la mediana del campione con la mediana della poopolazione(a questo punto, ipoteticamente uguale alla media): H0: mean(peso)=3300g H1: mean(peso)!=3300g
wilcox.test(data$Peso, mu=3300)##
## Wilcoxon signed rank test with continuity correction
##
## data: data$Peso
## V = 1495594, p-value = 0.9612
## alternative hypothesis: true location is not equal to 3300Il p-value>0.05 quindi non rifiuto l’ipotesi nulla.
shapiro.test(data$Lunghezza)##
## Shapiro-Wilk normality test
##
## data: data$Lunghezza
## W = 0.90941, p-value < 2.2e-16Neanche le lunghezze sono distribuite normalmente, come peraltro appariva gia’ dal grafico di densita’ di distribuzione. H0: mean(lunghezza)=500mm H1: mean(lunghezza)!=500mm
wilcox.test(data$Lunghezza, mu=500)##
## Wilcoxon signed rank test with continuity correction
##
## data: data$Lunghezza
## V = 877236, p-value < 2.2e-16
## alternative hypothesis: true location is not equal to 500p-value<0.05 quindi rifiuto H0 e concludo che la media delle lunghezze nel campione e’ significativamente diversa da quella della popolazione
- Per verificare differenze significative tra i due sessi, considero i due sessi come gruppi indipendenti e verifico per ognuno la normalita’ di ogni variabile e l’omogeneita’ delle varianze tra i gruppi; solo infine svolgo la verifica dell’uguaglianza fra medie (per le variabili peso, lunghezza e diametro del cranio):
shapiro.test(data$Peso[data$Sesso=="M"])##
## Shapiro-Wilk normality test
##
## data: data$Peso[data$Sesso == "M"]
## W = 0.96647, p-value = 2.321e-16shapiro.test(data$Peso[data$Sesso=="F"])##
## Shapiro-Wilk normality test
##
## data: data$Peso[data$Sesso == "F"]
## W = 0.96285, p-value < 2.2e-16library(car)## Caricamento del pacchetto richiesto: carDataleveneTest(data$Peso~data$Sesso)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.7931 0.3732
## 2498I dati del peso tra i gruppi non rispettano l’assunzione di normalita’ per il t-test (per l’omoschedasticita’ ho usato il Levene perche’ ho letto che e’ piu’ indicato in casi di non-normalita’ delle distribuzioni, ma anche con il test di Bartlett l’esito sarebbe stato lo stesso), quindi procedo con il non-parametrico. Il sistema di ipotesi e’ analogo per tutti, ossia l’H0 e’ che le medie dei due gruppi sono significativamente uguali, mentre l’H1 e’ che sono significativamente diverse:
wilcox.test(data$Peso~data$Sesso)##
## Wilcoxon rank sum test with continuity correction
##
## data: data$Peso by data$Sesso
## W = 538641, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0Siccome p-value<0.05, rifiuto l’ipotesi nulla di uguaglianza e concludo che ci sono differenze significative di peso tra i due sessi.
shapiro.test(data$Lunghezza[data$Sesso=="M"])##
## Shapiro-Wilk normality test
##
## data: data$Lunghezza[data$Sesso == "M"]
## W = 0.92028, p-value < 2.2e-16shapiro.test(data$Lunghezza[data$Sesso=="F"])##
## Shapiro-Wilk normality test
##
## data: data$Lunghezza[data$Sesso == "F"]
## W = 0.89953, p-value < 2.2e-16Data la non-normalita’ procedo direttamente col test non parametrico:
wilcox.test(data$Lunghezza~data$Sesso)##
## Wilcoxon rank sum test with continuity correction
##
## data: data$Lunghezza by data$Sesso
## W = 594455, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0Tra i due sessi ci sono differenze significative anche di lunghezza.
shapiro.test(data$Cranio[data$Sesso=="M"])##
## Shapiro-Wilk normality test
##
## data: data$Cranio[data$Sesso == "M"]
## W = 0.97046, p-value = 3.006e-15shapiro.test(data$Cranio[data$Sesso=="F"])##
## Shapiro-Wilk normality test
##
## data: data$Cranio[data$Sesso == "F"]
## W = 0.95543, p-value < 2.2e-16wilcox.test(data$Cranio~data$Sesso)##
## Wilcoxon rank sum test with continuity correction
##
## data: data$Cranio by data$Sesso
## W = 641638, p-value = 9.633e-15
## alternative hypothesis: true location shift is not equal to 0E ci sono differenze anche per il diametro del cranio.
- Per saggiare l’ipotesi che il numero di parti cesarei sia sproporzionato tra i tre ospedali, eseguo un test chi-quadrato:
tabella_contingenza <- table(data$Tipo.parto,data$Ospedale)
chisq.test(tabella_contingenza)##
## Pearson's Chi-squared test
##
## data: tabella_contingenza
## X-squared = 1.0972, df = 2, p-value = 0.5778Il p-value>0.05 non mi permette di rifiutare l’ipotesi nulla, secondo cui non ci sono differenze significative nel numero di parti cesarei tra i tre ospedali.
Analisi multidimensionale: 1. E’ gia’ emersa in precedenza una significativa differenza di peso tra i due sessi, che posso visualizzare tramite boxplot:
boxplot(Peso~Sesso, data=data,
ylab="peso del neonato (g)")
Quindi procedo ad effettuare qualche controllo sulle relazioni tra le
altre caratteristiche del neonato:
plot(data$Lunghezza, data$Peso,
main="Scatterplot peso-lunghezza",
xlab="Lunghezza (mm)",
ylab="Peso (g)")
plot(data$Cranio, data$Peso,
main="Scatterplot peso-diametro del cranio",
xlab="Diametro del cranio (mm)",
ylab="Peso (g)")
plot(data$Lunghezza, data$Cranio,
main="Scatterplot lunghezza-diametro del cranio",
xlab="Lunghezza (mm)",
ylab="Diametro del cranio (mm)")
cor(data$Lunghezza, data$Peso)## [1] 0.7960368cor(data$Cranio, data$Peso)## [1] 0.7048015cor(data$Lunghezza, data$Cranio)## [1] 0.603341Come era immaginabile, il peso (variabile risposta) e’ positivamente correlato sia con la lunghezza che con il diametro del cranio e c’e’ una correlazione lineare positiva anche tra lunghezza e cranio; tuttavia, non penso che siano queste le variabili piu’ utili per prevedere il peso dei futuri neonati in quanto saranno valori noti solo dopo il parto, percio’ mi concentrerei sulle relazioni con altre variabili note a priori:
plot(data$Gestazione, data$Peso, xlab="Settimane di gestazione", ylab="Peso")
plot(data$Anni.madre, data$Peso, xlab="Anni della madre", ylab="Peso")
plot(data$N.gravidanze, data$Peso, xlab="Numero di gravidanze sostenute", ylab="Peso")
Dagli scatterplot sembra esserci una moderata correlazione positiva tra
peso e settimane di gestazione, mentre non sembra esserci una
particolare correlazione tra peso e anni della madre e tra peso e numero
di gravidanze. Ma cio’ puo’ essere verificato ulteriormente:
cor(data$Peso,data$Gestazione)## [1] 0.5917687cor(data$Peso,data$Anni.madre)## [1] -0.02377346cor(data$Peso,data$N.gravidanze)## [1] 0.0024073Si conferma una moderata correlazione lineare positiva tra peso e settimane di gestazione e assenza di correlazione lineare con gli anni della madre e col numero di gravidanze. Per visualizzare gli effetti delle altre variabili sul peso, per prima cosa li visualizzo tramite boxplot (lo faccio anche per l’ospedale nonostante non sospetti che possa influire, ma per motivi clinici vari non e’ comunque da escludere):
boxplot(Peso~Tipo.parto, data=data,
xlab="tipo di parto",
ylab="peso del neonato (g)")
boxplot(Peso~Fumatrici, data=data,
xlab="madre fumatrice (0=no / 1=si)",
ylab="peso del neonato (g)")
boxplot(Peso~Ospedale, data=data,
xlab="ospedale",
ylab="peso del neonato (g)")
I boxplot non mostrano particolari differenze nei confronti tra tipo di
parto e tra diversi ospedali; si osserva una diversa variabilita’ nei
gruppi per quanto riguarda il fattore fumatrice. Dopodiche’ eseguo una
serie di t-test di confronto (previa verifica delle assunzioni
necessarie, altrimenti mi avvalgo di test non parametrici) quando ho 2
gruppi:
shapiro.test(data$Peso[data$Fumatrici==0])##
## Shapiro-Wilk normality test
##
## data: data$Peso[data$Fumatrici == 0]
## W = 0.96912, p-value < 2.2e-16shapiro.test(data$Peso[data$Fumatrici==1])##
## Shapiro-Wilk normality test
##
## data: data$Peso[data$Fumatrici == 1]
## W = 0.9701, p-value = 0.01867shapiro.test(data$Peso[data$Tipo.parto=="Nat"])##
## Shapiro-Wilk normality test
##
## data: data$Peso[data$Tipo.parto == "Nat"]
## W = 0.96791, p-value < 2.2e-16shapiro.test(data$Peso[data$Tipo.parto=="Ces"])##
## Shapiro-Wilk normality test
##
## data: data$Peso[data$Tipo.parto == "Ces"]
## W = 0.9783, p-value = 6.428e-09wilcox.test(data$Peso~data$Fumatrici)##
## Wilcoxon rank sum test with continuity correction
##
## data: data$Peso by data$Fumatrici
## W = 138162, p-value = 0.05971
## alternative hypothesis: true location shift is not equal to 0wilcox.test(data$Peso~data$Tipo.parto)##
## Wilcoxon rank sum test with continuity correction
##
## data: data$Peso by data$Tipo.parto
## W = 634748, p-value = 0.5315
## alternative hypothesis: true location shift is not equal to 0e anova quando ho piu’ di 2 gruppi (previa verifica delle assunzioni):
shapiro.test(data$Peso[data$Ospedale=="osp1"])##
## Shapiro-Wilk normality test
##
## data: data$Peso[data$Ospedale == "osp1"]
## W = 0.9623, p-value = 1.177e-13shapiro.test(data$Peso[data$Ospedale=="osp2"])##
## Shapiro-Wilk normality test
##
## data: data$Peso[data$Ospedale == "osp2"]
## W = 0.98001, p-value = 2.22e-09shapiro.test(data$Peso[data$Ospedale=="osp3"])##
## Shapiro-Wilk normality test
##
## data: data$Peso[data$Ospedale == "osp3"]
## W = 0.96439, p-value = 2.16e-13siccome la distribuzione non e’ normale, procedo con il corrispettivo non parametrico:
kruskal.test(Peso~Ospedale, data=data)##
## Kruskal-Wallis rank sum test
##
## data: Peso by Ospedale
## Kruskal-Wallis chi-squared = 3.0052, df = 2, p-value = 0.2225Questi test suggeriscono di non rifiutare l’ipotesi nulla secondo cui non ci sono differenze significative tra i gruppi delle 3 variabili considerate, il che mi farebbe pensare che il peso dei neonati non dipende da queste variabili (fatta forse eccezione per la variabile Fumatrici, che e’ un po’ al limite).
- Il modello di regressione lineare multipla chiarira’ la situazione:
attach(data)
mod1 <- lm(Peso ~ . ,data=data)
summary(mod1)##
## Call:
## lm(formula = Peso ~ ., data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1123.3 -181.2 -14.6 160.7 2612.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6735.1677 141.3977 -47.633 < 2e-16 ***
## Anni.madre 0.7983 1.1463 0.696 0.4862
## N.gravidanze 11.4118 4.6665 2.445 0.0145 *
## Fumatrici -30.1567 27.5396 -1.095 0.2736
## Gestazione 32.5265 3.8179 8.520 < 2e-16 ***
## Lunghezza 10.2951 0.3007 34.237 < 2e-16 ***
## Cranio 10.4725 0.4261 24.580 < 2e-16 ***
## Tipo.partoNat 29.5027 12.0848 2.441 0.0147 *
## Ospedaleosp2 -11.2216 13.4388 -0.835 0.4038
## Ospedaleosp3 28.0984 13.4972 2.082 0.0375 *
## SessoM 77.5473 11.1779 6.938 5.07e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 669.1 on 10 and 2489 DF, p-value: < 2.2e-16Questo primo modello mette in evidenza la significativita’ delle seguenti variabili nello spiegare l’andamento della variabile risposta: Gestazione, Lunghezza, Cranio e Sesso; le variabili Tipo.parto, Ospedale e N.gravidanze sono meno significative; le restanti variabili non sono significative. La bonta’ del modello e’ del 72.78%
- Il modello si puo’ migliorare, ed innanzitutto rimuovo le variabili non significative (nonostante l’eta’ della madre sia una variabile di controllo che generalmente e’ meglio tenere, in questo caso specifico non penso che sia poi cosi’ rilevante ed avendo molte variabili che invece sono significative preferisco rimuoverla):
mod2 <- update(mod1, ~.-Anni.madre -Fumatrici)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1113.18 -181.16 -16.58 161.01 2620.19
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.4293 135.9438 -49.340 < 2e-16 ***
## N.gravidanze 12.3619 4.3325 2.853 0.00436 **
## Gestazione 31.9909 3.7896 8.442 < 2e-16 ***
## Lunghezza 10.3086 0.3004 34.316 < 2e-16 ***
## Cranio 10.4922 0.4254 24.661 < 2e-16 ***
## Tipo.partoNat 29.2803 12.0817 2.424 0.01544 *
## Ospedaleosp2 -11.0227 13.4363 -0.820 0.41209
## Ospedaleosp3 28.6408 13.4886 2.123 0.03382 *
## SessoM 77.4412 11.1756 6.930 5.36e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared: 0.7287, Adjusted R-squared: 0.7278
## F-statistic: 836.3 on 8 and 2491 DF, p-value: < 2.2e-16Il mod2 mostra un aumento di significativita’ della variabile N.gravidanze, mentre il resto e’ invariato. Siccome per la variabile Ospedale ho che, rispetto a osp1, uno dei due livelli e’ significativo(osp3) e l’altro no(osp2), provo a cambiare il livello di riferimento:
data$Ospedale <- relevel(as.factor(data$Ospedale), ref="osp2")
mod3 <- update(mod2)
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1113.18 -181.16 -16.58 161.01 2620.19
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6718.4520 135.9956 -49.402 < 2e-16 ***
## N.gravidanze 12.3619 4.3325 2.853 0.00436 **
## Gestazione 31.9909 3.7896 8.442 < 2e-16 ***
## Lunghezza 10.3086 0.3004 34.316 < 2e-16 ***
## Cranio 10.4922 0.4254 24.661 < 2e-16 ***
## Tipo.partoNat 29.2803 12.0817 2.424 0.01544 *
## Ospedaleosp1 11.0227 13.4363 0.820 0.41209
## Ospedaleosp3 39.6635 13.3677 2.967 0.00303 **
## SessoM 77.4412 11.1756 6.930 5.36e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared: 0.7287, Adjusted R-squared: 0.7278
## F-statistic: 836.3 on 8 and 2491 DF, p-value: < 2.2e-16Il livello osp3 e’ significativo anche rispetto ad osp2, quindi provo a racchiudere insieme osp1 e osp2:
data$Ospedale_modificato <- ifelse(data$Ospedale %in% c("osp1","osp2"), "osp1_osp2", "osp3")
data$Ospedale_modificato <- factor(data$Ospedale_modificato)
mod4 <- lm(Peso ~.-Ospedale -Anni.madre -Fumatrici, data=data)
summary(mod4)##
## Call:
## lm(formula = Peso ~ . - Ospedale - Anni.madre - Fumatrici, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1118.5 -182.1 -16.6 162.8 2619.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6712.5106 135.7937 -49.432 < 2e-16 ***
## N.gravidanze 12.4340 4.3313 2.871 0.00413 **
## Gestazione 32.0128 3.7893 8.448 < 2e-16 ***
## Lunghezza 10.3028 0.3003 34.309 < 2e-16 ***
## Cranio 10.4962 0.4254 24.674 < 2e-16 ***
## Tipo.partoNat 29.2959 12.0809 2.425 0.01538 *
## SessoM 77.5462 11.1741 6.940 4.99e-12 ***
## Ospedale_modificatoosp3 34.2512 11.6255 2.946 0.00325 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2492 degrees of freedom
## Multiple R-squared: 0.7286, Adjusted R-squared: 0.7279
## F-statistic: 955.8 on 7 and 2492 DF, p-value: < 2.2e-16Ora tutte le variabili del modello sono significative, ma la sua bonta’ non e’ ancora migliorata,percio’ provo comunque a rimuovere la variabile ospedale perche’ ai fini di una generalizzazione non la reputo una variabile interessante:
mod5 <- update(mod4, ~.-Ospedale_modificato)
summary(mod5)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1129.31 -181.70 -16.31 161.07 2638.85
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.2971 135.9911 -49.322 < 2e-16 ***
## N.gravidanze 12.7558 4.3366 2.941 0.0033 **
## Gestazione 32.2713 3.7941 8.506 < 2e-16 ***
## Lunghezza 10.2864 0.3007 34.207 < 2e-16 ***
## Cranio 10.5057 0.4260 24.659 < 2e-16 ***
## Tipo.partoNat 30.0342 12.0969 2.483 0.0131 *
## SessoM 77.9285 11.1905 6.964 4.22e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2493 degrees of freedom
## Multiple R-squared: 0.7277, Adjusted R-squared: 0.727
## F-statistic: 1110 on 6 and 2493 DF, p-value: < 2.2e-16La bonta’ del modello non e’ molto cambiata quindi, a favore della semplicita’ non dovrebbe essere stata una scelta sbagliata. Inoltre, per logica, anche il tipo di parto non dovrebbe essere importante per prevedere un peso, quindi provo a semplificare ulteriormente:
mod6 <- update(mod5, ~.-Tipo.parto)
summary(mod6)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.44 -180.81 -15.58 163.64 2639.72
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.1445 135.7229 -49.226 < 2e-16 ***
## N.gravidanze 12.4750 4.3396 2.875 0.00408 **
## Gestazione 32.3321 3.7980 8.513 < 2e-16 ***
## Lunghezza 10.2486 0.3006 34.090 < 2e-16 ***
## Cranio 10.5402 0.4262 24.728 < 2e-16 ***
## SessoM 77.9927 11.2021 6.962 4.26e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1328 on 5 and 2494 DF, p-value: < 2.2e-16Anche stavolta la bonta’ del modello non e’ andata piu’ di tanto persa, percio’ tutto sommato il mod6 mi sembra il migliore compromesso, ma per sicurezza faccio un check:
BIC(mod1, mod2, mod3, mod4, mod5, mod6)## df BIC
## mod1 12 35241.97
## mod2 10 35228.03
## mod3 10 35228.03
## mod4 9 35220.88
## mod5 8 35221.75
## mod6 7 35220.10Anche se con lievi differenze, i valori di BIC confermano il mod6 come migliore. Verifico la multicollinearita’:
car::vif(mod6)## N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1.023475 1.669189 2.074689 1.624465 1.040054Non c’e’ problema di multicollinearita’.
- Possibili effetti non lineari:
mod7 <- update(mod6, ~.+I(Lunghezza^2)+I(Cranio^2)+I(Gestazione^2)+I(N.gravidanze^2))
summary(mod7)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) +
## I(N.gravidanze^2), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1195.58 -181.99 -11.69 162.33 1469.96
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.075e+03 1.162e+03 -0.925 0.354814
## N.gravidanze 2.822e+01 7.897e+00 3.574 0.000359 ***
## Gestazione 3.692e+02 6.740e+01 5.478 4.74e-08 ***
## Lunghezza -2.912e+01 4.434e+00 -6.567 6.21e-11 ***
## Cranio -5.376e+00 9.800e+00 -0.549 0.583360
## SessoM 7.240e+01 1.099e+01 6.588 5.41e-11 ***
## I(Lunghezza^2) 4.061e-02 4.542e-03 8.942 < 2e-16 ***
## I(Cranio^2) 2.328e-02 1.444e-02 1.612 0.107101
## I(Gestazione^2) -4.282e+00 8.839e-01 -4.845 1.35e-06 ***
## I(N.gravidanze^2) -2.646e+00 1.276e+00 -2.074 0.038228 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.2 on 2490 degrees of freedom
## Multiple R-squared: 0.7399, Adjusted R-squared: 0.739
## F-statistic: 787.2 on 9 and 2490 DF, p-value: < 2.2e-16BIC(mod6, mod7)## df BIC
## mod6 7 35220.10
## mod7 11 35129.99Il modello e’ leggermente migliorato, ma l’effetto quadratico del N.gravidanze pesa meno del suo effetto lineare, quindi semplifico:
mod8 <- update(mod7, ~.-I(N.gravidanze^2))
summary(mod8)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1187.03 -182.50 -12.03 162.60 1469.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.217e+03 1.160e+03 -1.049 0.2945
## N.gravidanze 1.441e+01 4.246e+00 3.394 0.0007 ***
## Gestazione 3.754e+02 6.738e+01 5.572 2.79e-08 ***
## Lunghezza -2.925e+01 4.436e+00 -6.592 5.28e-11 ***
## Cranio -4.979e+00 9.804e+00 -0.508 0.6116
## SessoM 7.232e+01 1.100e+01 6.577 5.83e-11 ***
## I(Lunghezza^2) 4.075e-02 4.544e-03 8.967 < 2e-16 ***
## I(Cranio^2) 2.276e-02 1.445e-02 1.575 0.1154
## I(Gestazione^2) -4.374e+00 8.834e-01 -4.951 7.86e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.4 on 2491 degrees of freedom
## Multiple R-squared: 0.7395, Adjusted R-squared: 0.7387
## F-statistic: 883.9 on 8 and 2491 DF, p-value: < 2.2e-16BIC(mod7,mod8)## df BIC
## mod7 11 35129.99
## mod8 10 35126.48Inoltre, l’aggiunta dei due effetti non lineari mi potrebbe permettere di rimuovere la variabile Cranio:
mod9 <- update(mod8, ~. -Cranio -I(Cranio^2))
summary(mod9)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Sesso +
## I(Lunghezza^2) + I(Gestazione^2), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1211.40 -195.64 -15.68 186.18 2212.03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.913e+03 1.011e+03 -2.880 0.00401 **
## N.gravidanze 2.575e+01 4.719e+00 5.457 5.33e-08 ***
## Gestazione 5.109e+02 6.963e+01 7.338 2.92e-13 ***
## Lunghezza -3.327e+01 4.511e+00 -7.375 2.22e-13 ***
## SessoM 8.429e+01 1.228e+01 6.865 8.35e-12 ***
## I(Lunghezza^2) 4.813e-02 4.622e-03 10.413 < 2e-16 ***
## I(Gestazione^2) -6.001e+00 9.164e-01 -6.549 7.00e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 300 on 2493 degrees of freedom
## Multiple R-squared: 0.6742, Adjusted R-squared: 0.6735
## F-statistic: 860 on 6 and 2493 DF, p-value: < 2.2e-16BIC(mod8, mod9)## df BIC
## mod8 10 35126.48
## mod9 8 35669.66Quest’ultima non e’ stata una buona scelta (e potevo aspettarmelo visto che logicamente il parametro Cranio pesa molto sul Peso) allora la reinserisco, provando anche l’interazione con la Lunghezza:
mod10 <- update(mod8, ~.+(Lunghezza*Cranio))
summary(mod10)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) +
## Lunghezza:Cranio, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1180.84 -181.86 -12.29 164.69 1332.10
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.596e+03 1.162e+03 -1.374 0.169703
## N.gravidanze 1.434e+01 4.235e+00 3.387 0.000718 ***
## Gestazione 3.263e+02 6.846e+01 4.767 1.98e-06 ***
## Lunghezza -2.058e+01 4.991e+00 -4.123 3.86e-05 ***
## Cranio -9.700e+00 9.859e+00 -0.984 0.325271
## SessoM 7.339e+01 1.097e+01 6.690 2.75e-11 ***
## I(Lunghezza^2) 4.918e-02 5.058e-03 9.723 < 2e-16 ***
## I(Cranio^2) 6.541e-02 1.835e-02 3.564 0.000372 ***
## I(Gestazione^2) -3.744e+00 8.969e-01 -4.174 3.09e-05 ***
## Lunghezza:Cranio -4.963e-02 1.322e-02 -3.754 0.000178 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.7 on 2490 degrees of freedom
## Multiple R-squared: 0.741, Adjusted R-squared: 0.74
## F-statistic: 791.4 on 9 and 2490 DF, p-value: < 2.2e-16BIC(mod8, mod10)## df BIC
## mod8 10 35126.48
## mod10 11 35120.19Il mod10 e’ ancora leggermente migliorato, ma siccome ci sono molti termini provo a verificare se posso evitare il fattore quadratico a favore della sola interazione:
mod11 <- update(mod10, ~.-I(Lunghezza^2) -I(Cranio^2))
summary(mod11)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2) + Lunghezza:Cranio, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1154.60 -181.54 -14.63 164.19 2898.29
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.145e+03 1.060e+03 -2.025 0.04300 *
## N.gravidanze 1.299e+01 4.321e+00 3.005 0.00268 **
## Gestazione 1.156e+02 6.662e+01 1.736 0.08269 .
## Lunghezza -2.589e+00 2.947e+00 -0.878 0.37979
## Cranio -8.032e+00 4.250e+00 -1.890 0.05887 .
## SessoM 7.376e+01 1.120e+01 6.585 5.54e-11 ***
## I(Gestazione^2) -1.019e+00 8.739e-01 -1.166 0.24367
## Lunghezza:Cranio 3.820e-02 8.652e-03 4.415 1.05e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.4 on 2492 degrees of freedom
## Multiple R-squared: 0.7297, Adjusted R-squared: 0.7289
## F-statistic: 961 on 7 and 2492 DF, p-value: < 2.2e-16BIC(mod10, mod11)## df BIC
## mod10 11 35120.19
## mod11 9 35211.03Questo tentativo e’ da scartare, il mod10 rimane il miglior candidato. Un ulteriore tentativo potrebbe essere valutare l’interazione anche con Gestazione:
mod12 <- update(mod10, ~.+(Gestazione*Lunghezza)+(Gestazione*Cranio))
summary(mod12)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) +
## Lunghezza:Cranio + Gestazione:Lunghezza + Gestazione:Cranio,
## data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1181.34 -181.80 -10.71 168.09 1320.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.334e+02 1.191e+03 -0.112 0.910856
## N.gravidanze 1.480e+01 4.216e+00 3.511 0.000454 ***
## Gestazione 1.650e+02 7.505e+01 2.199 0.027967 *
## Lunghezza -4.767e+00 5.916e+00 -0.806 0.420415
## Cranio -2.287e+01 1.018e+01 -2.245 0.024830 *
## SessoM 7.327e+01 1.092e+01 6.707 2.45e-11 ***
## I(Lunghezza^2) 5.633e-02 6.507e-03 8.656 < 2e-16 ***
## I(Cranio^2) 4.980e-02 1.872e-02 2.660 0.007873 **
## I(Gestazione^2) -4.535e+00 1.613e+00 -2.811 0.004970 **
## Lunghezza:Cranio -8.489e-02 1.672e-02 -5.078 4.09e-07 ***
## Gestazione:Lunghezza -2.798e-01 1.969e-01 -1.421 0.155420
## Gestazione:Cranio 1.061e+00 2.083e-01 5.097 3.72e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.4 on 2488 degrees of freedom
## Multiple R-squared: 0.7437, Adjusted R-squared: 0.7425
## F-statistic: 656.2 on 11 and 2488 DF, p-value: < 2.2e-16BIC(mod10,mod12)## df BIC
## mod10 11 35120.19
## mod12 13 35109.50L’interazione tra Gestazione e Lunghezza non e’ significativa, ma lo e’ quella con Cranio:
mod13 <- lm(Peso ~ Gestazione+N.gravidanze+Lunghezza+Cranio+Sesso+(Cranio*Lunghezza)+(Gestazione*Cranio)+I(Lunghezza^2)+I(Cranio^2)+I(Gestazione^2), data=data )
summary(mod13)##
## Call:
## lm(formula = Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio +
## Sesso + (Cranio * Lunghezza) + (Gestazione * Cranio) + I(Lunghezza^2) +
## I(Cranio^2) + I(Gestazione^2), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1167.37 -182.19 -11.71 168.33 1318.88
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.032e+02 1.190e+03 -0.171 0.864479
## Gestazione 1.802e+02 7.431e+01 2.425 0.015395 *
## N.gravidanze 1.483e+01 4.216e+00 3.516 0.000445 ***
## Lunghezza -7.234e+00 5.657e+00 -1.279 0.201045
## Cranio -2.058e+01 1.006e+01 -2.046 0.040877 *
## SessoM 7.275e+01 1.092e+01 6.662 3.30e-11 ***
## I(Lunghezza^2) 5.048e-02 5.042e-03 10.012 < 2e-16 ***
## I(Cranio^2) 5.483e-02 1.839e-02 2.981 0.002901 **
## I(Gestazione^2) -6.298e+00 1.032e+00 -6.102 1.21e-09 ***
## Lunghezza:Cranio -9.269e-02 1.579e-02 -5.870 4.94e-09 ***
## Gestazione:Cranio 1.014e+00 2.056e-01 4.932 8.69e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.5 on 2489 degrees of freedom
## Multiple R-squared: 0.7435, Adjusted R-squared: 0.7424
## F-statistic: 721.3 on 10 and 2489 DF, p-value: < 2.2e-16BIC(mod12,mod13)## df BIC
## mod12 13 35109.5
## mod13 12 35103.7In definitiva, il mod13 sembra il migliore tra quelli che considerano effetti non lineari ed interazioni: e’ pur vero che sono stati introdotti tanti parametri in piu’, pero’ ognuno di questi ha una significativita’ alta e maggiore dello stesso parametro preso con grado 1. Per ulteriore sicurezza confronto il mod13 con il mod6, i migliori modelli trovati, anche tramite un altro criterio di informazione:
AIC(mod6, mod13)## df AIC
## mod6 7 35179.33
## mod13 12 35033.82Sia il BIC che l’AIC mi permettono di scegliere il mod13, ma per sicurezza controllo che non ci sia multicollinearita’:
car::vif(mod13)## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif## Gestazione N.gravidanze Lunghezza Cranio
## 678.640428 1.026122 780.105639 960.672823
## Sesso I(Lunghezza^2) I(Cranio^2) I(Gestazione^2)
## 1.049694 561.399700 1433.526168 721.544338
## Lunghezza:Cranio Gestazione:Cranio
## 1932.504954 1638.162089C’e’ un’elevata multicollinearita’ quindi, considerato che le bonta’ dei modelli sono comunque simili, ritorno alle origini preferendo il modello piu’ semplice: il mod6.
- Diagnostica sui residui:
par(mfrow=c(2,2))
plot(mod6)
grafico1 (Residual vs fitted): i punti sono perlopiu’ distribuiti
casualmente intorno allo 0, tranne che per una deviazione nella parte
sinistra del grafico; grafico2 (QQ residuals): i punti sembrano seguire
una distribuzione normale poiche’ perlopiu’ giacciono sulla bisettrice
del grafico, anche se ci sono deviazioni alle estremita’ grafico3
(scale-location): la varianza sembra costante, fatta eccezione che per
la parte sinistra grafico4 (Residuals vs Leverage): c’e un’osservazione
(1551) che si trova oltre la distanza di Cook
Vado a verificare quanto osservato graficamente, partendo dal test di normalita’ dei residui:
shapiro.test(residuals(mod6))##
## Shapiro-Wilk normality test
##
## data: residuals(mod6)
## W = 0.97408, p-value < 2.2e-16I residui non sono distribuiti normalmente (evidentemente, le deviazioni all’estremita’ pesano di piu’ di quanto sperassi). Test di incorrelazione:
library(lmtest)## Caricamento del pacchetto richiesto: zoo##
## Caricamento pacchetto: 'zoo'## I seguenti oggetti sono mascherati da 'package:base':
##
## as.Date, as.Date.numericdwtest(mod6)##
## Durbin-Watson test
##
## data: mod6
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0I residui non sono autocorrelati. Test di omoschedasticita’:
bptest(mod6)##
## studentized Breusch-Pagan test
##
## data: mod6
## BP = 90.253, df = 5, p-value < 2.2e-16La varianza non e’ costante.
Ricerca dei valori leverage:
lev <- hatvalues(mod6)
plot(lev)
p <- sum(lev)
n <- nrow(data)
soglia <- 2*p/n
abline(h=soglia, col=2)
length(lev[lev>soglia])## [1] 152Ho 152 valori leverage
cook <- cooks.distance(mod6)
plot(cook)
Tra cui in realta’ e’ particolarmente influente sulle stime di
regressione solamente l’osservazione 1551, la quale e’ oltre la soglia
di 0.5 (ma, almeno, entro 1).
Ricerca degli outliers:
outliers <- rstudent(mod6)
plot(outliers)
abline(h=c(-2,2), col=2)
sum(outliers>2|outliers<(-2))## [1] 101E anche 101 potenziali outliers…
outlierTest(mod6)## rstudent unadjusted p-value Bonferroni p
## 1551 10.051908 2.4906e-23 6.2265e-20
## 155 5.027798 5.3138e-07 1.3285e-03
## 1306 4.827238 1.4681e-06 3.6702e-03…di cui pero’ solo 3 sono effettivamente tali.
Insomma questo modello non e’ proprio il massimo…quindi provo a rifare tutto sul mod13, il mio secondo miglior candidato:
par(mfrow=c(2,2))
plot(mod13)
Gia’ dai grafici la situazione non mi sembra migliore, anzi,
tutt’altro…
shapiro.test(residuals(mod13))##
## Shapiro-Wilk normality test
##
## data: residuals(mod13)
## W = 0.99103, p-value = 2.306e-11dwtest(mod13)##
## Durbin-Watson test
##
## data: mod13
## DW = 1.9586, p-value = 0.1501
## alternative hypothesis: true autocorrelation is greater than 0bptest(mod13)##
## studentized Breusch-Pagan test
##
## data: mod13
## BP = 39.829, df = 10, p-value = 1.816e-05…e neanche il mod13 supera i test di normalita’ e omoschedasticita’ dei residui.
A questo punto provo a rimuovere l’osservazione n.1551, la piu’ problematica, per provare a capire se il problema e’ dovuto al modello o a questa forte influenza su di esso:
remove_row <- which(cooks.distance(mod6)>0.5)
data_filtered <- data[-remove_row,]
model_filtered <- update(mod6, data=data_filtered)
summary(model_filtered)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = data_filtered)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1165.74 -179.59 -12.74 162.89 1410.88
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6683.4142 133.0802 -50.221 < 2e-16 ***
## N.gravidanze 13.1652 4.2557 3.094 0.002 **
## Gestazione 29.5891 3.7340 7.924 3.43e-15 ***
## Lunghezza 10.8927 0.3017 36.109 < 2e-16 ***
## Cranio 9.9187 0.4225 23.476 < 2e-16 ***
## SessoM 78.1348 10.9840 7.114 1.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.3 on 2493 degrees of freedom
## Multiple R-squared: 0.7372, Adjusted R-squared: 0.7367
## F-statistic: 1399 on 5 and 2493 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(model_filtered)
shapiro.test(residuals(model_filtered))##
## Shapiro-Wilk normality test
##
## data: residuals(model_filtered)
## W = 0.98886, p-value = 4.764e-13dwtest(model_filtered)##
## Durbin-Watson test
##
## data: model_filtered
## DW = 1.954, p-value = 0.1251
## alternative hypothesis: true autocorrelation is greater than 0bptest(model_filtered)##
## studentized Breusch-Pagan test
##
## data: model_filtered
## BP = 11.393, df = 5, p-value = 0.04411In realta’, neanche questa rimozione migliora il modello (peraltro probabilmente non e’ in generale una scelta del tutto idonea, ma volevo fare un tentativo). Direi che un modello di regressione lineare non puo’ spiegare abbastanza bene questo dataset e che un GLM sarebbe piu’ adatto, magari attraverso una link function logaritmica.
Considerato il valore dell’Adjusted-R2, la bonta’ del mod6 e’ del 72,65%.
Per fare una previsione in base alle 3 variabili disponibili dovrei utilizzare un modello costruito solo su di esse, ma non fornirebbe un risultato per nulla accurato, quindi utilizzo il mio “miglior” modello (mod6) imputando i valori mancanti con le rispettive medie condizionate al sesso, dato che tra i sessi ci sono differenze significative:
predict(mod6, data.frame(Gestazione=39, Sesso="F", N.gravidanze=3, Lunghezza=mean(data$Lunghezza[data$Sesso=="F"]), Cranio=mean(data$Cranio[data$Sesso=="F"])))## 1
## 3195.331- Grafici finali del mod6:
ggplot(data)+
geom_point(aes(x=Gestazione,
y=Peso,
col=Sesso), position="jitter")+
geom_smooth(aes(x=Gestazione,
y=Peso,
col=Sesso), se=F, method="lm")+
labs(title="Modello di previsione del Peso dei neonati",
x="N. di settimane di gestazione",
y="Peso")## `geom_smooth()` using formula = 'y ~ x'
ggplot(data)+
geom_point(aes(x=Lunghezza,
y=Peso,
col=Sesso), position="jitter")+
geom_smooth(aes(x=Lunghezza,
y=Peso,
col=Sesso), se=F, method="lm")+
labs(title="Modello di previsione del Peso dei neonati",
x="Lunghezza(mm)",
y="Peso")## `geom_smooth()` using formula = 'y ~ x'
ggplot(data)+
geom_point(aes(x=Cranio,
y=Peso,
col=Sesso), position="jitter")+
geom_smooth(aes(x=Cranio,
y=Peso,
col=Sesso), se=F, method="lm")+
labs(title="Modello di previsione del Peso dei neonati",
x="Diametro del cranio(mm)",
y="Peso")## `geom_smooth()` using formula = 'y ~ x'
ggplot(data)+
geom_point(aes(x=N.gravidanze,
y=Peso,
col=Sesso), position="jitter")+
geom_smooth(aes(x=N.gravidanze,
y=Peso,
col=Sesso), se=F, method="lm")+
labs(title="Modello di previsione del Peso dei neonati",
x="N. di gravidanze sostenute",
y="Peso")## `geom_smooth()` using formula = 'y ~ x'
Da questi grafici sembra che effettivamente il N.gravidanze non
influisca particolarmente sul peso (pendenza ~ 0), mentre le altre 3
variabili sono molto piu’ significative.