U ovoj lekciji ćemo se baviti jednostavnom linearnom regresijom (Simple Linear Regression). Njen cilj je opisivanje i kvantifikovanje veze između dva obeležja \(X\) - prediktor, promenljiva objašnjenja, nezavisna promenljiva \(Y\) - odgovor, ishod, zavisna promenljiva

Ako nam je dat dvodimenzioni uzorak \((x,y)\) sa grafika raspršenja (scatter plot) dobijamo ideju o postojanju i obliku zavisnosti. Učitavamo bazu u kojoj se nalaze podaci o broju umrlih (na 10 miliona ljudi) od raka kože u 49 država SAD-a. Kao prediktor koristimo kolonu Lat u kojoj se nalaze vrednosti geografske širine centra svake od država.

skinc=read.table("E:\\STAT501_Lesson01\\skincancer.txt",header=TRUE)
head(skinc)
##         State  Lat Mort Ocean  Long
## 1     Alabama 33.0  219     1  87.0
## 2     Arizona 34.5  160     0 112.0
## 3    Arkansas 35.0  170     0  92.5
## 4  California 37.5  182     1 119.5
## 5    Colorado 39.0  149     0 105.5
## 6 Connecticut 41.8  159     1  72.8
attach(skinc)
plot(Lat,Mort)

Možemo pretpostaviti linearnu vezu između ova dva obeležja \[y=b_0+b_1x\]

Uočavamo negativnu korelaciju (kako jedna promenljiva raste druga opada, i obrnuto) između ove dve promenljive, što znači negativno \(b_1\). Objašnjenje za to je da što je manja geografska širina, država je bliža ekvatoru, pa su njeni stanovnici izloženiji negativnom uticaju Sunca. Prava koja “najbolje” opisuje ovu vezu data je na sledećem grafiku.

Uvedimo prvo potrebne oznake za definisanje linearnog modela.

Pošto se ne bavimo determinističkim vezama, naš uzorak neće pripadati tačno pomenutoj pravoj, ali dovoljno je da bude blizu sa nekim manjim odstupanjem, pa da linearna veza bude dovoljno objašnjenje. Dakle, naš model je \[y_i=b_0+b_1x_i+e_i\] Parameter \(b_0\) zovemo i intercept (odsečak) a \(b_1\) slope (nagib). Ovaj model zovemo prosta linearna regresija (Simple Linear Regression). Kada se bavimo linearnom zavisnošću jedne zavisne sa više nezavisnih promenljivih radi se o višestrukoj linearnoj regresiji (Multiple Linear Regression). Iako koristimo nazive zavisna i nezavisna promenljiva, vrednosti \(x_i\) uopšte ne smatramo slučajnim u našem modelu (uslovno očekivanje u odnosu na prediktor). Slučajan je niz grešaka \(e_i\), pa odatle i niz \(y_i\).

Niz grešaka predviđanja ili reziduala se može izraziti i ovako \[e_i=y_i-\hat{y}_i.\] Sada ćemo objasniti po kom kriterijumu se određuje “najbolja prava” i kako se dobija njena jednačina. Kao specijalan slučaj razmatramo prostu linearnu regresiju gde je \(b_0=0\), ona se zove i regresija kroz koordinatni početak (Regression Through Origin). U tom slučaju, drugačija je ocena za \(b_1\), kao i neke druge osobine modela. Nju ćemo kasnije pomenuti. Dakle pretpostavljamo nadalje da je \(b_0\neq 0\) Najbolja prava je ona kojom pravimo najmanju grešku; zovemo je i regresiona prava. Kao ukupnu grešku predviđanja \(Q\) uzimamo sumu kvadrata reziduala. \[\begin{align} Q &=\sum_{i=1}^{n}(y_i-\hat{y}_i)^2 \\ &= \sum_{i=1}^{n}(y_i-b_0-b_1x_i)^2 \end{align}\]

Traženjem minimuma funkcije \(Q(b_0,b_1)\) dobijamo izraze za ocene parametara \(b_0\) i \(b_1\) metodom najmanjih kvadrata

\[\begin{align} \hat{b}_0=& \bar{y}-\hat{b}_1\bar{x} \\ \hat{b}_1=& \frac{\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^{n}(x_i-\bar{x})^2} \end{align}\]

Iz prethodnih izraza se može primetiti da tačka \((\bar{x},\bar{y})\) pripada regresionoj pravoj. Takođe, primetimo da znak koeficijenta pravca regresione prave zavisi samo od brojioca, koji ima oblik uzoračke kovarijacije. Kratko objašnjenje znaka kovarijacije i oblika linearne zavisnosti: Tačka \((\bar{x},\bar{y})\) deli ravan na četiri kvadranta. Ukoliko je većina tačaka u drugom i četvrtom kvadrantu (kao što je to bilo u našem primeru) onda je većina proizvoda oblika \((x_i-\bar{x})(y_i-\bar{y})\) negativna pa će i njihova suma biti negativna, kao i koeficijent \(b_1\). Ukoliko je većina tačaka u prvom i trećem kvadrantu, onda se dobija pozitivna vrednost za \(b_1\).

Lako se dobija da kod SLR važi i \(\sum e_i=0\) (ovo ne važi ako je \(b_0=0\)!)

Opisaćemo ukratko značenje parametara. Vrednost \(b_0\) nam govori koja je u proseku vrednost zavisne promenljive ako je vrednost prediktora 0 (ovo nema smisla ukoliko vrednost 0 nije u opsegu prediktora). Vrednost \(b_1\) nam govori za koliko se menja zavisna promenljiva ako se prediktor poveća za 1.

Regresionu pravu uvek možemo naći, nalaženjem minimuma funkcije \(Q\). Ipak, nekad linearna veza nije odgovarajuća. Teorijski, postoje neki uslovi koji trebaju biti ispunjeni (uslovi Gaus-Markova) a to su

Ako važe ovi uslovi, onda su ocene dobijene metodom najmanjih kvadrata dobre.

Treba nam neka informacija i o tome koliko je uopšte linearni model odgovarajući. Funkcija lm vraća objekat klase lm. Generička funkcija summary primenjena na taj objekat daje sve potrebne informacije o linearnom modelu.

model=lm(Mort~Lat,data=skinc)
summary(model)
## 
## Call:
## lm(formula = Mort ~ Lat, data = skinc)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.972 -13.185   0.972  12.006  43.938 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 389.1894    23.8123   16.34  < 2e-16 ***
## Lat          -5.9776     0.5984   -9.99 3.31e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.12 on 47 degrees of freedom
## Multiple R-squared:  0.6798, Adjusted R-squared:  0.673 
## F-statistic:  99.8 on 1 and 47 DF,  p-value: 3.309e-13

Na izlazu prvo dobijamo informacije o rezidualima. Vrednosti minimuma, maksimuma i kvartila tog niza. Sledeće dobijamo informacije o koeficijentima linearnog modela. Za dva parametra imamo dva reda. Prvi je slobodan član \(b_0\) a drugi koeficijent pravca, to jest parametar modela uz prediktor Lat. Za svaki od parametara modela dobijamo ocenu (dobijenu opisanom metodom najmanjih kvadrata) i standardnu grešku te ocene.

Pretpostavimo sada dodatno da greške imaju \(N(0,\sigma^2)\) raspodelu. Može se izvesti da \[\begin{align} E(\hat{b}_0)& =b_0, \quad D(\hat{b}_0)=\sigma^2\left[n^{-1}+\frac{\bar{x}_1^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}\right] \\ E(\hat{b}_1)& =b_1, \quad D(\hat{b}_1)=\sigma^2/\sum_{i=1}^{n}(x_i-\bar{x})^2) \end{align}\] Dakle, ocene parametara su nepristrasne. Odatle imamo da važi \[E(y_i)=b_0+b_1x_i,\] što znači da su regresionom pravom ocenjene vrednosti koje u proseku odgovaraju raznim vrednostima prediktora. Primetimo da disperzije ocena zavise od nepoznate disperzije grešaka. Jedna nepristrasna ocena za \(\sigma^2\) je \[s^2=\frac{1}{n-2}\sum_{i=1}^{n}e_i^2\] Često se označava i sa MSE. Faktor \(1/(n-2)\) javlja se zbog ocenjena dva parametra u SLR. Ova vrednost data je u R-u kao Residual standard error; broj stepeni slobode koji se navodi pored je upravo \(n-2\). Sada, zamenom \(\sigma^2\) sa \(s^2\) dobijamo ocene disperzija ocena. U drugoj koloni, Std Error nalaze se koreni ocena disperzija ocena, tj. \[\begin{align} s.e.(\hat{b}_0)=s\left[n^{-1}+\frac{\bar{x}^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}\right]^{1/2} s.e.(\hat{b}_1)=s/\left(\sum_{i=1}^{n}(x_i-\bar{x})^2)\right)^{1/2} \end{align}\]

Iz pretpostavke o normalnosti reziduala dobijamo i \(y_i \in N(b_0+b_1x_i,\sigma^2)\). Pošto su ocene \(\hat{b}_0\) i \(\hat{b}_1\) linearna kombinacija \(y_i\)-eva, i one imaju normalnu raspodelu. Odatle se može naslutiti Studentova raspodela pivota koji se koristi za testiranje i nalaženje intervala poverenja \[(\hat{b}_j-b_j)/s.e.(\hat{b}_j)\in t_{n-2}.\] Interval poverenja za vrednost parametra tada je \[(\hat{b}_1-s.e.(\hat{b}_1)t_{n-2,\alpha/2},\hat{b}_1+s.e.(\hat{b}_1)t_{n-2,\alpha/2}),\] gde je \(t_{n-2,\alpha/2}\) gornji \(\alpha/2\)-ti kvantil \(t_{n-2}\) raspodele. Poznata je sledeća veza između intervala poverenja i testiranja hipoteza. \((1-\alpha)100\%\)-ni interval poverenja za vrednost nekog parametra sadrži one vrednosti parametra za koje bi se pri testiranju odgovarajuće proste hipoteze (\(H_0: \hat{\theta}=\theta\), \(H_1: \hat{\theta}\neq\theta\) ) nulta hipoteza prihvatila na nivou značajnosti \(\alpha\). Zbog ovog, ako na primer interval poverenja za \(b_0\) sadrži nulu, to znači da se prihvata nulta hipoteza \(H_0: \hat{b}_0=b_0\). Treća i četvrta kolona se odnose upravo na pomenuto testiranje. Odgovarajuće test statistike imaju Studentovu t raspodelu, otud naziv t value. U poslednjoj koloni date su p-vrednosti testova. Kritične su velike apsolutne vrednosti test statistike, pa je p-vrednost testa upravo verovatnoća intervala u nazivu te kolone. Broj zvezdica pored vrednosti u poslednjoj koloni ne predstavlja ništa do oznake u kom se intervalu nalazi odgovarajuća vrednost. U našem primeru imamo oba puta tri zvezdice. To jednostavno naznačava da je broj sa leve strane u intervalu (0,0.001). Da smo, na primer, imali p-vrednost 0.06 pored nje bi stajala tačka.

Mi smo odabrali parametre modela tako da je suma kvadrata reziduala najmanja. Pitanje je i koliko je ta minimalna suma mala, jer ona govori o tome u kojoj meri je linearni model odgovarajući, koliko je disperzije početnog skupa \(1/n\sum_{i=1}^{n}(y_t-\bar{y})^2\) objašnjeno uzimanjem u obzir linearne veze sa prediktorom \(1/n\sum_{i=1}^{n}(y_t-\hat{y}_t)\). Jedna mera koja govori o tome koliko je smanjeno srednje kvadratno odstupanje je koeficijent determinacije \[R^2=\frac{\sum_{i=1}^{n}(\hat{y}_t-\bar{y})}{\sum_{i=1}^{n}(y_t-\bar{y})^2}=1-\frac{\sum_{i=1}^{n}(y_t-\hat{y}_t)}{\sum_{i=1}^{n}(y_t-\bar{y})^2}\] Ove sume su poznate pod nazivima SSR (Regression Sum of Squares), SSTO (Total Sum of Squares) i SSE (Error Sum of Squares), pa se \(R^2\) može lakše pamtiti ovako \[R^2=SSR/SSTO=1-SSE/SSTO.\] Ako je \(R^2=1\) to znači da sve tačke \(y_i\) pripadaju regresionoj pravoj. Ako je \(R^2=0\) onda je ocenjena linija horizontalna, pa prediktor ne objašnjava varijabilnost u zavisnoj promenljivoj ništa bolje od srednje vrednosti. Naziv \(R^2\) potiče od toga što je ovako definisana vrednost jednaka \(+/-\)korenu koeficijentu korelacije između prediktora i zavisne promenljive.

\[ \frac{\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum_{i=1}^{n}(x_i-\bar{x})^2\sum_{i=1}^{n}(y_i-\bar{y})^2}} \]

U to se možemo uveriti našim primerom. Pozitivna vrednost uzima se kada je znak \(b_1>0\), inače negativne

cor(Mort,Lat)
## [1] -0.8245178
cor(Mort,Lat)^2
## [1] 0.6798296

Može se dokazati da \(R^2\) uzima vrednosti u intervalu \([0,1]\). Ovim smo opisali Multiple R squared. Vrednost \(R^2\) se može često zloupotrebiti i pogrešno shvatiti. Navešćemo neke primere kada treba biti oprezan. \(R^2\) meri linearnu vezu između datih promenljivih. Ako se dobije da je \(R^2=0\) to znači da nema linearne zavisnosti, ali moŽda postoji neka nelinearna zavisnost. Velika vrednost \(R^2\) ne mora da znači da je linearna funkcija najbolja veza. Ovo se može videti sa grafika. Jedna neobična vrednost u uzorku moŽe puno da utiče na vrednost \(R^2\) Korelacije na implicira kauzalnost *Velika vrednost \(R^2\) ne znači prediktivnu moć

Osim \(t\)-testa kojim se kod SLR ocenjuje da li postoji lienarna zavisnost (\(b_1\neq 0\)), imamo i ANOVA \(F\)-test. Ideja je već rečena, razbiti ukupnu disperziju na deo koji je uhvaćen modelom i deo koji ostaje neobjašnjen modelom. Poznato je već \[SSTO=SSE+SSR\] Broj stepeni slobode za SSTO je \(n-1\) (kao što je imenilac kod ocene disperzije), a broj stepeni slobode za SSE je \(n-2\). Neka su MSE - Mean Square Error, MSR - Regression Mean Square dobijeni deljenjem SSE odnosno SSR sa svojim brojem stepeni slobode. Sledeća statistika ima \(F\) raspodelu \[\frac{MSR}{MSE}=\frac{SSR/1}{SSE/(n-2)}\] Vrednost ove statistike i pomenuti stepeni slobode prikazani su u poslednjem redu izlaza funkcije lm. Ona ima smisla u slučaju testiranja \(H_0: b_1=0\) protiv \(H_1: b_1\neq 0\) zato što je dobijeno \[\begin{align} E(MSR)&=\sigma^2+b_1^2\sum(x_i-\bar{x}) \\ E(MSE)&=\sigma^2 \end{align}\]

U prilog \(H_0\) idu vrednosti količnika bliske jedinici, a kritične su velike vrednosti. U ANOVA tabeli su detaljnije opisane vrednosti u vezi sa ovim testom

anova(model)
## Analysis of Variance Table
## 
## Response: Mort
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Lat        1  36464   36464  99.797 3.309e-13 ***
## Residuals 47  17173     365                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Važi da je F-statistika kvadrat t-statistike.

Regresija kroz koordinatni početak

U retkim slučajevima kada to ima smisla, možemo pretpostaviti da je \(\beta_0=0\). Model postaje \[y_i=\beta_1x_i+e_i.\] Ocena parametra je \[b_{RTO}=\frac{\sum_{i=1}^{n}x_iy_i}{\sum_{i=1}^{n}x_i^2}\] U ovom slučaju ne mora da važi da tačka koja je srednja vrednost uzorka pripada regresionoj pravoj. Uopšte, ne važi ni \(\sum_{i=1}^{n}e_i=0\), a \(R^2\) može da bude manje od nule. Zato se ono drugačije i definiše \[R^2=1-\frac{\sum_{i=1}^{n}e_i^2}{\sum_{i=1}^{n}y_t^2}\] Zbog ove drugačije definicije, pošto je imenilac u prvom slučaju dosta manji nego u drugom, ne možemo na osnovu \(R^2\) porediti modele koji imaju i one koji nemaju intercept.

Predviđanje

Najvažnija primena regresije je radi predviđanja. Mogu se postaviti dva pitanja: 1.Nalaženje srednje vrednosti zavisne promenljive \(\mu_Y\) za datu vrednost prediktora \(x_0\). 2.Predviđanje zavisne promenljive neke nove opservacije \(y_{new}\) sa datom vrednošću prediktora \(x_0\). Iako ova dva problema izgledaju isto, oni to nisu. Ova dva pitanja se u našem primeru prevode ovako 1. Koja je srednja vrednost stope smrtnosti za sve lokacije na 40 stepeni severne geografske širine? 2. Koja je predviđena vrednost stope smrtnosti jedne lokacije koja je na 40 stepeni severne geografske širine?

Što se tiče tačkaste ocene, ova dva problema imaju isti odgovor - jedino što imamo je regresiona prava, a greške ne možemo da uključimo u predviđanje zbog njihove nezavisnosti. U oba slučaja ocena je \[\hat{y}_0=\hat{b}_0+\hat{b}_1x_0\] Razlika je u intervalnom ocenjivanju. U prvom problemu tražimo interval poverenja za srednju vrednost zavisne promenljive kada je nezavisna jednaka \(x_h\). Očekivanje i disperzija predviđanja su \[\begin{align} E(\hat{y}_0)&=b_0+b_1x_0 \\ D(\hat{y}_0)&=\sigma^2\left[n^{-1}+\frac{(x_0-\bar{x})^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}\right] \end{align}\]

Ispostavlja se da i predviđanja, kada se skaliraju, imaju \(t_{n-2}\) raspodelu, pa se slično kao za parametre konstruiše interval poverenja. On je dat sa \[ \hat{y}_0 \pm t_{n-2,\alpha/2}\sqrt{s^2\left[n^{-1}+\frac{(x_0-\bar{x})^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}\right]}. \] Neka je \(x_0=40\). Ovu vrednost dobijamo pomoću funkcije predict

predict(model, newdata=data.frame(Lat=40),interval="confidence")
##        fit      lwr      upr
## 1 150.0839 144.5617 155.6061

Interval predviđanja koji se dobija za drugi problem dat je sa \[ \hat{y}_0 \pm t_{n-2,\alpha/2}\sqrt{s^2\left[1+n^{-1}+\frac{(x_0-\bar{x})^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}\right]}. \] Vrednost intervala predviđanja za novu opservaciju \(x_0=40\) u R-u dobijamo na sličan način

predict(model, newdata=data.frame(Lat=40),interval="predict")
##        fit     lwr      upr
## 1 150.0839 111.235 188.9329

Iz datih izraza za intervale poverenja vidimo da u drugom intervalu pod korenom imamo jedno dodatno \(s^2\). Pretpostavimo da su poznate stvarne vrednosti, srednja vrednost broja ljudi umrlih od raka kože na \(x_0=40\) je 150 po milionu ljudi a da je disperzija 400. Zanima nas predviđena vrednost smrtnosti u gradu Kolumbus, Ohajo (ovaj grad se nalazi tačno na 40 stepeni geografske širine i u sredini države Ohajo). Zbog poznate srednje vrednosti i disperzije, napravili bismo standardni interval poverenja za srednju vrednost normalne raspodele \(\mu_{Y}\pm2\sigma\). Problem je što nam srednja vrednost i disperzija nisu poznate. Ima smisla oceniti srednju vrednost \(\mu_{Y}\) odgovarajućom tačkom na regresionoj pravoj. Međutim različiti uzorci daju različite ocene za \(\mu_{Y}\). Moramo da uzmemo u obzir grešku koju pravimo ovom ocenom. Takođe, disperzija uglavnom nije poznata, nju ocenjujemo sa \(s^2\). Zbog ocenjivanja ovih nepoznatih, disperzija pri predviđanju nove zavisne promenljive sadrži dva dela. Prvi je varijabilnost kao posledica ocenjivanja \(\mu_Y\) sa \(\hat{y}_0\), a drugi je varijabilnost u vrednostima zavisne promenljive oko te srednje vrednosti. Prvi deo je onaj koji se nalazi pod korenom u izrazu za interval poverenja, drugi deo ocenjujemo se \(s^2\). Sabiranjem ove dve komponente dobijamo upravo ono što je pod korenom u izrazu za interval predviđanja.