# 8 Bivarijantno modeliranje: dvije kontinuirane varijable, 2. dio ## Linearna regresija ```{r, include=FALSE} # setwd("rossiter/scripts/") setwd("~/Dropbox/math/probability/R/rossiter/scripts") mhw = read.csv("../mhw.csv") ``` **Zadatak** Odredi pravac linearne regresija dobiven pomoću metode najmanjih kvadrata kada prinos slame ovisi o prinosu zrna. Odredimo model linearne regresije u slučaju kada je $y=ax+b$, koristeći R. ```{r, tidy=FALSE} x = mhw$grain y = mhw$straw model.straw.grain = lm(y ~ x) model.straw.grain ``` Ako nas zanimaju samo koeficijenti: ```{r} model.straw.grain$coefficients ``` ili, na drugačiji način ```{r} coefficients(model.straw.grain) ``` Znači, jednadžba pravca regresije je $y = 1.4305\cdot x + 0.8663$ Na temelju dobivenog modela `model.straw.grain` možemo, koisteći funkciju `predict()` predviđati vrijednos slame za proizvolju vrijednost zrna. Npr. koliki će biti prinos slame ako je prinos zrna 3.14 ili 4.13? ```{r} # Pažnja, u data.frame varijabla se mora zvati x, jer je # tako definiran model.straw.grain predict(model.straw.grain, data.frame(x = c(3.14, 4.13))) ``` **Zadatak** Iz podataka o prinosima zrnevlja i slame grafički prikažite točke i slici dodajte pravac regresije ```{r fig.width=7, fig.height=6, tidy=FALSE} plot(mhw$straw ~ mhw$grain) title("Prinos slame predvidjen prema prinosu zrna") abline(model.straw.grain) grid() text(4.5, 4.5, paste("koef. smjera: ", round(coefficients(model.straw.grain)[2], 2))) ``` Ista slika pomocu *ggplot2*. ```{r fig.width=7, fig.height=6, tidy=FALSE} ggplot(mhw, aes(grain, straw)) + geom_point(shape=1, size = 3) + # Use hollow circles geom_smooth(method=lm, # Add linear regression line se=FALSE) + # Don't add shaded confidence region labs(title="Prinos slame predvidjen prema prinosu zrna", x = "Prinos zrna", y = "Prinos slame") + annotate("text", x = 4.5, y = 4.5, label = paste("koef. smjera: ", round(coefficients(model.straw.grain)[2], 2))) + theme_bw() ``` Pomoću biblioteke *ggplot2* lako je dobiti i grafički prikaz intervala povjerenja za točke na pravcu regresije (podrazumijevanja vrijednost pouzdanosti je 95%). ```{r fig.width=7, fig.height=6, tidy=FALSE} ggplot(mhw, aes(grain, straw)) + geom_point(shape=1, size = 3) + # Use hollow circles geom_smooth(method=lm, # Add linear regression line ) + # Add shaded confidence region labs(title="Prinos slame predvidjen prema prinosu zrna", x = "Prinos zrna", y = "Prinos slame") + annotate("text", x = 4.5, y = 4.5, label = paste("koef. smjera: ", round(coefficients(model.straw.grain)[2], 2))) + theme_bw() ``` # Kako dobiti pdf file ``` list.files() ``` ``` system("pandoc --latex-engine=xelatex mhw-08-2.md -o mhw-08-2.pdf -H margins.sty") ``` Napomena. Style file margins.sty ima samo jedan redak i to `\usepackage[vmargin=1in,hmargin=1in]{geometry}` radi podrazumijevano prevelikih margina koje radi `pandoc` ### Analiza pravca regresije Podsjetnik: reziduali u regresijskoj analizi određeni su razlikom mjerenih i funkcijskih vrijednosti (modelskih vrijednosti) i trebali bi biti normalno distribuirani. **Zadatak** Prikažte histogram i sliku kvantil - kvantil za regresijske reziduale. Sumarno, numerički, prikažite osnovne podatke o rezidualima. ```{r, tidy=FALSE} re = residuals(model.straw.grain) hist(re, main = "Reziduali iz linearnog modela slama ~ zrna ") ``` Na prvi pogled čini se da reziduali nisu normalno distribuirani, jer je histogram previše strm. Provjerimo dojam na grafičkom prikazu kvantil - kvantil-a. ```{r, tidy=FALSE} qqnorm(re, main = "Q-Q plot reziduala iz linearnog modela slama ~ zrna") qqline(re) summary(re) ``` Iz slike se vidi da su najmanji reziduali premali a najveći preveliki da bi bili distribuirani po normalnoj razdiobi. Isto tako, iz sumarnog se prikaza vidi da je rang veći od 5 što je dosta velik broj. **Zadatak** Koristeći test Shapiro-Wilks testirajte normalnost reziduala. ```{r} shapiro.test(re) ``` "p-value" je manja od 0.05. Na temelju histograma, Q-Q plota i vrijednosti od "p-value" hipotezu da su reziduali distribuirani po normalnoj razdiobi treba odbaciti. **Zadatak** Nacrtaj Q-Q plot reziduala iz linearnog modela slama~zrno pomoću biblioteke *gglib2*. Na adresi http://stackoverflow.com/questions/4357031/qqnorm-and-qqline-in-ggplot2 piše > The ggplot package doesn't seem to contain code for calculating the > parameters of the qqline, so I don't know if it's possible to achieve > such a plot in a (comprehensible) one-liner. Znači, morati ćemo se koristiti posebnom funkcijom. Na istoj adresi je navedena funkcija: ```{r, tidy=FALSE} qqplot.data <- function (vec) # argument: vector of numbers { # following four lines from base R's qqline() y <- quantile(vec[!is.na(vec)], c(0.25, 0.75)) x <- qnorm(c(0.25, 0.75)) slope <- diff(y)/diff(x) int <- y[1L] - slope * x[1L] d <- data.frame(resids = vec) ggplot(d, aes(sample = resids)) + stat_qq() + geom_abline(slope = slope, intercept = int) } ``` Slika bez `qqline()` ne treba posebnu funkciju: ```{r, tidy=FALSE} redf=data.frame(re) ggplot(redf, aes(sample=re)) + stat_qq() ``` ali sliku sa `qqline()` trebamo nacrtati sa gore definiranom funkcijom. ```{r} qqplot.data(re) ``` **Zadatak** Usporedite modelske vrijednosti (vrijednosti dobivene pomoću pravca regresije) i mjerene vrijednosti slame. Bojom istaknite one točke koje su (vertikalno) od pravca regresije (modelirane vrijednosti) udaljene za više od jedan. ```{r, tidy=FALSE} fmsg = fitted(model.straw.grain) # lim = c(min(fmsg, mhw$straw), max(fmsg, mhw$straw)) df.fmsg = data.frame(fmsg, mstraw=mhw$straw) ggplot(df.fmsg, aes(fmsg, mstraw)) + geom_point(color = ifelse( (abs(fmsg - mhw$straw)) < 1 , "gray" , ifelse( (fmsg < mhw$straw), "blue", "red")), size = ifelse(abs(fmsg - mhw$straw) < 1, 3, 4)) + geom_abline(intercept = 0, slope = 1) + labs(title = "Mjerene vrijednosti koje su (vertikalno) za vise od 1 udaljene od pravca regresije", x = "Modelske vrijednosti (vrijednosti sa pravca regresije)", y = "Slama") ``` **Zadatak** Neki prinosi slame nisu dobro aproksimirani pomoću linearnog modela. Identificirajte one prinose slame kojima je apsolutna vrijednost reziduala veća od 1.8. Prikažite te retke iz `mhw` sortirane po veličini omjera zrna/slama (grain/straw). ```{r} re = residuals(model.straw.grain) which.high.res = which(abs(re) > 1.8) sort(re[which.high.res]) ``` Izdvojimo iz tablice `mhw` retke određene sa `which.high.res` pa toj tablici dodajmo stupac `gsr` omjera `grain/straw`. ```{r} high.res = mhw[which.high.res, ] high.res$gsr = high.res$grain/high.res$straw high.res ``` Prikažimo retke (rekorde) iz data framea `high.res` sortirane po omjeru `grain/straw`. ```{r} high.res[order(high.res$gsr), ] ``` ### Daljnja dijagonistika linearne regresije Ako je model linearne regresije odgovara podacima tada ne bi smjelo biti nikakve relacije (veze) između reziduala i modelskih vrijednosti (vrijednosti na pravcu linearne regresije). Funkcija `plot` omogućava grafički prikaz, u kojem se prikazuju mjerene točke i krivulja. Ako nema nikakve relacije između reziduala i vrijednosti na pravcu regresije tada bi nacrtana krivulja trebala biti horizontali pravac i trebao bi prolaziti kroz ishodište (tj. trebao bi ležati na x-osi). Točke na ovakvoj krivulji dobijemo pomoću funkcije `lowess()`. ```{r fig.width=4, fig.height=3, tidy=FALSE} plot(re ~ fmsg) title("Reziduali vs. modelirane vrijednosti") tocke = lowess(fmsg, re) lines(tocke$x, tocke$y, type="l", col="red") ``` Ovakvu sliku možemo dobiti jednostavnije opcijom `which = 1` funkcije `plot()`. ```{r} plot(model.straw.grain, which = 1) ``` Kako dobivena krivulja nije pravac, to ne možemo, u potpunosti, biti zadovoljni gornjom linearnom regresijom. Nacrtajmo "istu" sliku pomoću biblioteke `ggplot2`. ```{r fig.width=7, fig.height=6, tidy=FALSE} privdf = data.frame(x=fmsg, y=re) lowessdf = data.frame(lowess(fmsg, re)) ggplot(data = privdf, aes(x, y)) + geom_point(size = 2) + geom_line(data=lowessdf, aes(x, y), color = "red") + labs( title = "Krivulja dobivena lokalnom tezinskom lin. regresijom", x = "Modelirane vrijednosti", y = "Reziduali") + theme_bw() ``` ### Malo o funkciji `lowess()`. Upotreba ``` lowess(x, y, f=2/3, iter=3, delta=.01*diff(range(x))) ``` Kako rezultat funkcije ovisi o parametru `f`? - f: the smoother span. This gives the proportion of points in the plot which influence the smooth at each value. Larger values give more smoothness. ```{r} data(cars) plot(cars, main ="lowess(cars)") lines(lowess(cars), col=2) lines(lowess(cars, f=.2), col=3) lines(lowess(cars, f=.3), col=4) legend(5, 120, c(paste("f = ", c("2/3", ".2", ".3"))), lty = 1, col = 2:4) ``` ### Predviđanja vrijednosti Ako je određen pravac regresije i želimo predvidjeti vrijednost na temelju prediktora moguće su dvije greške: 1. Nepouzdanost u određivanju pravca regresije 2. Nepouzdanost u određivanju vrijednosti iz prediktora. Stoga je dobro odrediti intervale povjerenja za vrijednosti pravca regresije (confidence interval) i intervale povjerenja za predviđenu vrijednost (prediction interval). **Zadatak** Izračunajte vrijednost i 95% confidence interval i prediction interval prinosa slame za srednju vrijednost $\mu$ prinosa zrna, te za $\mu \pm k\cdot\sigma$ gdje je $k = 1, 2$. ```{r, tidy=FALSE} xpotez = mean(mhw$grain); sdev = sd(mhw$grain) # Pažnja, u p.frame varijabla se mora zvati x, jer je # tako definiran model.straw.grain p.frame = data.frame(x = seq(xpotez - 2*sdev, xpotez + 2*sdev, by = sdev)) predict(model.straw.grain, newdata = p.frame, interval="confidence", level=0.95) ``` ```{r, tidy=FALSE} predict(model.straw.grain, newdata = p.frame, interval="prediction", level=0.95) ``` Prediction interval je, prirodno, širi od confidence interval jer je zadnji ocjena procječne vrijednosti a prvi ocjena bilo kakave vrijednosti. **Zadatak** Nacrtajte točke koje prikazuju slamu vs. žita, pravac regresije i oba intervala povjerenja za prinos žita od 2 do 6 sa korakom 0.1. ```{r, tidy=FALSE} p.frame = data.frame(x = seq(from = 2, to = 6, by = 0.1)) pred.c = predict(model.straw.grain, newdata = p.frame, interval = "confidence", level = 0.95) pred.p = predict(model.straw.grain, newdata = p.frame, interval = "prediction", level = 0.95) plot (straw ~ grain, data = mhw, pch = 20) title(main = "Prinos slame prema predvidjanju prinosa zita", sub = "Prediction (blue) i confidence (red) intervali") abline(model.straw.grain) grid() lines(p.frame$x, pred.c[, "lwr"], col = "red", lwd = 1.5) lines(p.frame$x, pred.c[, "upr"], col = "red", lwd = 1.5) lines(p.frame$x, pred.p[, "lwr"], col = "blue", lwd = 1.5) lines(p.frame$x, pred.p[, "upr"], col = "blue", lwd = 1.5) points(mean(x), mean(y), pch=23, cex=2, bg="red") abline(h = mean(y), lty = 2) abline(v = mean(x), lty = 2) ``` **Zadatak** Usporedite regresiju slama vs. žito sa regresijom žito vs. slama. ```{r} model.grain.straw = lm (grain ~ straw, data = mhw) summary (model.grain.straw) ``` ```{r} summary (model.straw.grain) ``` U rezultatima se navodi "multiple R" (točnije njegov kvadrat). Značenje tog koeficijenta (bez kvadrata) je koeficijent korelacije između predviđene vrijednosti(na pravcu regresije) i izmjerene vrijednosti. Odredite koeficijent korelacije R za obje regresije ```{r} sqrt(summary(model.grain.straw)$r.squared) sqrt(summary(model.straw.grain)$r.squared) ``` Dobiveni rezultati su isti. Koja je posljedica te činjenice? Neka Y ovisi o X (tj. model je Y ~ X); slučajna vrijednost Y se sastoji od vrijednosti pravca regresije (a vrijednost na pravcu regresije je određena sa X) i slučajnog odstupanja od pravca regresije. Statistička ovisnost odstupanje od pravca regresije je opisano koeficijentom korelacije R. Kako su u oba modela koeficijenti isti to slijedi da poznavanje slame u prvom modelu (žito ~ slama) jednako određuje žito kao što poznavanje žita u drugom modelu (slama ~ žito) određuje slamu. **Zadatak** Može se činiti da je pravac regresije u modelu (žito ~ slama) inverzna funkcija pravca regresije u modelu (slama ~ žito). Provjerite da li je to istina? Ako je zadan pravac $y=a\cdot x + b$, njemu inverzan se dobije ako se po $y$ riješi izraz $x = a\cdot y + b$, tj. dobije se $$ y = \frac{1}{a} x - \frac{b}{a} $$ Pravac regresije u modelu (slama ~ žito) je $y = 1.431\cdot x + 0.866$. Inverzni pravac je $$ \begin{aligned} y &= \frac{1}{1.431}x - \frac{0.866}{1.431} \\ y &= 0.699\cdot x - 0.605 \end{aligned} $$ Budući je jednadžba pravca regresije u modelu (žito ~ slama) jednaka $y = 0.3723 \cdot x + 1.523$ to inverzni pravac prvog modela nema niti jedan koeficijent zajednički sa pravcem drugog modela. Zaljučujemo da regresija između dvije varijable nije simetrična već rezultat ovisi o tome koja je varijabla input a koja output (koja je nezavisna a koja zavisna).