Druga domaca naloga

Djenita Mustafoska

Prikaz pripravljenih podatkov za analizo:

library(readxl)
stanovanja <- read_excel("Apartments.xlsx")
head(stanovanja)

## # A tibble: 6 × 5
##     Age Distance Price Parking Balcony
##   <dbl>    <dbl> <dbl>   <dbl>   <dbl>
## 1     7       28  1640       0       1
## 2    18        1  2800       1       0
## 3     7       28  1660       0       0
## 4    28       29  1850       0       1
## 5    18       18  1640       1       1
## 6    28       12  1770       0       1

Razlaga podatkov:

Torej imamo podatke za:

Age, kar predstavlja koliko je stanovanje staro v letih.
Distance, kar predstavlja razdaljo stanovanja od centra v kilometrih.
Price, kar predstavlja ceno stanovanja na kvadratni meter v evrih, to bo tudi naša preučevana spremenljivka.
Parking, predstavlja če stanovanje ima partking (kar označimo z 1) ali nima parkinga (označimo z 0).
Balcony, predstavlja če stanovanje ima balkon (kar označimo z 1) ali nima balkona (označimo z 0).

Na podlagi naših podatkov je edina smiselna enota preučevanje cena stanovanja, predvidevamo da bo ta odvisna od starosti stanovanja, razdalje od centra in od tega ali ima uporabnik stanovanja dostop do balkona ali parkinga (to bomo tudi ugotovili z nadaljno analizo).

Vir:

Podatki pridobljeni od profesorja.

Cilj analize:

Napraviti smiselen linearen regresijski model za ceno stanovanja na kvadratni meter in pojasniti učinek regresijskih koeficientov v modelu.

Opredelitev modela multiple regresije:

Torej naša odvisna spremenljivka y_i modela multiple regresije bo cena stanovanja, pojasnjevalne spremenljivke x_ij bodo starost stanovanja, razdalja stanovanja od centra, in indikatorski spremenljivki parking in balkon, ki predstavljata če ima stanovanje parking in/ali če ima stanovanje balkon. Regresijska enačba, ki bo določala naš model:

Price = β₀ + β₁ * Age + β₂ * Distance + β₃ * Balcony + β₄ * Parking + ε

Grafični prikaz med izbranimi spremenljivkami

Napravimo razsevni diagram za vse pojasnjevalne spremenljivke z ocenjevano spremenljivko:

library(ggplot2)


library(ggpubr)

plot_colors <- c("darkgoldenrod1", "deepskyblue3", "seagreen3", "sienna3")

plot1 <- ggplot(stanovanja, aes(x = Age, y = Price)) +
  geom_point(color = plot_colors[1]) +
  labs(x = "Leta", y = "Cena") +
  ggtitle("Razsevni diagram starosti stanovanja in cene")

plot2 <- ggplot(stanovanja, aes(x = Distance, y = Price)) +
  geom_point(color = plot_colors[2]) +
  labs(x = "Razdalja od centra", y = "Cena") +
  ggtitle("Razsevni diagram razdalje od centra in cene")

plot3 <- ggplot(stanovanja, aes(x = Balcony, y = Price)) +
  geom_point(color = plot_colors[3]) +
  labs(x = "Balkon", y = "Cena") +
  ggtitle("Razsevni diagram balkona in cene")

plot4 <- ggplot(stanovanja, aes(x = Parking, y = Price)) +
  geom_point(color = plot_colors[4]) +
  labs(x = "Parking", y = "Cena") +
  ggtitle("Razsevni diagram parkinga in cene")

ggarrange(plot1, plot2, plot3, plot4, ncol = 2, nrow = 2)

Nekako iz grafov lahko sklepamo da bližje kot smo centru bo v povprečju višja cena stanovanja (negativna korelacija), lahko tudi sklepamo, da če stanovanje ima parking potem bo imelo v povprečju višjo ceno (pozitivna korelacija).

Ocena regresijske funkcije po metodi najmanjših kvadratov

lm(Price ~ Age + Distance + Balcony + Parking, data = stanovanja)

## 
## Call:
## lm(formula = Price ~ Age + Distance + Balcony + Parking, data = stanovanja)
## 
## Coefficients:
## (Intercept)          Age     Distance      Balcony      Parking  
##    2301.667       -6.799      -18.045        1.935      196.168

model <- lm(Price ~ Age + Distance + Balcony + Parking, data = stanovanja)
summary(model)

## 
## Call:
## lm(formula = Price ~ Age + Distance + Balcony + Parking, data = stanovanja)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -459.92 -200.66  -57.48  260.08  594.37 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2301.667     94.271  24.415  < 2e-16 ***
## Age           -6.799      3.110  -2.186  0.03172 *  
## Distance     -18.045      2.758  -6.543 5.28e-09 ***
## Balcony        1.935     60.014   0.032  0.97436    
## Parking      196.168     62.868   3.120  0.00251 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.7 on 80 degrees of freedom
## Multiple R-squared:  0.5004, Adjusted R-squared:  0.4754 
## F-statistic: 20.03 on 4 and 80 DF,  p-value: 1.849e-11

Torej, če bodo vse pojasnjevalne spremenljivke 0 (novo stanovanje v centru, vendar brez balkona in parkinga) bo po oceni z metodo najmanjših kvadratov stanovanje stalo 2301.667 evrov na kvadratni meter, če je stanovanje starejše za eno leto se bo v povprečju cena stanovanja zmanjšala za 6.799 evrov na kvadratni meter, za vsak dodaten kilometer od centra bo stanovanje v povprečju cenejše za 18.045 evrov na kvadratni meter, če ima stanovanje balkon bo v povprečju dražje za 1.935 evrov na kvadratni meter in če ima stanovanje parking bo dražje za 196.168 evrov na kvadratni meter.

Če pogledamo p-vrednost naših pojasnjevalnih spremenljivk, vidimo da so vse manjše od 0.05 razen pojasnjevalne spremenljivke za balkon, ki je 0.97436, kar pomeni da ta spremenljivka ni statistično značilna za naš model. Torej jo lahko odstranimo.

Izboljšan regresijski model za preučevano slučajno spremenljivko

lm(Price ~ Age + Distance + Parking, data = stanovanja)

## 
## Call:
## lm(formula = Price ~ Age + Distance + Parking, data = stanovanja)
## 
## Coefficients:
## (Intercept)          Age     Distance      Parking  
##    2302.675       -6.805      -18.046      196.083

model2 <- lm(Price ~ Age + Distance + Parking, data = stanovanja)

summary(model2)

## 
## Call:
## lm(formula = Price ~ Age + Distance + Parking, data = stanovanja)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -458.95 -199.49  -58.34  261.20  593.42 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2302.675     88.380  26.054  < 2e-16 ***
## Age           -6.805      3.085  -2.206  0.03024 *  
## Distance     -18.046      2.741  -6.584 4.22e-09 ***
## Parking      196.083     62.423   3.141  0.00235 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 272 on 81 degrees of freedom
## Multiple R-squared:  0.5004, Adjusted R-squared:  0.4819 
## F-statistic: 27.04 on 3 and 81 DF,  p-value: 3.227e-12

Predpostavke regresijskega modela

Poglejmo če so izpolnjene predpostavke obeh regresijskih modelov:

Linearnost regresijskega modela, je izpolnjena.
Zahteva da so vrednosti pojasnjevalnih spremenljivk fiksne (nestohastične) pri ponovitvah vzorcev, je izpolnjena.
Predpostavka o ničelni povprečni vrednosti residualov, bi morala biti izpolnjena, poglejmo.

mean(resid(model))

## [1] -9.11787e-15

mean(resid(model2))

## [1] 2.712863e-15

Vrednost je zelo blizu 0.

Preden naredimo preizkus ali ima residual upanje res enako 0, preverimo če je residual porazdeljen normalno. Postavimo naslednjo hipotezo za oba modela:

H₀ : Residual je porazdeljen normalno.
H₁ : Residual ni porazdeljen normalno.

Hipotezo preverimo s Shapirov-Wilkovim testom.

shapiro.test(resid(model))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.95749, p-value = 0.006869

shapiro.test(resid(model2))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model2)
## W = 0.95738, p-value = 0.006762

Pri obeh preizkusih dobimo p-vrednost zelo majhno (<0.05), hipotezo H₀ zavrnemo in nemoremo sklepati, da sta residuala porazdeljena normalno.

Poglejmo si še kako izgledata histogram residualov:

hist(resid(model),
     xlab = "Vrednost residuala",
     ylab = "Frekvenca",
     main = "Histogram residualov",
     col = "sienna1",
     ylim = range(0,15),
     xlim = range(-500,700))

hist(resid(model2),
     xlab = "Vrednost residuala",
     ylab = "Frekvenca",
     main = "Histogram residualov",
     col = "turquoise4",
     ylim = range(0,15),
     xlim = range(-500,700))

Histograma izgledata približno normalno porazdeljena.

Naredimo t-test za oba modela za naslednjo hipotezo:

H₀ : Residual ima pričakovano vrednosto enako nič, E(ε) \(=\) 0.
H₁ : E(ε) \(\neq\) 0.

t.test(resid(model), mu = 0)

## 
##  One Sample t-test
## 
## data:  resid(model)
## t = -3.1477e-16, df = 84, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -57.60409  57.60409
## sample estimates:
##    mean of x 
## -9.11787e-15

t.test(resid(model2), mu = 0)

## 
##  One Sample t-test
## 
## data:  resid(model2)
## t = 9.3653e-17, df = 84, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -57.60446  57.60446
## sample estimates:
##    mean of x 
## 2.712863e-15

Pri obeh preizkusih dobimo p-vrednost = 1 (>0.05), kar implicira da je upanje res 0, hipotezo H₀ ne zavrnemo. Predpostavimo, da je tretja predpostavka izpolnjena.

Predpostavka o homoskedastičnost

Preverimo, če imamo v našem modelu homoskedastičnost, postavimo ničelno hipotezo:

H₀ : Varianca je konstantna (model je homoskedastičen).
H₁ : Varianca ni konstantna (model je heteroskedastičen).

# Naložimo install.packages("lmtest")
library(lmtest)

bptest(model)

## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 6.4758, df = 4, p-value = 0.1663

bptest(model2)

## 
##  studentized Breusch-Pagan test
## 
## data:  model2
## BP = 6.4334, df = 3, p-value = 0.09232

V obeh primerih dobimo p-vrednost > 0.05 in ničelno hipotezo ne zavrnemo. Predostavimo, da je naš model homoskedastičen.

Odsotnost avtokorelacije

Narišimo razsevni diagram residualov, na podlagi katerega bomo videli avtokorealcijo.

plot(residuals(model), xlab = "opazovanje", ylab = "residual", main = "Razsevni diagram residualov")

plot(residuals(model2), xlab = "opazovanje", ylab = "residual", main = "Razsevni diagram residualov")

Razsevni diagram nakazuje, da ni avtokorelacije.

Nekoreliranost med pojasnjevalnimi spremenljivkami in slučajno spremenljivko residuala

residuals <- resid(model)

residuals2 <- resid(model2)

cov(residuals, stanovanja[, c("Age", "Balcony", "Distance", "Parking")])

##               Age      Balcony      Distance       Parking
## [1,] 3.762467e-14 4.180436e-15 -1.450559e-13 -7.092128e-15

cov(residuals2, stanovanja[, c("Age", "Distance", "Parking")])

##                Age      Distance      Parking
## [1,] -3.456653e-13 -4.397408e-14 2.998593e-16

Vrednosti kovariance so zelo majhne (blizu 0), zato sklepamo da je tudi ta predpostavka izpolnjena.

Število opazovanj presega število pojasnjevalnih spremenljivk, torej tudi ta predpostavka je izpolnjena
Variabilnost vrednosti pojasnjevalnih spremenljivk mora biti končno pozitivno število.

variance_Age <- var(stanovanja$Age)
variance_Distance <- var(stanovanja$Distance)
variance_Balcony <- var(stanovanja$Balcony)
variance_Parking <- var(stanovanja$Parking)

print(c(variance_Age,variance_Balcony,variance_Distance,variance_Parking))

## [1]  93.9644258   0.2487395 129.4375350   0.2529412

Tudi ta predpostavka je izpolnjena.

Predpostavimo tudi, da je naš model pravilno specificiran.
Preverimo odsotnost popolne multikolinearnosti.

cor(stanovanja[, c("Age", "Distance", "Balcony" , "Parking")])

##                  Age     Distance      Balcony     Parking
## Age       1.00000000  0.042908128 -0.055303626 -0.12399145
## Distance  0.04290813  1.000000000  0.001530351 -0.30711555
## Balcony  -0.05530363  0.001530351  1.000000000 -0.03406038
## Parking  -0.12399145 -0.307115547 -0.034060379  1.00000000

library(car)

vif(model)

##      Age Distance  Balcony  Parking 
## 1.019308 1.104263 1.004865 1.121335

vif(model2)

##      Age Distance  Parking 
## 1.015640 1.104172 1.119347

Na podlagi statistike VIF ugotavljamo, da v konkretnem primeru ni težave s premočno multikolinearnostjo, saj so vse vrednosti pod kritično mejo 5.

Vidimo, da je korelacija med prisotnostjo parkirišča in razdaljo od centra dokaj nizka (-0.307), kar implicira negativno korelacijo med razdaljo od centra in prisotnostjo parkirišča to pomeni, da bolj, ko smo oddaljeni od centra večja je verjetnost, da imamo na razpolago parkirišče, vendar je ta vrednost še vedno zelo majhna, zato lahko sklepamo, da je ta predpostavka tudi izpolnjena.

Iz podatkov odstranimo potencialne osamelce z visokim vplivom na ocenjevano regresijsko funkcijo.

ostanki <- round(rstandard(model), 3)

cookd <- round(cooks.distance(model), 3)

print(cookd)

##     1     2     3     4     5     6     7     8     9    10    11    12    13 
## 0.003 0.024 0.002 0.013 0.022 0.002 0.000 0.011 0.000 0.011 0.016 0.008 0.014 
##    14    15    16    17    18    19    20    21    22    23    24    25    26 
## 0.000 0.007 0.004 0.011 0.001 0.005 0.011 0.005 0.039 0.000 0.005 0.024 0.002 
##    27    28    29    30    31    32    33    34    35    36    37    38    39 
## 0.018 0.019 0.005 0.000 0.037 0.000 0.050 0.001 0.008 0.020 0.013 0.194 0.067 
##    40    41    42    43    44    45    46    47    48    49    50    51    52 
## 0.005 0.014 0.002 0.000 0.006 0.001 0.014 0.003 0.011 0.016 0.007 0.005 0.003 
##    53    54    55    56    57    58    59    60    61    62    63    64    65 
## 0.055 0.007 0.078 0.000 0.084 0.026 0.000 0.006 0.028 0.002 0.011 0.022 0.002 
##    66    67    68    69    70    71    72    73    74    75    76    77    78 
## 0.000 0.011 0.000 0.011 0.018 0.008 0.013 0.000 0.007 0.004 0.009 0.001 0.016 
##    79    80    81    82    83    84    85 
## 0.003 0.009 0.016 0.007 0.005 0.003 0.003

hist(ostanki,
     xlab = "Standardizirani ostanki",
     ylab = "Frekvenca",
     main = "Histogram standardiziranih ostankov",
     col = "slateblue3",
     xlim = range(-3,3))

ostanki2 <- round(rstandard(model2), 3)

cookd2 <- round(cooks.distance(model2), 3)

print(cookd2)

##     1     2     3     4     5     6     7     8     9    10    11    12    13 
## 0.003 0.023 0.002 0.011 0.017 0.002 0.000 0.008 0.000 0.011 0.015 0.009 0.012 
##    14    15    16    17    18    19    20    21    22    23    24    25    26 
## 0.000 0.008 0.004 0.009 0.001 0.006 0.010 0.005 0.038 0.000 0.004 0.023 0.002 
##    27    28    29    30    31    32    33    34    35    36    37    38    39 
## 0.017 0.019 0.005 0.000 0.039 0.000 0.051 0.001 0.007 0.020 0.013 0.217 0.070 
##    40    41    42    43    44    45    46    47    48    49    50    51    52 
## 0.005 0.013 0.001 0.000 0.006 0.001 0.015 0.004 0.011 0.015 0.007 0.006 0.002 
##    53    54    55    56    57    58    59    60    61    62    63    64    65 
## 0.057 0.007 0.092 0.000 0.087 0.026 0.000 0.006 0.023 0.002 0.011 0.017 0.002 
##    66    67    68    69    70    71    72    73    74    75    76    77    78 
## 0.000 0.008 0.000 0.011 0.015 0.009 0.012 0.000 0.008 0.004 0.009 0.001 0.015 
##    79    80    81    82    83    84    85 
## 0.004 0.011 0.015 0.007 0.006 0.002 0.004

hist(ostanki2,
     xlab = "Standardizirani ostanki",
     ylab = "Frekvenca",
     main = "Histogram standardiziranih ostankov",
     col = "slateblue3",
     xlim = range(-3,3))

Ostanki obeh modelov se porazdeljujejo približno normalno, vsi se nahajajo znotraj intervala [-3,3]. V konkretnem primeru nimamo težav z osamelci.

hist(cookd,
     xlab = "Cookove razdalje",
     ylab = "Frekvenca",
     main = "Histogram Cookovih razdalj",
     col = "slateblue3")

hist(cookd2,
     xlab = "Cookove razdalje",
     ylab = "Frekvenca",
     main = "Histogram Cookovih razdalj",
     col = "slateblue3")

V analizi na podlagi Cookovih razdalj smo se odločili, da ne bomo odstranili enot iz modelov, saj so vse razdalje pod mejno vrednostjo 1.

Vse predpostavke pri obeh modelih so izpolnjene.

Odločimo se, da je drugi model bolj primeren model za naše podatke na podlagi sledečega:

Drugi model ima višjo vrednost prilagojenega R-kvadrata kot prvi model, prilagojeni R-kvadrat upošteva število pojasnjevalnih spremenljivk v modelu in je manjši glede na vključitev manj pomembnih spremenljivk v model.
Drugi model ima manjšo standardno napako residuala, kar pomeni da se drugi model bolje prilega podatkom.
Drugi model ima večjo vrednost F-statistike kot prvi model in p-vrednost je v drugem modelu manjša kot v prvem. To pomeni da regresijski koeficienti dajejo boljšo napoved modela glede na podatke.
p-vrednost za pojasnjevalno spremenljivko balkon je v prvem modelu >0.05, kar pomeni da glede na drugi model ni značilne statistične razlike.

Prikaz ustreznega modela

model2 <- lm(Price ~ Age + Distance + Parking, data = stanovanja)

summary(model2)

## 
## Call:
## lm(formula = Price ~ Age + Distance + Parking, data = stanovanja)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -458.95 -199.49  -58.34  261.20  593.42 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2302.675     88.380  26.054  < 2e-16 ***
## Age           -6.805      3.085  -2.206  0.03024 *  
## Distance     -18.046      2.741  -6.584 4.22e-09 ***
## Parking      196.083     62.423   3.141  0.00235 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 272 on 81 degrees of freedom
## Multiple R-squared:  0.5004, Adjusted R-squared:  0.4819 
## F-statistic: 27.04 on 3 and 81 DF,  p-value: 3.227e-12

sqrt(summary(model2)$r.squared)

## [1] 0.7073978

Torej za naš model veljajo naslednje stvari. Model predpostavlja, da ko so vse pojasnjevalne spremenljivke 0 (oddaljenost od centra 0, stanovanje staro 0 let, brez balkona in parkirišča), da cena stanovanje v povprečju stane 2302.675 evrov na kvadratni meter, za vsako leto ko je stanovanje starejše cena v povprečju pade za 6.805 evrov na kvadratni meter, za vsak dodaten kilometer od centra cena v povprečju stanovanja na kvadrtni meter pade za 18.046 evrov in če ima stanovanje možnost uporabe parkirišča cena stanovanja v povprečju zraste za 196.083 evrov.

Opisna statistika podatkov uporabljenih v končnem modelu

stanovanja$Parking <- as.factor(stanovanja$Parking)
stanovanja$Balcony <- as.factor(stanovanja$Balcony)
summary(stanovanja)

##       Age           Distance         Price      Parking Balcony
##  Min.   : 1.00   Min.   : 1.00   Min.   :1400   0:42    0:48   
##  1st Qu.:12.00   1st Qu.: 4.00   1st Qu.:1710   1:43    1:37   
##  Median :18.00   Median :12.00   Median :1950                  
##  Mean   :18.55   Mean   :14.22   Mean   :2019                  
##  3rd Qu.:24.00   3rd Qu.:20.00   3rd Qu.:2290                  
##  Max.   :45.00   Max.   :45.00   Max.   :2820

Povprečna cena stanovanja je 2019 evrov, najcenejše stanovanje je 1400 in najdražje je 2820 evorv, dokaj majhen razpon v ceni. Povprečna razdalja od centra je 14.22 kilometrov najmanjša razdalja je 1 kilometer in največja 45 kilometrov oddaljenosti od centra. Razmerje parkirišč in balkonov, ki ga stanovanja imajo je priblišno enako (v prvem primeru 42:43 in v drugem primeru 48:37).

Razlaga koeficientov in analize

Parcialne koeficiente za naš model smo že razložili. Multipli determinacijski koeficient R-kvadrat, ki je 0.5004 nam pove koliko variabilnosti lahko razložimo v preučevani spremenljivki s pomočjo pojasnjevalnih spremenljivk, višja kot je vrednost boljša je razlaga variablinosti z modelom (boljši model imamo), v našem primeru je ta koeficient zadovoljivo visok (razloži približno 50.04% variance cene v modelu). F-statistika je v našem primeru 27.04, ta nam pove natančnost regersijskega modela. Večja vrednost F-statistike in manjša p-vrednost implicira, da vsaj ena pojansjevalna spremenljivka v modelu pomembno pripomore k razlagi variance. V našem primeru je F-statistika relativno velika in p-vrednost povezana z njo precej majhna, kar implicira, da je model statistično značilen.