House sales

# htmltools::includeHTML("https://kepler.gl/demo/map?mapUrl=https://dl.dropboxusercontent.com/s/nq8hdpw5h851axg/keplergl_vh7lis1e.json")

#  !pip install keplergl
#  !pip install matplot
#  
#  $ jupyter nbextension install --py --sys-prefix keplergl # can be skipped for notebook 5.3 and above
# $ jupyter nbextension enable --py --sys-prefix keplergl # can be skipped for notebook 5.3 and above
 
# from keplergl import KeplerGl
# with open('data/kc_house_data_numeric.csv', 'r') as f:
#     csvData = f.read()
# map_1 = KeplerGl(height=400)
# map_1.add_data(data=csvData, name='data_2')
# map_1
# 
# import matplotlib.pyplot as plt
# plt.show()

1 Introduzione

1.1 Consegna

Esaminiamo i dati kc_house_data relativi alle transazioni immobiliari per scoprire come il prezzo sia influenzato dalla metratura delle abitazioni, considerando eventualmente anche altre caratteristiche degli immobili. L’obiettivo di questa analisi è quindi quello di individuare eventuali connessioni tra le diverse variabili tramite la creazione di un modello di regressione.

Nel grafico (immagine seguente) sono presenti degli esagoni con una scala di colori dal blu al rosso, se la media delle case all’interno dell’esagono è bassa sarà più tendente al blu, viceversa, di colore rosso. L’altezza delle barre rappresenta la densità di case nell’area coperta dall’esagono. La mappa interattiva è disponibile cliccando qui.

knitr::include_graphics("image/kepler.png")

1.2 Importazione dati

Per cominciare, dobbiamo importare i dati \(kc_house_data\) e rappresentarli in modo appropriato. Utilizzando i pacchetti giusti e le funzioni adeguate, possiamo caricare i dati in un dataframe e visualizzarli in modo da avere una prima idea delle informazioni a disposizione e di come siano organizzati. Successivamente, possiamo utilizzare diverse tecniche di rappresentazione dei dati, come grafici e tabelle, per esplorare ulteriormente i dati e comprendere meglio la relazione tra prezzo e metratura delle abitazioni.

# dati presi da kaggle
df = "data/kc_house_data.csv" %>% 
  read.csv() %>% 
  tibble()

# creo nuovi file csv per la mappa Kepler.gl
"data/kc_house_data.csv" %>% 
  read.csv() %>% 
  tibble() %>% 
  mutate_at("price", as.numeric) %>% 
  write.csv(file = "data/kc_house_data_numeric.csv", quote = T)

"data/kc_house_data.csv" %>% 
  read.csv() %>% 
  tibble() %>% 
  mutate_at("price", as.numeric) %>% 
  slice(which(price > quantile(df$price,.999) %>% as.numeric)) %>% 
  write.csv(file = "data/lusso999.csv", quote = T)

df

pal = with(df, colorFactor(brewer.pal(10,"RdYlGn"), price))
dfPopup = df %>% 
  mutate(popup_info = paste("Prezzo della casa: ", price, " $", "</br>",
                            "Superficie del soggiorno: ", sqft_living, "</br>", 
                            "Superficie del seminterrato: ", sqft_basement, "</br>", 
                            "Superficie del lotto: ", sqft_lot, "</br>",
                            "Numero di piani: ", floors, "</br>", 
                            "Numero di bagni: ", bathrooms, "</br>",
                            "Numero di camere: ", bedrooms, "</br>",
                            "Vista sul mare: ", waterfront, "</br>",
                            "Data di vendita dell'immobile: ", date, "</br>",
                            "Numero di visite: ", view, "</br>", 
                            "Valutazione: ", grade, " ", condition, "</br>",
                            "Coordinate: ", lat, " ", long, "</br>",
                            "Anno di costruzione: ", yr_built, "</br>", 
                            "Anno del restauro: ", yr_renovated, "</br>"))
leaflet() %>% 
  addTiles() %>% 
  addCircleMarkers(data = dfPopup,
                   lat = ~ lat,
                   lng = ~ long,
                   radius = ~ 1,
                   color = ~ pal(price),
                   popup = ~ popup_info)

1.3 Descrizione dataset

Il dataset che ci troveremo ad analizzare è composto da 21597 esempi di case vendute nella contea di King, nello stato di Washington (USA), tra Maggio 2014 e Maggio 2015. Il dataset include informazioni su 21 variabili differenti, come prezzo, metratura, numero di camere da letto, anno di costruzione, e così via. Si tenga presente che questi dati coprono solo parzialmente la città di Seattle. L’obiettivo di questa analisi sarà quello di esplorare le relazioni tra le variabili presenti nel dataset e capire come influiscano sul prezzo delle case, utilizzando una vasta gamma di tecniche statistiche e di visualizzazione.

Breve descrizione delle 21 variabili (Data type):

• id: numero identificativo associato all’i-esima casa venduta (Numeric);

• date: la data in cui la casa i-esima è stata venduta (String);

• price: la variabile risposta (Numeric);

• bedrooms: n° di camere da letto per casa (Numeric);

• bathrooms: n° di bagni per camera (Numeric);

• sqftliving: n° di piedi quadrati del soggiorno (Numeric);

• sqftlot: n° di piedi quadrati del lotto (Numeric);

• floors: n° di piani (Numeric);

• waterfront: se la casa ha vista sul mare/fiume/lago (Numeric);

• view: n° di visite effettuate (Numeric);

• condition: valutazione complessiva delle condizioni della casa, dove 1 indica usurata e 5 indica eccellente. Per maggiori informazioni: (http://info.kingcounty.gov/assessor/esales/Glossary.aspx?type=r#g) (Numeric);

• grade: grado complessivo assegnato all’abitazione, basato sul King County grading system: 1 povera, 13 eccellent. (Numeric);

• sqftabove: piedi quadri della casa escluso il seminterrato (Numeric);

• sqftbasement: piedri quadri del seminterrato (Numeric);

• yrbuilt: anno di costruzione (Numeric);

• yrrenovated: anno di restauro della casa (Numeric);

• zipcode: codice zip (Numeric);

• lat: Latitude coordinate (Numeric);

• long: Longitude coordinate (Numeric);

• sqftliving15: zona giorno nel 2015(dato implicato in caso di restauri dell’area del lotto), si noti che questo dato potrebbe non aver avuto effetto sui piedri quadri del lotto (Numeric);

• sqftlot15: piedri quadri del lotto nel 2015 (dato implicato in caso di restauri dell’area del lotto) (Numeric).

La presenza della variabile grade nel dataset, basata su un sistema di classificazione specifico per la contea di King, suggerisce che i dati potrebbero essere stati raccolti da una fonte ufficiale come ad esempio l’ufficio delle tasse del contea. Inoltre, la presenza della variabile view, che indica il numero di visite effettuate ad una casa, suggerisce che i dati potrebbero provenire da un’agenzia immobiliare o un’altra fonte simile. Anche se non ci sono ragioni per mettere in discussione l’accuratezza generale dei dati, è sempre opportuno considerare con cautela l’applicabilità dei modelli basati su questo dataset a casi più generali, tenendo in considerazione la possibile limitatezza geografica e temporale dei dati e l’eventuale influenza di variabili non incluse nel dataset.

2 Analisi preliminare dei dati

2.1 Variabile risposta: Price

Per la scelta della variabile risposta il percorso è imposto, va utilizzata la variabile price, ma vanno considerati i seguenti accorgimenti per l’assunzione sul modello parametrico da utilizzare: - Per la costruzione di un modello di regressione lineare, condizione necessaria è che la variabile risposta sia una variabile quantitativa continua con distribuzione normale; - Per un modello glm gamma la variabile risposta deve essere quantitativa continua con supporto in R+ e deve assumere la forma funzionale per il glm gamma.

Per arrivare ai seguenti risultati ci arriveremo tramite delle trasformazioni.

ggplot(df, aes(x = price, y = after_stat(density))) +
  geom_histogram(color = "orchid",
               fill = "orchid",
               alpha = 0.5,
               bins = 100) + 
  geom_density(color = "paleturquoise",
               fill = "paleturquoise",
               alpha = 0.3, # densità del colore 
               kernel = "gaussian",
               adjust = 1) +
  geom_line(aes(y = dnorm(price, mean(price), sd(price))),
            color = "aquamarine") +
  geom_line(aes(y = dgamma(price,
                           fitdistrplus::fitdist(df$price, "gamma", method = "mme")$estimate[1],
                           fitdistrplus::fitdist(df$price, "gamma", method = "mme")$estimate[2])),
            color = "blueviolet") +
  xlim(c(0,3E6)) 

2.1.1 Traformazione della variabili risposta

I dati della variabile price non si distribuiscono perfettamente:

shapiro.test(df$price[sample(1:nrow(df), 5000)])

## 
##  Shapiro-Wilk normality test
## 
## data:  df$price[sample(1:nrow(df), 5000)]
## W = 0.72161, p-value < 2.2e-16

jarque.bera.test(df$price)

## 
##  Jarque Bera Test
## 
## data:  df$price
## X-squared = 1131390, df = 2, p-value < 2.2e-16

tibble(kurtosis = kurtosis(df$price), skewness = skewness(df$price))

Nei due grafici seguenti si può osservare la distribuzione dei quantili teorici di una distribuzione normale con i quantili teorici osservati. I quantili empirici del logaritmo di price si distribuiscono maggiormente sulla retta che divide i quadranti del grafico. Si può constatare che, nonostante la trasformazione logaritmica, la distribuzione continua a possedere un’assimetria con una coda di destra più pesante.

df %>% 
  dplyr::select(price) %>% 
  add_column("type" = "price") %>% 
  add_row(price = df$price %>% log, type = "log(price)") %>% 
  ggplot(aes(sample = price)) +
  geom_qq(color = "orchid",
          alpha = 0.8) +
  geom_qq_line(color = "aquamarine",
          alpha = .8,
          linewidth = 1) +
  facet_wrap(vars(type), 
             scale = "free_y")

Provo a effetuare una trasformazione logaritmica per provare a vedere se con questa trasformazione i miei dati si distribuiscono meglio secondo una normale o una distribuzione gamma.

df %>% 
  mutate_at("price",log) %>% 
  ggplot(aes(x = price, y = after_stat(density))) +
  geom_histogram(color = "orchid",
               fill = "orchid",
               alpha = 0.5,
               bins = 100) + 
  geom_density(color = "paleturquoise",
               fill = "paleturquoise",
               alpha = 0.3,
               kernel = "gaussian",
               adjust = 1) +
  geom_line(aes(y = dnorm(price, mean(price), sd(price))),
            color = "aquamarine") +
  geom_line(aes(y = dgamma(price,
                           fitdistrplus::fitdist(df$price %>% log, "gamma", method = "mme")$estimate[1],
                           fitdistrplus::fitdist(df$price %>% log, "gamma", method = "mme")$estimate[2])),
            color = "blueviolet")

Mentre prima la distribuzione gamma si adattava molto meglio ai dati della normale, ora le due distribuzioni appaiono molto simili.

df = df %>% mutate("lprice" = log(price))
paste("Skewness:", skewness(df$lprice),
      "Kurtosis: ", kurtosis(df$lprice))

## [1] "Skewness: 0.430944309171445 Kurtosis:  0.690268972854811"

shapiro.test(df$price[sample(1:nrow(df), 5000)])

## 
##  Shapiro-Wilk normality test
## 
## data:  df$price[sample(1:nrow(df), 5000)]
## W = 0.72231, p-value < 2.2e-16

shapiro.test(df$lprice[sample(1:nrow(df), 5000)])

## 
##  Shapiro-Wilk normality test
## 
## data:  df$lprice[sample(1:nrow(df), 5000)]
## W = 0.99297, p-value = 5.304e-15

2.2 Altre variabili

2.2.1 Zipcode

Osservando come sono formati i dati si nota che zipcode è considerata come una variabile numerica quando, invece, una sua classificazione adeguata sarebbe come variabile qualitativa nominale (non ordinale). Inoltre, si pensa che la variabile possa molto influire nel prezzo, in quanto la zona della città influenza di molto il prezzo delle abitazioni che vi sono dentro.

df = df %>% mutate_at("zipcode", as.factor)
df$zipcode = reorder(df$zipcode,df$lprice)

ggplot(df, aes(zipcode, lprice)) +
  geom_boxplot(outlier.colour = "orchid",
               outlier.alpha = 0.7,
               outlier.size = 0.7) +
  theme(axis.text.x=element_text(vjust = 1,size = 8, angle = 90))

Il grafico geospaziale seguente è utile per comprendere meglio come si distribuiscano le zone della città e per capire eventuali patterns che ci sono tra zone contigue.

# con brewer.pal creo una palette di 10 colori dal verde alrosso
# con colorRampPalette amplio i colori da 10 al numero di zipcode presenti
# con colorFactor assegno ogni colore al valore della media dei prezzi delle case nel zipcode pres in considerazione 
pal = with(df,colorFactor(colorRampPalette(brewer.pal(10,"RdYlGn"))(df$zipcode %>% levels() %>% length()),as.data.frame(tapply(df$price,df$zipcode,mean))[zipcode,]))
                        
leaflet() %>% 
  addTiles() %>% 
  addCircleMarkers(data = df,
                   lat = ~ lat,
                   lng = ~ long,
                   radius = ~ 1,
                   color = ~ pal(as.data.frame(tapply(df$price,df$zipcode,mean))[zipcode,]))

2.2.2 Grade

Provo a stimare un GAM tramite una spline naturale con 5 nodi.

ggplot(df,aes(grade,lprice)) + 
  geom_jitter(width = .2,
              color = "gray25",
              size = 0.7) +
  geom_smooth(method = 'gam',
              formula = y ~ ns(x, knots = 5),
              color = "orchid",
              fill = "orchid1",
              alpha = 0.3)

2.2.3 Bathrooms

Anche in questo caso viene stimata una spline naturale in modo da garantire buone stime sui valori estremi, in quanto la funzione non assumerà un comportamento polinomiale. Una valida alternativa poteva essere usare una B-spline.

ggplot(df,aes(bathrooms,lprice)) + 
  geom_jitter(width = .2,
              color = "gray25",
              size = 0.7) +
  geom_smooth(method = 'gam',
              formula = y ~ ns(x, knots = 5),
              color = "orchid",
              fill = "orchid1",
              alpha = 0.3)

2.2.4 Sqft_basement

Interessante notare come la variabile sqft_basement contenga tantissimi valori pari a zero che indicano una mancanza d’informazione. Quindi, nel caso sia molto correlata con altre variabili con meno o informazioni mancanti, converrà usare queste ultime. Le procedure appena illustrate andrebbero fatte con tutte le altre variabili indipendenti possibili.

ggplot(df,aes(sqft_basement,log(price))) + 
  geom_jitter(width = .2,
              color = "gray25",
              size = 0.7) +
  geom_smooth(method = 'gam',
              formula = y ~ ns(x, knots = 25),
              color = "orchid",
              fill = "orchid1",
              alpha = 0.3)

2.3 Collinearità

Bisogna evitare la multicollinearità nel modello. Per far ciò si va inanzitutto a visualizzare un correlogramma, ovvero un grafico che mostra le correlazioni tra le variabili.

Per evitare la collinearità si potrà fare un analisi delle componenti principali (PCA), evitare alcune variabili, o applicare una penalizzazione.

df %>% 
  select_if(is.numeric) %>% 
  dplyr::select(-c(id, lat, long, lprice)) %>% 
  cor %>% 
  corrplot::corrplot(method = "number",
                     hclust.method = "ward.D2",
                     diag = F,
                     type = "upper",
                     order = "hclust",
                     number.cex = .6)

3 Modello

3.1 Modello iniziale

Analizzando la correlazione tra variabili e osservando come le variabili si distribuiscano al variare della variabile risposta (lprice), si opteranno per delle determinate scelte: - sqft_living come una B-spline con 30 nodi - zipcode in quanto si è visto che come variabile categorica ha un’influenza decisa sulla variabile risposta - viene inserita anche l’interazione tra grade e sqft_living a causa della loro forte correlazione.

Già con questo modello iniziale si ottiene una devianza spiegata del 87.1%.

parametri = list(formula = lprice ~ zipcode + bs(sqft_living, knots = 30)*grade + bs(sqft_basement, knots = 30) + view + waterfront,
            data = df)
modIniziale = do.call(gam, append(parametri, list(family = gaussian(link = identity), method = "ML")))
modIniziale %>% summary

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## lprice ~ zipcode + bs(sqft_living, knots = 30) * grade + bs(sqft_basement, 
##     knots = 30) + view + waterfront
## 
## Parametric coefficients:
##                                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         1.417e+01  2.376e+00   5.965 2.49e-09 ***
## zipcode98168                        7.854e-02  1.778e-02   4.417 1.01e-05 ***
## zipcode98032                        2.740e-03  2.170e-02   0.126 0.899511    
## zipcode98001                        2.268e-02  1.679e-02   1.351 0.176569    
## zipcode98023                       -5.526e-03  1.598e-02  -0.346 0.729560    
## zipcode98148                        1.674e-01  2.852e-02   5.869 4.45e-09 ***
## zipcode98188                        1.082e-01  2.114e-02   5.119 3.09e-07 ***
## zipcode98198                        8.106e-02  1.766e-02   4.591 4.44e-06 ***
## zipcode98003                        4.102e-02  1.766e-02   2.323 0.020193 *  
## zipcode98178                        1.476e-01  1.792e-02   8.238  < 2e-16 ***
## zipcode98030                        7.197e-02  1.798e-02   4.003 6.26e-05 ***
## zipcode98055                        1.582e-01  1.778e-02   8.900  < 2e-16 ***
## zipcode98031                        1.054e-01  1.773e-02   5.943 2.84e-09 ***
## zipcode98042                        1.017e-01  1.575e-02   6.456 1.10e-10 ***
## zipcode98022                        1.252e-01  1.837e-02   6.816 9.64e-12 ***
## zipcode98106                        3.313e-01  1.708e-02  19.399  < 2e-16 ***
## zipcode98146                        2.776e-01  1.754e-02  15.829  < 2e-16 ***
## zipcode98092                        6.180e-02  1.692e-02   3.651 0.000262 ***
## zipcode98058                        1.990e-01  1.619e-02  12.294  < 2e-16 ***
## zipcode98108                        3.597e-01  1.942e-02  18.526  < 2e-16 ***
## zipcode98038                        1.960e-01  1.563e-02  12.546  < 2e-16 ***
## zipcode98133                        4.770e-01  1.600e-02  29.811  < 2e-16 ***
## zipcode98118                        4.641e-01  1.594e-02  29.106  < 2e-16 ***
## zipcode98010                        3.106e-01  2.329e-02  13.338  < 2e-16 ***
## zipcode98056                        3.494e-01  1.646e-02  21.234  < 2e-16 ***
## zipcode98155                        4.550e-01  1.623e-02  28.039  < 2e-16 ***
## zipcode98126                        5.474e-01  1.692e-02  32.346  < 2e-16 ***
## zipcode98014                        3.658e-01  2.174e-02  16.828  < 2e-16 ***
## zipcode98166                        3.361e-01  1.804e-02  18.625  < 2e-16 ***
## zipcode98019                        3.570e-01  1.930e-02  18.500  < 2e-16 ***
## zipcode98045                        3.656e-01  1.862e-02  19.636  < 2e-16 ***
## zipcode98125                        5.769e-01  1.647e-02  35.034  < 2e-16 ***
## zipcode98028                        4.444e-01  1.763e-02  25.206  < 2e-16 ***
## zipcode98059                        3.710e-01  1.615e-02  22.969  < 2e-16 ***
## zipcode98070                        3.954e-01  2.236e-02  17.682  < 2e-16 ***
## zipcode98011                        4.743e-01  1.920e-02  24.702  < 2e-16 ***
## zipcode98034                        5.636e-01  1.579e-02  35.696  < 2e-16 ***
## zipcode98024                        5.192e-01  2.517e-02  20.624  < 2e-16 ***
## zipcode98065                        4.151e-01  1.737e-02  23.898  < 2e-16 ***
## zipcode98136                        6.891e-01  1.799e-02  38.310  < 2e-16 ***
## zipcode98144                        6.757e-01  1.705e-02  39.629  < 2e-16 ***
## zipcode98072                        5.394e-01  1.780e-02  30.299  < 2e-16 ***
## zipcode98117                        8.194e-01  1.580e-02  51.848  < 2e-16 ***
## zipcode98103                        8.260e-01  1.564e-02  52.801  < 2e-16 ***
## zipcode98107                        8.479e-01  1.794e-02  47.273  < 2e-16 ***
## zipcode98027                        5.543e-01  1.654e-02  33.522  < 2e-16 ***
## zipcode98116                        7.687e-01  1.720e-02  44.686  < 2e-16 ***
## zipcode98115                        8.327e-01  1.570e-02  53.024  < 2e-16 ***
## zipcode98008                        6.900e-01  1.765e-02  39.095  < 2e-16 ***
## zipcode98029                        6.098e-01  1.730e-02  35.255  < 2e-16 ***
## zipcode98122                        8.161e-01  1.765e-02  46.251  < 2e-16 ***
## zipcode98007                        6.930e-01  2.098e-02  33.027  < 2e-16 ***
## zipcode98177                        6.280e-01  1.811e-02  34.680  < 2e-16 ***
## zipcode98052                        6.648e-01  1.577e-02  42.159  < 2e-16 ***
## zipcode98053                        6.243e-01  1.660e-02  37.607  < 2e-16 ***
## zipcode98077                        5.283e-01  1.926e-02  27.425  < 2e-16 ***
## zipcode98074                        5.807e-01  1.643e-02  35.331  < 2e-16 ***
## zipcode98033                        8.152e-01  1.637e-02  49.803  < 2e-16 ***
## zipcode98199                        8.777e-01  1.739e-02  50.460  < 2e-16 ***
## zipcode98075                        5.936e-01  1.706e-02  34.793  < 2e-16 ***
## zipcode98105                        9.787e-01  1.854e-02  52.780  < 2e-16 ***
## zipcode98119                        9.876e-01  1.963e-02  50.321  < 2e-16 ***
## zipcode98005                        8.045e-01  2.006e-02  40.104  < 2e-16 ***
## zipcode98102                        9.766e-01  2.324e-02  42.025  < 2e-16 ***
## zipcode98006                        6.937e-01  1.619e-02  42.853  < 2e-16 ***
## zipcode98109                        1.006e+00  2.277e-02  44.190  < 2e-16 ***
## zipcode98112                        1.078e+00  1.798e-02  59.961  < 2e-16 ***
## zipcode98040                        9.452e-01  1.787e-02  52.901  < 2e-16 ***
## zipcode98004                        1.175e+00  1.738e-02  67.584  < 2e-16 ***
## zipcode98039                        1.356e+00  3.042e-02  44.586  < 2e-16 ***
## bs(sqft_living, knots = 30)1       -2.661e+00  2.379e+00  -1.119 0.263323    
## bs(sqft_living, knots = 30)2       -9.354e-01  2.320e+00  -0.403 0.686820    
## bs(sqft_living, knots = 30)3       -4.597e+00  2.871e+00  -1.601 0.109350    
## bs(sqft_living, knots = 30)4        0.000e+00  0.000e+00     NaN      NaN    
## grade                               3.016e+08  1.057e+09   0.285 0.775429    
## bs(sqft_basement, knots = 30)1      7.132e-03  8.500e-03   0.839 0.401466    
## bs(sqft_basement, knots = 30)2     -6.349e-02  3.261e-02  -1.947 0.051521 .  
## bs(sqft_basement, knots = 30)3     -2.529e-01  9.085e-02  -2.784 0.005380 ** 
## bs(sqft_basement, knots = 30)4     -4.414e-01  1.902e-01  -2.321 0.020280 *  
## view                                7.438e-02  2.038e-03  36.499  < 2e-16 ***
## waterfront                          4.707e-01  1.668e-02  28.217  < 2e-16 ***
## bs(sqft_living, knots = 30)1:grade -3.016e+08  1.057e+09  -0.285 0.775429    
## bs(sqft_living, knots = 30)2:grade -3.016e+08  1.057e+09  -0.285 0.775429    
## bs(sqft_living, knots = 30)3:grade -3.016e+08  1.057e+09  -0.285 0.775429    
## bs(sqft_living, knots = 30)4:grade -3.016e+08  1.057e+09  -0.285 0.775429    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Rank: 84/85
## R-sq.(adj) =   0.87   Deviance explained = 87.1%
## -ML = -5286.5  Scale est. = 0.036025  n = 21597

L’interazione tra sqft_living e grade causa ha una S.E. molto elevato, così anche le rispettive variabili prese da sole. Provo a rimuovere l’interazione per vedere se le stime dei coefficienti migliorano.

parametri = list(formula = lprice ~ zipcode + bs(sqft_living, knots = 30) + bs(sqft_basement, knots = 30) + grade + view + waterfront,
            data = df)
modIniziale = do.call(gam, append(parametri, list(family = gaussian(link = identity), method = "ML")))
modIniziale %>% summary

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## lprice ~ zipcode + bs(sqft_living, knots = 30) + bs(sqft_basement, 
##     knots = 30) + grade + view + waterfront
## 
## Parametric coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    -1.471e+09  4.411e+09  -0.333 0.738806    
## zipcode98168                    8.041e-02  1.783e-02   4.510 6.51e-06 ***
## zipcode98032                   -3.447e-03  2.176e-02  -0.158 0.874119    
## zipcode98001                    1.815e-02  1.683e-02   1.078 0.281017    
## zipcode98023                   -1.070e-02  1.603e-02  -0.668 0.504259    
## zipcode98148                    1.618e-01  2.861e-02   5.654 1.59e-08 ***
## zipcode98188                    1.034e-01  2.121e-02   4.876 1.09e-06 ***
## zipcode98198                    7.452e-02  1.770e-02   4.210 2.57e-05 ***
## zipcode98003                    3.435e-02  1.771e-02   1.940 0.052441 .  
## zipcode98178                    1.447e-01  1.797e-02   8.050 8.74e-16 ***
## zipcode98030                    6.750e-02  1.803e-02   3.744 0.000182 ***
## zipcode98055                    1.552e-01  1.783e-02   8.702  < 2e-16 ***
## zipcode98031                    1.008e-01  1.778e-02   5.669 1.45e-08 ***
## zipcode98042                    9.591e-02  1.580e-02   6.072 1.28e-09 ***
## zipcode98022                    1.201e-01  1.843e-02   6.521 7.16e-11 ***
## zipcode98106                    3.256e-01  1.711e-02  19.024  < 2e-16 ***
## zipcode98146                    2.762e-01  1.758e-02  15.710  < 2e-16 ***
## zipcode98092                    5.545e-02  1.697e-02   3.268 0.001086 ** 
## zipcode98058                    1.941e-01  1.623e-02  11.959  < 2e-16 ***
## zipcode98108                    3.555e-01  1.947e-02  18.253  < 2e-16 ***
## zipcode98038                    1.913e-01  1.567e-02  12.209  < 2e-16 ***
## zipcode98133                    4.695e-01  1.604e-02  29.271  < 2e-16 ***
## zipcode98118                    4.597e-01  1.598e-02  28.759  < 2e-16 ***
## zipcode98010                    3.105e-01  2.336e-02  13.292  < 2e-16 ***
## zipcode98056                    3.452e-01  1.650e-02  20.919  < 2e-16 ***
## zipcode98155                    4.491e-01  1.627e-02  27.604  < 2e-16 ***
## zipcode98126                    5.399e-01  1.696e-02  31.843  < 2e-16 ***
## zipcode98014                    3.651e-01  2.181e-02  16.744  < 2e-16 ***
## zipcode98166                    3.317e-01  1.810e-02  18.328  < 2e-16 ***
## zipcode98019                    3.511e-01  1.935e-02  18.143  < 2e-16 ***
## zipcode98045                    3.590e-01  1.867e-02  19.228  < 2e-16 ***
## zipcode98125                    5.691e-01  1.651e-02  34.480  < 2e-16 ***
## zipcode98028                    4.370e-01  1.767e-02  24.724  < 2e-16 ***
## zipcode98059                    3.671e-01  1.620e-02  22.660  < 2e-16 ***
## zipcode98070                    3.909e-01  2.243e-02  17.430  < 2e-16 ***
## zipcode98011                    4.675e-01  1.925e-02  24.281  < 2e-16 ***
## zipcode98034                    5.579e-01  1.583e-02  35.235  < 2e-16 ***
## zipcode98024                    5.132e-01  2.525e-02  20.327  < 2e-16 ***
## zipcode98065                    4.079e-01  1.741e-02  23.425  < 2e-16 ***
## zipcode98136                    6.794e-01  1.802e-02  37.705  < 2e-16 ***
## zipcode98144                    6.649e-01  1.708e-02  38.934  < 2e-16 ***
## zipcode98072                    5.338e-01  1.785e-02  29.900  < 2e-16 ***
## zipcode98117                    8.112e-01  1.583e-02  51.229  < 2e-16 ***
## zipcode98103                    8.133e-01  1.566e-02  51.950  < 2e-16 ***
## zipcode98107                    8.342e-01  1.795e-02  46.458  < 2e-16 ***
## zipcode98027                    5.496e-01  1.658e-02  33.151  < 2e-16 ***
## zipcode98116                    7.595e-01  1.724e-02  44.063  < 2e-16 ***
## zipcode98115                    8.239e-01  1.573e-02  52.362  < 2e-16 ***
## zipcode98008                    6.825e-01  1.769e-02  38.572  < 2e-16 ***
## zipcode98029                    6.019e-01  1.734e-02  34.711  < 2e-16 ***
## zipcode98122                    8.022e-01  1.766e-02  45.423  < 2e-16 ***
## zipcode98007                    6.874e-01  2.105e-02  32.663  < 2e-16 ***
## zipcode98177                    6.217e-01  1.816e-02  34.242  < 2e-16 ***
## zipcode98052                    6.581e-01  1.580e-02  41.641  < 2e-16 ***
## zipcode98053                    6.160e-01  1.664e-02  37.024  < 2e-16 ***
## zipcode98077                    5.254e-01  1.931e-02  27.205  < 2e-16 ***
## zipcode98074                    5.764e-01  1.647e-02  34.993  < 2e-16 ***
## zipcode98033                    8.120e-01  1.642e-02  49.461  < 2e-16 ***
## zipcode98199                    8.707e-01  1.743e-02  49.947  < 2e-16 ***
## zipcode98075                    5.926e-01  1.710e-02  34.664  < 2e-16 ***
## zipcode98105                    9.686e-01  1.858e-02  52.132  < 2e-16 ***
## zipcode98119                    9.744e-01  1.966e-02  49.574  < 2e-16 ***
## zipcode98005                    8.013e-01  2.012e-02  39.832  < 2e-16 ***
## zipcode98102                    9.616e-01  2.323e-02  41.389  < 2e-16 ***
## zipcode98006                    6.937e-01  1.624e-02  42.723  < 2e-16 ***
## zipcode98109                    9.953e-01  2.282e-02  43.606  < 2e-16 ***
## zipcode98112                    1.069e+00  1.802e-02  59.353  < 2e-16 ***
## zipcode98040                    9.431e-01  1.792e-02  52.635  < 2e-16 ***
## zipcode98004                    1.173e+00  1.743e-02  67.298  < 2e-16 ***
## zipcode98039                    1.363e+00  3.044e-02  44.760  < 2e-16 ***
## bs(sqft_living, knots = 30)1    1.471e+09  4.411e+09   0.333 0.738806    
## bs(sqft_living, knots = 30)2    1.471e+09  4.411e+09   0.333 0.738806    
## bs(sqft_living, knots = 30)3    1.471e+09  4.411e+09   0.333 0.738806    
## bs(sqft_living, knots = 30)4    1.471e+09  4.411e+09   0.333 0.738806    
## bs(sqft_basement, knots = 30)1 -1.559e-03  8.458e-03  -0.184 0.853795    
## bs(sqft_basement, knots = 30)2 -1.087e-02  3.241e-02  -0.335 0.737259    
## bs(sqft_basement, knots = 30)3 -4.260e-01  8.994e-02  -4.737 2.18e-06 ***
## bs(sqft_basement, knots = 30)4 -4.492e-01  1.904e-01  -2.359 0.018354 *  
## grade                           8.754e-02  1.950e-03  44.894  < 2e-16 ***
## view                            7.601e-02  2.036e-03  37.330  < 2e-16 ***
## waterfront                      4.697e-01  1.671e-02  28.103  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.869   Deviance explained =   87%
## -ML = -5215.7  Scale est. = 0.036257  n = 21597

3.2 Famiglie parametriche

A scopo illustrativo vengono creati altri modelli con altre distribuzioni e di conseguenza altre famiglie distributive e funzioni legame. All’interno della tabella seguente ci sono alcuni parametri e test dei diversi modelli.

Alcune distribuzioni non potevano essere inserite, come la binomiale o la poisson. Ad esempio, considerando la poisson, essa è una distribuzione di probabilità discreta a differenza del logaritmo del prezzo delle case. Inoltre, essa esprime le probabilità per il numero di eventi che si verificano successivamente, quindi non compatibile con i nostri dati. Tra i modelli è stato inserito anche quello di quasipoisson che è diverso dalla distribuzione di poisson. Sul quasipoisson ci sono solo le assunzioni del secondo ordine e quindi non c’è l’ipotesi distributiva. Quindi, avrà in comune con la distribuzione di poisson la funzione legame. Ciò viene fatto solo a scopo illustrativo, in quanto è meglio utilizzare, se possibile, un’ipotesi distributiva che ha proprietà più forti; nel nostro caso la gamma sostituisce molto meglio la quasipoisson.

Partendo dalla prima colonna si nota che la quasipoisson non possiede AIC, questo è dovuto al fatto che l’AIC si calcola a partire dalla funzione di verosimiglianze. La quasipoisson non possiede una funzione di verosimiglianza in quanto non ha un’ipotesi distributiva, possiede una funzione di quasiverosimiglianza.

Per le distribuzioni Gamma e quasipoisson è stato stimato il fattore di scala, \(\psi\). Nella quasipoisson \(\psi\) è stimato.

La devianza, così come il criterio di validazione incrociata generalizzata (GCV), risulta nettamente minore per il modello costruito con la distribuzione gamma.

famiglie = c("quasipoisson", "Gamma", "gaussian")
modelli = lapply(1:length(famiglie), function(i) do.call(gam, append(parametri, list(family =  get(famiglie[i])))))
tibble(famiglie,
       aic = sapply(1:length(famiglie), function(i) modelli[[i]]$aic),
       scale = sapply(1:length(famiglie), function(i) modelli[[i]]$scale),
       sigma2 = sapply(1:length(famiglie), function(i) modelli[[i]]$sig2),
       deviance = sapply(1:length(famiglie), function(i) modelli[[i]]$deviance),
       gcv.ubre = sapply(1:length(famiglie), function(i) modelli[[i]]$gcv.ubre),
       )

4 Modelli annidati

4.1 Modello pieno

Modello con tutte le variabili d’interesse.

modPieno = glm(formula = paste("lprice ~",paste(df %>% dplyr::select(-lprice, -id, -lat, -long, -date, -price) %>% names, collapse = "+")) %>% as.formula,
               data = df,
               family = gaussian)
plot(modPieno)

4.2 Metodi iterativi

modAIC = stepAIC(modPieno, 
                 steps = 100,
                 direction = "backward")

## Start:  AIC=-11502.83
## lprice ~ bedrooms + bathrooms + sqft_living + sqft_lot + floors + 
##     waterfront + view + condition + grade + sqft_above + sqft_basement + 
##     yr_built + yr_renovated + zipcode + sqft_living15 + sqft_lot15
## 
## 
## Step:  AIC=-11502.83
## lprice ~ bedrooms + bathrooms + sqft_living + sqft_lot + floors + 
##     waterfront + view + condition + grade + sqft_above + yr_built + 
##     yr_renovated + zipcode + sqft_living15 + sqft_lot15
## 
##                 Df Deviance      AIC
## - sqft_lot15     1   736.55 -11504.2
## <none>               736.53 -11502.8
## - bedrooms       1   736.69 -11500.1
## - yr_built       1   737.41 -11478.9
## - floors         1   739.23 -11425.8
## - yr_renovated   1   741.54 -11358.4
## - bathrooms      1   742.35 -11334.7
## - sqft_lot       1   743.00 -11316.0
## - sqft_above     1   746.86 -11204.0
## - sqft_living15  1   759.31 -10846.9
## - condition      1   760.11 -10824.2
## - waterfront     1   764.06 -10712.3
## - view           1   764.14 -10710.2
## - sqft_living    1   768.07 -10599.3
## - grade          1   804.86  -9588.7
## - zipcode       69  2067.10  10646.2
## 
## Step:  AIC=-11504.22
## lprice ~ bedrooms + bathrooms + sqft_living + sqft_lot + floors + 
##     waterfront + view + condition + grade + sqft_above + yr_built + 
##     yr_renovated + zipcode + sqft_living15
## 
##                 Df Deviance      AIC
## <none>               736.55 -11504.2
## - bedrooms       1   736.71 -11501.6
## - yr_built       1   737.43 -11480.3
## - floors         1   739.26 -11426.8
## - yr_renovated   1   741.57 -11359.6
## - bathrooms      1   742.36 -11336.5
## - sqft_above     1   746.89 -11205.2
## - sqft_lot       1   748.63 -11154.8
## - sqft_living15  1   759.50 -10843.4
## - condition      1   760.16 -10824.9
## - waterfront     1   764.07 -10713.9
## - view           1   764.15 -10711.7
## - sqft_living    1   768.17 -10598.5
## - grade          1   804.86  -9590.7
## - zipcode       69  2069.30  10667.3

5 Qualità del modello

Il grafico del qq-plot ci mostra che l’assunzione di normalità dei residui non è rispettata, notiamo la presenza di code pesanti.

gam.check(modIniziale)

## 
## Method: ML   Optimizer: outer newton
## full convergence after 3 iterations.
## Gradient range [-0.001403008,-0.001403008]
## (score -5215.683 & scale 0.0362571).
## Hessian positive definite, eigenvalue range [10798.5,10798.5].
## Model rank =  81 / 81

5.1 Analisi dei residui

# standardizzo i residui e faccio un grafico
with(modIniziale,(residuals - mean(residuals))/sd(residuals)) %>% 
  as.data.frame() %>% 
  ggplot(aes(x = ., y = after_stat(density))) + 
  geom_histogram(color = "orchid",
               fill = "orchid",
               alpha = 0.5,
               bins = 40) + 
  geom_density(color = "paleturquoise",
               fill = "paleturquoise",
               alpha = 0.3, # densità del colore 
               kernel = "gaussian",
               adjust = 1) +
  geom_line(aes(y = dnorm(., 0, 1)),
            color = "aquamarine")

ks.test(with(modIniziale,(residuals - mean(residuals))/sd(residuals)),
        rnorm(length(df$price),0,1))

## 
##  Asymptotic two-sample Kolmogorov-Smirnov test
## 
## data:  with(modIniziale, (residuals - mean(residuals))/sd(residuals)) and rnorm(length(df$price), 0, 1)
## D = 0.052831, p-value < 2.2e-16
## alternative hypothesis: two-sided

6 Penalizzazione

Interessante vedere come tramite penalizzazione la stima del coefficiente della variabile indipendente sqft_living, che era la variabile maggiormente correlata con la variabile dipendente, questo dovuto al fatto che essa era fortemente correlata con altre variabili indipendenti.

modLasso = glmnet(makeX(df[, !names(df) %in% c("lprice","id","lat","long", "price", "date")]),
                       df$lprice,
                       alpha = 1)
coefplot(modLasso,
         intercept = F,
         interactive = T,cex = 0.5,
         coefficients = names(df)[!names(df) %in% c("lprice","id","lat","long", "price", "date")])

coefplot(modLasso,
         intercept = F,
         interactive = T,
         cex = 0.5,
         coefficients = names(modIniziale$coefficients)[names(modIniziale$coefficients) %like% "zipcode"])

plot(modLasso)