Jūs esate čia

Tiesinė regresija

Statistikos ir ekonomikos tyrimuose neretai siekiama ištirti reiškinių priklausomybę, pavyzdžiui:
* Kaip parduodamų ledų kiekis priklauso nuo oro temperatūros?
* Ar bakalauro darbo įvertinimas turi įtakos darbo užmokesčiui?
* …

Šias priklausomybes galima nustatyti regresinės analizės pagalba. Sudaromas regresinis modelis - statistinis modelis, padedantis nustatyti vieno kintamojo priklausomybę nuo kito. Tarkime, kad mūsų duomenis sudaro k porų : (xi, yi),i = 1, ..,k. Bendriausias regresijos modelis, siejantis kiekvieną x ir y:

Y=α + β X + ϵi,

čia Y = (y1, y2, ..,yk)T - priklausomo kintamasojo reikšmių vektorius, **T -transponuota,
X = (x1, x2, …, xk)T - nepriklausomo kintamamojo reikšmių vektorius,
α ir β - konstantos.
ϵi, i = 1, …, k - atsitiktinė paklaida.

Regresijos tikslas yra surasti lygtį $\hat{Y}=\hat{\alpha}+\hat{\beta}X$, kuri „geriausiai“ atitiktų turimus duomenis. Vienas iš būdų tai padarytiyra naudojantis Mažiausių kvadratų metodu: išbrėžiama tiesė, kuri minimizuoja paklaidų kvadratų sumą nuo tiesinio regresinio modelio. α yra konstanta, o β - krypties koeficientu.

Tiesinė daugialypė regresija tiriama, kai vieną priklausomąjį kintamąjį Y = (y1, y2, ..,yn)T sieja tiesinė priklausomybė su k nepriklausomų kintamųjų X = (X1, X2, ..,Xk)T. Tuomet tiesinė daugialypė regresija išreiškiama:

$$Y={\alpha}+{\beta}_1X_1+{\beta}_2X_2+…+{\beta}_kX_k+e={\alpha}+\vec{\beta}X+{\epsilon_i} ,$$

čia X nepriklausomų kintamųjų matrica, kurią sudaro:
$\vec{X_1}=(x_{11},x_{12},…,x_{1n})$,
$\vec{X_2}=(x_{21},x_{22},…,x_{2n})$,
… ,
$\vec{X_k}=(x_{k1},x_{k2},…,x_{kn})$;
α, $\vec{\beta}=({\beta}_1,{\beta}_2,…,{\beta}_k)^T$ yra nežinomos konstantos,
ϵi, i = 1, …, n - atsitiktinė paklaida.

Tikslas - regresijos tiesės lygtis, t.y. rasti tokius α ir β įverčius $\hat{\alpha}$, $\hat{\beta}$, kad

$$\hat{Y}=\hat{\alpha}+\hat{\beta}_1X_1+\hat{\beta}_2X_2+…+\hat{\beta}_kX_k ,$$

kad liekamosioji paklaida ϵi visur būtų mažiausia: ${\epsilon_i}=Y-\hat{Y}=Y-(\hat{\alpha}+\hat{\beta}X)$ , kiekvienam i = 1, 2, …, n.

Čia taip pat įvertinimas daromas naudojantis mažiausių kvadratų metodu.

Viena pagrindinė tiesinių regresinių modelių nauda - priklausomojo kintamojo prognozė. T.y. atradę tiesinę priklausomybę, galime nustatyti, kaip keisis Y reikšmė, pakeitus nepriklausomąjį kintamąjį X.

Šiame įraše pateikiamas pavyzdys kaip atliekama išsami tiesinė regresija naudojantis statistine programa R.

Pavyzdys

Turimi duomenys - nekilnojamojo turto (gyvenamųjų butų) kainos ir keletas jų charakteristikų. Tikslas - sumodeliuoti tiesinę regresiją, kuri padėtų prekybos agentams geriau suprasti kokią įtaką skirtingos charakteristikos daro galutinei būsto kainai.

1.Duomenys:

1.1. Susipažįstama su duomenimis:

#instaliuojamos analizei reikalingos bibliotekos
library(car)
library(knitr) 
datafull<-read.csv2("data.csv") 
head(datafull) #atspausdinami duomenys
##   kaina plotas aukstas garsoIzoliacija silumosLaidumas
## 1 55800     73       8        7.815161        7.606684
## 2 52700     73       1       11.176124       10.990996
## 3 53400     66       6        8.001061        8.460651
## 4 71600     86       2       10.421210       10.475897
## 5 46700     54       5        7.003640        6.687059
## 6 61200     75       4       11.187534       10.287033
##   atstumasIkiPrekybosCentro
## 1                0.80853061
## 2                0.08585002
## 3                0.68803387
## 4                0.72644151
## 5                1.14492192
## 6                0.88719906

Pateikiamos charakteristikos: plotas, aukštas, garsoIzoliacija, silmosLaidumas, atstumasIkiPrekybosCentro. Svarbu patikrinti nepriklausomų kintamųjų, šiuo atveju - charakteristikų, koreliaciją. Koreliacijos atveju, tie patys kintamieji, turi labai panašią įtaką kainai. Tad pirmiausia stebime kaip atrodo turimi duomenys, ar nėra akivaizdžių koreliacijų:

1

1.2. Pastebima, kad 2 rodikliai - garsoIzoliacija ir silumosLaidumas - tiesiškai susiję, todėl tikrinama galima duomenų koreliacija.

print(cor(datafull), digits = 3)
##                            kaina plotas aukstas garsoIzoliacija
## kaina                     1.0000 0.8567  0.1209          0.2517
## plotas                    0.8567 1.0000  0.0443          0.0635
## aukstas                   0.1209 0.0443  1.0000         -0.1377
## garsoIzoliacija           0.2517 0.0635 -0.1377          1.0000
## silumosLaidumas           0.2390 0.0585 -0.1257          0.9536
## atstumasIkiPrekybosCentro 0.0899 0.0732 -0.0115          0.0338
##                           silumosLaidumas atstumasIkiPrekybosCentro
## kaina                              0.2390                    0.0899
## plotas                             0.0585                    0.0732
## aukstas                           -0.1257                   -0.0115
## garsoIzoliacija                    0.9536                    0.0338
## silumosLaidumas                    1.0000                    0.0336
## atstumasIkiPrekybosCentro          0.0336                    1.0000

Kadangi koreliacijos koeficientas tarp aptartų rodiklių labai didelis (arti 1) - rodikliai labai susiję, t.y. abu apibūdina buto izoliaciją. Tai reikškia, kad koeficientų įverčiai tikriausiai bus paslinkti. Siekiant panaikinti šią problemą reikia arba panaikinti vieną iš kintamųjų arba juos kažkaip sujungti.
Sujungti (vienas iš būdų - skaičiuoti vidurkį) nepasirenkama, nes nenustatyta ar abiejų koeficienų duomenis apibūdina tie patys matavimo vienetai. Išvada: ištrinami garsoIzoliacijos duomenys. SilumosLaidumas yra svarbesnis, daro įtaką gyventojų išlaidoms, t.y mokesčiams už šildymą.

datafull$garsoIzoliacija<-NULL

1.3. Tikrinamos išskirtys
Viena iš priežasčių, kodėl kuriami modeliai nėra tikslūs - duomenyse yra išskirčių. Jos iškreipia pagrindines duomenų charakteristikas, tokias kaip vidurkis, dispersija. Siekiant sukurti kuo tikslesnę regresiją, išskirtis rekomenduojama pašalinti.
Sukuriamas modelis modIsskirtims ir braižomas kvantilių grafikas Q-Q plot, kuris parodo reikšmių pasiskirstymą:

  1. modIsskirtims<-lm(kaina~+plotas+aukstas+silumosLaidumas+atstumasIkiPrekybosCentro, datafull)
  2. qqPlot(modIsskirtims, id.n=2)

2

## 254 253 
##   1 254

Pastebima, kad duomenyse yra išskirčių, tad duomenys tikrinami pagal outlierTest diganostiką:

outlierTest(modIsskirtims)
##      rstudent unadjusted p-value Bonferonni p
## 253 14.486648         7.1195e-35   1.8084e-32
## 254 -6.882738         4.7778e-11   1.2136e-08

Rezultatas: pagal Bonferonni p korekciją (p-value<0.05) - eilutėse 253 ir 254 stebimos išskirtys. Eilutes su duomenimis, kuriuose yra išskirtys - pašaliname. Sukuriami nauji duomenys pavadinimu „data“, kurie bus naudojami tolimesnei analizei.

data<-datafull[-c(253,254),] 

2.Modelio kūrimas

2.1. Sukuriamas tiesinis modelis, kuriame panaudojami visi turimi kintamieji

mod1<-lm(kaina~+plotas+aukstas+silumosLaidumas+atstumasIkiPrekybosCentro, data)
kable(summary(mod1)$coef, digits=3)
Estimate Std. Error t value Pr(> t )
(Intercept) 7966.669 977.862 8.147 0.000
plotas 599.371 12.130 49.410 0.000
aukstas 318.658 63.493 5.019 0.000
silumosLaidumas 527.454 47.522 11.099 0.000
atstumasIkiPrekybosCentro 85.454 81.381 1.050 0.295

Iš summary lentelėje pateiktos p-value galime teigti, kad koeficientas prie kintamojo „atstumasIkiPrekybosCentro“ yra nereikšmingas (H0: b4=0 - priimame), todėl jis pašalinamas.

2.2. Atnaujinamas mod2 ir tikrinamas jo kintamųjų reikšmingumas

mod2<-update(mod1, kaina~plotas+aukstas+silumosLaidumas)
kable(summary(mod2)$coef, digits=3)
Estimate Std. Error t value Pr(> t )
(Intercept) 8035.795 975.845 8.235 0
plotas 600.334 12.098 49.622 0
aukstas 318.225 63.504 5.011 0
silumosLaidumas 528.818 47.514 11.130 0

Kadangi likę įverčiai reikšmingi - tikriname modelio tikslumą. Toliau tikrinamos prielaidos, kurias turi atitikti „geras“ modelis.

3.Modelio patikimumo parametrai

3.1. „R-Squared“:

(summary(mod2))$r.squared
## [1] 0.915869

R2 yra statistinis matas, kuris parodo, kaip duomenys yra „arti“ regresinės tiesės. R2 reikšmė yra intervale [0,1], tad kuo matas arčiau 1 - tuo regresinė tiesė labiau prisitaikiusi prie duomenų. Be to, šis rodiklis taip pat labai naudingas, kai siekiama palyginti keletą modelių.

3.2.Multikolinearumo tikrinimas
Statistikoje multikolinearumas yra fenomenas, kuris nurodo, kad regresijoje 2 ar daugiau kintamųjų yra koreliuoti. Šios problemos egzistavimą tikriname naudojant Variance Inflation Factor:

vif(mod2)
##          plotas         aukstas silumosLaidumas 
##        1.005730        1.017733        1.019520

Multikolinearumo problemos nėra. Taip teigiama, nes egzistuojančią vidinę koreliaciją parodo koeficientai, kurių reikšmė yra daugiau nei 10 (VIF sąlyga).

  • Jeigu modelio kintamieji turi multikolinearumo problemą, galimos regresijos pasekmės:
  1. sulaukiame, kad kintamieji yra nereikšmingi, nors žinome, kad jų bendras poveikis tikrai nėra nulinis;
  2. koreliuotų kintamųjų įgyti koeficientai - priešingų reikšmių, t.y. jie vienas kitą kompensuoja. Todėl prarandama aiški interpretacija.
  3. modelio interpretacija tampa neaiški: jei koreliuotą kintamąjį padidiname vientetu, nebegalime tiksliai nusakyti, kaip keičiasi priklausomas kintamasis.

3.3. Homoskedastiškumo tikrinimas
Homoskedastiškumas parodo, kad regresijos atsitiktinės paklaidos yra visur vienodos. Išbrėžiamos sukurto modelio paklaidos:

  1. par(mfrow=c(1,1))
  2. plot(mod2$res~mod2$fitted, main="Paklaidų išsibarstymas pagal reikšmes", col=2, ylab="mod2 paklaidos", xlab="Butų kainos")
  3. abline(0,0)

3

Kaip matome, paklaidos pasiskirsčiusios pagal vidurkį 0, daugiausia intervale (-5000,5000). Todėl spėjama, kad modelis - homoskedastiškas. Tai dar tikrinama pagal Breuch-Pagan testą (ncvTest)

ncvTest(mod2)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 3.218568    Df = 1     p = 0.0728073

Kadangi p>0.05, vadinasi, priimame testo H0 hipotezę, kuri teigia, kad modelis homoskedastiškas - visų stebėjimų paklaidų dispersija yra konstanta.

  • Jei homoskedastiškumas nėra tenkinamas, tai formaliai didelės blogybės nėra - įverčiai ir toliau liktų nepaslinkti ir suderinti. Visgi reiktų būti atidiem ir papildomai pasvarstyti tris svarbius niuansus:
  1. Nehomogeniškumas gali indikuoti modelio neadekvatumą. Pvz. gal liko nepastebėta netiesinė sąveiką? O gal liko didelių išskirčių?
  2. Statistinio reikšmingumo tikrinimui naudojama kovariacijų matrica gali būti neteisinga.
  3. Nors įverčiai nepaslinkti, bet dabar negalime teigti, kad jie yra efektyvūs - tokiu atveju gali būti ir geresnių vertinimo būdų.

3.4. Tikrinamas liekanų normalumas

4

Mėlyna linija šioje histogramoje parodo liekanų tankio pasiskirstymą. Ji primena varpo formą, tad spėjama, kad liekanos pasiskirsčiusios normaliai. Prielaida dar tikrinama Shapiro testu:

shapiro.test(mod2$res) 
## 
##  Shapiro-Wilk normality test
## 
## data:  mod2$res
## W = 0.9898, p-value = 0.07272

Pagal shapiro.test p-value>0.05 priimame H0: liekanų paklaidos yra normalios.

  • Liekanų normalumą svarbu patikrinti tam, kad nustatytume, ar iš duomenų ištraukta visa naudinga informacija ir panaudota modelyje. Pavyzdžiui, ar atrasta duomenų tendencija įvertinta modelyje. Jeigu atsakymas teigiamas, paklaidos yra baltasis triukšmas (išsidėsčiusios nepriklausomai ir jų vidurkis yra 0):
  1. plot(residuals(mod2), type="l", main="Modelio mod2 paklaidos", ylab="mod2 paklaidos", xlab="Stebėjimas", col=100)
  2. abline(0,0)

5

Jeigu sąlyga netenkinama - modelis dar nėra pakankamai „geras“, tad verta pasvarstyti jo tobulinimą.

3.5. Tikinama liekanų autokoreliacija
Geram modeliui taip pat svarbu, kad liekanos nebūtų susiję tarpusavyje. Tai tikrina Durbin-Watson testas:

durbinWatsonTest(mod2)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.05475525      1.885685   0.358
##  Alternative hypothesis: rho != 0

Kadangi testo p-value>0.05, H0: nėra koreliacijos tarp liekanų, priimame.

  • Jeigu testo nulinė hipotezė atmetama - įverčiai vis dar yra nepaslinkti ir suderinti, tačiau tampa nebeefektyvūs. Pagrindinė liekanų autokoreliacijos problema - netikslumai prognozuojant.

Vadinasi, kadangi modelis atitinka šiuos kriterijus, galime teigti, kad teisingas modelis kaina=plotas+aukstas+silumosLaidumas*. Atitinkamai įverčiai parodo kaip pakitus vienam vienetui iš šių charakteristikų, pakinta kaina: ***kaina=8035.8 + 600.33 * plotas + 3.18.23 * aukstas + 528.82 * silumosLaidumas + *e*** Interpretacija: 1 kv.m padidėjimas padidina kainą 600,33eurais, kiekvienas papildomas aukštas kainą padidina 318,23 eurais, 1 vnt. šilumos laidumo koeficiento padidėjimas padidina kainą 528,82 eurais. Laisvasis narys - 8035,80 eurai gali būti papildomas išlaidos pvz.:notaro paslaugos, buto lokacija, viešojo transporto prieinamumas ir pan.

Trumpai apie modelio „gerinimą“

Sudarytas modelis mod2 - geras, ir tinka pradedantiems mokytis regresinės analizės. Tačiau, duomenys neretai gali būti ranginiai, kas pastebima ir atidžiau išanalizavus mūsų turimus duomenis. Galime pastebėti, kad kintamasis aukstas - yra ranginis kintamasis. Tai galime matyti iš liekanų boxplot grafiko:

  1. boxplot(mod2$residuals ~ data$aukstas, col = "lightgray", main="mod2 liekanų boxplot pagal aukštą", ylab="Liekanos", xlab="Aukštas")
  2. abline(0,0)

6

Iš grafiko matome, kad pirmo aukšto kaina neatitinka regresinės tiesės lygties, t.y. paklaidos neigiamos, jų vidurkis nėra lygus 0. Tai reiškia, kad pirmo aukšto kainos yra gerokai mažesnės ir neatitinka sukurto modelio. Siekiant į tai atsižvelgti, galima kurti modelį, kuriame išskiriamas pirmas aukštas:

aukstas1<-data$aukstas==1
mod3<-lm(kaina~plotas +aukstas1+silumosLaidumas,data=data)
kable(summary(mod3)$coef,digits=2)
Estimate Std. Error t value Pr(> t )
(Intercept) 10235.63 748.04 13.68 0
plotas 598.84 9.83 60.94 0
aukstas1TRUE -5440.16 422.84 -12.87 0
silumosLaidumas 531.42 38.39 13.84 0

Gautame modelyje aukstas1TRUE yra kintamasis, kuris galioja tik tuo atveju, jei pasitinkto buvo aukštas yra pirmas. Kitu atveju - jis modelyje įgyja reikšmę 0. Matoma, kad visi kintamieji reikšmingi. Tikrinama ar šis mod3 yra tikslesnis pagal R^2 ir AIC kriterijus:

aic<-c(AIC(mod2), AIC(mod3))
r.sq<-c((summary(mod2))$r.squared, (summary(mod3))$r.squared)
pav<-c("mod2", "mod3")
tab<-data.frame(pav,aic,r.sq)
tab
##    pav      aic      r.sq
## 1 mod2 4673.994 0.9158690
## 2 mod3 4569.454 0.9444361

mod3 AIC mažesnis, R^2 koeficientas - didesnis, lyginant su mod2. Tad tikėtina, kad modelis mod3 yra tikslesnis už mod2.

Tokiu atveju, gauta lygtis yra kaina = 10235.63 + 598.84*plotas − 5440.16 * 1aukštas + silumosLaiduma**s*.

Komentuoti

Basic HTML

  • Web puslapių adresai ir el. pašto adresai automatiškai tampa nuorodomis.
  • Tags allowed: a, em, strong, u, s, cite, code, blockquote, ol, ul, li, dl, dt, dd, pre, p, br
  • Mathematics inside the configured delimiters is rendered by MathJax. The default math delimiters are $$...$$ and \[...\] for displayed mathematics, and $...$ and \(...\) for in-line mathematics.
  • Syntax highlight code surrounded by the <pre class="brush: lang">...</pre> tags, where lang is one of the following language brushes: php, python, r, sass, sql, vb.
  • Linijos ir paragrafai atskiriami automatiškai

Markdown

  • You can enable syntax highlighting of source code with the following tags: <code>, <blockcode>, <c>, <cpp>, <drupal5>, <drupal6>, <java>, <javascript>, <php>, <python>, <r>, <ruby>. The supported tag styles are: <foo>, [foo].
  • Quick Tips:
    • Two or more spaces at a line's end = Line break
    • Double returns = Paragraph
    • *Single asterisks* or _single underscores_ = Emphasis
    • **Double** or __double__ = Strong
    • This is [a link](http://the.link.example.com "The optional title text")
    For complete details on the Markdown syntax, see the Markdown documentation and Markdown Extra documentation for tables, footnotes, and more.
  • Typographic refinements will be added.
  • Mathematics inside the configured delimiters is rendered by MathJax. The default math delimiters are $$...$$ and \[...\] for displayed mathematics, and $...$ and \(...\) for in-line mathematics.

Plain text

  • HTML žymės neleidžiamos.
  • Web puslapių adresai ir el. pašto adresai automatiškai tampa nuorodomis.
  • Linijos ir paragrafai atskiriami automatiškai
CAPTCHA
This question is for testing whether or not you are a human visitor and to prevent automated spam submissions.
Target Image

Theme by Danetsoft and Danang Probo Sayekti inspired by Maksimer