Základy R

RStudio http://www.rstudio.org
Instalace: http://www.rstudio.com/ide/download/
Jedna z mnoha knih pro zacátecníky https://cran.r-project.org/doc/contrib/Paradis-rdebuts_en.pdf

Jednoduché výpocty

\(\textbf{Sedivý odstín predstavuje kod}\)

\(\textbf{Sedivý odstín s ## predstavuje provedení kodu (to co se vám zobrazí v konzoli)}\)

3+3

## [1] 6

Priradíme souctu promennou

soucet=3+3

Promenná se nevytiskne

Pokud chceme znát její hodnotu(y) napíseme její název

soucet

## [1] 6

Funkce

V programu R je nadefinována celá rada funkcí.

\(y =f(x)\) obecný zápis funkce

\(y =NázevFunkce(argumenty)\)

odmocnina=sqrt(10)  # zde je název funkce sqrt() a argument je hodnota 10.

odmocnina2=sqrt(soucet) # zde je argumentem funkce jiz název jiné promenné, prípadne jiné funkce

V totmo kurzu budeme casto generovat náhodná císla. K tomu nám pomùze napríklad funkce \(rnorm\). Chceme vygenerovat 1000 hodnot normálního rozdelení \(\epsilon \sim N(10,400)\)

epsilon=rnorm(1000,10,20)

Takto jsme vygenerovali 1000 hodnot, které se rídí normálním rozdelením se strední hodnotou 10 a rozptylem 20²!!

Nápovedu o argumentech funkce a co umí, získáme pomocí otazníku pred názvem funkce.

?rnorm

Vygenerujte 1000 hodnot ze Studentova rozdelení s 100 stupni volnosti.

Rko má výhodu v tom, ze témer vse najdete na internetu (via)

package

Program R obsahuje tzv. Core function. Pokud chceme vyuzít nejakou funkci mimo “core” musíme si jí bud naprogramovat a nebo vyuzijeme jiz naprogramované funkce, které nekdo vytvoril za nás. Tyto funkce jsou ve forme balíkù \(packages\), které musíme stáhnout.

install.packages("ggplot2")   # stáhneme a nainstalujeme balík funkcí, který se jmenuje $ggplot2$. Jeho vlastnosti najdete na internetu

Dále musíme vyvolat daný balík

library(ggplot2)  #pri kazdém novém spustení, je treba vyvolat balíky, které budeme potrebovat

## Warning: package 'ggplot2' was built under R version 3.2.5

Opakování statistiky

Ve statistice, ekonometrii deláme úsudky o chování základního souboru. Základní soubor:

vsichni zamestannci v cR
Vsichni studenti
atd.

Úsudky, které mùzeme napríklad delat:

prùmerná mzda je vetsí jak 20 000Kc
je rozptyl výsky u zen vetsí jak 5cm?

V ekonometrii nás budou zajímat prevázne vztahy mezi velicinami tzv. data generation process DGP

závisí mzda na vzdelání
jsou peníze neutrální
existuje vztah mezi inflací a nezamestnaností

Základní soubor (populace)

Neznáme základní soubor, ale pro pochopení pochopení, si jej vytvoríme:

zakladni_soubor=rnorm(1000,180,20)   #vytvoríme si základní soubor o 1000 jedincích
                                    #se strední hodnotou 180 a rozptylem 400

Prestoze jsme si soubory nazvali jako základní, tak kazdý z nás bude mít jiné hodnoty. Nekdy chceme,aby se hodnoty nemenily. Napr. já pri kontrole úkolù.

set.seed(353)        # stejné hodnoty 
zakladni_soubor2=rnorm(1000,180,20)
head(zakladni_soubor2)    # uvidíme nekolik prvních hodnot, uzitecné u velkých souborù

## [1] 172.6415 158.8621 174.5699 201.8297 162.7445 191.6623

set.seed(357)        # stejné hodnoty 
vyber1=sample(zakladni_soubor2,50)  # provedeme výber o velikosti 50 ze základního souboru

hist(vyber1)      # tento výber si zobrazíme pomocí histogramu.

Vyuzití balíku ggplot2

qplot(vyber1,               # název souboru obsahující vstupní data. 
      geom="histogram",
      binwidth = 10,  
      main = "Histogram výberového souboru", 
      xlab = "hodnoty",  
      fill=I("blue"), 
      col=I("red"))        # zbylé argumenty funkce zjistíme pomocí ?qplot

mean(vyber1)

## [1] 183.0523

summary(vyber1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   140.6   170.9   180.3   183.1   196.9   225.4

Simulace

V reálném svete neznáme celou populaci, tedy ani parametry a náhodnou slozku regresní rovnice. Nyní vsak budeme predpokládat, ze známe jak parametry, tak rozdelení náhodné slozky. Predpokládejme, ze máme jednoduchý regresní model:

\(y=\beta_0+\beta_1x+\epsilon\)

kdy \(\beta_0=0.5\) \(\beta_1=3\) a \(\epsilon \sim N(0,900)\)

Nasím úkolem je vygenerovat hodnoty \(y\).

set.seed(22)       # díky této fci budeme mít vsichni stejné hodnoty 
n=100           # velikost simulované populace
x=seq(from=1, to=n)   # vytvoríme nezávisle promennou x jako sekvenci od 1 do n
e = rnorm(n, 0, 30)             #vygenerujeme n-hodnot s normálním rozdelením, E(e)=0 a Var(e)=400
y = 0.5 + 3 * x + e            #nakonec sestavíme populacní regresní funkci
plot(x, y)                    #data si zobrazíme ve grafu 
abline(a=0.5,b=3)             #do grafu vykreslíme prímku se zadanými parametry

Co predstavuje prímka v grafu?
Proc nelezí body na prímce?
Vertikální vzdálenost mezi bodem a prímkou predstavuje?
Zmente rozptyl náhodné slozky, k cemu dochází?

populace =data.frame(y,x) #hodnoty z vygenerovane "populace" si ulozíme, k tomu slouzí príkaz data.frame()

Regrese

Nyní je nasím cílem odhadnout “neznámé” parametry \(\beta\)

Provedeme odhad \(y=\beta_0+\beta_1x+\epsilon\) pomocí metody nejmensích ctvercu

regp=lm(y~x,populace) # odhad lineárního regresního modelu pomocí OLS (MNC)
regp

## 
## Call:
## lm(formula = y ~ x, data = populace)
## 
## Coefficients:
## (Intercept)            x  
##       6.193        2.937

Co ríká výsledek? Chybí nám neco?

Výsledek je zapsán v tzv. listu, coz je datový formát

summary(regp)

## 
## Call:
## lm(formula = y ~ x, data = populace)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -52.904 -21.189  -5.778  17.854  94.807 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.1930     5.9741   1.037    0.302    
## x             2.9370     0.1027  28.596   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 29.65 on 98 degrees of freedom
## Multiple R-squared:  0.893,  Adjusted R-squared:  0.8919 
## F-statistic: 817.7 on 1 and 98 DF,  p-value: < 2.2e-16

Proc není výsledek stejný s populacními parametry? Vsak jsme odhadovali z dat celé populace.

Výberový soubor

Provedeme náhodný výber z populace.

set.seed(50)               # zajístí, ze vsichni budeme mít stejný náhodný výber 
vyber = sample(n, 50)     #funkce sample provádí náhodný výber 
xs = x[vyber]  # Vybraná pozorování ulozíme jako xs a ys.
ys = y[vyber]
vyber1=data.frame(y=ys,x=xs)   #data si ulozíme do data.frame

Provedeme odhad.

regs=lm(y~x,vyber1) # odhad lineárního regresního modelu pomocí OLS (MNc) z výberových dat
summary(regs)  # zobrazí tabulku výsledkù

## 
## Call:
## lm(formula = y ~ x, data = vyber1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -50.462 -22.779  -2.128  21.564  58.174 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.6727     7.4921   0.624    0.536    
## x             2.8909     0.1435  20.145   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 27.67 on 48 degrees of freedom
## Multiple R-squared:  0.8942, Adjusted R-squared:  0.892 
## F-statistic: 405.8 on 1 and 48 DF,  p-value: < 2.2e-16

V grafu si sestrojíme výberovou regresní prímku

plot(xs,ys)           # zobrazíme si data z výberového souboru 
abline(a=0.5,b=3)     # do grafu si zaneseme jak vypadá vztah mezi y a x v populaci 
abline(a=regs$coefficients[1],b=regs$coefficients[2],col = "red", lty = 3)  # V listu muzeme "listovat" pomocí $

abychom zjistili co vsechno umí daná funkce muzeme do konzole napsat \(?nazev funkce\), tedy napríklad \(?abline\)
co predstavuje vzdálenost mezi cervenou prímkou a bodem?

U výberu zmeníme set.seed(100) a odhad provedeme znovu. Novou regresní prímku si opet zaneseme do grafu.

Co se tedy mení s novým výberem? Parametry \(\beta\), nebo \(\hat{\beta}\) ?

Odhad parametru

Pro jednoduchou regresi platí, ze odhad parametru \(\beta_1\) je dán:

\(b_1=\dfrac{Cov(x,y)}{Var(x)}\)

b1=cov(xs,ys)/var(xs)

Co potrebujeme dále znát? Stací pouze hodnota odhadu?

Rozptyl odhadu b1 pro jednoduchou regresi je dán:

\(Var(b_1)=\dfrac{\hat{\sigma}^2}{\sum\limits_{i=1}^n (x_i-\bar{x})^2}\)

hsigma=var(regs$residuals)  # rozptyl residuí je odhadem rozptylu náhodné slozky 
suma=sum((xs-mean(xs))^2)
varb=hsigma/suma
sqrt(varb)

## [1] 0.142032

Vícerozmerná regrese

\(y=\beta_0+\beta_1x_1+\beta_2x_2+\beta_3 \log{x_3}+\epsilon\)

kdy \(\beta_0=0.5\) \(\beta_1=1.6\) \(\beta_2=0.5\) \(\beta_3=2.2\) a \(\epsilon \sim N(0,900)\)

Nasím úkolem je vygenerovat hodnoty \(y\).

set.seed(22)       # díky této fci budeme mít vsichni stejné hodnoty 
n=100           # velikost simulované populace
x1=seq(from=1, to=n)
x2=rnorm(n,100,10)
x3=runif(n,5,20)
e= rnorm(n, 0, 30)             #vygenerujeme n-hodnot s normálním rozdelením, E(e)=0 a Var(e)=400
y = 0.5 + 1.6 * x1+0.5*x2+2.2*x3 + e            #nakonec sestavíme populacní regresní funkci

Ulozte data a zobrazte si je.

data_regrese = data.frame(y,x1,x2,x3)

reg_vc=lm(y~x1+x2+x3,data_regrese) # odhad lineárního regresního modelu pomocí OLS (MNc) z výberových dat
summary(reg_vc)  # zobrazí tabulku výsledkù

## 
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = data_regrese)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -63.522 -17.761   3.404  17.351  72.425 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -40.60505   30.64647  -1.325  0.18833    
## x1            1.70016    0.09977  17.041  < 2e-16 ***
## x2            0.81545    0.29679   2.748  0.00717 ** 
## x3            2.23890    0.68033   3.291  0.00140 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.73 on 96 degrees of freedom
## Multiple R-squared:  0.7637, Adjusted R-squared:  0.7563 
## F-statistic: 103.4 on 3 and 96 DF,  p-value: < 2.2e-16

reg_vc2=lm(log(y)~x1+x2+x3,data_regrese) # odhad lineárního regresního modelu pomocí OLS (MNc) z výberových dat
summary(reg_vc2)  # zobrazí tabulku výsledkù

## 
## Call:
## lm(formula = log(y) ~ x1 + x2 + x3, data = data_regrese)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.76741 -0.13477  0.03836  0.14255  0.46462 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.4683639  0.2401388  14.443  < 2e-16 ***
## x1          0.0120292  0.0007818  15.387  < 2e-16 ***
## x2          0.0065823  0.0023256   2.830 0.005662 ** 
## x3          0.0183669  0.0053309   3.445 0.000847 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2251 on 96 degrees of freedom
## Multiple R-squared:  0.7291, Adjusted R-squared:  0.7206 
## F-statistic: 86.11 on 3 and 96 DF,  p-value: < 2.2e-16