R za novu generaciju IT stručnjaka

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 http://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  <NA>

2     2 526350040          20        RH           80    11622   Pave  <NA>

  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         <NA>      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    <NA>

2             0              0           0          120         0    <NA>

  Fence Misc.Feature Misc.Val Mo.Sold Yr.Sold Sale.Type Sale.Condition

1  <NA>         <NA>        0       5    2010       WD          Normal

2 MnPrv         <NA>        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 <- ames[ames$Yr.Sold==”2008”,]

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 <- lm(SalePrice ~ Gr.Liv.Area)

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  < 2e-16 ***

---

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: < 2.2e-16

Iz rezultata linearnog modela vidljivo je da je on statistički značajan ( p-value: < 2.2e-16), tj. cijena kuće je povezana s njezinom veličinom.

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 <- data.frame(Gr.Liv.Area  = 2000)

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 <- lm( SalePrice ~Gr.Liv.Area + Bedroom.AbvGr)

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   <2e-16 ***

Gr.Liv.Area      136.361      2.247   60.69   <2e-16 ***

Bedroom.AbvGr -29149.110   1372.135  -21.24   <2e-16 ***

---

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: < 2.2e-16

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 < 2.2e-16 ***

---

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 <- lm(SalePrice ~., data = ames[,10:25])

# Stepwise regression model - u oba smjera

step <- step(full, direction = “both”)

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

<none>                         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

<none>                         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).