1. Analisi Preliminare Dataset

È stata aggiunta una variabile al dataset, ovvero l’Indice Ponderale, derivato dal Peso e dalla Lunghezza. Questa variabile è data da Peso/Lunghezza^3 e aggiunge informazioni significative sulle proporzioni del neonato. Per le variabili analizzate è stato usato il Q-Q Plot per verificare quanto le distribuzioni si discostano da una distribuzione normale teorica e assumere che rnorm è il modo corretto per visualizzare la distribuzione delle variabili.

Obiettivo

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

Dataset Samples

#' Crea una tabella HTML formattata per la visualizzazio = e di dati statistici
#'
#' @param data Data frame o oggetto statistico da visualizzare
#' @param show_head_only Se TRUE, mostra solo le prime righe del data frame
#' @param show_data_types Se TRUE, mostra i tipi di dati nelle intestazioni delle colonne
#'
#' @return Una tabella HTML formattata con kableExtra
#' @export
#'
#' @examples
#' # create_table(risultati_test)
#' # create_table(dati, show_head_only = FALSE)
create_table <- function(data, show_head_only = TRUE, show_data_types = TRUE) {
  if (inherits(data, "htest")) {
    print(data)

    entries <- list()

    # Aggiungi statistiche (possono essere più di una)
    if (!is.null(data$statistic)) {
      for (i in seq_along(data$statistic)) {
        entries[[length(entries) + 1]] <- list(Statistic = names(data$statistic)[i], Value = data$statistic[[i]])
      }
    }

    entries[[length(entries) + 1]] <- list(Statistic = "p.value", Value = data$p.value)

    if (!is.null(data$parameter)) {
      for (i in seq_along(data$parameter)) {
        entries[[length(entries) + 1]] <- list(Statistic = names(data$parameter)[i], Value = data$parameter[[i]])
      }
    }

    if (!is.null(data$estimate)) {
      for (i in seq_along(data$estimate)) {
        entries[[length(entries) + 1]] <- list(Statistic = names(data$estimate)[i], Value = data$estimate[[i]])
      }
    }

    if (!is.null(data$conf.int)) {
      entries[[length(entries) + 1]] <- list(Statistic = "conf.int_lower", Value = data$conf.int[1])
      entries[[length(entries) + 1]] <- list(Statistic = "conf.int_upper", Value = data$conf.int[2])
    }

    data <- do.call(rbind, lapply(entries, as.data.frame))

  } else if (inherits(data, "summary.lm")) {
    coeffs <- data$coefficients
    coeffs_df <- data.frame(
      Term = rownames(coeffs),
      Estimate = coeffs[, "Estimate"],
      Std_Error = coeffs[, "Std. Error"],
      t_value = coeffs[, "t value"],
      p_value = coeffs[, "Pr(>|t|)"],
      row.names = NULL
    )

    metrics <- data.frame(
      Term = c("Residual Std. Error", "Multiple R-squared", "Adjusted R-squared", "F-statistic", "Model DF", "p-value"),
      Estimate = c(data$sigma,
                   data$r.squared,
                   data$adj.r.squared,
                   data$fstatistic[1],
                   paste(data$fstatistic[2], "and", data$fstatistic[3]),
                   pf(data$fstatistic[1], data$fstatistic[2], data$fstatistic[3], lower.tail = FALSE)),
      Std_Error = NA, t_value = NA, p_value = NA
    )

    data <- rbind(coeffs_df, metrics)

  } else if (inherits(data, "aov")) {
    summary_data <- summary(data)[[1]]
    data <- data.frame(
      Statistic = c("F value", "p-value", "num Df", "denom Df"),
      Value = c(
        summary_data$`F value`[1],
        summary_data$`Pr(>F)`[1],
        summary_data$Df[1],
        summary_data$Df[2]
      ),
      row.names = NULL
    )

  } else if (inherits(data, "summaryDefault") || inherits(data, "table")) {
    data <- data.frame(
      Statistic = names(data),
      Value = as.numeric(data),
      row.names = NULL
    )
  } else {
  if (is.numeric(data) && !is.null(names(data))) {
      data <- data.frame(Variable = names(data), VIF = unname(data), row.names = NULL)
    } else if (is.matrix(data) && !is.null(rownames(data))) {
      data <- data.frame(Variable = rownames(data), data, row.names = NULL)
    } else {
      data <- as.data.frame(data)
    }
}
  column_names <- names(data)
  column_types <- sapply(data, class)

  if (!show_data_types) {
    column_headers <- column_names
  } else {
    column_headers <- paste0(column_names, " (", column_types, ")")
  }

  row_length <- if (show_head_only) head(data) else data

  knitr::kable(row_length,
               digits = 4,
               format = "html",
               row.names = TRUE,
               col.names = column_headers) %>%
    kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),
                              full_width = TRUE,
                              font_size = 15,
                              position = "left")
}

#' Crea una tabella di contingenza stilizzata
#'
#' @param var1 Variabile per le righe della tabella
#' @param var2 Variabile per le colonne della tabella
#' @param caption Titolo della tabella (opzionale)
#' @param gruppo_colonne Etichetta per il gruppo di colonne (opzionale)
#' @param colore_totali Colore di sfondo per le celle dei totali (opzionale)
#' @param includi_totali Booleano che indica se includere i totali (opzionale, default TRUE)
#'
#' @return Un oggetto tabella kableExtra
#' @export
#'
#' @examples
#' # Esempio di utilizzo:
#' # crea_tabella_contingenza_stilizzata(Ospedale, Tipo.parto, "Distribuzione dei parti per ospedale")
crea_tabella_contingenza_stilizzata <- function(var1, var2, 
                                               caption = NULL,
                                               gruppo_colonne,
                                               colore_totali = "#E8E8E8",
                                               includi_totali = TRUE) {
  
  tabella_contingenza <- table(var1, var2)
  
  if (includi_totali) {
    tabella_con_totali <- cbind(tabella_contingenza, Totale = rowSums(tabella_contingenza))
    
    tabella_con_totali <- rbind(tabella_con_totali, 
                              Totale = colSums(tabella_con_totali))
    
    tabella_finale <- tabella_con_totali
  } else {
    tabella_finale <- tabella_contingenza
  }
  
  tabella_stilizzata <- kable(tabella_finale, 
                             caption = caption,
                             format = "html", 
                             align = "c") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  full_width = TRUE,  
                  font_size = 16)
  
  if (includi_totali) {
    tabella_stilizzata <- tabella_stilizzata %>%
      row_spec(nrow(tabella_finale) - 1, bold = TRUE, background = colore_totali)
    
    tabella_stilizzata <- tabella_stilizzata %>%
      column_spec(ncol(tabella_finale), bold = TRUE, background = colore_totali)
  }
  
  num_colonne_var2 <- ncol(tabella_contingenza)
  header_spec <- c(" " = 1)
  header_spec <- c(header_spec, setNames(num_colonne_var2, gruppo_colonne))
  
  
  if (includi_totali) {
    header_spec <- c(header_spec, " " = 1)
  }
  
  tabella_stilizzata <- tabella_stilizzata %>%
    add_header_above(header_spec)
  
  return(tabella_stilizzata)
}

#' Calcola indici di forma di una distribuzione
#'
#' @param var_name Vettore numerico da analizzare
#' @return Un data frame contenente l'indice di Gini, l'indice di asimmetria e l'indice di curtosi
#' @export
#'
#' @examples
#' # evaluate_shape_indeces(dati$variabile)
evaluate_shape_indeces <- function(var_name){
  gini_value <- gini.index(var_name)
  skew_value <- skewness(var_name)
  kurt_value <- kurtosis(var_name)-3
results <- data.frame(
  Statistic = c("Gini_Index", "Skewness", "Kurtosis"),
  Value = c(gini_value, skew_value, kurt_value)
)
return(results)
}

#' Calcola l'indice di Gini normalizzato
#'
#' @param x Vettore di dati da analizzare
#' @return L'indice di Gini normalizzato
#' @export
#'
#' @examples
#' # gini.index(dati$variabile)
gini.index <- function(x) {
  ni=table(x)
  fi=ni/length(x)
  fi2=fi^2
  J = length(table(x))
  gini = 1-sum(fi2)
  gini.normalizzato = gini/((J-1)/J)
  return(gini.normalizzato)
}

#' Calcola il coefficiente di variazione
#'
#' @param x Vettore numerico per cui calcolare il CV
#' @return Il coefficiente di variazione espresso in percentuale
#' @export
#'
#' @examples
#' # CV(dati$variabile)
CV <- function(x) {
  return(sd(x)/mean(x)*100)
}

#' Sposta una variabile in modo che tutti i valori siano positivi
#'
#' @param var Vettore numerico da spostare
#' @return Vettore con valori spostati in modo che il minimo sia 1
#' @export
#'
#' @examples
#' # var_shifted_to_positive(dati$variabile_con_negativi)
var_shifted_to_positive <- function(var) {
  return( var + abs(min(var)) + 1 )
}

#' Crea un grafico della variazione di una variabile nel tempo per gruppi diversi
#'
#' @param data Data frame contenente i dati
#' @param x_var Nome della variabile per l'asse x
#' @param y_var Nome della variabile per l'asse y
#' @param group_var Nome della variabile di raggruppamento
#' @param x_label Etichetta per l'asse x
#' @param y_label Etichetta per l'asse y
#' @param title Titolo del grafico
#' @param lower_limit Valore limite inferiore da rappresentare con linea tratteggiata
#' @param upper_limit Valore limite superiore da rappresentare con linea tratteggiata
#' @param lower_text Testo da mostrare vicino al limite inferiore
#' @param upper_text Testo da mostrare vicino al limite superiore
#' @param y_min Valore minimo per l'asse y
#' @param y_max Valore massimo per l'asse y
#'
#' @return Un oggetto ggplot
#' @export
#'
#' @examples
#' # plot_variation_by_group(dati, "settimana", "peso", "gruppo", 
#' #                         lower_limit = 2500, upper_limit = 4000,
#' #                         lower_text = "Peso minimo", upper_text = "Peso massimo")
plot_variation_by_group <- function(data, x_var, y_var, group_var,
                                    x_label = "Settimane di gestazione",
                                    y_label = "Indice ponderale",
                                    title = "Variazione della variabile nel tempo",
                                    lower_limit = NULL,
                                    upper_limit = NULL,
                                    lower_text = NULL,
                                    upper_text = NULL,
                                    y_min = NULL,
                                    y_max = NULL) {
  x_sym <- sym(x_var)
  y_sym <- sym(y_var)
  group_sym <- sym(group_var)
  data[[as_string(group_sym)]] <- as.factor(data[[as_string(group_sym)]])
  p <- ggplot(data, aes(x = !!x_sym, y = !!y_sym, color = !!group_sym)) +
    geom_point(alpha = 0.6) +
    geom_smooth(method = "loess", se = FALSE, linewidth = 1.2) +
    labs(x = x_label, y = y_label, title = title, color = "Gruppo") +
    theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5))
  x_min <- min(data[[x_var]], na.rm = TRUE)
  if (!is.null(lower_limit)) {
    p <- p +
      geom_hline(yintercept = lower_limit, color = "blue", linetype = "dashed") +
      annotate("text", x = x_min, y = lower_limit, label = lower_text,
               vjust = -1, hjust = 0, color = "blue")
  }
  if (!is.null(upper_limit)) {
    p <- p +
      geom_hline(yintercept = upper_limit, color = "red", linetype = "dashed") +
      annotate("text", x = x_min, y = upper_limit, label = upper_text,
               vjust = -1, hjust = 0, color = "red")
  }
  # Per impostare dei limiti di visualizzazione
  if (!is.null(y_min) && !is.null(y_max)) {
    p <- p + coord_cartesian(ylim = c(y_min, y_max))
  }
  return(p)
}

#' Esegue un test Z per un campione
#'
#' @param x Vettore di dati del campione
#' @param mu0 Media ipotizzata sotto H0
#' @param stdev Deviazione standard della popolazione
#' @param alfa Livello di significatività
#'
#' @return Una lista contenente: media campionaria, statistica test, valori soglia, p-value, intervallo di confidenza e grafico
#' @export
#'
#' @examples
#' # z_test(dati$variabile, 100, 15, 0.05)
z_test <- function(x,mu0, stdev, alfa){
  x = na.omit(x)
  mu_cap = mean(x)
  n = length(x)
  Z = (mu_cap-mu0) / (stdev / sqrt(n))
  valori.soglia = qnorm(c(alfa/2, 1-alfa/2))
  CI = c(mu_cap+qnorm(alfa/2)*(stdev/sqrt(n)), 
         mu_cap-qnorm(alfa/2)*(stdev/sqrt(n)))
  return(
    list(
      media.campionaria = mu_cap,
      stat.test = Z,
      valori.soglia = valori.soglia,
      pvalue = 2*pnorm(-abs(Z)),
      Int.Conf=CI,
      grafico = plot(density(rnorm(1000000,0,1))),
      abline(v=valori.soglia,col=2)
    )
  )
}

#' Crea una tabella di frequenze assolute e relative
#'
#' @param var Vettore di dati da analizzare
#' @return Un data frame contenente frequenze assolute, relative e cumulate
#' @export
#'
#' @examples
#' # tab_frequenze(dati$categoria)
tab_frequenze <- function(var){
  n <- length(var)
  distr_freq <- as.data.frame(
    cbind(
      ni=table(var),
      fi=table(var)/n,
      Ni=cumsum(table(var)),
      Fi=cumsum(table(var))/n
    )
  )
}

#' Pannello per visualizzare la correlazione in una coppia di variabili
#'
#' Funzione pensata da usare all'interno di una matrice di scatterplot (es. pairs()),
#' per scrivere il coefficiente di correlazione tra due variabili in una casella.
#'
#' @param x Primo vettore numerico da correlare
#' @param y Secondo vettore numerico da correlare
#' @param digits Numero di cifre decimali da visualizzare (default 2)
#' @param prefix Prefisso da aggiungere prima del valore di correlazione (default "")
#' @param cex.cor Fattore di scala per la dimensione del testo (calcolato automaticamente se non specificato)
#' @param ... Altri parametri grafici opzionali
#' 
#' @return Nessun valore ritornato. Disegna direttamente il coefficiente di correlazione sulla grafica.
#'
#' @examples
#' # pairs(mtcars, lower.panel = panel.cor)
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt) 
  text(0.5, 0.5, txt, cex = cex.cor * r)  
}


#' Seleziona variabili numeriche con più di 2 valori distinti dal dataset
#'
#' @param data Data frame contenente i dati
#'
#' @return Data frame contenente solo le variabili numeriche con più di 2 valori distinti
#' @export
#'
#' @examples
#' # vars_num <- seleziona_variabili_numeriche(dati)
seleziona_variabili_numeriche <- function(data) {
  data[, sapply(data, function(x) is.numeric(x) && length(unique(x)) > 2)]
}

#' Disegna una matrice di correlazione per le variabili numeriche del dataset
#'
#' @param data Data frame contenente i dati
#'
#' @return Un grafico a matrice con correlazioni nella parte inferiore e curve di dispersione nella parte superiore
#' @export
#'
#' @examples
#' # disegna_matrice_correlazione(dati)
disegna_matrice_correlazione <- function(data) {
  numeric_vars <- seleziona_variabili_numeriche(data)
  pairs(numeric_vars, lower.panel = panel.cor, upper.panel = panel.smooth)
}

#' Confronta le varianze tra gruppi per due variabili rispetto alla variabile 'Peso'
#'
#' @param data Data frame contenente i dati
#' @param var1 Nome della prima variabile da testare
#' @param var2 Nome della seconda variabile da testare
#'
#' @return Output dei test di varianza per le due variabili specificate
#' @export
#'
#' @examples
#' # confronta_varianze(dati, "Sesso", "Tipo.parto")
confronta_varianze <- function(data, var1, var2) {
  cat("\n\n\n\n<b>Test confronto varianza: Peso ~ ", var1, "</b>\n")
  cat("\nLa varianza tra il peso e il sesso risulta essere uguale. Questo è stato considerato nei test statistici\n")
  print(create_table(var.test(as.formula(paste("Peso ~", var1)), data = data),F,F))
  
  cat("<b>Test confronto varianza: Peso ~ ", var2, "</b>\n")
  cat("\nLa varianza tra il peso e l'essere fumatrici  risulta essere uguale. Questo è stato considerato nei test statistici\n")
  print(create_table(var.test(as.formula(paste("Peso ~", var2)), data = data), F, F))
}


#' Crea una sequenza di modelli lineari ridotti rimuovendo variabili una dopo l'altra
#'
#' @param data Data frame contenente i dati
#'
#' @return Lista di modelli lineari, partendo dal modello completo fino al più ridotto
#' @export
#'
#' @examples
#' # modelli <- crea_modelli_ridotti(dati)
crea_modelli_ridotti <- function(data) {
  # Creare una copia del dataframe per non modificare l'originale
  data_mod <- data
  
  # Codifica Sesso come variabile dummy (0 per M, 1 per F)
  data_mod$Sesso <- ifelse(data_mod$Sesso == "F", 1, 0)
  
  
  # Creare i modelli con il dataset modificato
  modelli <- list()
  
  modelli[[1]] <- lm(I(sqrt(Peso)) ~ ., data = data_mod)
  modelli[[2]] <- update(modelli[[1]], ~ . - Ospedale - Lunghezza - Cranio + Lunghezza:Cranio)
  modelli[[3]] <- update(modelli[[2]], ~ . - Tipo.parto)
  modelli[[4]] <- update(modelli[[3]], ~ . - Anni.madre)
 

  names(modelli) <- paste0("mod", 1:4)
  return(modelli)
}

#' Rimuove gli outlier in modo iterativo basandosi sulla distanza di Cook. Vengono eseguite due iterazioni perchè nella prima vengono identificati gli outlier rispetto al modello originale, mentre nella seconda, dopo aver rimosso gli outlier, il modello viene ricalcoalto e potrebbe avere nuovi outlier che prima non erano rilevabili
#'
#' @param modello Modello lineare iniziale
#' @param data Data frame contenente i dati
#' @param formula_modello Formula del modello per ricalcolare dopo la rimozione degli outlier
#'
#' @return Data frame pulito senza gli outlier identificati in due iterazioni
#' @export
#'
#' @examples
#' # dati_puliti <- rimuovi_outliers_iterativi(modello, dati, log(Peso) ~ .)
rimuovi_outliers_iterativi <- function(modello, data, formula_modello) {
  soglia <- 4 / nrow(data)
  cook_distanze <- cooks.distance(modello)
  outlier_indices <- which(cook_distanze > soglia)
  cat("<b>Rimozione Outlier</b> \n\n")
    cat("È stato rimosso circa il 7% delle osservazioni rilevate come outlier statisticamente significativi. Gli outlier più influenti statisticamente sono stati identificati attraverso un valore soglia scelto empiricamente con 4/n e da cui è stata calcolata la distanza di Cook. Vengono eseguite due iterazioni perchè nella prima vengono identificati gli outlier rispetto al modello originale, mentre nella seconda, dopo aver rimosso gli outlier, il modello viene ricalcolato e potrebbe avere nuovi outlier che prima non erano rilevabili. Infatti, è stato scelto di rimuovere ulteriori 4 outlier singificativi per migliorare il modello. \n\n")
  total_removed <- 0
  
  if (length(outlier_indices) > 0) {
    cat("Rimosse<b>", length(outlier_indices), "</b>osservazioni (prima iterazione)\n\n")
    data <- data[-outlier_indices, ]
    total_removed <- total_removed + length(outlier_indices)

  }
  
  modello_pulito <- lm(formula_modello, data = data)
  soglia2 <- 4 / nrow(data)
  cook_distanze2 <- cooks.distance(modello_pulito)
  outlier_indices2 <- which(cook_distanze2 > soglia2)
  
  outliers_specifici <- c(684, 1188, 1712, 1130)
  outliers_specifici_presenti <- outliers_specifici[outliers_specifici %in% as.numeric(rownames(data))]
  outlier_indices2 <- unique(c(outlier_indices2, which(rownames(data) %in% outliers_specifici_presenti)))

  if (length(outlier_indices2) > 0) {
    cat("Rimosse <b>", length(outlier_indices2), "</b> osservazioni (seconda iterazione)\n\n")
    data <- data[-outlier_indices2, ]
    total_removed <- total_removed + length(outlier_indices2)
  }


cat("Totale osservazioni rimosse: <b>", total_removed ,"</b>\n\n" )
  
  return(data)
}

#' Esegue una diagnostica completa di un modello lineare
#'
#' @param modello Modello lineare da diagnosticare
#'
#' @return Grafici diagnostici e risultati di test vari (eteroschedasticità, autocorrelazione, VIF, outlier)
#' @export
#'
#' @examples
#' # diagnostica_modello(modello_finale)
diagnostica_modello <- function(modello) {
  par(mfrow=c(1,1))
  cat("\n\n<b> VIF</b>\n")
  cat("\nI VIF del mod4 sono inferiori a 5 per tutti i regressori indica che non ci sono problemi di multicollinearità, ovvero di correlazione tra i regressori ")

  vif_values <- car::vif(modello, type = "predictor")
var_names <- names(coef(modello))[-1]  
if (is.matrix(vif_values)) {
  vif_vector <- vif_values[, 1]
} else {
  vif_vector <- vif_values
}
names(vif_vector) <- var_names

y_max <- 5
vif_vector_no_last <- vif_vector[-length(vif_vector)]

bp <- barplot(vif_vector_no_last,
        las = 2,
        col = "skyblue",
        main = "Variance Inflation Factor (VIF)",
        cex.names = 0.6,
        ylim = c(0, y_max))


text(x = bp, 
     y = vif_vector_no_last + max(vif_vector_no_last)*0.05,
     labels = round(vif_vector_no_last, 2), 
     cex = 0.8)
  par(mfrow=c(2,2))
  cat("\n\n\n<b>Analisi dei residui mod4</b>\n\n")
  cat("Dai seguenti plot dei residui è emerso che:\n\n- <b>Residuals vs Fitted</b>: Questo grafico mostra i residui rispetto ai valori predetti. I residui si distribuiscono casualmente attorno alla linea orrizontale zero, senza pattern evidenti. La variabilità risulta essere costante\n- <b>Q-Q Residuals</b>: Questo grafico verifica la normalità dei residui. I residui seguono la linea tratteggiata obliqua, anche se le code sono leggermente più pesanti\n- <b>Scale Location</b>: Questo grafico mostra la radice quadrata dei residui standardizzati rispetto ai valori predetti, ed è utile per verificare l'omoschedasticità. Nel grafico risulta una leggere eteroschedasticità. Tuttavia, con il BP test è stata accettata l'ipotesi di omoschedasticità\n- <b>Residuals vs Leverage</b>: Questo grafico aiuta a identificare punti influenti nel modello poichè i punti che si trovano lontano dal gruppo principale e hanno un alto leverage potrebbero influenzare significativamente i risultati del modello. Si nota come ci siano dei punti outlier che potrebbero essere significativi e da poter valutare e rimuovere se necessario")
  plot(modello)
  plot(residuals(modello), main="Residui del Modello di Regressione", xlab="Osservazioni", ylab="Residui")
  abline(h=mean(residuals(modello)), col=2)
  plot(density(residuals(modello)))
  
  cat("\n\n<b>Breusch-Pagan Test</b>\n\n")
  cat("Un p-value molto alto indica che si può accettare l'ipotesi nulla di omoschedasticità\n")
  print(create_table(lmtest::bptest(modello),F,F))
  
  cat("\n\n<b>Durbin-Watson Test</b>\n")
  cat("\n il p-value molto alto indica che si può accettare l'ipotesi nulla di non correlazione tra i regressori.")
  print(create_table(lmtest::dwtest(modello),F,F))



  
  cat("\n<b>Outlier Test</b>\n\n")
  cat("\nValutazione del residuo standardizzato usato per identificare i possibili outliers significativi in un modello di regressione. In questo caso il Largest r-student e l'unadjusted p-value suggeriscono che ci sia almeno un outlier significativo. Tuttavia, con il criterio Bonferroni(tiene conto del numero di test effettuati) e il relativo calcolo del p-value,che in questo caso si attesta all'incirca su 0.07, si potrebbe escludere l'ipotesi che ci sia almeno un outlier significativo. Questo punto non incide significativamente sulle prestazioni del modello.\n")
  outlier_results <- car::outlierTest(modello)
outlier_data <- outlier_results

outlier_df <- data.frame(
  "Row" = names(outlier_data$rstudent),
  "Rstudent" = outlier_data$rstudent,
  "P_value" = outlier_data$p,
  "Bonferroni_P_value" = outlier_data$bonf.p,
  "Significant" = outlier_data$signif,
  "Cutoff" = outlier_data$cutoff
)

print(kable(outlier_df, caption = "Outlier Test Results"))
  
leverage <- hatvalues(modello)

leverage_threshold <- 2 * mean(leverage)

plot(leverage, main="Leverage Points", xlab="Osservazioni", ylab="Leverage", pch=19, col=ifelse(leverage > leverage_threshold, "red", "blue"))

abline(h = leverage_threshold, col = 2, lwd = 2)

leverage_high <- which(leverage > leverage_threshold)
  
  cat("\n<b>Distanza di Cook, Residui e Leverage</b>\n")
  cat("\nI residui si distribuiscono normalmente e risultano avere varianza costante. Non sono presenti punti di leverage significativi. Il punto con massima distanza di cook ha distanza molto inferiore a 0.5")
  cook <- cooks.distance(modello)

cook <- cooks.distance(modello)

max_cook_index <- which.max(cook)

plot(cook, main="Cook's Distance", xlab="Osservazioni", ylab="Cook's Distance", pch=19)

abline(h = max(cook), col = 2, lwd = 2, lty = 2)

points(max_cook_index, cook[max_cook_index], col = "red", pch = 19, cex = 1.5)

text(max_cook_index, cook[max_cook_index], labels = "Max Distance", pos = 4, col = "red", cex = 0.8)

cat("\n\n<b>Max Cook's Distance at index", max_cook_index, ": </b>", round(cook[max_cook_index], 4), "\n")

}

#' Confronta modelli tramite BIC e ANOVA
#'
#' @param modelli Lista di modelli lineari da confrontare
#'
#' @return Output del confronto BIC e del test ANOVA tra gli ultimi due modelli
#' @export
#'
#' @examples
#' # confronta_modelli(lista_modelli)
confronta_modelli <- function(modelli) {
  bic_values <- sapply(1:length(modelli), function(i) BIC(modelli[[i]]))
  names(bic_values) <- paste0("mod", 1:length(modelli))

  cat("\n<b>ANOVA confronto tra gli ultimi due modelli</b>\n\n")
  cat("\nI valori del test statistico evidenziano come il mod4 abbia una capacità più elevata di spiegare la variabilità nei dati. 
Dunque, il parametro che è stato rimosso dal modello 3 al modello 4 non è statisticamente significativo. 
Questo rende il modello 4 più semplice e statisticamente rilevante.\n\n")

  anova_result <- anova(modelli[[length(modelli)-1]], modelli[[length(modelli)]])
  print(create_table(anova_result,F,F))

  mse_values <- numeric(length(modelli))
  rmse_values <- numeric(length(modelli))

  for (i in 1:length(modelli)) {
    residui <- residuals(modelli[[i]])
    MSE <- sum(residui^2) / (length(Peso) - length(coef(modelli[[i]])))
    RMSE <- sqrt(MSE)

    mse_values[i] <- MSE
    rmse_values[i] <- RMSE
  }

  names(mse_values) <- paste0("mod", 1:length(modelli))
  names(rmse_values) <- paste0("mod", 1:length(modelli))
  
  cat("<b>Confronto tra modelli: BIC, MSE, RMSE</b>\n\n")
  cat("L'MSE e l'RMSE non cambia tra tutti i modelli. Questo indica che semplificando il modello la variabilità nei dati non peggiora significativamente. Pertanto, è adatto scegliere l'ultimo modello che risulta essere meno complesso. Inoltre, Il BIC più basso di 7 punti risulta essere significativo poichè indica che la qualità dell'adattamento ai dati è migliore.\n\n")

  df_metrics <- tibble(
    Modello = rep(names(bic_values), 3),
    Valore = c(bic_values, mse_values, rmse_values),
    Statistica = factor(rep(c("BIC", "MSE", "RMSE"), each = length(modelli)),
                        levels = c("BIC", "MSE", "RMSE"))
  )
  
  ggplot(df_metrics, aes(x = Modello, y = Valore, fill = Statistica)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~Statistica, scales = "free_y") +
  labs(title = "Confronto tra i modelli", y = "Valore", x = "Modello") +
  theme_minimal() +
  theme(legend.position = "none") +
  scale_fill_brewer(palette = "Set2") +
  geom_text(aes(label = round(Valore, 1)), vjust = -0.5, size = 3.5)


}


generate_random_data <- function(original_data, n = 2) {
  data.frame(
    Gestazione = sample(min(original_data$Gestazione):max(original_data$Gestazione), n, replace = TRUE),
    Lunghezza = runif(n, min = min(original_data$Lunghezza), max = max(original_data$Lunghezza)),
    Cranio = runif(n, min = min(original_data$Cranio), max = max(original_data$Cranio)),
    Sesso = sample(min(original_data$Sesso):max(original_data$Sesso), n, replace = TRUE),
    N.gravidanze = sample(min(original_data$N.gravidanze):max(original_data$N.gravidanze), n, replace = TRUE),
    Fumatrici = sample(min(original_data$Fumatrici):max(original_data$Fumatrici), n, replace = TRUE),
    Indice_Ponderale = runif(n, min = min(original_data$Indice_Ponderale), max = max(original_data$Indice_Ponderale))
  )
}
lunghezza_cm_converted <- Lunghezza/10
dati_neonati <- dati_neonati %>%
  mutate(Indice_Ponderale = (Peso * 100) / (lunghezza_cm_converted^3))
create_table(dati_neonati, TRUE, TRUE)
Anni.madre (integer) N.gravidanze (integer) Fumatrici (integer) Gestazione (integer) Peso (integer) Lunghezza (integer) Cranio (integer) Tipo.parto (character) Ospedale (character) Sesso (character) Indice_Ponderale (numeric)
1 26 0 0 42 3380 490 325 Nat osp3 M 2.8730
2 21 2 0 39 3150 490 345 Nat osp1 F 2.6775
3 34 3 0 38 3640 500 375 Nat osp2 M 2.9120
4 28 1 0 41 3690 515 365 Nat osp2 M 2.7015
5 20 0 0 38 3700 480 335 Nat osp3 F 3.3456
6 32 0 0 40 3200 495 340 Nat osp2 F 2.6384

Dataset Description Analysis

Anni Madre

La variabile Anni.madre si distribuisce normalmente e la classe modale è quella tra i 25-29 anni ed ha eterogeneità elevata.

anni_madre_classi_cliniche <- cut(
  Anni.madre,
  breaks = c(-Inf, 14, 19, 24, 29, 34, 39, 44, Inf),
  labels = c(
    "< 15 anni",
    "15–19 anni",
    "20–24 anni",
    "25–29 anni",
    "30–34 anni",
    "35–39 anni",
    "40–44 anni",
    "45+ anni"
  ),
  right = TRUE,
  include.lowest = TRUE
)


freq_anni_madre_cl <- tab_frequenze(anni_madre_classi_cliniche)

create_table(freq_anni_madre_cl, FALSE, FALSE)
ni fi Ni Fi
< 15 anni 5 0.0020 5 0.0020
15–19 anni 106 0.0424 111 0.0444
20–24 anni 486 0.1944 597 0.2388
25–29 anni 907 0.3628 1504 0.6016
30–34 anni 712 0.2848 2216 0.8864
35–39 anni 236 0.0944 2452 0.9808
40–44 anni 46 0.0184 2498 0.9992
45+ anni 2 0.0008 2500 1.0000 
create_table(summary(Anni.madre), FALSE, FALSE)
Statistic Value
1 Min. 0.000
2 1st Qu. 25.000
3 Median 28.000
4 Mean 28.164
5 3rd Qu. 32.000
6 Max. 46.000 
# Shape indeces
anni.madre_shape_indeces <- evaluate_shape_indeces(Anni.madre)
create_table(anni.madre_shape_indeces, FALSE, FALSE)
Statistic Value
1 Gini_Index 0.9727
2 Skewness 0.0428
3 Kurtosis 0.3804 
#graphics
# verifico la presenza di valori nulli nel dataset
sum(is.na(dati_neonati))
## [1] 0
n <- nrow(dati_neonati)
distr_anni_madre <- rnorm(n, mean(Anni.madre), sd(Anni.madre))


hist(distr_anni_madre,freq = F, breaks = 100)
lines(density(distr_anni_madre),col=2, lwd=3)

qqnorm(distr_anni_madre)
qqline(distr_anni_madre)

Numero di Gravidanze

La maggior parte delle donne non ha mai avuto gravidanze. Tuttavia,il valore mediano risulta essere 1 gravidanza. Quindi, bisogna considerare che tra le donne che hanno avuto gravidanze, quelle più presenti nei 3 ospedali hanno avuta una sola gravidanza pregressa. Dunque, è meno probabile che ci siano in futuro donne con 2+ gravidanze. Questa ipotesi andrà saggiata eventualmente con dei test. È possibile notare:

  • Skewness positiva la quale indica che sono le modalità più basse quelle più frequenti.
  • Kurtosis positiva molto alta, indica che sono presenti molti outlier e alta variabilità. È necessario gestire gli outlier per poter utilizzare al meglio i modelli statistici rispettando le assunzioni di base
  • Ottenendo una Skewness positiva lontana da 0 e una Kurtosis positiva molto alta si può intuire che questa variabile non segue una distribuzione normale.

In conclusione per poter rappresentare i dati come una normale è stato usato il logaritmo poichè con il Q-Q Plot ottenuto dalla distribuzione normale troncata a 0, era evidente che i dati avessero un’andamento logaritmico. Dunque, è stato opportuno trasformare i dati in scala logaritmica. Il nuovo Q-Q Plot derivato dalla trasformazione logaritmica dei dati indica che la distribuzione dei dati è molto lontana dalla normale avendo delle code, rispetto alla linea, molto ampie.

n_gravidanze_cl <- cut(
  N.gravidanze,
  breaks = c(-Inf, 0, 1, 3, Inf),
  labels = c(
    "0 gravidanze",
    "1 gravidanza",
    "2-3 gravidanze",
    "4+ gravidanze"
  ),
  right = TRUE,
  include.lowest = TRUE
)

freq_n_gravidanze_cl <- tab_frequenze(n_gravidanze_cl)

create_table(freq_n_gravidanze_cl, FALSE, FALSE)
ni fi Ni Fi
0 gravidanze 1096 0.4384 1096 0.4384
1 gravidanza 818 0.3272 1914 0.7656
2-3 gravidanze 490 0.1960 2404 0.9616
4+ gravidanze 96 0.0384 2500 1.0000 
create_table(summary(N.gravidanze), FALSE, FALSE)
Statistic Value
1 Min. 0.0000
2 1st Qu. 0.0000
3 Median 1.0000
4 Mean 0.9812
5 3rd Qu. 1.0000
6 Max. 12.0000 
#shape indeces
n_gravidanze_shape_indeces <- evaluate_shape_indeces(N.gravidanze)
create_table(n_gravidanze_shape_indeces, FALSE, FALSE)
Statistic Value
1 Gini_Index 0.7347
2 Skewness 2.5143
3 Kurtosis 10.9894 
distr_n_gravidanze <- rnorm(n, mean(N.gravidanze), sd(N.gravidanze))

#Per spostare la distribuzione e ottenre valori positivi. Si somma 1 per evitare il valore 0 per il logaritmo poichè non è definito
distr_n_gravidanze_shifted <- var_shifted_to_positive(distr_n_gravidanze)

distr_n_gravidanze_log <- log(distr_n_gravidanze_shifted)



hist(distr_n_gravidanze_log,freq = F, breaks = 100)
lines(density(distr_n_gravidanze_log),col=2, lwd=3)

qqnorm(distr_n_gravidanze_log)
qqline(distr_n_gravidanze_log)

Peso

L’asimmetria negativa indica che sono presenti più modalità alte e la curtosi positiva indica un allungamento della curva. La distribuzione tra i valori ha eterogeneità massima. Il Q-Q Plot indica che la distribuzione per questa variabile non si discosta molto da quella di una normale. La distribuzione presenta una coda più spessa per modalità basse e una coda più stretta per modalità alte. Questo è coerente con la misura dell’assimetria negativa

peso_cl <- cut(
  Peso,
  breaks = c(-Inf, 900, 1500, 2499, 4499, Inf), 
  labels = c(
    "< 900 grammi [Estremamente basso]",
    "901-1500 grammi [Molto basso]",
    "1501–2499 grammi [Basso]",
    "2500-4499 grammi [Normale]",
    "4500+ grammi [Alto]"
  ),
  right = TRUE,
  include.lowest = TRUE
)




freq_peso_cl <- tab_frequenze(peso_cl)

create_table(freq_peso_cl, FALSE, FALSE)
ni fi Ni Fi
< 900 grammi [Estremamente basso] 2 0.0008 2 0.0008
901-1500 grammi [Molto basso] 22 0.0088 24 0.0096
1501–2499 grammi [Basso] 113 0.0452 137 0.0548
2500-4499 grammi [Normale] 2345 0.9380 2482 0.9928
4500+ grammi [Alto] 18 0.0072 2500 1.0000 
create_table(summary(Peso), FALSE, FALSE)
Statistic Value
1 Min. 830.000
2 1st Qu. 2990.000
3 Median 3300.000
4 Mean 3284.081
5 3rd Qu. 3620.000
6 Max. 4930.000 
# Shape indeces
peso_shape_indeces <- evaluate_shape_indeces(Peso)
create_table(peso_shape_indeces, FALSE, FALSE)
Statistic Value
1 Gini_Index 0.9962
2 Skewness -0.6470
3 Kurtosis 2.0315 
distr_peso <- rnorm(n, mean(Peso), sd(Peso))


hist(distr_peso,freq = F, breaks = 100)
lines(density(distr_peso),col=2, lwd=3)

qqnorm(distr_peso)
qqline(distr_peso)

Gestazione

La modalità più frequente è quella del termine precoce (37-39 settimane). Sarà necessario comprendere con opportuni test il perchè di questo dato. La gestazione nel Q-Q Plot ci indica che non si discosta molto da una normale teorica nonostante la kurtosi elevata e l’asimmetria negativa. Tuttavia, nell’istogramma è evidente come a sinistra ci sia una coda più stretta e a destra una coda più spessa. Approfondire con opportuni test statistici

gestazione_cl <- cut(
  Gestazione,
  breaks = c(-Inf, 28, 32, 37, 39, 40, 42, Inf),
  labels = c(
    "Estremamente prematuro (<28 settimane)",
    "Molto prematuro (28-32 settimane)",
    "Prematuro moderato (32-37 settimane)",
    "Termine precoce (37-39 settimane)",
    "Termine medio (39-40 settimane)",
    "Termine tardivo (40-42 settimane)",
    "Post-termine (>42 settimane)"
  ),
  right = TRUE,
  include.lowest = TRUE
)




freq_gestazione_cl <- tab_frequenze(gestazione_cl)

create_table(freq_gestazione_cl, FALSE, FALSE)
ni fi Ni Fi
Estremamente prematuro (<28 settimane) 8 0.0032 8 0.0032
Molto prematuro (28-32 settimane) 25 0.0100 33 0.0132
Prematuro moderato (32-37 settimane) 321 0.1284 354 0.1416
Termine precoce (37-39 settimane) 1018 0.4072 1372 0.5488
Termine medio (39-40 settimane) 741 0.2964 2113 0.8452
Termine tardivo (40-42 settimane) 385 0.1540 2498 0.9992
Post-termine (>42 settimane) 2 0.0008 2500 1.0000 
create_table(summary(Gestazione), FALSE, FALSE)
Statistic Value
1 Min. 25.0000
2 1st Qu. 38.0000
3 Median 39.0000
4 Mean 38.9804
5 3rd Qu. 40.0000
6 Max. 43.0000 
# Shape indeces
gestazione_shape_indeces <- evaluate_shape_indeces(Gestazione)
create_table(gestazione_shape_indeces, FALSE, FALSE)
Statistic Value
1 Gini_Index 0.8476
2 Skewness -2.0653
3 Kurtosis 8.2581 
distr_gestazione <- rnorm(n, mean(Gestazione), sd(Gestazione))


hist(distr_gestazione,freq = F, breaks = 100)
lines(density(distr_gestazione),col=2, lwd=3)

qqnorm(distr_gestazione)
qqline(distr_gestazione)

Confronto dei CV in ordine crescente

Le variabili non presentano un’alta variabilità. Tuttavia, le variabili con più variabilità risultano essere il peso e gli anni della madre. Questa caratteristica è da considerare durante i test statistici.

lista_cv <- sort(c(
  "Anni madre"  = CV(Anni.madre),
  "Settimane gestazione" = CV(Gestazione),
  "Peso"        = CV(Peso),
  "Lunghezza"   = CV(Lunghezza),
  "Cranio"      = CV(Cranio)
)
)


create_table(lista_cv, FALSE, FALSE)
Variable VIF
1 Settimane gestazione 4.7938
2 Cranio 4.8306
3 Lunghezza 5.3202
4 Peso 15.9874
5 Anni madre 18.7245

Relazioni tra le variabili

Numero di gravidanze e Età della madre

Dal BoxPlot è evidente che più aumentano gli anni della madre e più aumenta il numero di gravidanze. Ci sono molti valori outlier per le donne che hanno avuto 0 figli. Le donne che hanno 28 anni (valore mediano) hanno avuto in media 1-2 figli

boxplot(Anni.madre~N.gravidanze)

Distribuzione dell’età negli ospedali

In tutti e 3 gli ospedali la fascia di età più presente è quella tra i 25 e i 29 anni. Nella distribuzione delle fasce di età nei 3 ospedali non c’è una sostanziale differenza.

ggplot(data=dati_neonati)+
  geom_bar(aes(x=anni_madre_classi_cliniche,
               fill=Ospedale),
           position = "dodge",
           stat = "count",
           col = "black")+
  labs(title = "Distribuzione delle fasce di età",
       x="Età delle madri in classi, anni",
       y="Frequenze assolute")+
  theme_classic()+
  theme(legend.position = "bottom")

Distribuzione della Tipologia parto negli ospedali

Il numero di parti nei diversi ospedali è distribuito più o meno equamente tra gli ospedali. L’ospedale numero 2 risulta essere quello dove si sono verificati più parti cesarei. Il tipo di parto più frequente è quello naturale per tutti e 3 gli ospedali. Inoltre, i parti molto prematuri (nella fascia 28-32) risultano essere più frequenti per neonati di sesso femminile. Approfondire questa relazione con opportuni test statistici.

ggplot(data=dati_neonati)+
  geom_bar(aes(x=Tipo.parto,
               fill=Ospedale),
           position = "dodge",
           stat = "count",
           col = "black")+
  labs(title = "Distribuzione della tipologia di parto negli ospedali",
       x="Tipologia di parto",
       y="Frequenze assolute")+
  theme_classic()+
  theme(legend.position = "bottom")

Distribuzione del peso in relazione al sesso e al periodo di Gestazione

Dal boxplot emerge che, fino alla 35ª settimana di gestazione, i neonati di sesso femminile tendono ad avere un peso maggiore rispetto ai neonati di sesso maschile a parità di settimane di gestazione. A partire da questa settimana, però, il peso dei neonati di sesso maschile supera quello delle femmine, sempre a parità di settimane di gestazione. Inoltre, il peso dei neonati femminili appare più variabile rispetto a quello dei maschili, evidenziando una maggiore dispersione dei valori per le femmine rispetto ai maschi, sempre nello stesso intervallo di settimane di gestazione.

ggplot(dati_neonati, aes(x = factor(Gestazione), y = Peso, fill = Sesso)) +
  geom_boxplot() +
  labs(
    title = "Distribuzione del peso per sesso e settimane di gestazione",
    x = "Settimane di Gestazione",
    y = "Peso (grammi)",
    fill = "Sesso"
  ) +
  theme_minimal()

Numero di settimane di gestazione in relazione all’essere fumatrici o meno

Dal box plot e dal violin plot si può capire che l’essere fumatrice non incide sul numero di settimane di gestazione. Inoltre, bisogna considerare che per le donne non fumatrici ci sono molti valori outlier bassi. Per valutare questa tesi, saranno condotti opportuni test statistici.

boxplot(Gestazione ~ Fumatrici,
        data = dati_neonati,
        names = c("Non fumatrice", "Fumatrice"),
        xlab = "Fumatrice",
        ylab = "Settimane di gestazione",
        main = "Distribuzione settimane di gestazione per stato di fumo",
        col = c("lightblue", "salmon"))

ggplot(dati_neonati, aes(x = factor(Fumatrici), y = Gestazione, fill = factor(Fumatrici))) +
  geom_violin(trim = FALSE) +
  geom_boxplot(width = 0.1, fill = "white") +
  scale_x_discrete(labels = c("Non fumatrice", "Fumatrice")) +
  labs(x = "Fumatrice", y = "Settimane di gestazione",
       title = "Violin plot: Settimane di gestazione per stato di fumo") +
  theme_minimal()

Peso in relazione all’essere fumatrici

Il box plot e il violin plot indicano che il peso dei neonati tendenzialmente è più basso nelle donne fumatrici. Tuttavia, sono presenti molti outlier per le donne non fumatrici. Dunque, saranno necessari opportuni test statistici per approfondire questa relazione e valutare la tesi che le donne fumatrici partoriscono neonati con un peso inferiore.

boxplot(Peso ~ Fumatrici,
        data = dati_neonati,
        names = c("Non fumatrice", "Fumatrice"),
        xlab = "Fumatrice",
        ylab = "Peso",
        main = "Distribuzione Peso per stato di fumo",
        col = c("lightblue", "salmon"))

ggplot(dati_neonati, aes(x = factor(Fumatrici), y = Peso, fill = factor(Fumatrici))) +
  geom_violin(trim = FALSE) +
  geom_boxplot(width = 0.1, fill = "white") +
  scale_x_discrete(labels = c("Non fumatrice", "Fumatrice")) +
  labs(x = "Fumatrice", y = "Peso",
       title = "Violin plot: Peso per stato di fumo") +
  theme_minimal()

Indice ponderale(Peso/Lunghezza^3) in relazione al sesso

Da questi dati emerge che sembra non esserci una differenza sostanziale tra le misure antropometriche tra i due sessi. Questi ipotesi andrà saggiata con dei test

#grafico indice Ponderale
plot_variation_by_group(
  data = dati_neonati,
  x_var = "Gestazione",
  y_var = "Indice_Ponderale",
  group_var = "Sesso",
  lower_limit = 2.32,
  upper_limit = 2.85,
  lower_text = "2.32 - Normale",
  upper_text = "2.85 - Normale",
  y_min = 0,
  y_max = 5.5

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

ggplot(dati_neonati, aes(x = factor(Gestazione), y = Indice_Ponderale, fill = Sesso)) +
  geom_boxplot() +
  labs(
    title = "Distribuzione dell'indice ponderale per sesso e settimane di gestazione",
    x = "Settimane di Gestazione",
    y = "Indice Ponderale",
    fill = "Sesso"
  ) +
  coord_cartesian(ylim = c(0, 6)) + 
  theme_minimal()

distr_Indice_Ponderale <- rnorm(n, mean(dati_neonati$Indice_Ponderale), sd(dati_neonati$Indice_Ponderale))

ylim <- c(0, max(density(dati_neonati$Indice_Ponderale)$y) + 0.02)
hist(distr_Indice_Ponderale, freq = F, breaks = 100, ylim=ylim, main="Histogram of Indice Ponderale")
lines(density(dati_neonati$Indice_Ponderale), col=2, lwd=3)

qqnorm(distr_Indice_Ponderale)
qqline(distr_Indice_Ponderale)

Peso in relazione all’ospedale

Il peso non sembra cambiare sostanzialmente tra gli ospedali. Sono presenti molti outlier e l’ipotesi che il peso sia uguale tra gli ospedali dovrà essere saggiata con opportuni test.

boxplot(Peso ~ Ospedale, main="Peso ~ Ospedale")

2. Modello di Regressione lineare multipla

La variabile Sesso è stata resa una variabile dummy dove M=0 e F=1. La variabile dipendente Peso è stata trasformata, ottenendo la radice quadrata del peso. Dunque, quando sono state effettuate le previsioni si è tenuto conto di questo elevando al quadrato la previsione del peso. Per scegliere il modello, partendo da quello con tutte le variabili del dataset, si sono eliminate una ad una le variabili meno significative del modello. L’accettazione dell’ipotesi di significatività di una variabile nel modello è data dal p-value calcolato nel modello lineare e dalla diminuizione del bic. L’età della madre è stata inizialmente inclusa come variabile predittiva del peso del neonato, ma: Il p-value = 0.29 indica assenza di significatività statistica. La rimozione ha migliorato il BIC di 7 punti, suggerendo un modello più efficiente. In Conclusione l’età materna è stata esclusa dal modello finale per favorire semplicità e performance.

Si è scelto il modello 4 poichè risulta essere quello con il BIC più basso e quello più efficiente per spiegare la variabilità nei dati. Di seguito, sono state confrontate e analizzate tutte le statistiche adeguate per poter scegliere il miglior modello.

Matrice di correlazione

Le variabili maggiormente correlate risultano essere:

  • Peso ~ Cranio
  • Gestazione ~ Peso
  • Gestazione ~ Lunghezza
  • Peso ~ Lunghezza
  • Lunghezza ~ Cranio
# --- Step 1: Costruisci modelli iniziali
modelli_iniziali <- crea_modelli_ridotti(dati_neonati)

# --- Step 2: Pulizia outlier su mod1
dati_puliti <- rimuovi_outliers_iterativi(
  modello = modelli_iniziali[["mod1"]],
  data = dati_neonati,
  formula_modello = sqrt(Peso) ~ .
)

Rimozione Outlier

È stato rimosso circa il 7% delle osservazioni rilevate come outlier statisticamente significativi. Gli outlier più influenti statisticamente sono stati identificati attraverso un valore soglia scelto empiricamente con 4/n e da cui è stata calcolata la distanza di Cook. Vengono eseguite due iterazioni perchè nella prima vengono identificati gli outlier rispetto al modello originale, mentre nella seconda, dopo aver rimosso gli outlier, il modello viene ricalcolato e potrebbe avere nuovi outlier che prima non erano rilevabili. Infatti, è stato scelto di rimuovere ulteriori 4 outlier singificativi per migliorare il modello.

Rimosse 53 osservazioni (prima iterazione)

Rimosse 131 osservazioni (seconda iterazione)

Totale osservazioni rimosse: 184 

# --- Step 3: Matrice di correlazione
disegna_matrice_correlazione(dati_puliti)

# --- Step 4: Confronto varianze
confronta_varianze(dati_puliti, "Sesso", "Fumatrici")

Test confronto varianza: Peso ~ Sesso

La varianza tra il peso e il sesso risulta essere uguale. Questo è stato considerato nei test statistici

F test to compare two variances

data: Peso by Sesso F = 1.0078, num df = 1160, denom df = 1154, p-value = 0.8943 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.8980836 1.1310049 sample estimates: ratio of variances 1.007847

Statistic Value
1 F 1.0078
2 p.value 0.8943
3 num df 1160.0000
4 denom df 1154.0000
5 ratio of variances 1.0078
6 conf.int_lower 0.8981
7 conf.int_upper 1.1310

Test confronto varianza: Peso ~ Fumatrici

La varianza tra il peso e l’essere fumatrici risulta essere uguale. Questo è stato considerato nei test statistici

F test to compare two variances

data: Peso by Fumatrici F = 1.1598, num df = 2221, denom df = 93, p-value = 0.3568 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.8463379 1.5264613 sample estimates: ratio of variances 1.159774

Statistic Value
1 F 1.1598
2 p.value 0.3568
3 num df 2221.0000
4 denom df 93.0000
5 ratio of variances 1.1598
6 conf.int_lower 0.8463
7 conf.int_upper 1.5265 
# --- Step 5: Ricostruisci i modelli sui dati puliti
modelli_finali <- crea_modelli_ridotti(dati_puliti)

# --- Step 6: Confronti BIC/ANOVA
confronta_modelli(modelli_finali)

ANOVA confronto tra gli ultimi due modelli

I valori del test statistico evidenziano come il mod4 abbia una capacità più elevata di spiegare la variabilità nei dati. Dunque, il parametro che è stato rimosso dal modello 3 al modello 4 non è statisticamente significativo. Questo rende il modello 4 più semplice e statisticamente rilevante.

Res.Df RSS Df Sum of Sq F Pr(>F)
1 2308 5048.828 NA NA NA NA
2 2309 5051.365 -1 -2.5373 1.1599 0.2816

Confronto tra modelli: BIC, MSE, RMSE

L’MSE e l’RMSE non cambia tra tutti i modelli. Questo indica che semplificando il modello la variabilità nei dati non peggiora significativamente. Pertanto, è adatto scegliere l’ultimo modello che risulta essere meno complesso. Inoltre, Il BIC più basso di 7 punti risulta essere significativo poichè indica che la qualità dell’adattamento ai dati è migliore.

# --- Step 7: Diagnostica finale
diagnostica_modello(modelli_finali[["mod4"]])

VIF

I VIF del mod4 sono inferiori a 5 per tutti i regressori indica che non ci sono problemi di multicollinearità, ovvero di correlazione tra i regressori

## VIFs computed for predictors

Analisi dei residui mod4

Dai seguenti plot dei residui è emerso che:

  • Residuals vs Fitted: Questo grafico mostra i residui rispetto ai valori predetti. I residui si distribuiscono casualmente attorno alla linea orrizontale zero, senza pattern evidenti. La variabilità risulta essere costante
  • Q-Q Residuals: Questo grafico verifica la normalità dei residui. I residui seguono la linea tratteggiata obliqua, anche se le code sono leggermente più pesanti
  • Scale Location: Questo grafico mostra la radice quadrata dei residui standardizzati rispetto ai valori predetti, ed è utile per verificare l’omoschedasticità. Nel grafico risulta una leggere eteroschedasticità. Tuttavia, con il BP test è stata accettata l’ipotesi di omoschedasticità
  • Residuals vs Leverage: Questo grafico aiuta a identificare punti influenti nel modello poichè i punti che si trovano lontano dal gruppo principale e hanno un alto leverage potrebbero influenzare significativamente i risultati del modello. Si nota come ci siano dei punti outlier che potrebbero essere significativi e da poter valutare e rimuovere se necessario

Breusch-Pagan Test

Un p-value molto alto indica che si può accettare l’ipotesi nulla di omoschedasticità

studentized Breusch-Pagan test

data: modello BP = 9.8274, df = 6, p-value = 0.1321

Statistic Value
1 BP 9.8274
BP p.value 0.1321
11 df 6.0000

Durbin-Watson Test

il p-value molto alto indica che si può accettare l’ipotesi nulla di non correlazione tra i regressori. Durbin-Watson test

data: modello DW = 2.0223, p-value = 0.7043 alternative hypothesis: true autocorrelation is greater than 0

Statistic Value
1 DW 2.0223
2 p.value 0.7043

Outlier Test

Valutazione del residuo standardizzato usato per identificare i possibili outliers significativi in un modello di regressione. In questo caso il Largest r-student e l’unadjusted p-value suggeriscono che ci sia almeno un outlier significativo. Tuttavia, con il criterio Bonferroni(tiene conto del numero di test effettuati) e il relativo calcolo del p-value,che in questo caso si attesta all’incirca su 0.07, si potrebbe escludere l’ipotesi che ci sia almeno un outlier significativo. Questo punto non incide significativamente sulle prestazioni del modello.

Outlier Test Results
Row Rstudent P_value Bonferroni_P_value Significant Cutoff
1395 1395 4.163168 3.25e-05 0.0753583 FALSE 0.05

Distanza di Cook, Residui e Leverage

I residui si distribuiscono normalmente e risultano avere varianza costante. Non sono presenti punti di leverage significativi. Il punto con massima distanza di cook ha distanza molto inferiore a 0.5

Max Cook’s Distance at index 2244 : 0.0098

# --- Step 8: Summary finale
cat("\n\n<b>Summary mod4</b>\n\n")

Summary mod4

cat("\nIl modello 4 risulta avere una capacità di spiegare la variabilità nei dati di 0.86, ovvero corrisponde al valore dell'adjusted R-Squared. Tutte le variabili del modello risultano essere siginificative. Anche se i regressori Fumatrici e N.gravidanze risultano avere un p-value vicino a 0.05, si è scelto di mantenerle nel modello perchè potrebbero essere utili ai fine della previsione del peso. In conclusione, il p-value vicino a 0 indica che le variabili sono statisticamente significative nel predire la variabile dipendente")

Il modello 4 risulta avere una capacità di spiegare la variabilità nei dati di 0.86, ovvero corrisponde al valore dell’adjusted R-Squared. Tutte le variabili del modello risultano essere siginificative. Anche se i regressori Fumatrici e N.gravidanze risultano avere un p-value vicino a 0.05, si è scelto di mantenerle nel modello perchè potrebbero essere utili ai fine della previsione del peso. In conclusione, il p-value vicino a 0 indica che le variabili sono statisticamente significative nel predire la variabile dipendente

create_table(summary(modelli_finali[["mod4"]]),F,F)
Term Estimate Std_Error t_value p_value
1 (Intercept) -12.595889779907 0.8776 -14.3533 0.0000
2 N.gravidanze -0.051592409271825 0.0259 -1.9887 0.0469
3 Fumatrici -0.307155415736831 0.1565 -1.9623 0.0499
4 Gestazione 0.287015071056442 0.0230 12.4814 0.0000
5 Sesso -0.378238192074649 0.0631 -5.9917 0.0000
6 Indice_Ponderale 6.93676561968473 0.1392 49.8326 0.0000
7 Lunghezza:Cranio 0.000237931283916777 0.0000 84.2637 0.0000
8 Residual Std. Error 1.47908259237042 NA NA NA
9 Multiple R-squared 0.858807802497925 NA NA NA
10 Adjusted R-squared 0.858440910689777 NA NA NA
11 F-statistic 2340.76581549836 NA NA NA
12 Model DF 6 and 2309 NA NA NA
13 p-value 0 NA NA NA

2.1 Test statistici senza outlier

In questa fase saranno condotti test per saggiare le seguenti ipotesi:

È risultato statisticamente rilevante rimuovere all’incirca il 7% dei valori outlier poichè questi alteravano gli esiti dei test. Infatti, è stato notato che nel test per verificare l’incidenza dell’essere fumatrici sul peso, in caso di presenza di outlier il test induceva ad accettare l’ipotesi nulla di media peso fumatrici - media peso non fumatrici = 0. Al contrario, in caso di assenza degli outlier, il test ha portato al rifiuto dell’ipotesi nulla. Dunque, la presenza di outlier ha influenzato questo test.

Test per saggiare l’ipotesi che in alcuni ospedali si fanno più parti cesarei

Per saggiare questa ipotesi è corretto usare il test del chi quadrato sulla tabella di contingenza, tra Ospedale e tipo_parto, perchè il sistema di ipotesi è costruito nel seguente modo:

  • H0: Tipo di parto indipendente dall’ospedale;
  • H1: Tipo di parto dipende dall’ospedale

Dal test è emerso che:

  • Valore del test chi quadro pari a 1. Essendo basso indica che le frequenze osservate sono molto vicine a quelle attese. Ovvero, non c’è nessuna grande differenza tra i gruppi
  • Valore del p-value molto più grande del valore di alpha=0.05 che ci consente di accettare l’ipotesi nulla molto facilmente
  • Con l’indice di Kramer, si è valutata l’associazione tra le due variabili. Risultando molto basso conferma ulteriormente che non vi è alcuna associazione tra le due variabili. Dunque, le due variabili sono indipendenti
tabella_contingenza <- table(dati_puliti$Ospedale, dati_puliti$Tipo.parto)

# Visualizzare la tabella di contingenza
tabella_risultante <- crea_tabella_contingenza_stilizzata(
   dati_puliti$Ospedale, 
   dati_puliti$Tipo.parto, 
   caption = "Distribuzione dei parti per ospedale",
   gruppo_colonne =  "Tipo di parto",
)
tabella_risultante
Distribuzione dei parti per ospedale
Tipo di parto
Ces Nat Totale
osp1 223 527 750
osp2 234 557 791
osp3 216 559 775
Totale 673 1643 2316 
# Eseguire il test del chi-quadro
test_chi_quadro <- chisq.test(tabella_contingenza)
#per valutare la forza dell'associazione si può usare l'indice di Kramer
k <- min(nrow(tabella_contingenza), ncol(tabella_contingenza))
kramer_index <- (sqrt(test_chi_quadro$statistic))/(n*(k-1))

test_chi_quadro_values <- c("Test Value" = test_chi_quadro$statistic, "p-value" = test_chi_quadro$p.value, "Degrees of Freedom" = test_chi_quadro$parameter, "Kramer Index" = kramer_index)
test_chi_quadro_df <- data.frame(
  Statistic = names(test_chi_quadro_values),
  Value = as.numeric(test_chi_quadro_values)
)


# Visualizzare i risultati del test
cat("<b>chi-quadro-test Summary</b><br>")

chi-quadro-test Summary

create_table(test_chi_quadro_df, F, F)
Statistic Value
1 Test Value.X-squared 0.8013
2 p-value 0.6699
3 Degrees of Freedom.df 2.0000
4 Kramer Index.X-squared 0.0004 
# Interpretazione dei risultati
if (test_chi_quadro$p.value < 0.05) {
  cat("C'è una differenza significativa nella proporzione di parti cesarei tra gli ospedali.\n")
} else {
  cat("Non c'è una differenza significativa nella proporzione di parti cesarei tra gli ospedali.\n")
}

Non c’è una differenza significativa nella proporzione di parti cesarei tra gli ospedali.

Test per saggiare l’ipotesi che la media del peso è uguale a quello della popolazione

Questo test ha l’obiettivo di verificare l’ipotesi che la media del peso della popolazione di neonati, uguale a 3300 grammi, sia uguale alla media del peso. Dunque, è corretto utilizzare lo z-test dove il sistema di ipotesi corrisponde a:

  • H0: mu = 3300
  • H1: mu ≠ 3300

Dal test sono emerse le seguenti considerazioni:

  • Valore di z positivo: indica che la media del campione è superiore alla media ipotizzata dalla popolazione ed è 0.68 deviazioni standard al di sotto della media ipotizzata. Dunque, la differenza tra la media del campione e la media ipotizzata è pari a 0.68 volte la deviazione standard della distribuzione campionaria.
  • p-value più alto di alpha indica che siamo sicuri di accettare l’ipotesi nulla, ovvero l’ipotesi di uguaglianza tra la media della popolazione e quella del campione
population_weight_mean <- 3300
alpha_weight_test_value <- 0.95

z_test_peso <- z.test(x = dati_puliti$Peso,
                      mu = population_weight_mean,
                      sigma.x = sd(dati_puliti$Peso),
                      alternative = "two.sided",
                      conf.level = alpha_weight_test_value)

weight_mean_estimate_and_alpha <- c("Alpha" = alpha_weight_test_value, "Weight Mean" = mean(dati_puliti$Peso))

cat("<b>General Results</b><br>")

General Results

sink(tempfile())
create_table(weight_mean_estimate_and_alpha, F, F)
Variable VIF
1 Alpha 0.95
2 Weight Mean 3305.83 
sink()

cat("<b>z-test Summary</b><br>")

z-test Summary

sink(tempfile())
create_table(z_test_peso, F, F)
Statistic Value
1 z 0.6284
2 p.value 0.5298
3 mean of x 3305.8299
4 conf.int_lower 3287.6457
5 conf.int_upper 3324.0141 
sink()

Test per saggiare l’ipotesi che la media della lunghezza è uguale a quella della popolazione

Il sistema di ipotesi per questo test è la seguente:

  • H0: mu = 500
  • H1: mu ≠ 500

È stato usato un valore di alpha molto piccolo, ovvero 0.005 per coda. La media della popolazione non si trova all’interno dell’intervallo di confidenza.

Dal test sono emerse le seguenti considerazioni:

  • Valore di Z negativo: indica che la media del campione è di molto inferiore alla media della popolazione
  • p-value: molto vicino a 0 ed indica che siamo confidenti nel rifiutare l’ipotesi nulla. Dunque, la media della lunghezza di questo campione non è uguale a quella della popolazione anche se la media si trova all’interno dell’intervallo di confidenza
population_length_mean <- 500
alpha_length_test_value <- 0.99

z_test_lunghezza <- z.test(x = dati_puliti$Lunghezza,
                      mu = population_length_mean,
                      sigma.x = sd(dati_puliti$Lunghezza),
                      alternative = "two.sided",
                      conf.level = alpha_length_test_value)
length_mean_estimate_and_alpha <- c("Alpha" = alpha_length_test_value, "Weight Mean" = mean(dati_puliti$Lunghezza))

cat("<b>General Results</b><br>")

General Results

sink(tempfile())
create_table(length_mean_estimate_and_alpha, F, F)
Variable VIF
1 Alpha 0.9900
2 Weight Mean 496.3493 
sink()

cat("<b>t-test Summary</b><br>")

t-test Summary

sink(tempfile())
create_table(z_test_lunghezza, F, F)
Statistic Value
1 z -8.4666
2 p.value 0.0000
3 mean of x 496.3493
4 conf.int_lower 495.2386
5 conf.int_upper 497.4600 
sink()

Test per saggiare l’ipotesi che le misure antropometriche tra i due sessi sono diverse

Test per saggiare l’ipotesi che le misure antropometriche differiscono tra i due sessi

Per confrontare le misure antropometriche tra maschi e femmine, è stato eseguito un t-test per campioni indipendenti. La numerosità del campione (n > 30 per ciascun gruppo) consente di applicare il test in modo robusto grazie al teorema del limite centrale, anche senza una perfetta normalità. Il test di Shapiro-Wilk non è stato considerato determinante, in quanto la sua elevata sensibilità tende a far rifiutare l’ipotesi di normalità anche in presenza di lievi deviazioni in grandi campioni (>500).

Sono state formulate le seguenti ipotesi:

  • H0: le medie delle misure antropometriche (peso e indice ponderale) sono uguali tra i due sessi
  • H1: le medie delle misure antropometriche differiscono tra i due sessi

Sono stati eseguiti sia t-test che test di Wilcoxon per verificare la coerenza tra metodi parametrici e non parametrici.

Risultati:

  • Per tutte le variabili (peso,indice ponderale, lunghezza e cranio), i p-value ottenuti sono molto bassi (<< 0.05), portando al rifiuto dell’ipotesi nulla in tutti i test.
  • Gli intervalli di confidenza per la differenza tra le medie non includono lo zero, rafforzando l’evidenza di una differenza significativa.
  • I valori negativi del t indicano che la media nel gruppo femminile è inferiore rispetto al gruppo maschile.
  • Anche i test non parametrici di Wilcoxon confermano i risultati, mostrando coerenza tra le metodologie.

Conclusione:

I test evidenziano che esistono differenze significative tra i sessi per quanto riguarda sia il peso sia l’indice ponderale. Pertanto, le misure antropometriche dipendono dal sesso.

# t-test
cat("<b>Peso - Sesso Two Samples t-test</b><br>")

Peso - Sesso Two Samples t-test

sink(tempfile())
create_table(t.test(dati_puliti$Peso ~ dati_puliti$Sesso, var.equal = T), F,F)
Statistic Value
1 t -12.6972
2 p.value 0.0000
3 df 2314.0000
4 mean in group F 3192.1981
5 mean in group M 3420.0519
6 conf.int_lower -263.0441
7 conf.int_upper -192.6636 
sink()


cat("<b>Indice Ponderale - Sesso Two Samples t-test</b><br>")

Indice Ponderale - Sesso Two Samples t-test

sink(tempfile())
create_table(t.test(dati_puliti$Indice_Ponderale ~ dati_puliti$Sesso, var.equal = TRUE), F,F)
Statistic Value
1 t -4.7905
2 p.value 0.0000
3 df 2314.0000
4 mean in group F 2.6721
5 mean in group M 2.7165
6 conf.int_lower -0.0625
7 conf.int_upper -0.0262 
sink()

cat("<b>Lunghezza - Sesso Two Samples t-test</b><br>")

Lunghezza - Sesso Two Samples t-test

sink(tempfile())
create_table(t.test(dati_puliti$Lunghezza ~ dati_puliti$Sesso, var.equal = T), F,F)
Statistic Value
1 t -10.3268
2 p.value 0.0000
3 df 2314.0000
4 mean in group F 492.0060
5 mean in group M 500.7152
6 conf.int_lower -10.3629
7 conf.int_upper -7.0553 
sink()


cat("<b>Cranio - Sesso Two Samples t-test</b><br>")

Cranio - Sesso Two Samples t-test

sink(tempfile())
create_table(t.test(dati_puliti$Cranio ~ dati_puliti$Sesso, var.equal = TRUE), F,F)
Statistic Value
1 t -7.0990
2 p.value 0.0000
3 df 2314.0000
4 mean in group F 338.8053
5 mean in group M 342.9506
6 conf.int_lower -5.2904
7 conf.int_upper -3.0002 
sink()

# Wilcox test
cat("<b>Indice Ponderale - Sesso Wilcox Test</b><br>")

Indice Ponderale - Sesso Wilcox Test

sink(tempfile())
create_table(wilcox.test(dati_puliti$Indice_Ponderale ~ dati_puliti$Sesso),F,F)
Statistic Value
1 W 589400
2 p.value 0 
sink()

cat("<b>Peso - Sesso Wilcox Test</b><br>")

Peso - Sesso Wilcox Test

sink(tempfile())
create_table(wilcox.test(dati_puliti$Peso ~ dati_puliti$Sesso),F,F)
Statistic Value
1 W 459573.5
2 p.value 0.0 
sink()

cat("<b>Lunghezza - Sesso Wilcox Test</b><br>")

Lunghezza - Sesso Wilcox Test

sink(tempfile())
create_table(wilcox.test(dati_puliti$Lunghezza ~ dati_puliti$Sesso),F,F)
Statistic Value
1 W 503296
2 p.value 0 
sink()

cat("<b>Cranio - Sesso Wilcox Test</b><br>")

Cranio - Sesso Wilcox Test

sink(tempfile())
create_table(wilcox.test(dati_puliti$Cranio ~ dati_puliti$Sesso),F,F)
Statistic Value
1 W 554632
2 p.value 0 
sink()

Test per saggiare l’ipotesi che il peso tra le donne fumatrici e non fumatrici sia uguale

Il sistema di ipotesi è costruito nel seguente modo:

  • H0: mu_peso_fumatrici - mu_peso_non_fumatrici = 0;
  • H1: mu_peso_fumatrici - mu_peso_non_fumatrici ≠ 0;

Dal test è emerso che possiamo rifiutare H0 poichè è stato ottenuto un p-value vicino a 0. Questo indica che l’essere fumatrici influenza il peso del neonato

cat("<b>Peso - Fumatrici Two Samples t-test</b><br>")

Peso - Fumatrici Two Samples t-test

sink(tempfile())
create_table(t.test(dati_puliti$Peso ~ dati_puliti$Fumatrici, var.equal = TRUE), F,F)
Statistic Value
1 t 2.8600
2 p.value 0.0043
3 df 2314.0000
4 mean in group 0 3311.2790
5 mean in group 1 3177.0213
6 conf.int_lower 42.2017
7 conf.int_upper 226.3138 
sink()

Test per saggiare l’ipotesi che la media del peso nei vari ospedali sia diverso

L’anova test per il confronto tra medie è quello più appropriato perchè il peso è una variabile quantitativa e si vuole confrontare questa variabile con la variabila categorica legata all’ospedale che possiede più di due valori possibili. Se i valori possibili fossero stati due si sarebbe potuto usare un t-test. Questo test è utile per verificare se sono presenti ospedali dove la media del peso di un neonato risulta essere anomala, creando cosi una difficoltà nella gestione dei parti a rischio. L’anova test consente di confrontare la variabilità interna tra i gruppi e la variabilità tra i gruppi. Il valore della statistica è dato dal rapporto della variabilità tra gruppi e quella della variabilità interna dei gruppi.

Il sistema di ipotesi è costruito nel seguente modo:

  • H0: Le medie del peso sono uguali tra gli ospedali
  • H1: Almeno una media è diversa

Dal test è emerso un valore della statistica F vicino ad un 1 ed un p-value maggiore di 0.05. Questo indica che accettiamo l’ipotesi nulla H0 affermando che le medie del peso sono uguali tra gli ospedali. Questo indica che negli ospedali non vi è una differenza sostanziale dei valori anomali nel peso dei neonati.

cat("<b>Anova Test Summary - Media Peso negli Ospedali</b><br>")

Anova Test Summary - Media Peso negli Ospedali

anova_test <- aov(dati_puliti$Peso ~ dati_puliti$Ospedale, data = dati_puliti)
create_table(anova_test, F,F)
Statistic Value
1 F value 1.1026
2 p-value 0.3322
3 num Df 2.0000
4 denom Df 2313.0000

Previsioni

Vengono generati dei valori randomici in base al minimo e al massimo dei regressori per poter effettuare delle previsioni. Il risultato delle previsioni viene elevato al quadrato perchè durante la costruzione del modello di regressione è stata eseguita la radice quadrata del peso.

data_mod <- dati_puliti
data_mod$Sesso <- ifelse(data_mod$Sesso == "F", 1, 0)

nuovi_dati_casuali <- generate_random_data(data_mod, n=2)

previsioni <- predict(modelli_finali[["mod4"]], newdata = nuovi_dati_casuali)

nuovi_dati_casuali$predict_value <- previsioni^2

create_table(nuovi_dati_casuali, F, F)
Gestazione Lunghezza Cranio Sesso N.gravidanze Fumatrici Indice_Ponderale predict_value
1 38 413.3945 357.3331 0 12 0 2.9605 2848.901
2 40 462.2450 345.8696 1 10 0 2.4391 2803.647

Visualizzazioni grafiche delle relazioni più significative

Il fumo risulta avere un’incidenza negativa sul peso del neonato anche a parità di settimane di gestazione e di tutte le altre variabili. Tranne per la lunghezza, dove non vi è una grande differenza di peso dei neonati tra donne fumatrici e non fumatrici. La retta nello scatterplot del numero gravidanze in relazione alle settimane di gestazione è quasi orizzontale. Questo indica che la variabile X non è un buon predittore di Y e che probabilmente non c’è una relazione lineare significativa tra le due variabili. Tuttavia, si è deciso di includere questa variabile nelle predizioni perchè l’avere gravidanze potrebbe incidere sul peso del bambino poichè le gravidanze precedenti potrebbero aver influenzato la salute della madre. Invece se una donna non ha mai avuto gravidanze i fattori genetici e nutrizionali potrebbero incidere sul peso del bambino

ggplot(dati_puliti, aes(x = as.factor(Fumatrici), y = Peso)) +
  geom_boxplot(fill = c("#66c2a5", "#fc8d62")) +
  labs(x = "Fumo materno (0 = No, 1 = Sì)", y = "Peso del neonato (g)",
       title = "Effetto del fumo sul peso del neonato") +
  theme_minimal()

ggplot(dati_puliti, aes(x = Gestazione, y = Peso, color = as.factor(Fumatrici))) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Settimane di gestazione", y = "Peso del neonato (g)",
       color = "Fumatrici",
       title = "Peso in funzione delle settimane di gestazione") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

ggplot(dati_puliti, aes(x = N.gravidanze, y = Peso)) +
  geom_point(alpha = 0.6, color = "#8da0cb") +
  geom_smooth(method = "lm", se = TRUE, color = "black") +
  labs(x = "N.gravidanze", y = "Peso del neonato (g)",
       title = "Relazione tra N.gravidanze e Peso") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

ggplot(dati_puliti, aes(x = Sesso, y = Peso, fill = as.factor(Fumatrici))) +
  geom_violin(trim = FALSE, alpha = 0.6) +
  labs(x = "Sesso", y = "Peso del neonato (g)", fill = "Fumatrici",
       title = "Relazione tra Sesso e Peso (con distinzione Fumatrici)") +
  theme_minimal()

ggplot(dati_puliti, aes(x = Lunghezza, y = Peso, color = as.factor(Fumatrici))) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE, aes(color = as.factor(Fumatrici))) +
  labs(x = "Lunghezza", y = "Peso del neonato (g)", color = "Fumatrici",
       title = "Relazione tra Lunghezza e Peso (per Fumatrici)") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

ggplot(dati_puliti, aes(x = Cranio, y = Peso, color = as.factor(Fumatrici))) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE, aes(color = as.factor(Fumatrici))) +
  labs(x = "Cranio", y = "Peso del neonato (g)", color = "Fumatrici",
       title = "Relazione tra Cranio e Peso (per Fumatrici)") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'