0 - Preparazione dell’ambiente di lavoro

In questa sezione si procede alla preparazione dell’ambiente di lavoro, quindi all’installazione dei pacchetti necessari per svolgere l’analisi sui dati. 

options(repos = c(CRAN = "https://cran.rstudio.com/")) # Per creazione del documento con knit
install.packages("dplyr") # Per manipolazioni sui dati e creazione di grafici
## 
## The downloaded binary packages are in
##  /var/folders/xd/gg9fwq69289958mc26w3vmbw0000gn/T//RtmpB3BHCA/downloaded_packages
install.packages("ggplot2") # Per creazione di grafici
## 
## The downloaded binary packages are in
##  /var/folders/xd/gg9fwq69289958mc26w3vmbw0000gn/T//RtmpB3BHCA/downloaded_packages
install.packages("moments") # Per calcolo di asimmetria e curtosi
## 
## The downloaded binary packages are in
##  /var/folders/xd/gg9fwq69289958mc26w3vmbw0000gn/T//RtmpB3BHCA/downloaded_packages
install.packages("gridExtra") # Per personalizzazione della finestra di output
## 
## The downloaded binary packages are in
##  /var/folders/xd/gg9fwq69289958mc26w3vmbw0000gn/T//RtmpB3BHCA/downloaded_packages
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(moments)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine


1 - Preparazione dei dati

In questa sezione viene effettuato l’import dei dati da github (repository Profession.AI) e vengono mostrate le prime 10 righe (osservazioni) del dataset. Infine si procede all’estrazione delle colonne come variabili e vengono descritte. 

url<-"https://raw.githubusercontent.com/ProfAI/statistica-descrittiva-per-datascience/refs/heads/main/progetto%20finale/Real%20Estate%20Texas.csv" # Link dataset in repository github
dati <- read.csv(url, sep=",") # Import dei dati
head(dati,10) # Visualizzazione prime 10 osservazioni
##        city year month sales volume median_price listings months_inventory
## 1  Beaumont 2010     1    83 14.162       163800     1533              9.5
## 2  Beaumont 2010     2   108 17.690       138200     1586             10.0
## 3  Beaumont 2010     3   182 28.701       122400     1689             10.6
## 4  Beaumont 2010     4   200 26.819       123200     1708             10.6
## 5  Beaumont 2010     5   202 28.833       123100     1771             10.9
## 6  Beaumont 2010     6   189 27.219       122800     1803             11.1
## 7  Beaumont 2010     7   164 22.706       124300     1857             11.7
## 8  Beaumont 2010     8   174 25.237       136800     1830             11.6
## 9  Beaumont 2010     9   124 17.233       121100     1829             11.7
## 10 Beaumont 2010    10   150 23.904       138500     1779             11.5
attach(dati) # Estrazione colonne

Commento
Le variabili presenti nel dataset sono:
-city: città.
-year: anno di riferimento.
-month: mese di riferimento.
-sales: numero totale di vendite.
-volume: valore totale delle vendite (espresso in milioni di dollari).
-median_price: prezzo mediano di vendita in dollari.
-listings: numero totale di annunci attivi.
-months_inventory: quantità di tempo necessaria per vendere tutte le inserzioni correnti al ritmo attuale delle vendite (espresso in mesi).

2 - Tipo di variabili

In questa sezione si esplora per linee generali la struttura del dataset e vengono identificati i tipi di variabili. Inoltre viene controllata l’eventuale presenza di valori mancanti.

str(dati) # Struttura dataset
## 'data.frame':    240 obs. of  8 variables:
##  $ city            : chr  "Beaumont" "Beaumont" "Beaumont" "Beaumont" ...
##  $ year            : int  2010 2010 2010 2010 2010 2010 2010 2010 2010 2010 ...
##  $ month           : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ sales           : int  83 108 182 200 202 189 164 174 124 150 ...
##  $ volume          : num  14.2 17.7 28.7 26.8 28.8 ...
##  $ median_price    : num  163800 138200 122400 123200 123100 ...
##  $ listings        : int  1533 1586 1689 1708 1771 1803 1857 1830 1829 1779 ...
##  $ months_inventory: num  9.5 10 10.6 10.6 10.9 11.1 11.7 11.6 11.7 11.5 ...
print(paste("Esito controllo valori mancanti: ",any(is.na(dati)))) # Controllo presenza na
## [1] "Esito controllo valori mancanti:  FALSE"

Commento
Il dataset è interamente bilanciato (valori mancanti assenti) ed è composto da 8 variabili e 240 osservazioni (241 se si considera anche l’intestazione iniziale).
-city: variabile qualitativa nominale che può assumere come valori Beaumont, Bryan-College Station, Tyler, Wichita Falls.
-year: variabile quantitativa discreta o anche qualitativa ordinale che può assumere valori 2010, 2011, 2012, 2013, 2014.
-month: variabile qualitativa nominale ciclica che può assumere valori 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 (1=Gennaio, 2=Febbraio,…).
-sales: variabile quantitativa discreta.
-volume: variabile quantitativa continua.
-median_price: variabile quantitativa continua.
-listings: variabile quantitativa discreta.
-months_inventory: variabile quantitativa continua
Osservazione: Tutte le variabili quantitative continue all’interno del dataset sono espresse su scala di rapporti.

3 -Indici di posizione, variabilità e forma

In questa sezione vengono calcolati e mostrati (quando possibile e sensato) gli indici di posizione, di variabilità e di forma. In alcuni casi vengono mostrate le distribuzioni di frequenza e alcuni indicatori.
Per ragioni di legibilità e confrontabilità sono create tabelle e dataframe. Inoltre, piuttosto che fare un summary generale, vengono calcolati volta per volta i diversi indicatori e quantili scelti (0.1, 0.25, 0.5, 0.75, 0.90).

3.1 - Indici di posizione

# Funzione per calcolare la moda
calcola_moda <- function(x) { 
  freq <- table(x)
  moda <- as.numeric(names(freq)[freq == max(freq)])
  return(moda)
}
# Costruzione distribuzione di frequenze per city, year e month
# city
freq_ass <- table(city)
freq_rel <- table(city)/length(city) # Tabella distribuzione frequenze relative
distr_freq_city<-data.frame(cbind(freq_ass,freq_rel)) # Concatenazione e costruzione dataframe
colnames(distr_freq_city) <- c("Frequenza assoluta", "Frequenza relativa")
distr_freq_city
##                       Frequenza assoluta Frequenza relativa
## Beaumont                              60               0.25
## Bryan-College Station                 60               0.25
## Tyler                                 60               0.25
## Wichita Falls                         60               0.25
# year
freq_ass <- table(year) 
freq_rel <- table(year)/length(year)
distr_freq_year<-data.frame(cbind(freq_ass,freq_rel))
colnames(distr_freq_year) <- c("Frequenza assoluta", "Frequenza relativa")
distr_freq_year
##      Frequenza assoluta Frequenza relativa
## 2010                 48                0.2
## 2011                 48                0.2
## 2012                 48                0.2
## 2013                 48                0.2
## 2014                 48                0.2
# month
freq_ass <- table(month) 
freq_rel <- round(table(month)/length(month),3)
distr_freq_month<-data.frame(cbind(freq_ass,freq_rel))
colnames(distr_freq_month) <- c("Frequenza assoluta", "Frequenza relativa")
distr_freq_month
##    Frequenza assoluta Frequenza relativa
## 1                  20              0.083
## 2                  20              0.083
## 3                  20              0.083
## 4                  20              0.083
## 5                  20              0.083
## 6                  20              0.083
## 7                  20              0.083
## 8                  20              0.083
## 9                  20              0.083
## 10                 20              0.083
## 11                 20              0.083
## 12                 20              0.083

Commento
Dalle distribuzioni di frequenza per i caratteri qualitativi si nota l’equidistribuzione delle osservazioni nelle diverse modalità. Dunque, non c’è una moda precisa poichè, tutte le modalità di city, year e month potrebbero essere qualificate, a rigor di definizione, come moda.

# Creazione dataframe con indici di posizione
# Per le variabili qualitative ha senso calcolare alcuni indicatori solo per year che è ordinale

df_statistiche <- data.frame(
  Variabile = c("city","year", "month", "sales", "volume", "median_price", "listings", "months_inventory"),
  Moda = c( 
    NA, NA, NA, NA, NA, NA, NA, calcola_moda(months_inventory) # non calcolo moda per caratteri continui con molte modalità
  ),
  Media = c(
    NA, mean(year), NA, mean(sales), mean(volume), mean(median_price), mean(listings), mean(months_inventory) # calcolo media anche per year dato che è ordinale
  ),
  Minimo = c(
    NA, min(year), NA, min(sales), min(volume), min(median_price), min(listings), min(months_inventory)
  ),
  Massimo = c(
    NA, max(year), NA, max(sales), max(volume), max(median_price), max(listings), max(months_inventory)
  ),
  Q10 = c(
    NA, quantile(year, 0.1), NA, quantile(sales, 0.1), quantile(volume, 0.1),
    quantile(median_price, 0.1), quantile(listings, 0.1), quantile(months_inventory, 0.1)
  ),
  Q25 = c(
    NA, quantile(year, 0.25), NA, quantile(sales, 0.25), quantile(volume, 0.25), 
    quantile(median_price, 0.25), quantile(listings, 0.25), quantile(months_inventory, 0.25)
  ),
  Mediana = c(
    NA, quantile(year, 0.5), NA, quantile(sales, 0.5), quantile(volume, 0.5),
    quantile(median_price, 0.5), quantile(listings, 0.5), quantile(months_inventory, 0.5)
  ),
  Q75 = c(
    NA, quantile(year, 0.75), NA, quantile(sales, 0.75), quantile(volume, 0.75),
    quantile(median_price, 0.75), quantile(listings, 0.75), quantile(months_inventory, 0.75)
  ),
  Q90 = c(
    NA, quantile(year, 0.9), NA, quantile(sales, 0.9), quantile(volume, 0.9),
    quantile(median_price, 0.9), quantile(listings, 0.9), quantile(months_inventory, 0.9)
  )
)

# arrotondo tutti i valori del dataframe tranne prima colonna e poi la riaggiungo per la visualizzazione
print(cbind(df_statistiche[1],(round(df_statistiche[-1],3)))) 
##          Variabile Moda      Media    Minimo    Massimo       Q10       Q25
## 1             city   NA         NA        NA         NA        NA        NA
## 2             year   NA   2012.000  2010.000   2014.000  2010.000   2011.00
## 3            month   NA         NA        NA         NA        NA        NA
## 4            sales   NA    192.292    79.000    423.000   101.900    127.00
## 5           volume   NA     31.005     8.166     83.547    13.097     17.66
## 6     median_price   NA 132665.417 73800.000 180000.000 99960.000 117300.00
## 7         listings   NA   1738.021   743.000   3296.000   899.900   1026.50
## 8 months_inventory  8.1      9.193     3.400     14.900     6.690      7.80
##      Mediana        Q75        Q90
## 1         NA         NA         NA
## 2   2012.000   2013.000   2014.000
## 3         NA         NA         NA
## 4    175.500    247.000    302.100
## 5     27.062     40.893     53.739
## 6 134500.000 150050.000 158850.000
## 7   1618.500   2056.000   2946.700
## 8      8.950     10.950     12.210

Commento
Le valiabili qualitative (nominali e ordinali) sono tutte bilanciate ed equamente distribuite nelle diverse modalità. L’unica variabile qualitativa per cui ha senso calcolare i quantili e altri indicatori è year, che ovviamente presenta quantili equidistanti.
Media e mediana di tutte le variabili quantitative non coincidono, ma spesso sono molto vicine (se sono tanto o poco vicine tanto da incidere sulla forma della distribuzione sarà valutato nella sezione dedicata agli indici di forma).
Osservando i quantili per le variabili quantitative è possibile formulare prime ipotesi sulla forma delle distribuzioni.
Per esempio, nel caso di sales, il 10° quantile è 101.9, invece il 90° quantile è 302.1 (notevolmente più lontano dalla mediana), quindi il 50% della densità si trova sotto la mediana in un sottoinsieme più piccolo del supporto rispetto al 50% a destra della mediana che è diluito su un sottoinsieme del supporto più ampio.
Sempre riguardo i quantili, talvolta la mediana, talvolta il quantile di ordine 0.75 sono posizionati intorno al centro dei vari supporti. Per esempio, nel caso di volume, il quantile di ordine 0.75, piuttosto che la mediana, è più vicino al centro del supporto (circa 37.6). Ciò significa che il 75% della densità distributiva si trova nella prima metà del supporto. D’altra parte, sempre per volume la distanza dei quantili superiori dalla mediana sembra simile a quella dei quartili inferiori dalla mediana (anche in questo caso, se questa similitudine implica simmetria o meno, sarà valutato con gli indicatori di forma).
Si evidenzia che la distanza tra il valore minimo e massimo sembra particolarmente ampia nel caso di volume, median price e listings.

3.2 - Indici di variabilità

# Funzione per il calcolo del coefficiente di variazione
CV <-function(x){ 
  return( sd(x)/mean(x) )
}

# Funzione per il calcolo dell'indice di Gini
gini.index <- function(x){ 
  ni = table(x)
  fi = ni/length(x)
  fi2 = fi^2
  J = length(table(x))
  gini = 1-sum(fi2)
  gini.norm = gini/((J-1)/J)
  return(gini.norm)
}

# Calcolo separato degli indicatori e aggregazione in dataframe
range_values <- c(
  NA,
  range(year)[2] - range(year)[1],
  NA,
  range(sales)[2] - range(sales)[1],
  round(range(volume)[2] - range(volume)[1],3),
  range(median_price)[2] - range(median_price)[1],
  range(listings)[2] - range(listings)[1],
  range(months_inventory)[2] - range(months_inventory)[1]
)

IQR_values <- c(
  NA,
  IQR(year),
  NA,
  IQR(sales),
  round(IQR(volume),3),
  IQR(median_price),
  IQR(listings),
  IQR(months_inventory)
)

var_values <- c(
  NA,
  NA, # per year ha senso solo il range interquartile e non la varianza
  NA,
  format(var(sales),scientific=FALSE),
  format(var(volume),scientific=FALSE),
  format(var(median_price),scientific=FALSE),
  format(var(listings),scientific=FALSE),
  format(var(months_inventory),scientific=FALSE)
)

sd_values <- c(
  NA,
  NA,
  NA,
  round(sd(sales),3),
  round(sd(volume),3),
  round(sd(median_price),3),
  round(sd(listings),3),
  round(sd(months_inventory),3)
)

CV_values <- c(
  NA,
  NA,
  NA,
  round(CV(sales),3),
  round(CV(volume),3),
  round(CV(median_price),3),
  round(CV(listings),3),
  round(CV(months_inventory),3)
)

# È ovvio che l'indice di Gini per city, year e month è sicuramente 1 dato che assume tutte le modalità e sono equamente distribuite. 
# Per tutte le variabili quantitative l'indice assumerebbe valori vicinissimi a 1 proprio perchè sono variabili continue, dunque sarebbe più informativo e significativo se si calcolasse su una suddivisione in classi. 
gini_values <- c( 
  gini.index(city), 
  gini.index(year),
  gini.index(month),
  NA,
  NA,
  NA,
  NA,
  NA
)

# Creazione del dataframe
df <- data.frame(
  Variabile = c("city","year", "month", "sales", "volume", "median_price", "listings", "months_inventory"),
  Range = range_values,
  IQR = IQR_values,
  Varianza = var_values,
  Deviazione_Std = sd_values,
  CV = CV_values,
  Gini = gini_values
)

print(df)
##          Variabile      Range       IQR  Varianza Deviazione_Std    CV Gini
## 1             city         NA        NA      <NA>             NA    NA    1
## 2             year      4.000     2.000      <NA>             NA    NA    1
## 3            month         NA        NA      <NA>             NA    NA    1
## 4            sales    344.000   120.000    6344.3         79.651 0.414   NA
## 5           volume     75.381    23.233  277.2707         16.651 0.537   NA
## 6     median_price 106200.000 32750.000 513572983      22662.149 0.171   NA
## 7         listings   2553.000  1029.500    566569        752.708 0.433   NA
## 8 months_inventory     11.500     3.150  5.306889          2.304 0.251   NA

Commento
Nel caso delle variabili quantitative, per avere un indice di Gini realmente informativo è necessaria la suddivisione in classi. Invece, nel caso di caratteri qualitativi, è ovviamente 1 perchè le modalità, come già evidenziato, sono equamente distribuite.
I confronti di variabilità sono da fare sul coefficiente di variazione. La variabile che presenta maggiore variabilità è volume (con CV=0.537), però altre variabili con alta variabilità sono listings e sales. La variabile quantitativa che ha il CV più basso è median_price.

Per le variabili qualitative come city, year e month non ha senso calcolare le misure di variabilità.

3.3 - Indici di forma

# Calcolo asimmetria (skewness) e curtosi (kurtosis - 3) per ciascuna variabile. Per city, year e month non ha senso
skewness_values <- c(
  NA,
  NA,
  NA,
  skewness(sales),
  skewness(volume),
  skewness(median_price),
  skewness(listings),
  skewness(months_inventory)
)

kurtosis_values <- c(
  NA,
  NA,
  NA,
  kurtosis(sales) - 3,
  kurtosis(volume) - 3,
  kurtosis(median_price) - 3,
  kurtosis(listings) - 3,
  kurtosis(months_inventory) - 3
)

# Creazione del dataframe
df <- data.frame(
  Variabile = c("city","year", "month", "sales", "volume", "median_price", "listings", "months_inventory"),
  Asimmetria = skewness_values,
  Curtosi = kurtosis_values
)

print(df)
##          Variabile  Asimmetria    Curtosi
## 1             city          NA         NA
## 2             year          NA         NA
## 3            month          NA         NA
## 4            sales  0.71810402 -0.3131764
## 5           volume  0.88474203  0.1769870
## 6     median_price -0.36455288 -0.6229618
## 7         listings  0.64949823 -0.7917900
## 8 months_inventory  0.04097527 -0.1744475

Commento
Osservando i segni di asimmetria e curtosi si traggono le seguenti conclusioni.
L’unica variabile asimmetrica negativa è median_price (quindi la coda sinistra della distribuzione è più allungata), e coerentemente con gli indici di posizione la media è minore della mediana. Invece la variabile più asimmetrica positivamente è volume (anche sales ha un’asimmetria particolarmente intensa) che quindi ha una coda destra più allungata e la densità si concentra sulle modalità inferiori (come ipotizzato nella discussione dei quantili negli indici di posizione). Anche nel caso di asimmetria positiva si riscontra coerenza con gli indici di posizione dato che la media è maggiore della mediana. Inoltre, si trova conferma dell’ipotesi avanzata precendentemente, ovvero per sales, volume, ecc. (tipico dell’asimmetria positiva) la densità si concentra sulle modalità con valori inferiori del supporto.
L’unica variabile leptocurtica è volume (quindi la distribuzione presenta un picco più alto rispetto ad una normale, di conseguenza la densità è più concentrata in corrispondenza di mezzi locali), tutte le altre distribuzioni sono platicurtiche (quindi la densità è maggiormente diluita nelle code delle distribuzioni ed il picco è più basso).
Volume e months_inventory sono leggermente tendenti alla mesocurtosi.

4 - Suddivisione in classi della variabile volume e calcolo dell’indice di Gini

In questa sezione la variabile volume viene suddivisione in classi usando alcuni valori caratteristici per definire l’estremo superiore e inferiore di ogni classe (minimo, mediana, media, massimo). Questa prima proposta è poi confrontata con la suddivisione in classi secondo sottointervalli (più o meno equidistanti) del supporto.

# PRIMA PROPOSTA
# Creazione delle classi di volume in base ai valori caratteristici
dati$volume_class <- cut(volume, breaks = c(8.166, 27.063, 31.005, 85.547)) 
# Frequenze per le classi di volume (per costruzione grafico)
freq_volume <- table(dati$volume_class)

# Creazione del grafico delle frequenze per le classi di volume
ggplot(as.data.frame(freq_volume), aes(x = Var1, y = Freq/length(dati$volume_class))) + # var1 e freq vengono create di default con il dataframe
  geom_col(fill = "steelblue") +
  labs(title = "Distribuzione di frequenza per classi di 'volume'",
       x = "Classi di volume", y = "Frequenza") +
  theme_minimal()

# Calcolo dell'indice di Gini per la distribuzione di frequenza delle classi di volume
gini_volume <- gini.index(dati$volume_class) # Osservazione: se si passa freq_volume allora viene applicato 2 volte table (1 tramite funzione e 1 per def di freq_volume)
gini_volume
## [1] 0.8561719
# SECONDA PROPOSTA
dati$volume_class <- cut(volume, breaks = c(seq(1,100,6))) # Suddivisione del supporto in sottointervalli

freq_volume <- table(dati$volume_class)

ggplot(as.data.frame(freq_volume), aes(x = Var1, y = Freq/length(dati$volume_class))) +
  geom_col(fill = "steelblue") +
  labs(title = "Distribuzione di frequenza per classi di 'volume'",
       x = "Classi di volume", y = "Frequenza") +
  theme_minimal()

gini_volume <- gini.index(dati$volume_class)
gini_volume
## [1] 0.9368889

Commento La conclusione più evidente è che, inevitabilmente la costruzione delle classi impatta sull’indice di Gini quindi è preferibile scegliere le classi in base a valori del supporto.
Nella prima proposta, la prima e l’ultima classe sono le più numerose, invece la classe centrale è notevolmente sbilanciata verso il basso (scarsissima numerosità, inferiore al 10%). In questa distribuzione in classi l’indice di Gini è molto elevato, oltre 0.85, ed indica una quasi equiripartizione (in media).
Nella seconda proposta è scelta una suddivisione in classi per sottointervalli del supporto. Da ciò si evidenzia una distribuzione coerente con gli indici di forma e di posizione.

Approfondimento
La distribuzione empirica sui dati grezzi di volume è coerente con la seconda proposta di suddivisione in classi e con gli indicatori calcolati precedentemente. Infatti, la forma che segue la distribuzione in classi (seconda proposta) è molto simile alla stima del kernel (approssimativa) mostrata di seguito.

hist(volume, probability = TRUE, main = "Istogramma con Densità Kernel", xlab = "Volume", col = "lightblue", border = "black")
lines(density(volume), col = "red", lwd = 2)

# La distribuzione è asimmetrica positiva (come ci si aspettava)

Osservazione
La variabile volume ha supporto 0,+inf. La stima della densità con kernel gaussiano (di default in density) potrebbe essere distorto in prossimità dell’estremo inferiore del supporto. Per dati unidimensionali con supporto limitato sarebbe meglio applicare uno stimatore BGK (Botev,2010). Ad ogni modo la reale distribuzione della variabile, in questo caso, è da indagare tra Dagum, log-Normale, F, Chi-square (o loro simili).

5 - Calcolo dell’indice di Gini per la variabile city

In questa sezione viene calcolato l’indice di Gini per la variabile city, che ci si aspetta sia 1 perchè le unità statistiche si distribuiscono equamente nelle diverse modalità.

# Dimostrazione
gini_city <- gini.index(city)
gini_city
## [1] 1


6 - Calcolo delle probabilità

In questa sezione viene calcolata la probabilità che si verifichino alcuni eventi. Ovvero, la probabilità che, presa una riga casualmente: a) essa riporti la città “Beaumont”, b) che essa riporti il mese di Luglio, c) che essa riporti il mese di dicembre 2012.

totale <- nrow(dati)
prob_beaumont <- sum(city == "Beaumont") / totale
prob_luglio <- sum(month == "7") / totale
prob_dec2012 <- sum(month == "12" & year == 2012) / totale

print(paste("Probabilità di una osservazione casuale con città Beaumont: ", prob_beaumont))
## [1] "Probabilità di una osservazione casuale con città Beaumont:  0.25"
print(paste("Probabilità di una osservazione casuale con mese Luglio: ", round(prob_luglio,3)))
## [1] "Probabilità di una osservazione casuale con mese Luglio:  0.083"
print(paste("Probabilità di una osservazione casuale con mese Dicembre e anno 2012: ", round(prob_dec2012,3))) # 1,67%
## [1] "Probabilità di una osservazione casuale con mese Dicembre e anno 2012:  0.017"


7 - Creazione di una colonna che rappresenti il prezzo medio

In questa sezione viene creata una colonna, aggiunta poi al dataset, che riporta il prezzo medio. Per la creazione della nuova variabile vengono coinvolte le volonne volume e sales.

dati$mean_price <- (dati$volume * 1e6) / dati$sales # moltiplicazione per un milione dato che nel dataset è scalata e l'ordine di grandezza si riduce dividendo per sales
head(dati[c("volume", "sales", "mean_price")])
##   volume sales mean_price
## 1 14.162    83   170626.5
## 2 17.690   108   163796.3
## 3 28.701   182   157697.8
## 4 26.819   200   134095.0
## 5 28.833   202   142737.6
## 6 27.219   189   144015.9


8 - Creazione di una colonna che rappresenti l’efficacia

In questa sezione viene creata una colonna, similmente alla sezione precedente, che rappresenti l’efficacia degli annunci di vendita. Vengono coinvolte le variabili sales e listings così da indicare quante vendite sono state generate da ogni annuncio attivo.

dati$efficacia <- dati$sales / dati$listings
head(dati[, c("sales", "listings", "efficacia")])
##   sales listings  efficacia
## 1    83     1533 0.05414220
## 2   108     1586 0.06809584
## 3   182     1689 0.10775607
## 4   200     1708 0.11709602
## 5   202     1771 0.11405985
## 6   189     1803 0.10482529


9 - Analisi condizionata

In questa sezione viene presentata un’analisi condizionata, utilizzando dplyr, delle variabili sales, volume e listings rispetto a city, year e month.
Per ciascun box presentato, in ogni grafico le diverse modalità hanno un colore unicamente associato. Quindi, le legende saranno omesse altrimenti causerebbero troppa ridondanza.

# Condizionamento a city
sintesi <- dati %>%
  group_by(city) %>%
  summarise(
    media_sales = mean(sales),
    sd_sales = sd(sales),
    media_volume = mean(volume),
    sd_volume = sd(volume),
    media_listings = mean(listings),
    sd_listings = sd(listings)
  )
sintesi
## # A tibble: 4 × 7
##   city    media_sales sd_sales media_volume sd_volume media_listings sd_listings
##   <chr>         <dbl>    <dbl>        <dbl>     <dbl>          <dbl>       <dbl>
## 1 Beaumo…        177.     41.5         26.1      6.97          1679.        91.1
## 2 Bryan-…        206.     85.0         38.2     17.2           1458.       253. 
## 3 Tyler          270.     62.0         45.8     13.1           2905.       227. 
## 4 Wichit…        116.     22.2         13.9      3.24           910.        73.8
g1<-ggplot(sintesi, aes(x = city, y = media_sales, fill = city)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin = media_sales - sd_sales, ymax = media_sales + sd_sales), # bande di errore con ds
                position = position_dodge(width = 0.9), width = 0.25) +
  labs(title = "Media numero vendite (sales) per città con Deviazione Standard", 
       x = "City", y = "Sales") +
  theme_minimal()+
  theme(legend.position = "none")

g2<-ggplot(sintesi, aes(x = city, y = media_volume, fill = city)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin = media_volume - sd_volume, ymax = media_volume + sd_volume), 
                position = position_dodge(width = 0.9), width = 0.25) +
  labs(title = "Media valore delle vendite (volume) per città con Deviazione Standard", 
       x = "City", y = "Volume") +
  theme_minimal()+
  theme(legend.position = "none")

g3<-ggplot(sintesi, aes(x = city, y = media_listings, fill = city)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin = media_listings - sd_listings, ymax = media_listings + sd_listings), 
                position = position_dodge(width = 0.9), width = 0.25) +
  labs(title = "Media annunci attivi (listings) per città con Deviazione Standard", 
       x = "City", y = "Listings") +
  theme_minimal()+
  theme(legend.position = "none")

grid.arrange(g1, g2, g3, ncol = 1) # Output grafico combinato

Commento
Tyler ha la media delle vendite più alta, seguita da Bryan-College Station. Invece, Wichita Falls ha medie di vendita significativamente più basse. Tyler e Bryan-College Station mostrano anche una maggiore variabilità (deviazione standard) nelle vendite.
Tyler ha anche il volume medio più alto, simile a quanto osservato per le vendite, e Bryan-College Station ha un volume medio leggermente inferiore. Beaumont e Wichita Falls hanno volumi medi molto più bassi. La variabilità del volume segue un andamento simile a quella delle vendite.
Per quanti riguarda il numero medio di annunci, Tyler ha il numero medio di annunci più alto, superando di gran lunga tutte le altre città. In Tyler e Bryan-College Station si osserva una maggiore variabilità nel numero di annunci.
Sembra che Tyler e Bryan-College Station sono le città con le performance migliori in termini di vendite, volume e numero di annunci. 

# Condizionamento a year
sintesi2 <- dati %>%
  group_by(year) %>%
  summarise(
    media_sales = mean(sales),
    sd_sales = sd(sales),
    media_volume = mean(volume),
    sd_volume = sd(volume),
    media_listings = mean(listings),
    sd_listings = sd(listings)
  )
sintesi2
## # A tibble: 5 × 7
##    year media_sales sd_sales media_volume sd_volume media_listings sd_listings
##   <int>       <dbl>    <dbl>        <dbl>     <dbl>          <dbl>       <dbl>
## 1  2010        169.     60.5         25.7      10.8          1826         785.
## 2  2011        164.     63.9         25.2      12.2          1850.        780.
## 3  2012        186.     70.9         29.3      14.5          1777.        738.
## 4  2013        212.     84.0         35.2      17.9          1678.        744.
## 5  2014        231.     95.5         39.8      21.2          1560.        707.
colori <- c("2010" = "#F0F022", "2011" = "#66FF66", "2012" = "#6666FF", "2013" = "#C433FF", "2014" = "#33C4FF")
g1<-ggplot(sintesi2, aes(x = year, y = media_sales, fill = as.factor(year))) + # esplicito as.factor perchè altrimenti i colori seguono una gradazione che non è visivamente ottimale
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin = media_sales - sd_sales, ymax = media_sales + sd_sales), 
                position = position_dodge(width = 0.9), width = 0.25) +
  scale_fill_manual(values = colori) +
  labs(title = "Media numero vendite (sales) per anno con Deviazione Standard", 
       x = "Year", y = "Sales") +
  theme_minimal()+
  theme(legend.position = "none")

g2<-ggplot(sintesi2, aes(x = year, y = media_volume, fill = as.factor(year))) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin = media_volume - sd_volume, ymax = media_volume + sd_volume), 
                position = position_dodge(width = 0.9), width = 0.25) +
  scale_fill_manual(values = colori) +
  labs(title = "Media valore delle vendite (volume) per anno con Deviazione Standard", 
       x = "Year", y = "Volume") +
  theme_minimal()+
  theme(legend.position = "none")

g3<-ggplot(sintesi2, aes(x = year, y = media_listings, fill = as.factor(year))) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin = media_listings - sd_listings, ymax = media_listings + sd_listings), 
                position = position_dodge(width = 0.9), width = 0.25) +
  scale_fill_manual(values = colori) +
  labs(title = "Media annunci attivi (listings) per anno con Deviazione Standard", 
       x = "Year", y = "Listings") +
  theme_minimal()+
  theme(legend.position = "none")

grid.arrange(g1, g2, g3, ncol = 1)

Commento
Il valore delle vendite ed il numero di vendite presentano un trend crescente dal 2010 al 2014, con un’associato aumento della variabilità. Il numero di annunci attivo, invece, presenta un trend decrescente negli stessi anni ed una leggera flessione della variabilità.

# Condizionamento a month
sintesi3 <- dati %>%
  group_by(month) %>%
  summarise(
    media_sales = mean(sales),
    sd_sales = sd(sales),
    media_volume = mean(volume),
    sd_volume = sd(volume),
    media_listings = mean(listings),
    sd_listings = sd(listings)
  )
head(sintesi3,12)
## # A tibble: 12 × 7
##    month media_sales sd_sales media_volume sd_volume media_listings sd_listings
##    <int>       <dbl>    <dbl>        <dbl>     <dbl>          <dbl>       <dbl>
##  1     1        127.     43.4         19.0      8.37          1647.        705.
##  2     2        141.     51.1         21.7     10.1           1692.        711.
##  3     3        189.     59.2         29.4     12.0           1757.        727.
##  4     4        212.     65.4         33.3     14.5           1826.        770.
##  5     5        239.     83.1         39.7     19.0           1824.        790.
##  6     6        244.     95.0         41.3     21.1           1833.        812.
##  7     7        236.     96.3         39.1     21.4           1821.        827.
##  8     8        231.     79.2         38.0     18.0           1786.        816.
##  9     9        182.     72.5         29.6     15.2           1749.        803.
## 10    10        180.     75.0         29.1     15.1           1710.        779.
## 11    11        157.     55.5         24.8     11.2           1653.        741.
## 12    12        169.     60.7         27.1     12.6           1558.        693.
colori <- c(
  "1" = "red",   
  "2" = "#4682B4",   
  "3" = "#32CD32",  
  "4" = "#FFD700",  
  "5" = "#8A2BE2",  
  "6" = "#FF1493",  
  "7" = "#00FFFF",   
  "8" = "grey",   
  "9" = "#800080",   
  "10" = "#00FF00",  
  "11" = "orange",  
  "12" = "#0000FF"
)
g1<-ggplot(sintesi3, aes(x = month, y = media_sales, fill = as.factor(month))) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin = media_sales - sd_sales, ymax = media_sales + sd_sales), 
                position = position_dodge(width = 0.9), width = 0.25) +
  scale_fill_manual(values = colori) +
  labs(title = "Media numero vendite (sales) per mese con Deviazione Standard", 
       x = "Month", y = "Sales") +
  scale_x_continuous(breaks = 1:12) +
  theme_minimal()+
  theme(legend.position = "none")

g2<-ggplot(sintesi3, aes(x = month, y = media_volume, fill = as.factor(month))) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin = media_volume - sd_volume, ymax = media_volume + sd_volume), 
                position = position_dodge(width = 0.9), width = 0.25) +
  scale_fill_manual(values = colori) +
  labs(title = "Media valore delle vendite (volume) per mese con Deviazione Standard", 
       x = "Month", y = "Volume") +
  scale_x_continuous(breaks = 1:12) +
  theme_minimal()+
  theme(legend.position = "none")

g3<-ggplot(sintesi3, aes(x = month, y = media_listings, fill = as.factor(month))) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin = media_listings - sd_listings, ymax = media_listings + sd_listings), 
                position = position_dodge(width = 0.9), width = 0.25) +
  scale_fill_manual(values = colori) +
  labs(title = "Media annunci attivi (listings) per mese con Deviazione Standard", 
       x = "Month", y = "Listings") +
  scale_x_continuous(breaks = 1:12) +
  theme_minimal()+
theme(legend.position = "none")

grid.arrange(g1, g2, g3, ncol = 1)

Commento
L’andamento della media del numero di vendite e del valore totale delle vendite nei diversi mesi dell’anno è approssimativamente di tipo sinusoidale. Infatti il trend è crescente fino a Giugno (month=6) e decrescente fino a Novembre (month=11) con un finale rialzo in Dicembre (month=12).
Per quanto riguarda l’andamento della media degli annunci attivi, l’andamento è leggermente campanulare (quasi costante) con un picco nel mese di Giugno.
L’andamento della variabilità segue l’andamento della media quindi nei mesi in cui si registrano medie alte, anche la variabilità è elevata.

Approfondimento
Tralasciando varie questioni statistiche a livello di microdati (come ad esempio il matching e il tracciamento annuncio-realizzazione della vendita), si può avanzare l’ipotesi che gli annunci attivi generino con maggiore efficacia volume (ovvero valore) di vendite e numero di vendite solo per una parte dell’anno.
Con specifici modelli di durata e regressivi potrebbe essere utile comprendere quanto tempo passa tra l’inserimento di un annuncio e la realizzazione in vendita effettiva. Si può avanzare l’ipotesi più precisa che da Dicembre a Maggio si verifichi una espansione del mercato immobiliare (mantenendo un eccesso di offerta), invece nella seconda metà dell’anno una leggera contrazione.
Essendo movimenti di mercato di breve periodo e abbastanza “veloci”, per valutare i meccanismi di causa-effetto dei periodi di espansione si necessiterebbe di set di dati più ampi, ma soprattutto di più dettagliati in termini microeconomici e di microdati.

10 - Grafici con ggplot2

In questa sezione vengono mostrate ulteriori visualizzazioni grafiche tramite l’utilizzo di ggplot2.

10.1 - Boxplots

# Boxplot del prezzo mediano per città
dati$city <- as.factor(dati$city)
ggplot(dati, aes(x = city, y = median_price)) +
  geom_boxplot(fill = "lightblue", color = "darkblue") +
  labs(title = "Boxplot del prezzo mediano per città",
       x = "Città", y = "Prezzo Mediano") +
  theme_minimal()

Commento
La distribuzione del prezzo mediano in Bryan-College Station presenta la mediana più alta, invece Wichita Falls si osserva una mediana più bassa.
Per quanto riguarda la variabilità dei prezzi mediani, la maggiore estensione del boxplot (quindi il range che i valori possono assumere) si osserva in Wichita Falls e Beaumont.
I punti al di fuori dei “baffi” rappresentano valori anomali e si nota, quindi, la presenza di valori anomali (in questo caso più elevati) sia in Beaumont che in Wichita Falls e Bryan-College Station. La presenza di outliers indica che in alcune città esistono case con prezzi che si discostano significativamente dalla media.
Beaumont e Tyler si collocano in una fascia intermedia di prezzi mediani, con Tyler che mostra una leggera tendenza verso prezzi più alti.

# Boxplot sales per città affiancato nei diversi anni
ggplot(dati, aes(x = city, y = sales, fill = as.factor(year))) +
  geom_boxplot() +
  labs(title = "Boxplot del numero di vendite per città e anno",
       x = "Città", y = "Numero di vendite") +
  theme_minimal() +
  scale_fill_brewer(palette = "Set3", name = "Anno")

Commento
In tutte le città si osserva un aumento della mediana del numero di vendite nel corso dei diversi anni considerati. L’unica eccezione si registra per Wichita Falls, dove tale valore è oscillato nel corso degli anni ma mediamente è rimasto stabile.
Un fatto particolarmente rilevante è che in Bryan-College Station e Tyler la variabilità del numero di vendite è notevolmente elevato, in contrapposizione a Wichita Falls dove la variabilità è molto bassa.
L’unico valore anomalo si osserva per il numero di vendite in Beaumont nel 2012.

10.2 - Grafici a barre sovrapposte

# Grafico a barre sovrapposte per il totale delle vendite nei vari mesi e città
vendite_mensili <- dati %>%
  group_by(month, city, year) %>%
  summarise(tot_sales = sum(sales)) %>%
  ungroup() 
## `summarise()` has grouped output by 'month', 'city'. You can override using the
## `.groups` argument.
colori_citta <- c(
  "Beaumont" = "#1b9e77",
  "Bryan-College Station" = "#d95f02",
  "Tyler" = "#7570b3",
  "Wichita Falls" = "#e7298a"
)

g1<-ggplot(vendite_mensili, aes(x = month, y = tot_sales, fill = city)) +
  geom_col(position = "stack") +
  scale_fill_manual(values = colori_citta) +
  scale_x_continuous(breaks = 1:12) +
  labs(title = "Vendite mensili per città",
       x = "Mese", y = "Totale vendite") +
  theme_minimal()

g2<-ggplot(vendite_mensili, aes(x = month, y = tot_sales, fill = city)) +
  geom_col(position = "fill") + # grafico normalizzato
  scale_fill_manual(values = colori_citta) +
  scale_x_continuous(breaks = 1:12) +
  labs(title = "Vendite mensili normalizzate per città",
       x = "Mese", y = "Proporzione vendite") +
  theme_minimal()

grid.arrange(g1, g2, ncol = 1)

Commento
In termini assoluti, in quasi tutte le città il numero di vendite aumenta nella prima metà dell’anno, solamente in Wichita Falls restano più o meno costanti nel corso dei mesi. Come evidenziato precedentemente, il numero di vendite più grandi si concentra nei mesi tra Maggio e Agosto.
In termini relativi si nota chiaramente che la città che contribuisce principalmete al numero di vendite totali è Tyler, seguita da Bryan-College Station (con numero di vendite un pò variabile), Beaumont (con numero di vendite abbastanza consistente e costante nei diversi mesi) e infine Wichita Falls.

# Grafico a barre sovrapposte per vendite mensili per città, con differenziazione per anno
colori_citta <- c(
  "Beaumont" = "#1b9e77",
  "Bryan-College Station" = "#d95f02",
  "Tyler" = "#7570b3",
  "Wichita Falls" = "#e7298a"
)

g1 <- ggplot(vendite_mensili, aes(x = month, y = tot_sales, fill = city)) +
  geom_col(position = "stack") +
  scale_fill_manual(values = colori_citta) +
  scale_x_continuous(breaks = 1:12) +
  labs(
    title = "Vendite mensili per città e anni",
    x = "Mese",
    y = "Totale vendite",
    fill = "Città"
  ) +
  theme_minimal() +
  facet_wrap(~ as.factor(year))
g1

# grafico sales-month-city-year normalizzato
g2 <- ggplot(vendite_mensili, aes(x = month, y = tot_sales, fill = city)) +
  geom_col(position = "fill") +
  scale_fill_manual(values = colori_citta) +
  scale_x_continuous(breaks = 1:12) +
  labs(
    title = "Vendite mensili per città e anni",
    x = "Mese",
    y = "Totale vendite",
    fill = "Città"
  ) +
  theme_minimal() +
  facet_wrap(~ as.factor(year))
g2

Commento
Quanto detto per il grafico precedente è valido parzialmente anche se si considerano i diversi anni. Infatti, Wichita Falls presenta un trend di vendite mensili relativamente basso ma costante con una lieve tendenza al decremento, in Beaumont si osserva un numero si vendite abbastanza costistente e costante nel corso sia dei mesi che degli anni. Sei si considera Bryan-College Station il numero di vendite sul totale è circa il 25% e, pur subendo fluttuazioni di breve periodo, il trend è abbastanza costante. Infine, come detto precedentemente, in Tyler si concentra la più elevata percentuale di numero di vendire e tale qualifica si conferma in tutti gli anni considerati.

10.3 - Line chart

In questa sottosezione viene mostrato il line chart delle vendite (sales) nel tempo e per città congiuntamente. Per fare ciò viene implementata l’interazione tra mese e anno così da prendere una serie storica più lunga. Inoltre, è usata l’interazione per catturare più informazione, in accordo con le statistiche sulla variabilità dato che year aveva bassa variabilità (poche informazioni) e, invece, month era caratterizzata da un’elevata variabilità.

dati$data <- as.Date(paste(dati$year, sprintf("%02d", dati$month), "01", sep = "-"))

# Creazione del grafici
g1<-ggplot(dati, aes(x = data, y = sales, color = city)) +
  geom_line() +
  scale_x_date(date_breaks = "1 year", date_labels = "1/1/%Y")+
  labs(title = "Andamento delle vendite nel tempo per città (line chart)",
       x = "Data", y = "Numero di vendite") +
  theme_minimal()

g2<-ggplot(dati, aes(x = data, y = sales, color = city)) +
  geom_jitter() + # risalto dei punti
  geom_line() +
  scale_x_date(date_breaks = "1 year", date_labels = "1/1/%Y")+
  labs(title = "Andamento delle vendite nel tempo per città",
       x = "Data", y = "Numero di vendite") +
  theme_minimal()

grid.arrange(g1, g2, ncol = 1)

Commento
Il line chart riflette esattamente quanto già evidenziato in precedenza per quanto riguarda i confronti tra i numeri delle vendite in ogni città nel corso dei mesi e anni. Vengono confermati i diversi trend e i confronti relativi. In aggiunta, da questo grafico si vedono in modo marcato le componenti stagionali e cicliche per il numero delle vendite in Beaumont e Bryan-College Station, infatti ogni 12 mesi sembrano avere uno stesso andamento seppur con un trend crescente. Ciò è meno evidente per Tyler e Wichita Falls che hanno rispettivamente un trend crescente e costante, infatti non sembrano emergere pattern di andamento nel corso dei mesi e degli anni.