V tem delu domače naloge bom analiziral podatke o ceni nepremičnin ter vpliv določenih faktorjev nanje.

Za izbor spremenljivk sem se odločil predvsem zato, ker me je zanimalo, kako poleg kvadrature, še en “osnoven” (=število spalnic) ter, po lastnem mnenju, bolj luksuzen (=klet) faktor vplivata na samo ceno nepremičnine. Kakorkoli, neke vrste nenapisano pravilo je, da takoj ko nepremičnini dodamo nekaj z dodano vrednostjo (klet, prezračevanje, parkirišče,…), bo tudi cena nepremičnine narasla.

Vir podatkov: https://www.kaggle.com/datasets/yasserh/housing-prices-dataset

Opis spremenljivk:

Uvoz podatkov, faktoriranje kategorialnih spremenljivk “klet” in “spalnice”:

podatki <- read.table("/cloud/project/housing_csv.csv", header=TRUE, sep=";", dec=",")

#faktoriram kategorialne spremenljivke
podatki$spalniceF <- factor(podatki$spalnice,
                            levels = c("3", "2", "1", "4", "5", "6"),
                            labels = c("3", "2", "1", "4", "5", "6"))

podatki$kletF <- factor(podatki$klet,
                            levels = c("no", "yes"),
                            labels = c("0", "1"))

head(podatki, 10) #Prikažemo prvih 10 opazovanj
##    ID   cena kvadrati spalnice klet spalniceF kletF
## 1   1 133.00     7.42        4   no         4     0
## 2   2 122.50     8.96        4   no         4     0
## 3   3 122.50     9.96        3  yes         3     1
## 4   4 122.15     7.50        4  yes         4     1
## 5   5 114.10     7.42        4  yes         4     1
## 6   6 108.50     7.50        3  yes         3     1
## 7   7 101.50     8.58        4   no         4     0
## 8   8 101.50    16.20        5   no         5     0
## 9   9  98.70     8.10        4  yes         4     1
## 10 10  98.00     5.75        3   no         3     0

Opis številskih spremenljivk:

library(pastecs)
round(stat.desc(podatki[c(-1, -4, -5, -6, -7)], 2))
##              cena kvadrati
## median         43        5
## mean           48        5
## SE.mean         1        0
## CI.mean.0.95    2        0
## var           350        5
## std.dev        19        2
## coef.var        0        0

Matrika razsevnih grafikonov za številski spremenljivki:

library(car)
scatterplotMatrix(podatki[ , c(-1, -4, -5, -6, -7)],  
                  smooth = FALSE)

Korelacijska matrika za številski spremenljivki:

library(Hmisc)
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
rcorr(as.matrix(podatki[ , c(-1, -4, -5, -6, -7)])) #Izpišemo korelacijsko matriko
##          cena kvadrati
## cena     1.00     0.54
## kvadrati 0.54     1.00
## 
## n= 545 
## 
## 
## P
##          cena kvadrati
## cena           0      
## kvadrati  0

V regresijski model “fit1” vključim odvisno spremenljivko “cena” in pojasnjevalni spremenljivki “kvadrati” in “spalniceF”;

Preverim multikolinearnost med spremenljivkami z VIF-statistiko ter izračunam povprečni VIF;

Prikažem histogram standardiziranih ostankov ter izvedem Shapiro-Wilkov test normalnosti porazdelitve standardiziranih ostankov; ta podatek sicer lahko spustimo, saj je vzorec dovolj velik (n=545);

prikažem histogram Cookovih razdalij

fit1 <- lm(cena ~ kvadrati + spalniceF,
           data = podatki) 

vif(fit1) #Multikolinearnost
##               GVIF Df GVIF^(1/(2*Df))
## kvadrati  1.030344  1        1.015059
## spalniceF 1.030344  5        1.002994
mean(vif(fit1))
## [1] 1.67979
podatki$StdOstanki <- round(rstandard(fit1), 3) 
podatki$CooksD <- round(cooks.distance(fit1), 3)

hist(podatki$StdOstanki, 
     xlab = "Standardizirani ostanki", 
     ylab = "Frekvenca", 
     main = "Histrogram standardiziranih ostankov")

shapiro.test(podatki$StdOstanki)
## 
##  Shapiro-Wilk normality test
## 
## data:  podatki$StdOstanki
## W = 0.96374, p-value = 2.43e-10
hist(podatki$CooksD, 
     xlab = "Cookove razdalje", 
     ylab = "Frekvenca", 
     main = "Histrogram Cookovih razdalj")

VIF statistika za vse spremenljivke ne preseže vrednosti 1,1 in je s tem daleč od mejne vrednosti 5, s katero bi potencialno potrdili multikolinearnost; prav tako je povprečna vrednost VIF statistike dovolj nizka, da lahko to predpostavko pretirane multikolinearnosti zavrnemo.

Prikažemo 10 najnižjih in 10 najvišjih vrednosti standardiziranih ostankov ter 10 najvišjih vrednosti Cookovih razdalij:

head(podatki[order(podatki$StdOstanki), c("ID", "StdOstanki")], 10)
##      ID StdOstanki
## 404 404     -3.221
## 126 126     -2.331
## 453 453     -2.304
## 212 212     -2.247
## 402 402     -2.211
## 537 537     -1.900
## 488 488     -1.895
## 230 230     -1.881
## 532 532     -1.867
## 489 489     -1.850
head(podatki[order(-podatki$StdOstanki), c("ID", "StdOstanki")], 10)
##    ID StdOstanki
## 1   1      4.633
## 4   4      3.871
## 3   3      3.619
## 2   2      3.483
## 6   6      3.357
## 5   5      3.345
## 10 10      3.139
## 14 14      2.987
## 19 19      2.848
## 21 21      2.833
head(podatki[order(-podatki$CooksD), c("ID", "CooksD")], 10)
##      ID CooksD
## 113 113  0.172
## 396 396  0.172
## 537 537  0.059
## 404 404  0.042
## 446 446  0.039
## 529 529  0.039
## 126 126  0.038
## 1     1  0.037
## 35   35  0.035
## 29   29  0.029

Opazimo, da enote 1, 2, 3, 4, 5, 6, 10 in 404 pri standardiziranih ostankih padejo izven ranga [-3, 3], pri Cookovih razdaljah pa izstopajo enote 113, 396 in 537, ki niso dovolj blizu vrednosti razdalij drugih enot; vse omenjene enote tako izločimo iz testiranja s funkcijo “filter”:

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:Hmisc':
## 
##     src, summarize
## The following objects are masked from 'package:pastecs':
## 
##     first, last
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
podatki <- podatki %>%
  filter(!ID %in% c(1, 2, 3, 4, 5, 6, 10, 113, 396, 404, 537))

Preverimo še homoskedastičnost z Breusch-Paganovim testom:

library(olsrr)
## 
## Attaching package: 'olsrr'
## The following object is masked from 'package:datasets':
## 
##     rivers
ols_test_breusch_pagan(fit1)
## 
##  Breusch Pagan Test for Heteroskedasticity
##  -----------------------------------------
##  Ho: the variance is constant            
##  Ha: the variance is not constant        
## 
##               Data               
##  --------------------------------
##  Response : cena 
##  Variables: fitted values of cena 
## 
##          Test Summary           
##  -------------------------------
##  DF            =    1 
##  Chi2          =    134.9764 
##  Prob > Chi2   =    3.341104e-31

Zavrnemo H0 pri p<0,001, kar pomeni kršitev; potrebno bo uporabiti robustne standardne napake, ki so “odporne” na heteroskedastičnost;

Ponovimo regresijski model z Whitovimi robustnimi standardnimi napakami:

library(estimatr)
fit_robust1 <- lm_robust(cena ~ kvadrati + spalniceF,
                         se_type = "HC1",
                        data = podatki)

Multipli korelacijski koeficient:

sqrt(summary(fit_robust1)$r.squared)
## [1] 0.629363
fit_robust1 <- lm_robust(cena ~ kvadrati + spalniceF,
                         se_type = "HC1",
                        data = podatki)

Prikažemo rezultate multiple regresijske analize z uporabo robustnih standardnih napak:

summary(fit_robust1) #Rezultati regresijske analize
## 
## Call:
## lm_robust(formula = cena ~ kvadrati + spalniceF, data = podatki, 
##     se_type = "HC1")
## 
## Standard error type:  HC1 
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper  DF
## (Intercept)   28.693     1.5988  17.947 1.403e-56  25.5523   31.834 528
## kvadrati       3.921     0.3141  12.483 1.614e-31   3.3041    4.538 528
## spalniceF2   -10.553     1.0885  -9.694 1.460e-20 -12.6909   -8.414 528
## spalniceF1   -16.116     3.9211  -4.110 4.586e-05 -23.8186   -8.413 528
## spalniceF4     4.217     1.9312   2.184 2.942e-02   0.4235    8.011 528
## spalniceF5     7.872     3.7767   2.084 3.760e-02   0.4530   15.291 528
## 
## Multiple R-squared:  0.3961 ,    Adjusted R-squared:  0.3904 
## F-statistic: 70.88 on 5 and 528 DF,  p-value: < 2.2e-16
sqrt(summary(fit_robust1)$r.squared) #Multipli korelacijski koeficient
## [1] 0.629363
library(ppcor)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:olsrr':
## 
##     cement
## The following object is masked from 'package:dplyr':
## 
##     select
parc_kor <- pcor(podatki[ ,c("cena", "kvadrati")]) #Izračun parcialnih korelacijskih koeficientov
round(parc_kor$estimate, 2)
##          cena kvadrati
## cena     1.00     0.54
## kvadrati 0.54     1.00
round(parc_kor$p.value, 3)
##          cena kvadrati
## cena        0        0
## kvadrati    0        0

Nov regresijski model: dodamo spremenljivko “kletF”, uporaba Whitovih standardnih napak

fit_robust2 <- lm_robust(cena ~ kvadrati + spalniceF + kletF,
                         se_type = "HC1",
                        data = podatki)

summary(fit_robust2)
## 
## Call:
## lm_robust(formula = cena ~ kvadrati + spalniceF + kletF, data = podatki, 
##     se_type = "HC1")
## 
## Standard error type:  HC1 
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper  DF
## (Intercept)   27.081     1.5894  17.039 3.512e-52  23.9587   30.203 527
## kvadrati       3.902     0.3063  12.740 1.340e-32   3.3005    4.504 527
## spalniceF2    -9.953     1.0792  -9.222 6.974e-19 -12.0731   -7.833 527
## spalniceF1   -14.433     3.9523  -3.652 2.864e-04 -22.1976   -6.669 527
## spalniceF4     4.385     1.9265   2.276 2.322e-02   0.6010    8.170 527
## spalniceF5     6.675     3.5009   1.907 5.711e-02  -0.2025   13.552 527
## kletF1         4.402     1.2056   3.652 2.867e-04   2.0339    6.771 527
## 
## Multiple R-squared:  0.411 , Adjusted R-squared:  0.4043 
## F-statistic: 68.31 on 6 and 527 DF,  p-value: < 2.2e-16

Primerjava testov fit_robust1 in fit_robust2 z Waldovim testom:

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
waldtest(fit_robust1, fit_robust2, test = "F")
## Wald test
## 
## Model 1: cena ~ kvadrati + spalniceF
## Model 2: cena ~ kvadrati + spalniceF + kletF
##   Res.Df Df      F    Pr(>F)    
## 1    528                        
## 2    527  1 13.334 0.0002867 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Waldov test nam prikaže, da se je z dodano spremenljivko “kletF” delež variabilnosti povečal (p<0,001).

Pojasnila parcialnih regresijskih koeficientov: