R za novu generaciju IT stručnjaka
- Objavljeno u RAČUNALA
U današnjem ćemo tekstu prezentirati primjer jednostavne analize podataka unutar sustava R. Podaci koje ćemo ovom prilikom koristiti dolaze zajedno s instalacijom sustava i sadržavaju skup podataka o cijenama prodanih kuća u gradu Amesu, Iowa u razdoblju od 2006. do 2010.
Dodatno, skup podataka moguće je pronaći i na drugim web stranicama poput https://lib.stat.cmu.edu/datasets/boston. Skup podataka sadrži brojne varijable o kvaliteti i količini fizičkih atributa stambenih kuća u Iowi prodanih u razdoblju od 2006. do 2010. godine. To su najčešće varijable koje kupac želi znati o nekretnini (kvadratura, broj spavaćih soba i kupaonica, veličina parcele itd.).
Kako ovaj skup podataka dolazi instalacijom sustava R, dovoljno je u konzolu upisati ime objekta (skupa) kako bi se ispisao na konzoli. Također, sve funkcije sustava možemo odmah primijeniti na podatke.
Pogledajmo prva dva retka skupa podataka.
ames[1:2,]
Order PID MS.SubClass MS.Zoning Lot.Frontage Lot.Area Street Alley
1 1 526301100 20 RL 141 31770 Pave
2 2 526350040 20 RH 80 11622 Pave
Lot.Shape Land.Contour Utilities Lot.Config Land.Slope Neighborhood
1 IR1 Lvl AllPub Corner Gtl NAmes
2 Reg Lvl AllPub Inside Gtl NAmes
Condition.1 Condition.2 Bldg.Type House.Style Overall.Qual Overall.Cond
1 Norm Norm 1Fam 1Story 6 5
2 Feedr Norm 1Fam 1Story 5 6
Year.Built Year.Remod.Add Roof.Style Roof.Matl Exterior.1st Exterior.2nd
1 1960 1960 Hip CompShg BrkFace Plywood
2 1961 1961 Gable CompShg VinylSd VinylSd
Mas.Vnr.Type Mas.Vnr.Area Exter.Qual Exter.Cond Foundation Bsmt.Qual
1 Stone 112 TA TA CBlock TA
2 None 0 TA TA CBlock TA
Bsmt.Cond Bsmt.Exposure BsmtFin.Type.1 BsmtFin.SF.1 BsmtFin.Type.2
1 Gd Gd BLQ 639 Unf
2 TA No Rec 468 LwQ
BsmtFin.SF.2 Bsmt.Unf.SF Total.Bsmt.SF Heating Heating.QC Central.Air
1 0 441 1080 GasA Fa Y
2 144 270 882 GasA TA Y
Electrical X1st.Flr.SF X2nd.Flr.SF Low.Qual.Fin.SF Gr.Liv.Area
1 SBrkr 1656 0 0 1656
2 SBrkr 896 0 0 896
Bsmt.Full.Bath Bsmt.Half.Bath Full.Bath Half.Bath Bedroom.AbvGr
1 1 0 1 0 3
2 0 0 1 0 2
Kitchen.AbvGr Kitchen.Qual TotRms.AbvGrd Functional Fireplaces
1 1 TA 7 Typ 2
2 1 TA 5 Typ 0
Fireplace.Qu Garage.Type Garage.Yr.Blt Garage.Finish Garage.Cars
1 Gd Attchd 1960 Fin 2
2 Attchd 1961 Unf 1
Garage.Area Garage.Qual Garage.Cond Paved.Drive Wood.Deck.SF
1 528 TA TA P 210
2 730 TA TA Y 140
Open.Porch.SF Enclosed.Porch X3Ssn.Porch Screen.Porch Pool.Area Pool.QC
1 62 0 0 0 0
2 0 0 0 120 0
Fence Misc.Feature Misc.Val Mo.Sold Yr.Sold Sale.Type Sale.Condition
1 0 5 2010 WD Normal
2 MnPrv 0 6 2010 WD Normal
SalePrice
1 215000
2 105000
Tijekom analize postupno ćemo se upoznavati s podacima putem jednostavnih pitanja kao što je: “Koja je dimenzija podataka?”
dim(ames)
[1] 2930 82
Koliko varijabli se nalazi u skupu podataka?
dim(ames)[2] #ili ncol(ames)
[1] 82
Koliko opservacija (promatranja) sadrži skup podataka?
dim(ames)[1] # ili nrow(ames)
[1] 2930
Koja je klasa objekta?
class(ames)
[1] “data.frame”
Imena varijabli u skupu podataka.
names(ames)
[1] “Order” “PID” “MS.SubClass”
[4] “MS.Zoning” “Lot.Frontage” “Lot.Area”
[7] “Street” “Alley” “Lot.Shape”
[10] “Land.Contour” “Utilities” “Lot.Config”
[13] “Land.Slope” “Neighborhood” “Condition.1”
…
Sumarna statistika po varijablama - ispis za prvih 6 varijabli.
summary(ames)[,1:6]
Order PID MS.SubClass MS.Zoning
Min. : 1.0 Min. :5.263e+08 Min. : 20.00 A (agr): 2
1st Qu.: 733.2 1st Qu.:5.285e+08 1st Qu.: 20.00 C (all): 25
Median :1465.5 Median :5.355e+08 Median : 50.00 FV : 139
Mean :1465.5 Mean :7.145e+08 Mean : 57.39 I (all): 2
3rd Qu.:2197.8 3rd Qu.:9.072e+08 3rd Qu.: 70.00 RH : 27
Max. :2930.0 Max. :1.007e+09 Max. :190.00 RL :2273
RM : 462
Lot.Frontage Lot.Area
Min. : 21.00 Min. : 1300
1st Qu.: 58.00 1st Qu.: 7440
Median : 68.00 Median : 9436
Mean : 69.22 Mean : 10148
3rd Qu.: 80.00 3rd Qu.: 11555
Max. :313.00 Max. :215245
NA’s :490
Koliki je interkvartilni raspon varijable SalePrice - iznosa za koji je prodana nekretnina?
IQR(ames$SalePrice)
[1] 84000
Kolike su minimalne i maksimalne cijene nekretnina prodanih u gradu Amesu od 2006. do 2010. godine?
min(ames$SalePrice)
[1] 12789
max(ames$SalePrice)
[1] 755000
Koliko različitih vrijednosti nalazimo u varijabli Garage.Cond? Za koliko opservacija nemamo podatak u varijabli?
summary(ames$Garage.Cond)
Ex Fa Gd Po TA NA’s
1 3 74 15 14 2665 158
Iz ispisa sumarne statistike varijable vidljivo je da postoje 4 razine vrijednosti varijable te da za 158 opservacija u ovoj varijabli nema zapisa.
Koliki je raspon cijena nekretnina u ulicama koje imaju popločenje (Street == Pave)? Ulice su ili popločene ili pošljunčane (Pave/Grvl).
range(ames[ames$Street==”Pave”,]$SalePrice)
[1] 12789 755000
Nacrtat ćemo stem-and-leaf dijagram varijable SalePrice.
##
## The decimal point is 5 digit(s) to the right of the |
##
## 0 | 113444444
## 0 | 55555556666666666666666666666666777777777777777778888888888888888888+109
## 1 | 00000000000000000000000000000000000000000000000000000000000000000000+846
## 1 | 55555555555555555555555555555555555555555555555555555555555555555555+805
## 2 | 00000000000000000000000000000000000000000000000000000000000000000000+355
## 2 | 55555555555555555555555555555555555555555555555555555555555666666666+160
## 3 | 00000000000000000000000111111111111111111122222222222222222222222222+36
## 3 | 55555555555566666666666777777777777788888888888888889999999999
## 4 | 000000000111111112222222223334444
## 4 | 55556666777788899
## 5 | 00044
## 5 | 5566889
## 6 | 1123
## 6 |
## 7 |
## 7 | 56
Podatke ćemo prikazati i odgovarajućim grafikonima.
Cijena nekretnine – prikaz u 30 perioda podjele
Cijena nekretnine - prikaz prema godini prodaje
Koliko je kamina (Fireplaces) imala najskuplje prodana nekretnina u 2008. godini?
prodane_08
prodane_08[which.max(prodane_08$SalePrice),]$Fireplaces
[1] 2
Izračunat ćemo kovarijancu varijabli cijene kuće SalePrice i ukupne bruto površine Gr.Liv.Area. Kovarijanca je mjera linearne povezanosti dviju varijabli.
cov(ames$Gr.Liv.Area , ames$SalePrice)
[1] 28542200
Kako je kovarijanca mjera koja ovisi o mjernim jedinicama, ne govori nam mnogo tako dugo dok je ne standardiziramo, tj. izrazimo u obliku Pearson korelacijskog koeficijenta R.
cor( ames$SalePrice, ames$Gr.Liv.Area, method=”pearson”)
[1] 0.7067799
Približno 50 posto varijabilnosti cijene nekretnine (0.7 x 0.7) može se objasniti njezinom bruto površinom. Kako bismo jednostavnije odgovorili na ovakva pitanja, nacrtat ćemo dijagram raspršenja, engl. scatterplot, za varijable SalePrice (y os) i Gr.Liv.Area (x os). Izrada linearnog modela
Budući da smo zaključili da je veza naše varijable od interesa (cijena nekretnine) u linearnoj vezi s varijablom bruto površine, želimo ovu vezu i matematički opisati na cjelokupnom rasponu mogućih vrijednosti.
Kako ne bismo stalno pisali ime objekta (ames), napravit ćemo attach objekta i od sada pisati samo imena varijabli te izračunati jednadžbu pravca - linearni model.
attach(ames)
fit
summary(fit)
Call:
lm(formula = SalePrice ~ Gr.Liv.Area)
Residuals:
Min 1Q Median 3Q Max
-483467 -30219 -1966 22728 334323
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13289.634 3269.703 4.064 4.94e-05 ***
Gr.Liv.Area 111.694 2.066 54.061
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
Residual standard error: 56520 on 2928 degrees of freedom
Multiple R-squared: 0.4995, Adjusted R-squared: 0.4994
F-statistic: 2923 on 1 and 2928 DF, p-value:
Iz rezultata linearnog modela vidljivo je da je on statistički značajan ( p-value:
Interval pouzdanosti modela - predefinirano za 95 % interval ako ne tražimo drugačije.
confint(fit) #predefinirano za 0.95
2.5 % 97.5 %
(Intercept) 6878.4845 19700.7842
Gr.Liv.Area 107.6429 115.7451
Ako želimo procijeniti vrijednost neke nekretnine poznate bruto površine, trebali bismo izraditi, temeljem naših podataka i izmjerenoj kovarijanci, optimalan matematički model (pravac) kojim ćemo moći dati procjenu vrijednosti nekretnine za sve vrijednosti površina unutar kojih smo imali mjerenja u skupu podataka. Pravac se kroz skup točaka optimalno provlači metodom najmanjih kvadrata.
Ako želimo dobiti parametre pravca (linearnog modela), pogledamo objekt fit u koji smo spremili rezultat linearnog modeliranja, rezultat funkcije lm. Kako glasi jednadžba pravca kojim bruto površinom kuće objašnjavamo njezinu cijenu?
fit
Call:
lm(formula = SalePrice ~ Gr.Liv.Area)
Coefficients:
(Intercept) Gr.Liv.Area
13289.6 111.7
Vizualizirat ćemo linearni model (pravac) koji je optimalno provučen kroz točke te ga možemo koristiti za procjenu cijene nekretnine za bilo koju kvadraturu kuće. Procijenit ćemo cijenu nekretnine za 2000 kvadratnih feeta.
new
predict(fit, new,interval=”prediction”)
fit lwr upr
1 236677.6 125809 347546.2
Slično možemo napraviti i za neku drugu varijablu od interesa ili kombinacijom više varijabli. Primjerice, koliki postotak cijene kuće možemo objasniti površinom (Gr.Liv.Area) i brojem spavaćih soba (Bedroom.AbvGr).
fit2
summary(fit2)
Call:
lm(formula = SalePrice ~ Gr.Liv.Area + Bedroom.AbvGr)
Residuals:
Min 1Q Median 3Q Max
-581397 -27840 -862 24526 329875
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 59496.236 3741.249 15.90
Gr.Liv.Area 136.361 2.247 60.69
Bedroom.AbvGr -29149.110 1372.135 -21.24
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
Residual standard error: 52620 on 2927 degrees of freedom
Multiple R-squared: 0.5664, Adjusted R-squared: 0.5661
F-statistic: 1912 on 2 and 2927 DF, p-value:
Stavljanjem dvije varijable kao objašnjavajuće, uspjeli smo objasniti 56,64 % varijabiliteta cijene nekretnina.
Zanima nas je li drugi model statistički bolji od modela samo s jednom varijablom. Provjerit ćemo omjere varijabiliteta (napraviti analizu varijance - ANOVA) koje objašnjavaju jedan i drugi model. Prisjetimo se da smo model samo s jednom varijablom (bruto površina) zvali fit, dok smo model s 2 varijable zvali fit2.
anova(fit, fit2, test=”F”)
Analysis of Variance Table
Model 1: SalePrice ~ Gr.Liv.Area
Model 2: SalePrice ~ Gr.Liv.Area + Bedroom.AbvGr
Res.Df RSS Df Sum of Sq F Pr(>F)
1 2928 9.3549e+12
2 2927 8.1052e+12 1 1.2497e+12 451.29
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
Razlika u objašnjenom varijabilitetu dva modela je statistički veoma značajna - Pr(>F) 2.2e-16.
Dodavanjem većeg broja varijabli u model vrlo vjerojatno ćemo povećavati i stupanj objašnjenog varijabiliteta varijable od interesa, ali u nekom trenutku model će postati teško objašnjiv i kompleksan, ali i doveli do tzv. overfittinga. Zato je veoma važno od velikog skupa dostupnih varijabli (u našem slučaju 82) na optimalan način odabrati onaj podskup varijabli kojim ćemo dobiti što jednostavniji model, a koji pak objašnjava u najvećoj mjeri našu varijablu od interesa.
Metode kojima se radi selekcija podskupa varijabli kompleksan je dio statistike i obrade podataka. Radi demonstracije, napravit ćemo tzv. “stepwise” selekciju varijabli. U ovu proceduru trebali bismo staviti samo varijable koje imaju linearnu povezanost s varijablom od interesa, ali u ovoj demonstraciji stavljamo skup od 15 varijabli (redom od pozicije 10 do 25). Generalno, stepwise selekcijske metode mogu biti engl. forward, backward ili stepwise, a razlikuje ih način na koji stavljaju, odnosno izbacuju varijable koje ne doprinose značajno općoj kvaliteti modela i u ovome tekstu ih nećemo detaljno objašnjavati. Statistika koju u pozadini ove metode koriste je Akaike Information Criterion (AIC) - statistička mjera koja traži model s najmanjom Root Mean Squared Error (RMSE), ali pritom penalizira ukupan broj varijabli uključenih u model.
Postoje i druge metode selekcije podskupa varijabli u modelu koje koriste druge statistike kao indikatore optimalnog modela.
# model sa svim varijablama
full
# Stepwise regression model - u oba smjera
step
Start: AIC=61695.59
SalePrice ~ Land.Contour + Utilities + Lot.Config + Land.Slope +
Neighborhood + Condition.1 + Condition.2 + Bldg.Type + House.Style +
Overall.Qual + Overall.Cond + Year.Built + Year.Remod.Add +
Roof.Style + Roof.Matl + Exterior.1st
Df Sum of Sq RSS AIC
- Utilities 2 1.7908e+09 3.8313e+12 61693
3.8295e+12 61696
- Overall.Cond 1 4.6711e+09 3.8341e+12 61697
- Condition.2 7 3.1903e+10 3.8614e+12 61706
- Land.Contour 3 2.2329e+10 3.8518e+12 61707
- Condition.1 8 4.4214e+10 3.8737e+12 61713
- Lot.Config 4 4.1502e+10 3.8710e+12 61719
- Land.Slope 2 4.3752e+10 3.8732e+12 61725
- Year.Remod.Add 1 6.5731e+10 3.8952e+12 61743
- Year.Built 1 6.9731e+10 3.8992e+12 61746
- House.Style 7 8.8164e+10 3.9176e+12 61748
- Roof.Matl 7 1.0584e+11 3.9353e+12 61761
- Exterior.1st 15 1.3284e+11 3.9623e+12 61766
- Roof.Style 5 1.0921e+11 3.9387e+12 61768
- Bldg.Type 4 4.4784e+11 4.2773e+12 62012
- Neighborhood 27 1.1457e+12 4.9752e+12 62408
- Overall.Qual 1 1.1960e+12 5.0254e+12 62490
Step: AIC=61692.96
SalePrice ~ Land.Contour + Lot.Config + Land.Slope + Neighborhood +
Condition.1 + Condition.2 + Bldg.Type + House.Style + Overall.Qual +
Overall.Cond + Year.Built + Year.Remod.Add + Roof.Style +
Roof.Matl + Exterior.1st
Df Sum of Sq RSS AIC
3.8313e+12 61693
- Overall.Cond 1 4.5076e+09 3.8358e+12 61694
+ Utilities 2 1.7908e+09 3.8295e+12 61696
- Condition.2 7 3.2008e+10 3.8633e+12 61703
- Land.Contour 3 2.2298e+10 3.8535e+12 61704
- Condition.1 8 4.4153e+10 3.8754e+12 61711
- Lot.Config 4 4.0451e+10 3.8717e+12 61716
- Land.Slope 2 4.3649e+10 3.8749e+12 61722
- Year.Remod.Add 1 6.6441e+10 3.8977e+12 61741
- Year.Built 1 7.0486e+10 3.9017e+12 61744
- House.Style 7 8.8666e+10 3.9199e+12 61746
- Roof.Matl 7 1.0591e+11 3.9372e+12 61759
- Exterior.1st 15 1.3386e+11 3.9651e+12 61764
- Roof.Style 5 1.0995e+11 3.9412e+12 61766
- Bldg.Type 4 4.4838e+11 4.2796e+12 62009
- Neighborhood 27 1.1443e+12 4.9755e+12 62405
- Overall.Qual 1 1.1961e+12 5.0274e+12 62487
# odvajamo ames iz putanje za traženje
detach(ames)
Dijagram raspršenja za varijable.
Dijagram raspršenja varijabli s dodatnim linearnim modelom.
Primijetite da su kategorijske varijable razdijeljene u tzv. dummy varijable, pa se tako svaki nivo faktora pojavljuje kao prediktor u rezultatu.
Koristeći ovih 15 varijabli, objasnili smo 78,84 % varijabiliteta varijable od interesa (cijena nekretnine).