Funkcia pre výpočet Lineárnej regresie

Na odhad parametrov lineárnej regresnej funkcie (jednoduchá metóda najmenších štvorcov) môžme použiť viacero metód:

1) Výpočet vektora parametrov b (cez matice)

Postup:

1. načítanie vstupných vektorov:

#nacitanie obratu:

Yt = scan()

6.1 7.3 9.6 10.2 10.1 11.3 12.2 12.5 13.2#nacitanie cenovych indexov:

Pt = scan()

103 102 100 94 98 97 98 97 96#nacitanie prijmov obyv.:

It = scan()

110 114 130 135 141 152 160 165 170

 

2. postupné počítanie podľa vzorca:

#vytvorenie pomocneho vektora, bude to prvy stlpec matice A tvoreny jednotkami

#funkcia replicate(x=1,times=9): pom = rep(1,9)#vytvorenie matice A pomocou

#funkcia array(vektory_matice,rozmery_matice):

A = array( b(pom,Pt,It) , b(9,3) )#vypisanie matice A – len kvoli kontrole:

A#transponovanie matice A:

Atran = t(A)#inverzna matica k matici A’A:

AAinv = solve(Atran %*% A)

#funkcia solve(matica) vytvori inverznu maticu, matice sa nasobia vzdy pomocou %*%
#vypocet estimatora parametrov:
b = AAinv %*% Atran %*% Yt

#vypisanie:
b

Výpis vektoru parametrov modelu:

> b

odhad parametra b0	10.43922788
odhad parametra b1	-0.13838383
odhad parametra b2	0.09476636

koeficienty regresnej funkcie:

10.43923 … b0 – lokujúca konštanta/ aditívny člen/ absolútny člen
-0.13838 … b1
0.09477 … b2

 

2) Výpočet vektora parametrov b pomocou funkcie lm()

Postup:

1. načítanie vstupných vektorov:#nacitanie obratu:

Yt = scan()
6.1 7.3 9.6 10.2 10.1 11.3 12.2 12.5 13.2#nacitanie cenovych indexov:
Pt = scan()
103 102 100 94 98 97 98 97 96#nacitanie prijmov obyv:
It = scan()
110 114 130 135 141 152 160 165 170

2. zavolanie funkcie lm():

#lm() - funkcia na výpočet lineárneho modelu
#cbind() - funkcia, ktorá spája stĺpce do matice
#ulozenie vysledkov funkcie lm(zav_velicina ~ vektor/matica_nezav_veliciny/velicin):
regresia = lm(Yt ~ cbind(Pt,It))#vypisanie vektoru parametrov modelu
regresia#vypisanie podrobnejsieho vysledku funkcie lm()
summary(regresia)Výpis vektoru parametrov modelu:>regresiaCall:
lm(formula = Yt ~ cbind(Pt, It))

Coefficients:

(Intercept)	cbind(Pt, It)Pt	cbind(Pt, It)It
10.43923	-0.13838	0.09477

koeficienty regresnej funkcie:
10.43923 … b0 – lokujúca konštanta/ aditívny člen/ absolútny člen
-0.13838 … b1
0.09477 … b2

Výpis podrobnejšieho výsledku funkcie lm():

> summary(regresia)
Call:

lm(formula = Yt ~ cbind(Pt, It))
Residuals:

Min	1Q	Median	3Q	Max
-0.50999	-0.13967	-0.06466	0.15977	0.67953

 

Coefficients:

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	10.439228	7.450424	1.401	0.2107
cbind(Pt, It)Pt	-0.138384	0.066137	-2.092	0.0813 .
cbind(Pt, It)It	0.094766	0.008787	10.785	3.76e-05 ***

---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.374 on 6 degrees of freedom
Multiple R-squared: 0.9814, Adjusted R-squared: 0.9752
F-statistic: 158 on 2 and 6 DF, p-value: 6.47e-06

Estimate – odhad parametrov
Std. Error – odha štandardných odchýlok parametrov modelu
t value – hodnota t štatistiky, ak t < t0,05[v], kde v = n-(k+1) je počet voľnosti, parameter b je nevýznamný
Pr(> |t|) – ak je hodnota menšia ako zadaná hodnota hladiny významnosti, zamietame h0 a prijímame h1
cbind(Pt, It) – funkcia, ktorá spája stĺpce do matice

 

3) Vlastná funkcia

Táto funkcia vypočíta a vypíše odhadnuté hodnoty y, vykreslí graf hodnôt y a y^,  a vypíše detailný výpis z funkcie lm (parametre, R2…)

Postup:
1. nakopírovať funkciu do R-cran
2. načítať dáta
3. zavolať funkciu reg s dvomi parametrami (vektor hodnôt y, matica hodnôt x) a následne vypísanie výsledku

1.
reg = function (Yt, m)
  {
    regr = lm (Yt~m)
    nr = nrow(m)
    nc = ncol(m)
    y = c()
    b = c()   
    for (i in 1:(nc+1))
    {
      b[i]=regr$coef[i]
    }    
    for (i in 1:nr)
    {
      ypom = 0;
      for (j in 2:(nc+1))
      {
        ypom = ypom + b[j]*m[i,(j-1)]
      }
      y[i]=b[1]+ypom
    }
    plot (Yt, type = "o", col="green",xlab="cas t",ylab="hodnoty y", main = "Linearna regresia")
    lines (y, type = "o", col="red")
    cat("Linear regresion","\n",
      "--------------------------------------","\n",
      "Estimated values of Y:","\n","\n",
      " ",y,"\n",
      "--------------------------------------","\n",
      "Summary from lm function:","\n"
    )
    summary (regr)
  }
—————————————————————————
2.
gh = c(103,102,100,94,98,97,98,97,96)
    kh = c(110,114,130,135,141,152,160,165,170)
    yh = c(6.1,7.3,9.6,10.2,10.1,11.3,12.2,12.5,13.2)
3.
reg(yh,cbind (gh,kh))

Linear regresion
 --------------------------------------
 Estimated values of Y:
 
   6.609993 7.127442 8.920471 10.22461 10.23967 11.42048 12.04023 12.65245 13.26466
 --------------------------------------
 Summary from lm function:

Call:
lm(formula = Yt ~ m)

Residuals:
     Min       1Q   Median       3Q      Max
-0.50999 -0.13967 -0.06466  0.15977  0.67953

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.439228   7.450424   1.401   0.2107    
mgh         -0.138384   0.066137  -2.092   0.0813 .  
mkh          0.094766   0.008787  10.785 3.76e-05 ***
—
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.374 on 6 degrees of freedom
Multiple R-squared: 0.9814,     Adjusted R-squared: 0.9752
F-statistic:   158 on 2 and 6 DF,  p-value: 6.47e-06