INTRODUZIONE

Il presente report analizza il Boston Housing Dataset, un insieme di dati che raccoglie informazioni sulle caratteristiche strutturali, ambientali e socio-economiche delle abitazioni, nonché sul prezzo medio degli immobili in diverse aree della città di Boston. L’obiettivo principale dell’analisi è individuare i fattori che influenzano il valore delle case e sviluppare un modello predittivo affidabile in grado di stimarne il prezzo.

Il dataset è stato originariamente introdotto da Harrison e Rubinfeld (1978) ed è stato successivamente utilizzato da Belsley, Kuh e Welsch (1980) nel contesto dello studio dei modelli di regressione. Per il presente lavoro, i dati sono stati reperiti dalla piattaforma Kaggle (link: https://www.kaggle.com/datasets/mdraselsarker/boston-housing-dataset?resource=download).

L’analisi si articola in diverse fasi. In primo luogo, viene condotta un’esplorazione preliminare dei dati (Exploratory Data Analysis, EDA) al fine di comprendere la distribuzione delle variabili e le relazioni esistenti. Successivamente, si applica un’analisi delle componenti principali (PCA) per ridurre la dimensionalità del dataset, seguita da un’analisi dei cluster volta a individuare gruppi omogenei di aree urbane. Infine, viene implementato un modello di regressione lineare multipla per stimare il prezzo degli immobili e valutarne le performance predittive.

Librerie utilizzate:

library(ggplot2)
library(corrplot)
library(dplyr)
library(cluster)
library(factoextra) 
library(car)
library(reshape2)
library(sandwich)
library(lmtest)

ANALISI ESPLORATIVA (EDA)

L’analisi ha inizio con una Exploratory Data Analysis (EDA), finalizzata a una prima esplorazione dei dati. In questa fase preliminare è stata innanzitutto esaminata la struttura del dataset, verificando la tipologia delle variabili.

setwd("C:/Users/Sara/Desktop/Report")
df = read.csv("Boston.csv")

dim(df)
## [1] 506  14
names(df)
##  [1] "CRIM"    "ZN"      "INDUS"   "CHAS"    "NOX"     "RM"      "AGE"    
##  [8] "DIS"     "RAD"     "TAX"     "PTRATIO" "B"       "LSTAT"   "Price"

Il dataset è composto da 506 osservazioni e 14 variabili, di cui una è la variabile target, rappresentata da Price che indica il prezzo medio delle abitazioni occupate dai proprietari, espresso in migliaia di dollari. Le altre variabili sono:

sapply(df, class)
##      CRIM        ZN     INDUS      CHAS       NOX        RM       AGE       DIS 
## "numeric" "numeric" "numeric" "integer" "numeric" "numeric" "numeric" "numeric" 
##       RAD       TAX   PTRATIO         B     LSTAT     Price 
## "integer" "integer" "numeric" "numeric" "numeric" "numeric"

Il dataset è composto da 11 variabili numeriche continue e 3 variabili discrete (TAX, RAD), di cui una dicotomica (CHAS).

colSums((is.na(df))) #verifico se ci sono valori nulli
##    CRIM      ZN   INDUS    CHAS     NOX      RM     AGE     DIS     RAD     TAX 
##       0       0       0       0       0       0       0       0       0       0 
## PTRATIO       B   LSTAT   Price 
##       0       0       0       0
sum(duplicated((df))) #verifico se ci sono valori duplicati
## [1] 0

Contestualmente, è stata verificata la presenza di valori mancanti o osservazioni duplicate, che non sono risultati presenti nel dataset.

Successivamente, attraverso il calcolo delle principali statistiche descrittive, ottenute mediante la funzione summary(), è stato possibile analizzare le caratteristiche principali del dataset.

summary(df)
##       CRIM                ZN             INDUS            CHAS        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
##  1st Qu.: 0.08205   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
##       NOX               RM             AGE              DIS        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       RAD              TAX           PTRATIO            B         
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      LSTAT           Price      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00

Dall’analisi descrittiva emerge che la variabile CHAS presenta una distribuzione fortemente sbilanciata: la maggior parte delle osservazioni corrisponde ad abitazioni che non confinano con il fiume (0), mentre solo una piccola parte ha valore 1. Prima della visualizzazione, la variabile è stata trasformata in un fattore dicotomico con due livelli espliciti:

df$CHAS <- factor(df$CHAS, levels = c(0,1), labels = c("No", "Si"))

ggplot(df, aes(x=CHAS)) +
  geom_bar(fill = "skyblue", color = "black") +
  labs(title = "Distribuzione della variabile CHAS",
       x = "Confina con il fiume Charles",
       y = "Numero di abitazioni") +
  theme_minimal()

Il grafico a barre conferma la distribuzione sbilanciata, in cui la maggioranza delle abitazioni appartiene alla categoria “No”, mentre solo una piccola quota si affaccia sul fiume.

par(mfrow = c(1, 4))
boxplot(df$CRIM, main='CRIM',col='Sky Blue')
boxplot(df$ZN, main='ZN',col='Sky Blue')
boxplot(df$RM, main='RM',col='Sky Blue')
boxplot(df$B, main='B',col='Sky Blue')

L’analisi grafica mediane boxplot evidenzia caratteristiche distintive anche nelle altre principali variabili. La variabile CRIM presenta una distribuzione fortemente asimmetrica a destra, con una mediana bassa (0.25) e numerosi valori estremi (outlier) che arrivano quasi fino ad 89. Ciò indica che, pur essendoci molte zone con un basso tasso di criminalità, esistono alcune aree con valori eccezionalmente elevati.

La variabile ZN mostra invece una concentrazione la maggior parte dei valori su 0 (nel summary sia la mediana che il primo quartile sono uguali a 0), con pochi casi isolati di valori elevati, rappresentati come outlier nel boxplot. Ciò suggerisce che la gran parte delle abitazioni non si trova in zone con ampie aree residenziali.

Per quanto riguarda RM, la distribuzione risulta più simmetrica e compatta, con mediana e media intorno a 6.2 e pochi outlier. La concentrazione dei dati attorno a un valore centrale indica una relativa omogeneità nella dimensione degli alloggi.

Infine, la variabile B presenta valori medi elevati (~356) ma con ampia variabilità, come evidenziato da numerosi outlier verso il basso.

par(mfrow = c(2, 4), mar = c(4,4,2,1))

hist(df$INDUS, main = "INDUS", col = "lightblue", breaks = 20, border = "black")
hist(df$TAX, main = "TAX", col = "lightgreen", breaks = 20, border = "black")
hist(df$PTRATIO, main = "PTRATIO", col = "pink",breaks = 20, border = "black")
hist(df$NOX, main = "NOX", col = "lightyellow", breaks = 20, border = "black")
hist(df$AGE, main = "AGE", col = "lightgray", breaks = 20, border = "black")
hist(df$DIS, main = "DIS", col = "orange",breaks = 20, border = "black")
hist(df$RAD, main = "RAD", col = "lightcoral",breaks = 20, border = "black")
hist(df$LSTAT, main = "LSTAT", col = "lightcyan", breaks = 20, border = "black")

Successivamente ho esaminato la distribuzione delle altre variabili del dataset tramite istogrammi. I risultati evidenziano i seguenti aspetti:

df %>% 
  ggplot(aes(x = Price)) +
  geom_histogram(
    aes(y = after_stat(density)), bins = 30, fill = "coral", color = "black", alpha = 0.6) +
  geom_density(
    color = "darkred", linewidth = 1.2) +
  labs(
    title = "Distribuzione di Price",
    x = "Price",
    y = "Densità") +
  theme_bw()

Dal grafico della variabile Price si osserva una distribuzione unimodale con un picco intorno a valori di prezzo tra 20 e 25, indicando che la maggior parte degli immobili si colloca in questa fascia di prezzi. La presenza di una coda destra più estesa suggerisce una discreta asimmetria positiva, con alcuni immobili aventi valori di prezzo più elevati rispetto alla media.

Data la presenza d asimmetrie significative in diverse variabili del dataset, è stata condotta un’analisi preliminae degli outlier per identificare valori anomali che potessero influenzare l’interpretazione dei risultati. Per individuare i valori anomali nelle variabili numeriche del dataset, è stato utilizzato il metodo dell’intervallo interquartile (IQR), in cui si è estratto il sottoinsieme delle colonne numeriche del dataset originale. Per ogni variabile numerica x sono stati calcolati:

Per individuare i valori anomali (outlier), si estendono i limiti naturali dei dati oltre il 25° e 75° percentile, calcolando:

Infine si sono contati i valori al di fuori di questi limiti.

df_num <- df%>%select_if(is.numeric)

outliers <- sapply(df_num, function(x) {
  Q1 <- quantile(x, 0.25, na.rm= TRUE)
  Q3 <- quantile(x, 0.75, na.rm = TRUE)
  IQRV <- Q3 - Q1
  lower <- Q1- 1.5 * IQRV
  upper <- Q3 + 1.5 * IQRV
  sum(x<lower | x > upper, na.rm = TRUE)
  })

outliers
##    CRIM      ZN   INDUS     NOX      RM     AGE     DIS     RAD     TAX PTRATIO 
##      66      68       0       0      30       0       5       0       0      15 
##       B   LSTAT   Price 
##      77       7      40

L’individuazione dei valori anomali e stata effettuata sulle variabili numeriche utilizzando il metodo dell intervallo interquartile IQR. Per ciascuna variabile sono stati calcolati il primo quartile Q1 e il terzo quartile Q3 e l intervallo interquartile e stato definito come la differenza tra Q3 e Q1. I valori sono stati considerati outlier quando risultavano inferiori al limite inferiore o superiori al limite superiore determinati a partire dall IQR.

Per gestire i valori anomali nelle variabili numeriche, è stato applicato un procedimento di capping ai percentili 1° e 99°. Questo metodo consiste nel sostituire tutti i valori inferiori al 1° percentile con il valore stesso del 1° percentile e, analogamente, sostituire tutti i valori superiori al 99° percentile con il valore limite corrispondente.In questo modo si “limita” l’estensione dei valori estremi senza eliminarli, preservando la struttura generale dei dati ma riducendo l’impatto degli outlier sulle analisi e sui modelli.

df_num_capped <- as.data.frame(lapply(df_num, function(x) {
  low_limit <- quantile(x, probs = 0.01)
  high_limit <- quantile(x, probs = 0.99)
  x <- ifelse(x < low_limit, low_limit,
              ifelse(x > high_limit, high_limit, x))
  return(x)
}))

summary(df_num_capped)
##       CRIM                ZN           INDUS             NOX        
##  Min.   : 0.01361   Min.   : 0.0   Min.   : 1.254   Min.   :0.3980  
##  1st Qu.: 0.08205   1st Qu.: 0.0   1st Qu.: 5.190   1st Qu.:0.4490  
##  Median : 0.25651   Median : 0.0   Median : 9.690   Median :0.5380  
##  Mean   : 3.37518   Mean   :11.3   Mean   :11.119   Mean   :0.5548  
##  3rd Qu.: 3.67708   3rd Qu.:12.5   3rd Qu.:18.100   3rd Qu.:0.6240  
##  Max.   :41.37033   Max.   :90.0   Max.   :25.650   Max.   :0.8710  
##        RM             AGE              DIS             RAD        
##  Min.   :4.524   Min.   :  6.61   Min.   :1.207   Min.   : 1.000  
##  1st Qu.:5.886   1st Qu.: 45.02   1st Qu.:2.100   1st Qu.: 4.000  
##  Median :6.208   Median : 77.50   Median :3.207   Median : 5.000  
##  Mean   :6.287   Mean   : 68.58   Mean   :3.779   Mean   : 9.549  
##  3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.:5.188   3rd Qu.:24.000  
##  Max.   :8.335   Max.   :100.00   Max.   :9.223   Max.   :24.000  
##       TAX           PTRATIO            B              LSTAT       
##  Min.   :188.0   Min.   :13.00   Min.   :  6.73   Min.   : 2.883  
##  1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38   1st Qu.: 6.950  
##  Median :330.0   Median :19.05   Median :391.44   Median :11.360  
##  Mean   :407.8   Mean   :18.45   Mean   :356.72   Mean   :12.642  
##  3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23   3rd Qu.:16.955  
##  Max.   :666.0   Max.   :21.20   Max.   :396.90   Max.   :33.919  
##      Price      
##  Min.   : 7.01  
##  1st Qu.:17.02  
##  Median :21.20  
##  Mean   :22.54  
##  3rd Qu.:25.00  
##  Max.   :50.00
par(mfrow = c(2,3)) 
par(cex=0.6) 
vars <- c("CRIM", 
          "ZN", 
          "RM", 
          "Price",
          "B", 
          "DIS") 
for (v in vars){ 
  qqnorm(df_num_capped[[v]], 
           main = paste("Q-Q Plot per", v)) 
  qqline(df_num_capped[[v]], col="mediumpurple") }

Nonostante l’applicazione del capping per limitare l’effetto degli outlier estremi, diverse variabili continuano a mostrare una forte asimmetria nelle loro distribuzioni, come evidenziato dai QQ plot. Questa situazione è probabilmente dovuta alla natura intrinseca delle variabili, che presentano distribuzioni fortemente sbilanciate con code pesanti o valori concentrati in fasce ristrette. Il capping infatti riduce solo l’impatto degli outlier più estremi, senza modificare la forma complessiva della distribuzione.

In questo caso si considera l’esclusione di queste variabili particolarmente problematiche per migliorare la qualità complessiva del modello.

vars <- df_num_capped[, c("NOX", "RM", "LSTAT", "PTRATIO", "Price", "INDUS", "TAX", "AGE", "DIS")]
ggplot(df, aes(x = CHAS, y = Price)) +
  geom_boxplot(fill = "lightgreen") +
  labs(title = "Distribuzione del prezzo in base a CHAS",
       x = "Confina con fiume Charles",
       y = "Prezzo") +
  theme_minimal()

L’analisi bivariata tra il prezzo delle abitazioni e la variabile CHAS evidenzia differenze significative in funzione della vicinanza al fiume Charles. Il boxplot mostra come le abitazioni situate in prossimità del fiume presentino una mediana del prezzo più elevata rispetto a quelle che non confinano con il fiume. In particolare:

cor_matrix <- cor(vars)

corrplot(
  cor_matrix,
  method = "color",
  type = "upper",
  addCoef.col = "black",
  number.cex = 0.7,
  tl.cex = 0.8,
  diag = FALSE,
  col = colorRampPalette(c("steelblue", "white", "firebrick"))(200)
)

La matrice di correlazione evidenzia diverse relazioni significative tra le variabili analizzate. Si osservano correlazioni positive molto forti tra l’indice industriale (INDUS) e il livello di inquinamento da ossidi di azoto (NOX) (r = 0.77), nonché tra INDUS e l’imposta sulla proprietà (TAX) (r = 0.65). Questi risultati indicano che le aree con maggiore attività industriale tendono a presentare livelli di inquinamento e tassazioni più elevati.

Il numero medio di stanze per abitazione (RM) risulta fortemente correlato negativamente con NOX (r = -0.62) e con la percentuale di popolazione a basso status socioeconomico (LSTAT) (r = -0.36), mentre è positivamente correlato con il prezzo delle abitazioni (Price) (r = 0.71), suggerendo che abitazioni più spaziose sono associate a zone meno inquinate, con popolazioni meno svantaggiate e con valori immobiliari maggiori. Il prezzo delle abitazioni presenta inoltre una correlazione negativa marcata con LSTAT (r = -0.74), sottolineando come aree con popolazioni a basso status siano associate a prezzi inferiori.

Altre relazioni significative includono la correlazione negativa tra il rapporto studenti-insegnanti (PTRATIO) e il prezzo (r = -0.48) e la correlazione positiva tra l’età delle abitazioni (AGE) e le imposte (TAX) (r = 0.51).

Per migliorare l’interpretazione delle relazioni tra variabili nel matrix plot, è stata implementata una funzione personalizzata che sostituisce il classico scatter plot nel pannello superiore con la visualizzazione del coefficiente di correlazione di Pearson. La funzione calcola il valore della correlazione tra le due variabili considerando solo le osservazioni complete, formattando il risultato con due cifre decimali per una migliore leggibilità. Il valore numerico viene mostrato centralmente nel pannello, e assume colore rosso quando la correlazione è forte (valore assoluto superiore a 0.7), altrimenti viene visualizzato in nero, favorendo così un’immediata individuazione delle associazioni significative.

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor = 1.2, ...) {
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- cor(x, y, use = "complete.obs")
  txt <- formatC(r, format = "f", digits = digits)
  txt <- paste0(prefix, txt)
  col <- ifelse(abs(r) > 0.7, "red", "black")
  text(0.5, 0.5, txt, cex = cex.cor, col = col)
}

pairs(vars, upper.panel = panel.cor, lower.panel = panel.smooth)

Nel pairs plot sono rappresentate le principali relazioni bivariate tra le variabili del dataset. Nei pannelli inferiori sono riportati gli scatter plot con la curva di regressione locale (loess), che evidenziano la direzione e la forma delle relazioni tra ogni coppia di variabili. Ad esempio, si osserva una relazione inversa tra la percentuale di popolazione a basso status economico (LSTAT) e il prezzo dell’abitazione (Price). Al contrario, il numero medio di stanze per abitazione (RM) risulta correlato positivamente con il prezzo.

Nei pannelli superiori sono invece visualizzati i coefficienti di correlazione di Pearson fra le variabili, con valori numerici arrotondati a due cifre decimali. Le correlazioni più forti (valore assoluto superiore a 0.7) sono evidenziate in rosso per facilitarne l’individuazione. Tra queste si segnalano:

ANALISI DELLE COMPONENTI PRINCIPALI (PCA)

L’Analisi delle Componenti Principali (PCA) è una tecnica di riduzione della dimensionalità ampiamente utilizzata nell’analisi multivariata dei dati. Essa consente di sintetizzare l’informazione contenuta in un set di variabili correlate trasformandole in un numero minore di variabili, chiamate componenti principali, che sono combinazioni lineari ortogonali delle variabili originali.

Le componenti principali sono ordinate in modo tale da catturare progressivamente la massima varianza dei dati, facilitando così la visualizzazione, l’interpretazione e la successiva modellazione, soprattutto quando il dataset presenta elevate correlazioni tra le variabili o un numero molto grande di feature.

Un passaggio cruciale nella PCA è la standardizzazione delle variabili, ovvero la trasformazione delle variabili affinché abbiano media zero e varianza uno. Questo è particolarmente importante quando le variabili originali sono misurate su scale differenti, poiché la PCA è sensibile alla scala dei dati e tenderebbe a privilegiare le variabili con varianza maggiore, distorcendo i risultati.

La standardizzazione garantisce che ogni variabile contribuisca in modo equo alla determinazione delle componenti principali, evitando preponderanze dovute all’unità di misura o alla dispersione.

# Standardizza
vars_scaled <- scale(vars)
pca_result <- prcomp(vars_scaled, center = TRUE, scale. = TRUE)
summary(pca_result)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.2566 1.2234 0.89421 0.69859 0.60157 0.46882 0.45725
## Proportion of Variance 0.5658 0.1663 0.08885 0.05423 0.04021 0.02442 0.02323
## Cumulative Proportion  0.5658 0.7321 0.82098 0.87520 0.91541 0.93983 0.96307
##                           PC8     PC9
## Standard deviation     0.4399 0.37272
## Proportion of Variance 0.0215 0.01544
## Cumulative Proportion  0.9846 1.00000

Nel nostro caso, la prima componente principale (PC1) spiega circa il 56.6% della varianza totale, una quota significativa del dataset. La seconda componente (PC2) aggiunge un ulteriore 16.6%. Insieme, le prime due componenti spiegano oltre il 73% della varianza complessiva, indicando che gran parte dell’informazione multidimensionale è sintetizzata efficacemente in uno spazio bidimensionale. Le componenti dalla terza in poi spiegano porzioni di varianza sempre più piccole, con la terza componente che contribuisce con circa il 8.9%.

fviz_eig(pca_result, addlabels = TRUE, ylim = c(0, 60),
         main = "Varianza cumulativa spiegata")

Il grafico della varianza cumulativa spiegata (scree plot) mostra chiaramente come la prima componente (PC1) catturi la maggior parte dell’informazione (56.6%), seguita da un deciso calo nell’apporto della seconda componente (16.6%). Usando il metodo del gomito, si identifica il punto in cui l’incremento aggiuntivo di varianza spiegata si riduce sensibilmente: in questo caso, dopo la seconda componente la pendenza della curva si appiattisce notevolmente. Pertanto, si decide di mantenere le prime due componenti principali che insieme spiegano oltre il 73% della varianza totale, garantendo un buon compromesso tra la semplicità del modello e la quantità di informazione preservata.

pca_result$rotation
##                PC1         PC2         PC3         PC4        PC5          PC6
## NOX      0.3723236 -0.30201554  0.14179403 -0.19905540  0.0167699 -0.426349270
## RM      -0.2596499 -0.50730325 -0.36800742  0.04628507 -0.5730561 -0.215483672
## LSTAT    0.3685906  0.20284296  0.28412033  0.08272247 -0.4309880  0.526469680
## PTRATIO  0.2314461  0.34286204 -0.75118885  0.41425778  0.1147280 -0.026993201
## Price   -0.3172504 -0.47054010 -0.12578240  0.03947401  0.3043325  0.598369196
## INDUS    0.3841968 -0.15452358 -0.09661700 -0.25598099  0.3345003 -0.123414454
## TAX      0.3462726 -0.07960311 -0.39583576 -0.59830046 -0.2319013  0.323760840
## AGE      0.3470844 -0.29193883  0.12929442  0.52557646 -0.2843284 -0.007458728
## DIS     -0.3404459  0.39481632 -0.02994807 -0.28018496 -0.3650791 -0.125461522
##                 PC7         PC8          PC9
## NOX     -0.29598396 -0.21490151  0.628163232
## RM       0.33025503 -0.22737089 -0.045399789
## LSTAT    0.33903503 -0.29073126  0.270255324
## PTRATIO -0.02364551 -0.08064199  0.268385475
## Price   -0.03917119  0.14642605  0.432349206
## INDUS    0.71956168  0.32978628 -0.002142242
## TAX     -0.35617076  0.11082490 -0.247729999
## AGE     -0.19377186  0.61251574 -0.087949721
## DIS      0.06413051  0.54169619  0.449717000

I coefficienti delle variabili, noti come loadings, rappresentano i pesi con cui ciascuna variabile originaria contribuisce alle componenti principali derivate dalla PCA. Questi coefficienti indicano quanto e in che direzione una variabile influenza ciascuna componente.

Per esempio, nella prima componente (PC1), sono importanti variabili come NOX, LSTAT, INDUS e TAX con segno positivo, mentre RM e DIS hanno un valore negativo. Questo significa che PC1 rappresenta un insieme di caratteristiche che vanno da zone più inquinate e con condizioni peggiori a zone con case più grandi e più distanti dall’inquinamento.

La seconda componente (PC2) ha invece carichi negativi su RM e Price e positivi su DIS e PTRATIO, indicando una diversa combinazione di variabili legate alla qualità abitativa e alla distanza dai servizi.

# PC1 vs PC2 
fviz_pca_var(pca_result, 
             axes = c(1, 2), 
             col.var = "contrib", 
             gradient.cols = c("blue", "yellow", "red"), 
             repel = TRUE, 
             title = "Mappa variabili - PC1 vs PC2")

Il grafico rappresenta la proiezione delle variabili originali nello spazio definito dalle prime due componenti principali (PC1 e PC2) della PCA. Queste due componenti spiegano complessivamente circa il 73% della varianza totale, fornendo una sintesi efficace delle informazioni multidimensionali.

Le frecce indicano la direzione e l’intensità della correlazione di ciascuna variabile con le due componenti: frecce più lunghe corrispondono a variabili maggiormente correlate con la componente rappresentata sull’asse, mentre l’angolo tra due frecce segnala la relazione di correlazione tra due variabili (freccia parallela = correlazione positiva, freccia opposta = correlazione negativa, angolo retto = nessuna correlazione).

Dal grafico emergono diverse osservazioni rilevanti:

La colorazione delle frecce, basata sulla contribuzione delle variabili alle componenti, mette in evidenza come NOX, DIS e Price siano tra le variabili più influenti nella definizione delle prime due componenti.

ANALISI CLUSTER

Il clustering è una tecnica di analisi dei dati che permette di suddividere un insieme di osservazioni in gruppi (cluster) omogenei, in cui gli elementi all’interno di uno stesso gruppo sono più simili tra loro rispetto a quelli appartenenti ad altri gruppi. Questo approccio è utile per scoprire strutture nascoste nel dataset e per raggruppare dati con caratteristiche simili.

Nel contesto di dati ad alta dimensionalità, come quelli analizzati con la PCA, il clustering viene spesso applicato sulle componenti principali per semplificare la struttura e migliorare la qualità del raggruppamento.

Prima di eseguire il clustering, è importante determinare il numero ottimale di cluster da utilizzare. Per questo si utilizzano metodi come:

Per il metodo del gomito viene creato un vettore wss per memorizzare la somma delle distanze quadratiche all’interno dei cluster (Within-Cluster Sum of Squares) per valori di cluster da 1 a 10. Per ogni valore di k, il clustering k-means viene eseguito sui dati scalati (vars_scaled), con 25 inizializzazioni casuali per ridurre l’effetto di soluzioni subottimali. Il valore totale di somma delle distanze interne è salvato nel vettore wss.

Infine, viene tracciato un grafico che mostra come la WSS diminuisca all’aumentare del numero di cluster. Il grafico aiuta a identificare il punto in cui aggiungere ulteriori cluster non migliora significativamente la qualità del raggruppamento, ovvero il cosiddetto “gomito”. Questo punto è considerato il numero ottimale di cluster da utilizzare.

wss <- numeric(10)
set.seed(123)
for(k in 1:10){
  wss[k] <- kmeans(vars_scaled, centers = k, nstart = 25)$tot.withinss
}
plot(1:10, wss, type="b",
     xlab = "Numero di cluster",
     ylab = "WSS",
     main = "Metodo del gomito")

Nel grafico si osserva un forte calo di WSS passando da 1 a 2 cluster, seguito da una riduzione più graduale per valori maggiori di cluster. Il “gomito”, cioè il punto in cui la curva inizia a livellarsi, si presenta intorno a 2-3 cluster.

Successivamente si calcola la media dell’indice silhouette per un numero di cluster variabile da 2 a 10, al fine di valutare la qualità del clustering k-means.

avg_sil <- numeric(10)

for(k in 2:10){
  kmeans_k <- kmeans(vars_scaled, centers = k, nstart = 25)
  sil <- silhouette(kmeans_k$cluster, dist(vars_scaled))
  avg_sil[k] <- mean(sil[, 3])
}

plot(2:10, avg_sil[2:10], type = "b", 
     xlab = "Numero di cluster", ylab = "silhouette",
     main = "Metodo silhouette")

Nel grafico, il valore massimo della silhouette media si trova a 2 cluster, suggerendo che questa configurazione produce la miglior separazione e coesione dei gruppi. Aumentando il numero di cluster, l’indice tende a diminuire, indicando che ulteriori suddivisioni non migliorano la qualità del clustering.

Questo risultato conferma e rafforza l’indicazione già emersa dal metodo del gomito, evidenziando che 2 cluster rappresentano una scelta adeguata per i dati analizzati.

set.seed(123) # riproducibilità
kmeans_result2 <- kmeans(vars_scaled, centers = 2, nstart = 25)
kmeans_result2$size
## [1] 228 278

I valori indicano il numero di osservazioni presenti in ciascuno dei due cluster identificati dall’algoritmo di clustering. In questo caso, il primo cluster contiene 228 elementi, mentre il secondo cluster ne contiene 278. Questa distribuzione delle dimensioni suggerisce che i due gruppi sono abbastanza bilanciati in termini di numerosità, il che può facilitare ulteriori analisi e interpretazioni, evitando problemi legati a cluster troppo sbilanciati o piccoli. Per valutare la qualità del clustering effettuato con il metodo k-means, sono stati considerati due indicatori principali:

# calcolo della coesione totale
withinss2 <- kmeans_result2$tot.withinss

# variabilita
totss2 <- kmeans_result2$betweenss / kmeans_result2$totss * 100

withinss2
## [1] 2592.428
totss2
## [1] 42.96088

Per k=2, si ottiene una coesione interna (tot.withinss) pari a 2592.428 e una percentuale di variabilità spiegata tra i cluster del 42.96%. Questi valori indicano che i cluster non sono molto compatti internamente e la separazione tra i gruppi non è particolarmente marcata.

Per verificare se è possibile migliorare la qualità del clustering, si procede a calcolare il modello con k=3 cluster, confrontando i valori di coesione interna e di variabilità spiegata ottenuti.

set.seed(123)
kmeans_result3 <- kmeans(vars_scaled, centers = 3, nstart = 25)
kmeans_result3$size
## [1] 184 130 192

Il clustering con k=3 ha prodotto tre gruppi di dimensioni relativamente bilanciate: il primo cluster contiene 184 elementi, il secondo 130, e il terzo 192.

# calcolo della coesione totale
withinss3 <- kmeans_result3$tot.withinss

# variabilita
totss3 <- kmeans_result3$betweenss / kmeans_result3$totss * 100

withinss3
## [1] 2139.345
totss3
## [1] 52.92971

Per k=3, il valore di coesione interna (tot.withinss) è 2139.345, inferiore rispetto a valori ottenuti con un numero minore di cluster, il che indica che i punti all’interno di ciascun cluster sono più vicini tra loro e quindi i gruppi sono più compatti.

La percentuale di variabilità spiegata dal modello è circa il 52.93%, suggerendo che più della metà della variabilità totale dei dati è catturata dalla differenziazione tra i cluster. Ciò indica una buona separazione tra i gruppi e giustifica la scelta di k=3 come numero di cluster adeguato.

Per approfondire la comprensione dei gruppi individuati dal clustering con k=3, sono state integratoe le informazioni risultanti dall’assegnazione dei cluster direttamente nel dataset originale. In particolare, la variabile cluster è stata aggiunta al data frame contenente le variabili di interesse tramite la funzione mutate(), consentendo di associare ad ogni osservazione il proprio cluster di appartenenza come fattore categorico.

Successivamente, per ciascun cluster, sono state calcolate le medie di tutte le variabili numeriche mediante la funzione summarise() combinata con group_by(). Questa operazione ha prodotto una tabella riassuntiva che evidenzia le caratteristiche medie di ogni gruppo, facilitando l’identificazione delle differenze distintive tra i cluster.

vars <- vars %>%
  mutate(cluster = as.factor(kmeans_result3$cluster))

# calcola la media per cluster
summary_cluster <- vars %>%
  group_by(cluster) %>%
  summarise(across(where(is.numeric), mean))

summary_cluster
## # A tibble: 3 × 10
##   cluster   NOX    RM LSTAT PTRATIO Price INDUS   TAX   AGE   DIS
##   <fct>   <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1       0.677  5.98 18.2     19.5  16.8 19.1   593.  91.2  2.01
## 2 2       0.467  7.05  5.64    16.5  33.3  4.66  290.  46.8  5.12
## 3 3       0.497  6.06 12.0     18.8  20.8  7.84  310.  61.7  4.56

Dall’analisi delle medie delle variabili chiave per ciascuno dei tre cluster emergono profili distinti che aiutano a comprendere le differenze sottostanti tra i gruppi individuati dal clustering:

fviz_cluster(kmeans_result3, data = vars_scaled,
             geom = "point", 
             ellipse.type = "convex", 
             palette = "ginger",
             ggtheme = theme_minimal(),
             main = "Clustering con k=3")

La figura illustra la segmentazione ottenuta mediante il metodo di clustering con k=3, visualizzata sulle prime due componenti principali (Dim1 e Dim2) che cumulativamente spiegano il 73,2% della variabilità complessiva del dataset. L’analisi grafica evidenzia una netta separazione tra i tre cluster, a testimonianza dell’efficacia del modello nel differenziare gruppi omogenei. Si osserva una leggera sovrapposizione fra i cluster 1 e 3, indicativa di una zona di confine dove alcune osservazioni potrebbero avere caratteristiche comuni a entrambi i gruppi, possibile effetto di una riduzione dimensionale o di una vicinanza strutturale reale nei dati. Il cluster 2 appare ben separato e compatto rispetto agli altri due, suggerendo un gruppo con caratteristiche particolarmente omogenee.

REGRESSIONE LINEARE MULTIPLA

La regressione lineare multipla è una tecnica statistica utilizzata per modellare la relazione tra una variabile dipendente continua e più variabili indipendenti. Questo metodo permette di quantificare l’effetto di ciascuna variabile esplicativa sul valore della variabile di interesse, controllando contemporaneamente l’influenza delle altre. Nel presente studio, la regressione lineare multipla sarà applicata per analizzare la relazione tra il prezzo degli immobili e le variabili selezionate, al fine di fornire un modello predittivo e interpretativo efficace.

Per garantire una valutazione affidabile delle performance del modello, il dataset è stato suddiviso in due sottoinsiemi: un set di addestramento, utilizzato per stimare i parametri del modello, e un set di test, riservato alla verifica delle capacità predittive su dati nuovi e non utilizzati in fase di stima. La suddivisione è stata effettuata in modo casuale ma riproducibile, assicurando che la stessa divisione possa essere replicata per confronti futuri e per garantire la trasparenza dell’analisi. Infine, è stata controllata la dimensione dei due sottoinsiemi per confermare che la proporzione tra training e test sia rispettata, garantendo un’adeguata quantità di dati per entrambe le fasi.

set.seed(123)  # riproducibilità
n <- nrow(df)
train_indices <- sample(seq_len(n), size = floor(0.75 * n))  # 75% training
train_data <- df[train_indices, ]
test_data <- df[-train_indices, ]

dim(train_data)
## [1] 379  14
dim(test_data)
## [1] 127  14

La suddivisione del dataset originale ha prodotto due sottoinsiemi coerenti con le proporzioni previste: il training set contiene 379 osservazioni mentre il test set ne contiene 127, entrambi con 14 variabili.

Per ottimizzare il modello di regressione lineare multipla e migliorare l’equilibrio tra complessità e capacità predittiva, ho applicato la procedura di selezione automatica delle variabili tramite la funzione step(). Questa metodologia, nota come stepwise selection, consente di aggiungere o rimuovere iterativamente variabili dal modello in base a criteri statistici (come l’Akaike Information Criterion, AIC), al fine di individuare il sottoinsieme di predittori più rilevanti.

model1 <- lm(Price ~  ., data = train_data)
model_step <- step(model1)
## Start:  AIC=1192.59
## Price ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + 
##     TAX + PTRATIO + B + LSTAT
## 
##           Df Sum of Sq     RSS    AIC
## - INDUS    1      0.20  8187.4 1190.6
## - AGE      1      0.63  8187.8 1190.6
## <none>                  8187.2 1192.6
## - B        1     91.95  8279.1 1194.8
## - TAX      1    125.88  8313.1 1196.4
## - CRIM     1    197.96  8385.1 1199.6
## - ZN       1    265.67  8452.9 1202.7
## - NOX      1    305.50  8492.7 1204.5
## - CHAS     1    318.85  8506.0 1205.1
## - RAD      1    333.24  8520.4 1205.7
## - PTRATIO  1    765.69  8952.9 1224.5
## - DIS      1   1025.27  9212.5 1235.3
## - RM       1   1071.57  9258.8 1237.2
## - LSTAT    1   2448.24 10635.4 1289.7
## 
## Step:  AIC=1190.6
## Price ~ CRIM + ZN + CHAS + NOX + RM + AGE + DIS + RAD + TAX + 
##     PTRATIO + B + LSTAT
## 
##           Df Sum of Sq     RSS    AIC
## - AGE      1      0.67  8188.1 1188.6
## <none>                  8187.4 1190.6
## - B        1     92.68  8280.1 1192.9
## - TAX      1    177.84  8365.2 1196.7
## - CRIM     1    197.77  8385.2 1197.6
## - ZN       1    273.98  8461.4 1201.1
## - CHAS     1    322.96  8510.3 1203.3
## - NOX      1    336.58  8524.0 1203.9
## - RAD      1    378.24  8565.6 1205.7
## - PTRATIO  1    771.76  8959.1 1222.7
## - DIS      1   1070.51  9257.9 1235.2
## - RM       1   1090.39  9277.8 1236.0
## - LSTAT    1   2451.19 10638.6 1287.8
## 
## Step:  AIC=1188.63
## Price ~ CRIM + ZN + CHAS + NOX + RM + DIS + RAD + TAX + PTRATIO + 
##     B + LSTAT
## 
##           Df Sum of Sq     RSS    AIC
## <none>                  8188.1 1188.6
## - B        1     93.64  8281.7 1190.9
## - TAX      1    177.17  8365.2 1194.7
## - CRIM     1    197.49  8385.5 1195.7
## - ZN       1    274.77  8462.8 1199.1
## - CHAS     1    325.79  8513.8 1201.4
## - NOX      1    357.37  8545.4 1202.8
## - RAD      1    379.18  8567.2 1203.8
## - PTRATIO  1    772.20  8960.3 1220.8
## - RM       1   1118.87  9306.9 1235.2
## - DIS      1   1216.33  9404.4 1239.1
## - LSTAT    1   2636.59 10824.6 1292.4
summary(model_step)
## 
## Call:
## lm(formula = Price ~ CRIM + ZN + CHAS + NOX + RM + DIS + RAD + 
##     TAX + PTRATIO + B + LSTAT, data = train_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.6715  -2.6671  -0.5173   1.6211  24.6904 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  38.898906   5.849474   6.650 1.07e-10 ***
## CRIM         -0.103016   0.034625  -2.975 0.003122 ** 
## ZN            0.054315   0.015477   3.509 0.000505 ***
## CHASSi        3.660510   0.957928   3.821 0.000156 ***
## NOX         -16.331036   4.080491  -4.002 7.59e-05 ***
## RM            3.350842   0.473174   7.082 7.31e-12 ***
## DIS          -1.560359   0.211327  -7.384 1.04e-12 ***
## RAD           0.292992   0.071071   4.123 4.64e-05 ***
## TAX          -0.010568   0.003750  -2.818 0.005095 ** 
## PTRATIO      -0.889102   0.151127  -5.883 9.07e-09 ***
## B             0.006808   0.003323   2.049 0.041205 *  
## LSTAT        -0.596534   0.054875 -10.871  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.723 on 367 degrees of freedom
## Multiple R-squared:  0.7454, Adjusted R-squared:  0.7377 
## F-statistic: 97.67 on 11 and 367 DF,  p-value: < 2.2e-16

L’algoritmo ha iniziato con un modello completo che includeva tutte le 13 variabili predittive disponibili e ha valutato iterativamente l’impatto della rimozione di ciascuna variabile, confrontando il valore di AIC risultante.

Il processo ha identificato inizialmente due variabili (AGE e INDUS) la cui esclusione riduceva l’AIC (da 1192 a 1188), segnalando un miglioramento del compromesso tra bontà di adattamento e parsimonia del modello. Procedendo con la rimozione di queste variabili, la procedura ha fermato la selezione quando ulteriori esclusioni avrebbero peggiorato il modello (ovvero portato a un aumento dell’AIC).

Il modello finale quindi comprende 11 variabili: CRIM, ZN, CHAS, NOX, RM, DIS, RAD, TAX, PTRATIO, B e LSTAT. I coefficienti stimati mostrano che la maggioranza delle variabili mantiene un effetto significativo sul prezzo, confermando la rilevanza di tali predittori nella spiegazione della variabilità dei dati. Il modello finale presenta un Adjusted R-squared molto simile a quello iniziale (circa 0.7377), confermando che la semplificazione non ha sacrificato l’efficacia predittiva.

model2 <- lm(Price ~ CRIM + ZN + CHAS + NOX + RM + DIS + RAD + TAX + PTRATIO + B + LSTAT, data = train_data)
vif(model2)
##     CRIM       ZN     CHAS      NOX       RM      DIS      RAD      TAX 
## 1.773329 2.331141 1.066539 3.723866 1.834417 3.511882 6.401656 6.764313 
##  PTRATIO        B    LSTAT 
## 1.800462 1.347165 2.661182

Dopo aver identificato una forte correlazione tra la variabile TAX e altre variabili predittive del modello, si è scelto di escluderla per ridurre la multicollinearità.

model3 <- lm(Price ~ CRIM + ZN + CHAS + NOX + RM + DIS + RAD + PTRATIO + B + LSTAT, data = train_data)
summary(model3)
## 
## Call:
## lm(formula = Price ~ CRIM + ZN + CHAS + NOX + RM + DIS + RAD + 
##     PTRATIO + B + LSTAT, data = train_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.3916  -2.9432  -0.4002   1.7085  24.7055 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.168837   5.871770   6.330 7.13e-10 ***
## CRIM         -0.099736   0.034931  -2.855  0.00454 ** 
## ZN            0.046262   0.015354   3.013  0.00277 ** 
## CHASSi        3.919673   0.962453   4.073 5.70e-05 ***
## NOX         -19.020137   4.004579  -4.750 2.93e-06 ***
## RM            3.535863   0.472995   7.475 5.68e-13 ***
## DIS          -1.498626   0.212162  -7.064 8.16e-12 ***
## RAD           0.142739   0.047432   3.009  0.00280 ** 
## PTRATIO      -0.946749   0.151142  -6.264 1.05e-09 ***
## B             0.007200   0.003351   2.148  0.03235 *  
## LSTAT        -0.602904   0.055343 -10.894  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.768 on 368 degrees of freedom
## Multiple R-squared:  0.7399, Adjusted R-squared:  0.7328 
## F-statistic: 104.7 on 10 and 368 DF,  p-value: < 2.2e-16
vif(model3)
##     CRIM       ZN     CHAS      NOX       RM      DIS      RAD  PTRATIO 
## 1.771324 2.251669 1.056708 3.520204 1.799096 3.474144 2.798531 1.767473 
##        B    LSTAT 
## 1.344810 2.656666

Il modello di regressione lineare stimato spiega circa il 73.4% della variabilità della variabile dipendente Price (R² = 0.7328), suggerendo una buona capacità predittiva complessiva.

I coefficienti stimati indicano l’effetto marginale di ciascun predittore sul prezzo, tenendo costanti le altre variabili. La maggior parte dei coefficienti risulta statisticamente significativa (p-value < 0.01), in particolare variabili come RM, LSTAT, NOX e PTRATIO mostrano forti effetti con elevata significatività. Per quanto riguarda la multicollinearità, i valori di VIF riportati sono tutti inferiori a 5, con il valore massimo di circa 3.52 (per NOX), indicazione che la collinearità tra variabili è contenuta e non dovrebbe compromettere la stabilità delle stime.

Per valutare la bontà del modello di regressione e verificare l’ipotesi di normalità degli errori, si è calcolata la distribuzione dei residui del modello. In particolare, è stato estratto il vettore dei residui dal modello (model3) e ne sono stati calcolati la media e la deviazione standard. Successivamente, si è costruito un istogramma con densità normalizzata dei residui affiancato dalla curva teorica della distribuzione normale con media e deviazione standard corrispondenti.

residui <- residuals(model3)
residuals_df <- data.frame(residui = residui)

# Calcolo media e deviazione standard dei residui
mu <- mean(residui)
sigma <- sd(residui)

ggplot(residuals_df, aes(x = residui)) +
  geom_histogram(aes(y = ..density..), 
                 color = 'black', fill = 'skyblue', binwidth = 1) +
  stat_function(fun = dnorm, args = list(mean = mu, sd = sigma),
                color = "red", size = 1) +
  ggtitle("Istogramma dei residui") +
  ylab("Densità")

L’istogramma dei residui del modello mostra una distribuzione approssimativamente centrata attorno a zero, condizione necessaria per un modello adeguato. Tuttavia, si osserva una leggera asimmetria positiva, con una coda più estesa verso valori residui positivi rispetto a quanto atteso da una distribuzione normale teorica. Questa deviazione suggerisce una non perfetta normalità degli errori, che potrebbe indicare la presenza di eteroschedasticità, outlier o componenti sistematiche non adeguatamente catturate dal modello attuale.

Per approfondire queste osservazioni, si è analizzato il grafico Residui vs Valori Predetti, che consente di verificare l’ipotesi di omoscedasticità, ovvero che la varianza degli errori sia costante per ogni livello dei valori predetti. Infine, si è esaminato il Q-Q plot dei residui, utile per verificare la normalità della distribuzione degli errori confrontando quantili osservati con quantili teorici di una distribuzione normale.

par(mfrow = c(1, 2))

# Residui vs valori predetti
plot(model1$fitted.values, residuals(model3),
     xlab = "Valori Predetti",
     ylab = "Residui",
     main = "Residui vs valori predetti")
abline(h = 0, col = "red")

# QQ-plot dei residui
qqnorm(residuals(model3))
qqline(residuals(model3), col = "red")

par(mfrow = c(1, 1))

Il grafico Residui vs Valori Predetti mostra una distribuzione dei residui abbastanza uniforme lungo tutto l’intervallo dei valori predetti, senza variazioni evidenti nella dispersione. Questo indica che la varianza degli errori è approssimativamente costante, confermando l’assunzione di omoscedasticità necessaria per la validità del modello.

Nel Q-Q plot, i residui si allineano in modo soddisfacente lungo la linea di riferimento, suggerendo una distribuzione dei residui vicina alla normalità. Sono tuttavia presenti leggere deviazioni nelle code, in particolare nella parte superiore destra, coerenti con una lieve asimmetria già osservata nell’istogramma.

Per accertare la validità dell’assunzione di omoschedasticità nel modello, è stato applicato il test studentized di Breusch-Pagan.

bptest(model3)
## 
##  studentized Breusch-Pagan test
## 
## data:  model3
## BP = 43.423, df = 10, p-value = 4.177e-06

Il test ha restituito una statistica di 56.427 con 10 gradi di libertà e un p-value estremamente basso (4.177e-06), inferiore alla soglia convenzionale di significatività (0.05). Tale risultato permette di rifiutare l’ipotesi nulla di varianza costante degli errori, indicando la presenza di eteroschedasticità. Questa condizione implica che la varianza degli errori varia in modo significativo in relazione ai valori delle variabili predittive, il che potrebbe compromettere l’efficienza degli stimatori ordinari e l’affidabilità degli intervalli di confidenza e dei test statistici basati su errori standard classici.

Al fine di correggere la presenza di eteroschedasticità rilevata nel modello di regressione (come evidenziato dal test di Breusch-Pagan), è stata applicata una stima degli errori standard robusti utilizzando il metodo HC3 tramite la funzione coeftest() del pacchetto lmtest con la matrice di varianza-covarianza robusta.

Il metodo HC3 tratta gli errori in modo da ridurne il bias legato alla varianza non costante degli errori (eteroschedasticità), modificando la matrice di varianza-covarianza degli stimatori dei coefficienti per renderla più affidabile rispetto ai metodi più semplici come HC0 o HC1. Questa procedura non modifica le stime puntuali dei coefficienti, ma aggiorna gli errori standard, i valori t e i p-value tenendo conto della non costanza della varianza degli errori. Di conseguenza, le inferenze statistiche risultano più affidabili, riducendo il rischio di trarre conclusioni errate dovute a varianze non omogenee degli errori.

coeftest(model3, vcov = vcovHC(model3, type = "HC3"))
## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)  37.1688365   9.3219967  3.9872 8.065e-05 ***
## CRIM         -0.0997355   0.0390352 -2.5550 0.0110198 *  
## ZN            0.0462620   0.0161210  2.8697 0.0043461 ** 
## CHASSi        3.9196734   1.3463724  2.9113 0.0038187 ** 
## NOX         -19.0201371   4.1798612 -4.5504 7.283e-06 ***
## RM            3.5358626   0.9618044  3.6763 0.0002719 ***
## DIS          -1.4986263   0.2594870 -5.7753 1.635e-08 ***
## RAD           0.1427392   0.0574894  2.4829 0.0134775 *  
## PTRATIO      -0.9467489   0.1293975 -7.3166 1.608e-12 ***
## B             0.0071996   0.0039412  1.8267 0.0685480 .  
## LSTAT        -0.6029043   0.1074885 -5.6090 4.004e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Il modello di regressione lineare multipla, stimato con correzione dei coefficienti tramite errori standard robusti HC3, evidenzia effetti significativi di molteplici predittori sulla variabile risposta “Prezzo”.

Le variabili CRIM, ZN, CHAS, NOX, RM, DIS, RAD, PTRATIO e LSTAT risultano fortemente associate al prezzo con livelli di significatività elevati, mentre la variabile B si mostra marginalmente significativa (p ≈ 0.068).

I coefficienti stimati indicano che un aumento di NOX e LSTAT è associato a una diminuzione del prezzo, coerentemente con la natura di queste caratteristiche, mentre un incremento di RM o CHAS tende a far aumentare il prezzo mediamente. L’adozione degli errori standard robusti HC3 ha permesso di ottenere inferenze più affidabili correggendo gli effetti di eventuale eteroschedasticità, rafforzando la fiducia nei risultati.

predicted_test <- predict(model3, newdata = test_data)

ggplot(data = test_data, aes(x = Price, y = predicted_test)) +
  geom_point(color = "mediumpurple", alpha = 0.6, size = 2) +
  geom_abline(slope = 1, intercept = 0, color = "red", size = 1) + 
  labs(
    title = "Confronto tra valori osservati e predetti sul test set",
    x = "Prezzi osservati",
    y = "Prezzi predetti"
  ) +
  theme_minimal()

Lo scatterplot mostra un confronto tra i prezzi osservati e quelli predetti dal modello sul test set. In generale, si osserva che le predizioni seguono abbastanza bene la linea di identità, indicando una buona capacità del modello di stimare i valori reali. Tuttavia, si nota una certa dispersione attorno alla retta, soprattutto per i prezzi più elevati, dove il modello tende a sottostimare o sovrastimare in qualche caso. Questa variabilità suggerisce che il modello non riesce a catturare completamente tutte le dinamiche nei dati, soprattutto per valori estremi.

Per valutare la precisione delle predizioni del modello, sono state calcolate due metriche di errore comunemente utilizzate in ambito di regressione: il Root Mean Squared Error (RMSE) e il Mean Absolute Error (MAE). Il RMSE fornisce una misura della deviazione quadratica media tra i valori osservati e quelli predetti, attribuendo un peso maggiore agli errori più grandi, mentre il MAE rappresenta la media degli scarti assoluti, offrendo una valutazione più diretta e meno influenzata da outlier.

residuals <- test_data$Price - predicted_test  # residui

RMSE <- sqrt(mean(residuals^2))
MAE <- mean(abs(residuals))

RMSE
## [1] 4.962718
MAE
## [1] 3.698338

Le metriche di errore calcolate sul set di test indicano una buona capacità predittiva del modello. Il Root Mean Squared Error (RMSE) è pari a circa 4.96, mentre il Mean Absolute Error (MAE) risulta circa 3.70. Questi risultati indicano che, mediamente, le predizioni del modello si discostano dai valori reali rispettivamente di 4.96 e 3.70 unità di prezzo. Considerando la scala e la natura della variabile target, tali errori sono da ritenersi accettabili, confermando l’affidabilità del modello nel fornire stime precise.

CONCLUSIONI

L’analisi svolta sul Boston Housing Dataset ha permesso di identificare e quantificare i principali fattori che influenzano il prezzo medio delle abitazioni nell’area di Boston.

Dall’analisi esplorativa è emersa la presenza di distribuzioni asimmetriche e di valori anomali in alcune variabili, opportunamente gestiti con tecniche di capping per garantirne la robustezza nell’analisi successiva.

L’Analisi delle Componenti Principali ha sintetizzato efficacemente l’informazione multidimensionale riducendo il dataset a due componenti che spiegano complessivamente il 73% della varianza, facilitando la comprensione delle relazioni tra le variabili.

Il clustering ha identificato tre distinti gruppi di quartieri, caratterizzati da differenti condizioni ambientali, socio-economiche e abitative. In particolare, è stato possibile distinguere cluster rappresentanti aree svantaggiate con elevati livelli di inquinamento e minor qualità abitativa, cluster di qualità intermedia e cluster contraddistinti da caratteristiche ambientali favorevoli e prezzi elevati.

La regressione lineare multipla, sviluppata e raffinata attraverso procedure di selezione automatica e gestione della multicollinearità, ha mostrato un’elevata capacità predittiva, spiegando circa il 73% della variabilità del prezzo. Variabili quali il numero medio di stanze (RM), la percentuale di popolazione a basso status socio-economico (LSTAT), e la vicinanza al fiume Charles (CHAS) sono risultate fattori significativi con impatti attesi sul prezzo degli immobili.

La validazione del modello su un set di test indipendente ha confermato la sua buona capacità predittiva, con errori medi (RMSE ~4.96, MAE ~3.70) coerenti con l’ordine di grandezza della variabile target. L’analisi dei residui ha evidenziato una lieve eteroschedasticità, corretta mediante l’adozione di errori standard robusti, migliorando così l’affidabilità inferenziale del modello.