Podsjetimo se.
Područje uzgoja smo podijeliti na dva dijela, sjeverni i južni i istražili smo, da li, statistički, postoji razlika u prinosima.
Podjelu smo učinili tako da smo svakom polju pridružili logičku vrijednost TRUE ako je polje u sjevernom dijelu, a FALSE ako nije. Polja sa indeksom retka manjim od 11 su u sjevernom a ostala u južno dijelu.
Zadatak Promijenite data frame mhw tako da dodate stupac
(varijablu) in.north koji pokazuje da li je polje u sjevernom dijelu
ili ne.
mhw = cbind(mhw, in.north = (mhw$r < 11))
str(mhw)
## 'data.frame': 500 obs. of 5 variables:
## $ r : int 1 2 3 4 5 6 7 8 9 10 ...
## $ c : int 1 1 1 1 1 1 1 1 1 1 ...
## $ grain : num 3.63 4.07 4.51 3.9 3.63 3.16 3.18 3.42 3.97 3.4 ...
## $ straw : num 6.37 6.24 7.05 6.91 5.93 5.59 5.32 5.52 6.03 5.66 ...
## $ in.north: logi TRUE TRUE TRUE TRUE TRUE TRUE ...
Zadatak prikažite rezultate, iz predhodnih područja, regresije modela
(prinos slame) ~ (polovice polja) i
(prinos slame) ~ (prinos zrna).
model.straw.ns = lm(straw ~ in.north, data = mhw)
summary(model.straw.ns)
##
## Call:
## lm(formula = straw ~ in.north, data = mhw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.181 -0.608 -0.108 0.572 2.359
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.7484 0.0549 122.90 < 2e-16 ***
## in.northTRUE -0.4671 0.0777 -6.02 3.5e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.868 on 498 degrees of freedom
## Multiple R-squared: 0.0677, Adjusted R-squared: 0.0659
## F-statistic: 36.2 on 1 and 498 DF, p-value: 3.48e-09
Grafički prikažite dobivenu linearnu regresiju
plot(mhw$in.north, mhw$straw)
abline(model.straw.ns)
Koeficijenti pravca \( y = \alpha + \beta\cdot x \) koji se mogu očitati iz numeričkih rezultata 6.7484, -0.4671 kao i slika pokazuju da pravac regresije pada od sjevernog djela prema južnom što znači da je očekivana vrijednost prinosa na sjeveru veća nego na jugu.
Veličina \( R^2 \) je 0.0677. To je vrlo mali broj što znači da je dio varijance koja se može objasniti prelaskom sa sjevera na jug (utjecaj modela na promjenu varijance) vrlo mali.
model.straw.grain = lm(straw ~ grain, data = mhw)
summary(model.straw.grain)
##
## Call:
## lm(formula = straw ~ grain, data = mhw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0223 -0.3529 0.0104 0.3734 3.0342
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8663 0.2387 3.63 0.00031 ***
## grain 1.4305 0.0601 23.82 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.615 on 498 degrees of freedom
## Multiple R-squared: 0.533, Adjusted R-squared: 0.532
## F-statistic: 567 on 1 and 498 DF, p-value: <2e-16
Koeficijenti pravca regresije \( y = \alpha + \beta\cdot x \) u modelu (straw ~ grain) koji se mogu očitati iz numeričkih rezultata 0.8663, 1.4305 pokazuje da pravac regresije raste, pa znači da se porastom prinosa zrna očekuje porast prinosa slame i to za svaku jedinicu porasta prinosa zrna porast slame je 1.4305.
U ovom modelu veličina \( R^2 \) je 0.5326 što znači da se, približno, pola varijance prinosa slame može objasniti modelom (model utiče na promjenu varijance prinosa slame sa, približno, 50%).
Zadatak Modelirajte prinos slame kao posljedica kombiniranog utjecaja nezavisnih prediktora: polovice polja (sjevernog i južnog dijela) i prinosa zrna.
model.straw.ns.grain = lm (straw ~ in.north + grain, data = mhw)
summary (model.straw.ns.grain)
##
## Call:
## lm(formula = straw ~ in.north + grain, data = mhw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2548 -0.3189 -0.0276 0.3042 2.7871
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0461 0.2178 4.8 2.1e-06 ***
## in.northTRUE -0.5132 0.0500 -10.3 < 2e-16 ***
## grain 1.4499 0.0546 26.5 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.559 on 497 degrees of freedom
## Multiple R-squared: 0.614, Adjusted R-squared: 0.613
## F-statistic: 396 on 2 and 497 DF, p-value: <2e-16
Iz numeričkih rezultata (koeficijenti su 1.0461, -0.5132, 1.4499) slijedi da su oba prediktora potrebna da se opiše očekivanje prinosa slame (ni jedan koeficijent nije blizu nule). Prelazak sa južnog dijela na sjeverni znači pad prinosa slame (južni dio je opisan sa FALSE (0), sjeverni sa TRUE (1)), a porast prinosa zrna znači rast prinosa slame. Detaljnije, očekivni prinosi slame se mijenjaju ovako. Prelazak sa južnog dijela polja na sjeverni znači pad očekivanja prinosa slame za, približno, 0,5132 (u aplsolutnom iznosu), a jedinični porast prinosa zrna znači porast očekivanjog prinosa slame za 1,4499. Primjetite da koeficijent 1,4499 prinosa zrna nije isti kao koeficjent u modelu (straw ~ grain), kada se sjever i jug nisu uzimali u obzir; tada je taj koeficijetn iznosio 1,4305.
Daljni komentar rezultata je da je \( R^2 \) jednak 0.6142 što znači da se preko 60% varijance prinosa slame može objasinit modelom.
U zadnjem primjeru radi se o aditivnom modelu sa jednim neprekidnim i jednim diskretnim prediktorom. Zbog diskretnosti drugog prediktora postoji po jedna pravac regresije za svaku vrijednost diskretnog prediktora. U našem slučaju, budući diskretni prediktor ima dvije moguće vrijednosti, postoje dva pravca regresije, jedan za sjeveri a drugi za južni dio polja. Grafički ćemo prikazati ta dva pravca, zajedno sa pravcem iz modela (straw ~ grain) (tj. kada se sjever i jug ne uzimaju u obzir).
attach(mhw)
plot(straw ~ grain, col = ifelse(in.north, "blue", "slategray"),
pch = 20, xlab = "Prinos zrna", ylab = "Prinos slame")
title (main = "Prinos slame ovisno o prinosu zrna")
title (sub = "Sjeverni dio: plavo, Juzni dio: sivo, Cijelo polje: crveno")
abline (coefficients(model.straw.ns.grain)["(Intercept)"],
coefficients(model.straw.ns.grain)["grain"], col = "slategray")
abline (coefficients(model.straw.ns.grain)["(Intercept)"] +
coefficients(model.straw.ns.grain)["in.northTRUE"],
coefficients(model.straw.ns.grain)["grain"], col = "blue")
abline (model.straw.grain, lty = 2, col = "red")
Obratite pažnju da je odsječak pravca na y osi, za južni dio polja, jednak
coefficients(model.straw.ns.grain)["(Intercept)"]
a za sjeverni dio
coefficients(model.straw.ns.grain)["(Intercept)"] +
coefficients(model.straw.ns.grain)["in.northTRUE"].
Dalje, koeficijenti smjera su isti pa znači da su,
kod linarnog aditivnom modela pravci regresije paralelni.
Zadatak Usporedite \( R^2 \) tri modela
model.straw.ns
model.straw.grain
model.straw.ns.grain
i komentirajte rezultat.
Numerički podaci su
summary(model.straw.ns)$r.squared
## [1] 0.06774
summary(model.straw.grain)$r.squared
## [1] 0.5326
summary(model.straw.ns.grain)$r.squared
## [1] 0.6142
Od tri modela model model.straw.ns.grain najbolje objašnjava
varijabilnost prinosa slame — blizu dvije trećine (61%) varijabilnosti
je posljedica modela. Model model.straw.grain objašnjava oko
pola varijabilnosti (53%), dok model model.straw.ns vrlo malo utiče
na varijabilnost prinosa slame (6%).
Drugi način da usporedimo dva hijerarijska modela
(modeli su hijerarskijski ako složeniji ima sve prediktore jednostavnijeg
modela) je pomoću analize varijance. To je statistički test, tako da
možemo odrediti da li je složeniji model dokazano bolji. Za dobivanje
rezultata koristiti ćamo funkciju anova().
Funkcija anova(), u slučaju primjene na hijerarskijske modele
postavalja nul hipotezu da nema poboljšanja pri prelasku na
složeniji model. Uspoređuju se zbrojevi kvadrata reziduala
(Residual Sum of Squares — RSS) dva modela. To je iznos varijabiliteta
koji nije objašnjen modelom, tako da je manji RSS znači bolji rezultat.
Tada se računa F–omjer između obje varijance i vjerojatnost da se
dobije tako veliki F–omjer ako je nul hipoteza da nema
modelskog poboljšanja istinita. Vjerojatnost greške prvog tipa
(da ćemo odbaciti istinitu hipotezu) se prikazuje kao Pr(>F);
manja vrijednost je bolja.
anova(model.straw.ns.grain, model.straw.ns)
## Analysis of Variance Table
##
## Model 1: straw ~ in.north + grain
## Model 2: straw ~ in.north
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 497 155
## 2 498 375 -1 -220 704 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model.straw.ns.grain, model.straw.grain)
## Analysis of Variance Table
##
## Model 1: straw ~ in.north + grain
## Model 2: straw ~ grain
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 497 155
## 2 498 188 -1 -32.9 105 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Na temelja dobivenih rezultata možemo zaključiti
U oba slučaja je RSS složenijeg modela manji nego RSS jednostavnijeg pa složeniji model bolje objašnjava varijetet vrijednosti prinosa slame.
Vjerojatnost da ćemo učiniti grešku ako odbacimo nul hipotezu da poboljšanja nema pri prelasku na složeniji model je, u oba slučaja, praktički nula, pa nul hipotezu treba odbaciti.
Vidjeli smo da je linearni aditivni model superioran modelima sa jednim prediktorom. U tom modelu smatra se da su uticaji prediktora na prinos slame nezavisni, no tako ne mora biti. U našem slučaju to znači da je relacija između dvije polovice polja i prinosa zrna različita na te dvije polovice i da različito djeluje na prinos slame. Takav slučaj se zove interakcija. To daje dvije različite linearne regresije (neparalelni pravci) za razliku od dvije paralelne regresije u slučaju aditivnog linearnog modela.
Zadatak Modelirajte prinos slame kao kombinirani efekt dva prediktora čiji doprinos, međusobno, nije nezavisan.
Izraz u formuli za regresiju, koji daje doprinos interakcije dva prediktora, se
zadaje sa
in.north : straw
stoga formula glasi
straw ~ in.north + grain + in.north : grain
no ta forumla se zapisuje kraće
straw ~ in.north * grain
model.straw.ns.grain.i = lm (straw ~ in.north * grain, data = mhw)
summary (model.straw.ns.grain.i)
##
## Call:
## lm(formula = straw ~ in.north * grain, data = mhw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2242 -0.3169 -0.0398 0.3136 2.7574
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.2930 0.2979 4.34 1.7e-05 ***
## in.northTRUE -1.0375 0.4350 -2.39 0.017 *
## grain 1.3872 0.0752 18.44 < 2e-16 ***
## in.northTRUE:grain 0.1328 0.1094 1.21 0.225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.559 on 496 degrees of freedom
## Multiple R-squared: 0.615, Adjusted R-squared: 0.613
## F-statistic: 265 on 3 and 496 DF, p-value: <2e-16
Značenje koeficijenata iz predhodnog rezultata je slijedeće:
(Intercept) : procijenjeni prinos slame kada je prinos zrna jednak nuli u
južnoj polovici polja;
in.northTRUE : razlika u prosječnom prinosu slame na sjeveru u odnosu
na jug;
grain : porast prinosa slame za jedinični porast prinosa zrna u
južnom dijelu;
in.northTRUE:grain : razlika u koeficijentu smjera (razlika u porast
prinosa slame za jedinčni porast prinosa zrna) na sjeveru u odnosu
na jug.
Iz podataka se vidi da je koeficijent smjera na jugu 1.3872 a na sjeveru 1.3872 + 0.1328 što iznosi 1.5199. Razlika nije velika pa zaključujemo da ti koeficijenti nisu bitno različiti. To znači da interakcija zrna sa sjeverom i jugom ne utiče bitno na prinos slame.
Redak in.northTRUE:grain sadrži i test nul hipoteze da je utjecaj tog
retka na razliku jednak nuli i p-value je visok, 0.2255, pa ne treba
odbaciti hipotezu da je utjecaj zavisnosti zanemariv.
Zadatak Nacrtajte te dvije regresije i regresije aditivnog lineranog modela (tj. bez međusobnog utjecaja prinosa zrna i sjevera i juga).
Da bi mogli nasrtati te dvije regresije koristiti ćemo se opcijom
subset = funkcije lm() da izdvojimo observacije posebno za sjever
a posebno za jug.
plot (straw ~ grain, col = ifelse(in.north, "blue", "slategray"),
pch = 20, xlab = "Prinos zrna", ylab = "Prinos slame")
title (main = "Slama vs. zrna, prema polovicama polja")
title (sub = "Interakcija: neprekidne linije, aditivni model: prekinute linije")
abline (lm (straw ~ grain, subset = in.north), col = "blue")
abline (lm (straw ~ grain, subset = !in.north), col = "slategray")
# Isto je ako se zadnja dva retka napišu ovako
# abline (lm (straw ~ grain*in.north, subset = !in.north), col = "slategray")
# abline (lm (straw ~ grain*in.north, subset = in.north), col = "blue")
abline (coefficients(model.straw.ns.grain)["(Intercept)"],
coefficients(model.straw.ns.grain)["grain"], col = "slategray", lty = 2)
abline (coefficients(model.straw.ns.grain)["(Intercept)"] +
coefficients(model.straw.ns.grain)["in.northTRUE"],
coefficients(model.straw.ns.grain)["grain"], col = "blue", lty = 2)
Zadata Usporedite interaktivni i linearni aditivni model pomoću \( R^2 \) i pomoću analize varijance.
summary (model.straw.ns.grain)$r.squared
## [1] 0.6142
summary (model.straw.ns.grain.i)$r.squared
## [1] 0.6154
anova(model.straw.ns.grain.i, model.straw.ns.grain)
## Analysis of Variance Table
##
## Model 1: straw ~ in.north * grain
## Model 2: straw ~ in.north + grain
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 496 155
## 2 497 155 -1 -0.46 1.47 0.23
Iz numeričkih rezultata je vidljivo da se veličine \( R^2 \) jako malo razlikuju što znači da interaktivni model ne doprinosi većem postotku utjecaja na varijabilitet prinosa slame (Za podjsećanje: \( R^2 \) je postotak varijabilnosti varijable koji se može objasniti modelom ili drugačije \( R^2 \) = (Varijacija zbog modela) / (Ukupna varijacija varijable). Općenito: \( 0 \leq R^2 \leq 1 \) i što je \( R^2 \) veći model bolje odgovara podacima).
ANOVA pokazuje da hipotezu o jednakosti modela ne treba odbaciti jer je p-value = 0.23.
Kao i kod jednodimenzionalne regresije i višedimenzionalu treba istražiti da bi se utvrdilo da li odgovara modelskim predpostavkama.
Zadatak Prikažite dijagnostičke crteže za aditivni model:
reziduali vs. vrijednosti pravca regresije (fit),
normalni QQ–plot reziduala,
reziduale vs. leverage
Cook-ova udaljenost vs. leverage.
To su crteži koji nastaju pomoću funkcije plot.lm() dobiveni
kada opcija which= poprimi vrijednosti 1, 2, 5 i 6. Također ćemo
opcijom id.n specificirati da se istaknu ekstremalne točke.
Za razumijevanje slika korisno je znati slijedeće pojmove:
Standardizirani rezidual je rezidual podijeljen sa svojom standardnom devijacijom tako da standardizirani rezidual ima devijaciju jednaku jedan.
Snaga točke (Leverage value) je mjera utjecaja i-te točke na i-tu fitovanu vrijednost yhat_i, a Kukova udaljenost (Cook's distance) je normalizirana mjera utjecaja i-te točke na sve fitovane vrijednosti.
Napomena: Leverage value se dobije funkcijom hatvalues(model), a Cook's distance
funkcijom cooks.distance().
par (mfrow = c(2,2))
plot.lm (model.straw.ns.grain, which = c(1, 2, 5, 6), id.n = 10)
par (mfrow = c(1, 1))
Identificirajmo podpolja sa velikim rezidualima koristeći funkciju
rstandard() (“standardizirani rezidual”) i prikažimo samo te
retke podataka data frame-a.
(selected = which(abs(rstandard(model.straw.ns.grain)) > 3))
## 15 35 184 285 292 295 311 337 362
## 15 35 184 285 292 295 311 337 362
rstandard(model.straw.ns.grain)[selected]
## 15 35 184 285 292 295 311 337 362
## 5.001 -4.046 -3.193 3.178 -3.065 3.811 3.874 3.907 3.407
Dalje, formirajmo data frame sa svim informacijama za ta podpolja,
skupa sa rezidualima; koristiti ćemo funkciju cbind da dodamo stupac.
mhw.hires = cbind(mhw[selected, ], sres = rstandard(model.straw.ns.grain)[selected])
rm(selected)
str(mhw.hires)
## 'data.frame': 9 obs. of 6 variables:
## $ r : int 15 15 4 5 12 15 11 17 2
## $ c : int 1 2 10 15 15 15 16 17 19
## $ grain : num 3.46 4.42 4.59 3.7 4.86 3.73 3.74 3.05 4.26
## $ straw : num 8.85 5.2 5.41 7.67 6.39 8.58 8.63 7.64 8.61
## $ in.north: logi FALSE FALSE TRUE TRUE FALSE FALSE ...
## $ sres : num 5 -4.05 -3.19 3.18 -3.06 ...
Konačno poredajmo retke zadnjeg rezultata po rezidualima; koristiti
ćemo funkciju order():
mhw.hires[order(mhw.hires$sres), ]
## r c grain straw in.north sres
## 35 15 2 4.42 5.20 FALSE -4.046
## 184 4 10 4.59 5.41 TRUE -3.193
## 292 12 15 4.86 6.39 FALSE -3.065
## 285 5 15 3.70 7.67 TRUE 3.178
## 362 2 19 4.26 8.61 TRUE 3.407
## 295 15 15 3.73 8.58 FALSE 3.811
## 311 11 16 3.74 8.63 FALSE 3.874
## 337 17 17 3.05 7.64 FALSE 3.907
## 15 15 1 3.46 8.85 FALSE 5.001
Vizualizirati ćemo lokaciju tih podpolja; ona sa visokim pozitivnim rezidualima plavo, sa visokim negativnim rezidualima crveno:
plot (c, r, ylim = c(20,1), pch = 20, cex = 3,
col = ifelse(rstandard(model.straw.ns.grain) > 3, "darkblue",
ifelse (rstandard(model.straw.ns.grain) < -3 , "red" ,
ifelse (in.north, "lightblue", "gray"))),
xlab = "stupci", ylab ="reci")
abline (h = 10.5)
title (main = "Prikaz velikih reziduala")
title (sub = "Pozitivni plavi, negativni crveni")
U aditivnom modelu postoji jedan pravac regresije koji se translatira, posebno za svaku polovicu polja. U interakcijskom modelu postoje pravac regresije i njegova devijacija koja potječe od druge polovice polja.
Drugi način je da odustanemo od jednog pravca regresije i da formiramo dva nezavisna pravca svaki za svoju polovicu polja. Takav model se zove ugnježđeni model (nested model). Tu nam nije cilj naći vezu između slame i zrna na cijelom polju, već u svakoj polovici polja posebno.
Ugnježđeni model se specificira zanakom / u formuli (to nije matematičko
dijeljenje). Za naš data.frame fromula će biti
straw ~ in.north/grain
i čita se “straw ovisno o grain za svaku vrijednost varijable in.north
posebno”
model.straw.ns.grain.nest = lm (straw ~ in.north/grain, data = mhw)
summary (model.straw.ns.grain.nest)
##
## Call:
## lm(formula = straw ~ in.north/grain, data = mhw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2242 -0.3169 -0.0398 0.3136 2.7574
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.2930 0.2979 4.34 1.7e-05 ***
## in.northTRUE -1.0375 0.4350 -2.39 0.017 *
## in.northFALSE:grain 1.3872 0.0752 18.44 < 2e-16 ***
## in.northTRUE:grain 1.5199 0.0794 19.14 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.559 on 496 degrees of freedom
## Multiple R-squared: 0.615, Adjusted R-squared: 0.613
## F-statistic: 265 on 3 and 496 DF, p-value: <2e-16