Da bi ilustrirali robusne metode treba nam skup podataka koji sadrži
nekoliko podatka koji od razumnih odskaču po veličini (sadrži outlierse)
Zbog toga ćemo neke od podataka iz data framea mhw promijeniti.
Kontaminirati ćemo sjeverozapdni dio polja.
Polje je, u odnosu na strane svijeta, položano ovako:
N
1. 2. . . . 25.
-------------------
1. | | | ... | |
-------------------
2. | | | ... | |
-------------------
W . I
.
.
-------------------
20. | | | ... | |
-------------------
S
Prvo ćemo kopirati orginalne podatke
mhw.c = mhw
zatim ćemo odredite indekse sjeverozapadnog dijela
ix = (mhw.c$r < 5) & (mhw.c$c < 5)
rownames(mhw.c[ix, ])
## [1] "1" "2" "3" "4" "21" "22" "23" "24" "41" "42" "43" "44" "61" "62"
## [15] "63" "64"
subset(mhw.c, ix) # ili ovako mhw.c[ix, ]
## r c grain straw
## 1 1 1 3.63 6.37
## 2 2 1 4.07 6.24
## 3 3 1 4.51 7.05
## 4 4 1 3.90 6.91
## 21 1 2 4.15 6.85
## 22 2 2 4.21 7.29
## 23 3 2 4.29 7.71
## 24 4 2 4.64 8.23
## 41 1 3 4.06 7.19
## 42 2 3 4.15 7.41
## 43 3 3 4.40 7.35
## 44 4 3 4.05 7.89
## 61 1 4 5.13 7.99
## 62 2 4 4.64 7.80
## 63 3 4 4.69 7.50
## 64 4 4 4.04 6.66
i na kraju ćamo promijeniti podatke o prinosu zrna na tom sjeverozapadnom dijelu
mhw.c[ix, "grain"] = mhw.c[ix, "grain"]/4
Zadatak Sumarno prikažite numerički i grafički efekt kontaminacije
summary(mhw$grain)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.73 3.64 3.94 3.95 4.27 5.16
summary(mhw.c$grain)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.908 3.590 3.920 3.850 4.260 5.160
sd(mhw$grain)
## [1] 0.4583
sd(mhw.c$grain)
## [1] 0.6764
par (mfrow = c(1, 2))
hist (mhw$grain, xlab = " Prinos zrna", main = "Originalni podaci",
breaks = seq (0, 6, by = .25))
rug (mhw$grain)
hist (mhw.c$grain, xlab = " Prinos zrna", main = "Kontaminirani podaci",
breaks = seq (0, 6, by = .25))
rug (mhw.c$grain)
par (mfrow = c(1, 1))
Zadatak Izračunajte srednju vrijednost i standardnu devijaciju originalnog i kontaminiranog prinosa zrna. Također izračunajte medijan i interkvartilni razmak. Za sve ove podatke izračunajte postotak promjene.
# Sredine
mean (mhw$grain)
## [1] 3.949
mean (mhw.c$grain)
## [1] 3.846
# Standardne devijacije
sd (mhw$grain)
## [1] 0.4583
sd (mhw.c$grain)
## [1] 0.6764
# Mediani
median (mhw$grain)
## [1] 3.94
median (mhw.c$grain)
## [1] 3.915
# Interkvantilni razmaci
IQR (mhw$grain)
## [1] 0.6325
IQR (mhw.c$grain)
## [1] 0.67
# Postotak promjena sredine
(mean (mhw$grain) - mean (mhw.c$grain)) / mean (mhw$grain)
## [1] 0.02604
# Postotak promjena standardne devijacije
(sd(mhw$grain) - sd(mhw.c$grain)) / sd(mhw$grain)
## [1] -0.4759
# Postotak promjene mediana
(median (mhw$grain) - median (mhw.c$grain)) / median (mhw$grain)
## [1] 0.006345
# Postotak promjene interkvartilnog razmaka
(IQR (mhw$grain) - IQR (mhw.c$grain)) / IQR (mhw$grain)
## [1] -0.05929
Ako razdioba nije približno normalna tada nije uputno koristiti Pearsonov koeficijent korelacije; bolje je upotrebiti Spearmanov \( \rho \) koji je, zapravo Pearsonov koeficijent između ranka mjerenih veličina (rank svakom elementu pridruži njegov redni broj koji odgovara veličini tog elementa).
Zadatak Usporedite rank zrna i slame u prvih osam redaka orginalnog skupa podataka
head(cbind(mhw[, c("grain", "straw")], rank(mhw$grain), rank(mhw$straw)), n = 8)
## grain straw rank(mhw$grain) rank(mhw$straw)
## 1 3.63 6.37 123.0 254.5
## 2 4.07 6.24 299.0 219.5
## 3 4.51 7.05 445.5 356.5
## 4 3.90 6.91 228.0 329.0
## 5 3.63 5.93 123.0 136.0
## 6 3.16 5.59 23.5 70.5
## 7 3.18 5.32 26.0 36.0
## 8 3.42 5.52 62.5 59.0
Ako pogledamo, npr., prvi redak numeričkih rezultata vidimo da je za to mjerenje prinosa žita, po veličini, na 123. mjestu (do njih 500) što je ispod prvog kvantila a prinost slame na 254. mjestu što je oko mediana. Znači da je prinos zrna znatno niže rangiran nego prinos slame.
Ukupno, gledano, na svih osam redaka, rezultati ranga se (osim u prvom retku) ne razlikuju drastično, no nisu podjednaki.
Zadatak Grafički prikažite scatterplot ranka orginalnih podataka pored scatterplota orginalnih podataka, a ispod njih to isto ali za kontaminirane podatke.
par (mfrow = c(2, 2))
plot (rank(mhw$grain), rank(mhw$straw), xlab = "Rank zrna", ylab = "Rank slame",
pch = 1, main = "Orginalni podaci. Rank")
plot (mhw$grain, mhw$straw, xlab = "Prinos zrna", ylab = "Prinos slame",
pch = 20, main = "Orginalni podaci. Prinosi")
plot (rank(mhw.c$grain), rank(mhw.c$straw), xlab = "Rank zrna", ylab = "Rank slame",
pch = 1, main = "Kontaminirani podaci. Rank")
plot (mhw.c$grain, mhw.c$straw, xlab = "Prinos zrna", ylab = "Prinos slame",
pch = 20, main = "Kontaminirani podaci. Prinosi")
par (mfrow = c(1, 1))
Uspoređujući slike u svakom od redaka vidimo da su prikazi ranka više raspršeni nego prikazi prinosa. To se objašnjava činjenicom da se, pri prelasku sa numeričkih podataka na rank, dio informacija izgubio - izgubile su se numeričke vrijednosti prinosa a ostali su samo redni brojevi (po veličini) prinosa.
Dalje, može se vidjeti da je kontaminacija manje utjecala na prikaz ranka nege na prikaz prinosa (usporedba po stupcima); stoga očekujemo da će korelacija zasnovana na ranku biti manje podložna utjecaju kontaminacije.
Zadatak Izračunajte Spearmanov koeficijent korelacije između zrna i slame i usporedite ga sa Pearsonovim. Učinite to i za orginalne i za kontaminirane podatke.
cor(mhw$grain, mhw.c$straw, method = "pearson")
## [1] 0.7298
cor(mhw.c$grain, mhw.c$straw, method = "pearson")
## [1] 0.3597
cor(mhw$grain, mhw$straw, method = "spearman")
## [1] 0.7196
cor(mhw.c$grain, mhw.c$straw, method = "spearman")
## [1] 0.6168
Vidljivo je da je Pearsonov koeficijent korelacije bitno osjetljiviji na kontaminaciju podataka od Spearmanovog koeficijenta.
Slijedi da je u slučaju da podaci sadrže outlierse treba koristiti Spearmanov koeficijent korelacije.
Spearmanov koeficijent korelacije je statistička mjera monotonosti veze između uparenih podataka tj. taj koeficijent pokazuje u kojoj mjeri ako jedna varijabla naraste (padne) druga naraste (padne) i to bez zahtjeva linearnosti te veze. Ako jedna varijabla raste a druga pada tada je Spearmanov koeficijent negativan.
Rečeno na drugi način, taj koeficijetn opisuje koliko dobro se veza između dvaju varijabli može opisati monotonom funkcijom.
Zadatak Izračunajte Kendallov koeficijent korelacije i objasnite njegovo značenje
cor(mhw$grain, mhw$straw, method = "kendall")
## [1] 0.5401
cor(mhw.c$grain, mhw.c$straw, method = "kendall")
## [1] 0.4672
U statistici asocijacija je bilo koja veza između dvaju mjerenih veličina koja ih čini statistički zavisnim.
Kendallov korelacioni koeficijent (Kendallov tau), u statistici, se koristi kao mjera asocijativnosti dvije mjerene veličine.
Znači Kendallov koeficijent mjeri koliko su dvije varijable asocirane (povezane) da bi bile statistički zavisne.
Drugim rječima Kendallov koefijent mjeri veličinu veze dvaju skupa podataka koje te podatke čini statistički zavisnima.
Točnije rečeno tau je mjera korelacije ranka, tj. mjera sličnosti podataka kada su rangirani pomoću ranka.
Svojstva Kendellovog koeficijenta korelacije
Zadatak Prikažite numeričke rezultate primjene linearne regresije na orginalne i kontaminirane podatke.
( model.straw.grain = lm (straw ~ grain, data = mhw) )
##
## Call:
## lm(formula = straw ~ grain, data = mhw)
##
## Coefficients:
## (Intercept) grain
## 0.866 1.430
( model.straw.grain.c = lm (straw ~ grain, data = mhw.c) )
##
## Call:
## lm(formula = straw ~ grain, data = mhw.c)
##
## Coefficients:
## (Intercept) grain
## 4.678 0.478
Dobiveni rezultati ukazuju na drastičnu promjenu koeficijenata pravca regresije. Npr. za orginalne podatka je koeficjent smjera jednak 1.4305 što je iznad jedan (prinos slame raste brže nego prinos zrna i to za probližno 40%), a za kontaminirane je jednak 0.4777 što je znatno ispod jedan (prinos slame raste duplo sporije od prinosa zrna).
Zadatak Grafički prikažite dijagnostičke slike za regresiju orginalnih a zatim i kontaminiranih podataka.
par (mfrow = c(2, 2))
plot (model.straw.grain)
par (mfrow = c (1, 1))
par (mfrow = c(2, 2))
plot (model.straw.grain.c)
par (mfrow = c (1, 1))
Novi pravac regresije ćemo odrediti sa ciljem da umanjimo utjecaj kontaminiranih točaka pa ćemo dobiti “robusni” pravac regresije.
Postoji više načina da odredimo robusni pravac regresije. Ovdje ćemo
upotrebiti funkciju lqs() (least trimmed squares) iz bibilioteke MASS.
Ukratko, koeficijent \( \beta \) se određuje tako da se minimizira
zbroj kvadrata residuala za podskupove točaka iz orginalnog skupa
podataka (kardinalni broj podskupa ne smije biti manji od približno
polovice broje ukupnih točaka). Procedura nema analitičkog rješenje
već se mora izvoditi iterativno.
Točnije rečeno
(http://http://socserv.mcmaster.ca/jfox/Books/Companion/appendix/Appendix-Robust-Regression.pdf)
Odredimo početne koeficijente \( \alpha \) i \( \beta \) pravca regresije npr. po standardnoj metodi najmanjih kvadrata, uzimajući u obzir sve točke.
Izračunamo kvadrate rezidula u svim točkama i poredamo ih po veličini \[ e_1^2 \leq e_2^2 \leq \ldots \leq e_n^2 \]
Izaberemo prvih (najmanjih) \( m \) reziduala, oni određuju \( m, m < n \) točaka pa odredimo \( \alpha \) i \( \beta \) iz zahtjeva minimizacije rezidula u tim točkama. Kao što je već rečeno \( m \) je malo veće od polovice \( n \).
Procedure 2. i 3. ponavljamo do konvergencije koeficijenata.
require(MASS)
## Loading required package: MASS
(model.straw.grain.c.r <- lqs(straw ~ grain, data = mhw.c))
## Call:
## lqs.formula(formula = straw ~ grain, data = mhw.c)
##
## Coefficients:
## (Intercept) grain
## -0.0631 1.6857
##
## Scale estimates 0.548 0.531
Drugi način je Theil-Sen-ov pravac regresije. Pravac regresije ima oblik \( y = \alpha + \beta\cdot x \).
Sa svake dvije točke skupa podataka je određen koeficijent smjera. Tada je koeficijent smjera pravca regresije \( \beta \) jednak medijanu dobivenih koeficijenata smjera.
Kada je određen \( \beta \), za svaku točku \( (x_i, y_i) \) se izračuna odsječak pravca \( y_i - \beta x_i \) kroz tu točku i medijan dobivenih odsječaka je odsječak pravca regresije.
# Thiel-Sen
library("zyp")
## Loading required package: Kendall
(model.straw.grain.c.ts = zyp.sen(straw ~ grain, data = mhw.c))
##
## Call:
## NULL
##
## Coefficients:
## Intercept grain
## 1.25 1.34
Treći način je primjena funkcije rlm() iz bibilioteke MASS
(rlm je skraćenica od “robust linear model”)
(model.straw.grain.c.rlm = rlm(straw ~ grain, data = mhw.c))
## Call:
## rlm(formula = straw ~ grain, data = mhw.c)
## Converged in 12 iterations
##
## Coefficients:
## (Intercept) grain
## 2.6231 0.9874
##
## Degrees of freedom: 500 total; 498 residual
## Scale estimate: 0.679
Prikažimo na istoj slici kontaminirane točke i svih pet pravaca regresije.
with(mhw.c, plot(straw ~ grain))
abline(model.straw.grain.c.r)
abline(model.straw.grain.c, lty = 3, col = "magenta")
abline(1.25, 1.34, lty = 5, col = "darkgreen")
abline(model.straw.grain, lty = 2, col = "blue")
abline (model.straw.grain.c.rlm, lty = 6, col ="darkgray")
legend(1.5, 8.5, legend = c("robust fit, lqs()", "linear fit",
"Theil Sen fit",
"fit to uncontaminated", "rlm() fit"),
lty = c(1, 3, 5, 2, 6),
col = c("black", "magenta", "darkgreen", "blue", "darkgray"))
