Untitled
Ekonometrija kroz R II dio
Uklonimo sve objekte iz okruzenja (isto mozemo uraditi i klikom misa).
rm(list=ls())R ima mnostvo paketa. Jedan on njih ima na raspolaganju podatke koji su primjeri za anlizu. Prvo instaliramo, a kad je jednom instaliran paket potrebno ga je pozvati. Nazvacemo objekat stednja i spremiti ga u okruzenje. Kakvi su podaci? Da vidimo strukturu str(), odnosno klase pojedinih promjenljivih.Imamo 50 observacija za po 5 promjenljivih.Podatke mozemo sagledati po kolonama. - [,1] sr numeric agregatne licne stednje - [,2] pop15 numeric % of population under 15 - [,3] pop75 numeric % of population over 75 - [,4] dpi numeric real per-capita disposable income - [,5] ddpi numeric % growth rate of dpi
#install.packages("datasets.load")
library("datasets.load")
data(LifeCycleSavings)
stednja<- data.frame(LifeCycleSavings)str(LifeCycleSavings) ## 'data.frame': 50 obs. of 5 variables:
## $ sr : num 11.43 12.07 13.17 5.75 12.88 ...
## $ pop15: num 29.4 23.3 23.8 41.9 42.2 ...
## $ pop75: num 2.87 4.41 4.43 1.67 0.83 2.85 1.34 0.67 1.06 1.14 ...
## $ dpi : num 2330 1508 2108 189 728 ...
## $ ddpi : num 2.87 3.93 3.82 0.22 4.56 2.43 2.67 6.51 3.08 2.8 ...Kao sto smo rekli, deskriptivna statistika je obavezna. Posmatramo, trazimo da li imamo nesto neobicno, posmatramo minimalne i maksimalne vrijednosti. U konretnom slucaju, tesko je nesto zakljuciti u smislu postojanja manjkavosti.
summary(LifeCycleSavings)## sr pop15 pop75 dpi
## Min. : 0.600 Min. :21.44 Min. :0.560 Min. : 88.94
## 1st Qu.: 6.970 1st Qu.:26.21 1st Qu.:1.125 1st Qu.: 288.21
## Median :10.510 Median :32.58 Median :2.175 Median : 695.66
## Mean : 9.671 Mean :35.09 Mean :2.293 Mean :1106.76
## 3rd Qu.:12.617 3rd Qu.:44.06 3rd Qu.:3.325 3rd Qu.:1795.62
## Max. :21.100 Max. :47.64 Max. :4.700 Max. :4001.89
## ddpi
## Min. : 0.220
## 1st Qu.: 2.002
## Median : 3.000
## Mean : 3.758
## 3rd Qu.: 4.478
## Max. :16.710Dakle, pogledajmo par grafickih prikaza:
plot(stednja)
hist(stednja$sr)
hist(stednja$pop15)
hist(stednja$pop75)
hist(stednja$dpi)
hist(stednja$ddpi)
Pogledajmo sada model u kojem posmatramo stednju kao zavisnu promjenljivu, a sve ostale promjenljive kao regresore.
stednja.lm <- lm(sr~pop15+pop75+dpi+ddpi, data=stednja)
summary(stednja.lm)##
## Call:
## lm(formula = sr ~ pop15 + pop75 + dpi + ddpi, data = stednja)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.2422 -2.6857 -0.2488 2.4280 9.7509
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.5660865 7.3545161 3.884 0.000334 ***
## pop15 -0.4611931 0.1446422 -3.189 0.002603 **
## pop75 -1.6914977 1.0835989 -1.561 0.125530
## dpi -0.0003369 0.0009311 -0.362 0.719173
## ddpi 0.4096949 0.1961971 2.088 0.042471 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.803 on 45 degrees of freedom
## Multiple R-squared: 0.3385, Adjusted R-squared: 0.2797
## F-statistic: 5.756 on 4 and 45 DF, p-value: 0.0007904Mozemo testirati da li koeficijenti koji se odnose na regresore imaju statisticku znacajnost. Odnosno, imamo hipotezu da \(\beta_1=\beta_2=\beta_2=\beta_4=0\) posebno - da bi u potpunosti bilo jasno sta softver izbaci kao rezultat. Kako je p vrijednost mala, mozemo odbaciti navedenu hipotezu. Ovdje smo imali F test - zajednicki smo testirali sve koeficejnte. Pojedinacne p vrijednosti se odnose na pojedine koeficijente, te omogucavaju interpretaciju rezultata u tom smislu. Dakle, uradimo istu stvar direktno (dobijamo isti rezultat kao sto je u ispisu). U sustini, uslovno govoreci, testiramo da li je vrijednost \((1-R^2)\) statisticki znacajno razlicita od nule.
TSS<-sum((stednja$sr-mean(stednja$sr))^2)
RSS<-sum(stednja.lm$res^2)
((TSS-RSS)/4)/(RSS/45)## [1] 5.7556811-pf(5.75579, 4,45)## [1] 0.0007902729Zasto u imeniocu imamo 4, odnosno 45? (n-1 stepeni slobode)
Pogledajmo nesto drugaciji model, odnosno model u kojem imamo jedan manje regresor nego u prethodnom slucaju.
stednja.lm2<- lm(sr~pop75+dpi+ddpi, data=stednja)Kako mozemo “pjeske” ispitati statisticku znacajnost pojedino koeficijenata, odnosno konkretno onog vezanog za promjenljivu pop75. Izracunajmo RSS i statistiku F testa.
sum(stednja.lm2$res^2)## [1] 797.7249(797.72-650.71)/(650.71/45)## [1] 10.166511-pf(10.16651, 1, 45)## [1] 0.002603116Mozemo sada dalje odrediti statitiku t testa, odnosno p vrijednost.
sqrt(10.167)## [1] 3.1885732*(1-pt(3.188,45))## [1] 0.002606756Dakle, odbacujemo \(H_0\) da je pop75 jednako nuli, odnosno da nije statisticki znacajan.
Napomenimo da poredenje dva “ugnjezedena” (nested) modela u ovakvim slucajevima je mozda pogodnije putem ANOVA. Analiza varijacija, polazi od varijacija…
anova(stednja.lm2, stednja.lm)Pogledajmo neke grafikone.
plot(stednja.lm, scale=0, page=1, all.terms=T)
plot(stednja.lm2, scale=0, page=1, all.terms=T)
Sto se tice samo pristupa rezultatima:
stednja.lm <- lm(sr~pop15+pop75+dpi+ddpi, data=stednja)
output<-summary(stednja.lm)
SSR <- deviance(stednja.lm)
LL <- logLik(stednja.lm)
DegreesOfFreedom <- stednja.lm$df
Yhat <- stednja.lm$fitted.values
Coef <- stednja.lm$coefficients
Resid <- stednja.lm$residuals
s <- stednja.lm$sigma
RSquared <- stednja.lm$r.squared
CovMatrix <- s^2*stednja.lm$cov # ili vcov()
aic <- AIC(stednja.lm)