6.4. Analiza varijanse – ANOVA (klasična ANOVA i Welch ANOVA)

U slučaju testiranja hipoteza na osnovu više od dva uzorka, umjesto primjene t-testa na sve moguće kombinacije dva uzorka od ukupnog broja uzoraka i otežanog kontrolisanja greške tipa I, može se primijeniti analiza varijanse na osnovu koje testiramo hipotezu da li su razlike aritmetičkih sredina testiranog obilježja u svim uzorcima jednake ili da li postoji makar jedan par uzoraka čija je razlika aritmetičkih sredina statistički značajna. Pretpostavke koje su važile za t-test (normalnost i homogenost varijansi) važe i za klasičnu analizu varijanse uz dodatnu pretpostavku o nezavisnosti jedinica uzorka. Klasična analiza varijanse podrazumijeva postupak poređenja faktorske i rezidualne sume kvadrata odstupanja. Naime, ukupan varijabilitet nekog obilježja može se podijeliti na sumu kvadrata odstupanja unutar uzoraka (rezidualni varijabilitet) i između uzoraka (faktorski varijabilitet). Formule za računanje faktorske i rezidualne sume kvadrata odstupanja date su sljedećim izrazima:

$$V_f = \frac{\sum_{i=1}^{r}n_i(\bar{X_i} - \bar{X})^{2}}{r - 1}$$

gdje $V_f$ predstavlja faktorsku varijansu, $n_i$ predstavlja veličinu $i$-tog uzorka, $\bar{X_i}$ aritmetičku sredinu $i$-tog uzorka, $\bar{X}$ zajedničku aritmetičku sredinu svih uzoraka, dok broj uzoraka;

$$V_r = \frac{\sum_{i=1}^{r}\sum_{j=1}^{n}n_i(\bar{X_{ij}} - \bar{X_i})^{2}}{n - r}$$

gdje $V_r$ predstavlja rezidualnu varijansu, $X_{ij}$ predstavlja $j$-tu vrijednost obilježja u $i$-tom uzorku, $\bar{X}$ aritmetičku sredinu $i$-tog uzorka.

Odnos ova dva varijabiliteta, u slučaju ispunjenosti pretpostavki ANOVA, prati raspored sa stepenima slobode $r -1$ i $n - r$, gdje je broj uzoraka (modaliteti faktora varijabiliteta), a $n$ veličina uzoraka:

$$F = \frac{V_f}{V_r}$$

Primjer 45: U fajlu LGD.csv dati su podaci o realizovanim vrijednostima LGD-ija (kolona LGDR), kao i LGD segment (kolona segment.e) kome određeni kredit pripada. Pod pretpostavkom ispunjenosti uslova normalnosti, homogenosti varijansi i nezavisnosti opservacija, za nivo značajnosti od 5% ispitati da li postoji bar jedan par prosječnih vrijednosti LGD-ija datih segmenata koji su statistički različit jedan od drugog.

> db <- read.csv("LGD.csv", header = TRUE)
> alpha <- 0.05
> n <- nrow(db)
> n
[1] 500
> r <- length(unique(db$segment.e))
> r
[1] 4
> #faktorska varijansa (Vf)
> ng <- tapply(db$LGDR, db$segment.e, length)
> ng 
segment1 segment2 segment3 segment4 
      41      233      197       29 
> asg <- tapply(db$LGDR, db$segment.e, mean)
> asg
 segment1  segment2  segment3  segment4 
0.3027159 0.4555191 0.6006755 0.7529565 
> fsk <- (sum(ng * (asg - mean(db$LGDR))^2)) / (r - 1)
> fsk
[1] 1.91905
> #rezidualna varijansa (Vr)
> asg.vektor <- asg[db$segment.e]
> rsk <- sum((db$LGDR - asg.vektor)^2) / (n - r)
> rsk
[1] 0.002294377
> #F statistika
> f.stat <- fsk / rsk
> f.stat
[1] 836.4143
> #p vrijednost
> p.val <- 1 - pf(f.stat, r - 1, n - r)
> p.val < alpha
[1] TRUE
> #r funkcija
> anova.r <- aov(LGDR ~ segment.e, data = db)
> summary(anova.r)
             Df Sum Sq Mean Sq F value Pr(>F)    
segment.e     3  5.757  1.9190   836.4 <2e-16 ***
Residuals   496  1.138  0.0023                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Pošto je izračunata p vrijednost (p.val) manja od nivoa značajnosti 5%, zaključujemo da prosječni LGD bar jednog segmenta (segment.e) jeste statistički značajno različit u odnosu na prosječni LGD ostalih segmenata.

Da bi se utvrdilo kod kojeg je para segmenata razlika u prosječnom LGD-iju statistički značajna, možemo primijeniti t-test sa korekcijom p vrijednosti zbog višestrukog poređenja aritmetičkih sredina. Primjenju se različite metode korekcije p vrijednosti, a jedna od korišćenih je tzv. Bonferroni metoda, koja penalizuje p vrijednost t-testa množeći je sa brojem parova na koje je primijenjen t-test. Pored ovog metoda često korišćen metod kod klasične ANOVA jeste i tzv. Tukeyev metod razlike višestrukog t-testa (za dodatni primjer pogledati R funkciju ?TukeyHSD).

Primjer 46: Na osnovu podataka iz Primjera 45, i za isti nivo značajnosti, ispitati kod kojih parova segmenata prosječnog LGD-ija postoji statistički značajna razlika. Za korekciju izračunatih pojedinačnih p vrijednosti t-testa, primijeniti Bonferroni metodu.

> ug <- sort(unique(db$segment.e))
> ug
[1] "segment1" "segment2" "segment3" "segment4"
> #nezavisni t-test aritmetickih sredina segmenta 1 i segmenta 2
> par12 <- t.test(x = db$LGDR[db$segment.e%in%ug[1]], y = db$LGDR[db$segment.e%in%ug[2]], paired = FALSE)
> par12

        Welch Two Sample t-test

data:  db$LGDR[db$segment.e %in% ug[1]] and db$LGDR[db$segment.e %in% ug[2]]
t = -20.493, df = 58.006, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1677286 -0.1378776
sample estimates:
mean of x mean of y 
0.3027159 0.4555191 

> par12$p.val * ncol(combn(length(ug), 2))
[1] 1.815098e-27
> #korisnicka funkcija za generalizaciju p vrijednosti visestrukog t-testa i Bonferroni korekcije
> p.val.vtt <- function(x, y) { 
+ ug <- sort(unique(y))
+ tc <- combn(length(ug), 2)
+ tcn <- ncol(tc)
+ p.val.corr <- rep(NA, tcn)
+ for(i in 1:tcn) {
+ tc.l <- tc[, i]
+ seg.x <- tc.l[1]
+ seg.y <- tc.l[2]
+ p.val.itt <- t.test(x = x[y%in%ug[seg.x]], y = x[y%in%ug[seg.y]], paired = FALSE)$p.val 
+ p.val.corr[[i]] <- p.val.itt * tcn 
+ }
+ names(p.val.corr) <- apply(tc, 2, function(x) paste0(ug[x[1]], " VS ", ug[x[2]]))
+ return(p.val.corr)
+ }
> p.val.vtt(x = db$LGDR, y = db$segment.e)
segment1 VS segment2 segment1 VS segment3 segment1 VS segment4 segment2 VS segment3 
        1.815098e-27         2.956596e-45         1.069345e-41        9.848637e-109 
segment2 VS segment4 segment3 VS segment4 
        4.188637e-26         5.093785e-17 
> #r funkcija za visestruki t-test sa korekcijom p vrijednosti
> pw.tt <- pairwise.t.test(db$LGDR, db$segment.e, p.adjust.method = "bonferroni", pool.sd = FALSE)
> pw.tt$p.value
             segment1      segment2     segment3
segment2 1.815098e-27            NA           NA
segment3 2.956596e-45 9.848637e-109           NA
segment4 1.069345e-41  4.188637e-26 5.093785e-17 

U slučaju kada pretpostavka homogenosti varijanse između uzoraka nije ispunjena, veću moć otkrivanja statistički značajne razlike između aritmetičkih sredina ima tzv. Welch analiza varijanse. Test statistika Welch ANOVA računa se po sljedećoj fomuli:

$$F = \frac{\frac{\sum_{j=1}^{K}w_j(\bar{x_j} - \bar{x}^{'})}{K-1}}{1 + \frac{2(K-2)}{K^2-1}\sum_{j=1}^{K}(\frac{1}{n_j-1})(1 - \frac{w_j}{w})^{2}}$$

za izračunate:

$$w_j = \frac{n_j}{s_j^2}$$

$$w = \sum_{j=1}^{K}w_j$$

$$\bar{x}^{'} = \frac{\sum_{j=1}^{K}w_j\bar{x_j}}{w}$$

gdje $K$ predstavlja broj uzoraka za koje računamo razliku između aritmetičkih sredina, $s_{j}^{2}$ varijansu analiziranog obilježja za $j$-ti uzorak i $n_j$ broj opservacija analiziranog obilježja unutar $j$-tog uzorka. Izračunata statistika prati $F$ raspored sa $df1$ i $df2$ stepeni slobode definisane sljedećim formulama:

$$df1 = K - 1$$

$$df2 = \frac{K^2 - 1}{3\sum_{j=1}^{K}(\frac{1}{n_j-1})(1- \frac{w_j}{w})^2}$$

Primjer 47: Za podatke iz Primjera 45, pod pretpostavkom da uslov homogenosti varijanse nije ispunjen, ispitati da li postoji bar jedan par prosječnih vrijednosti realizovanih LGD-ija datih segmenata koji je statistički različit jedan od drugog.

> alpha <- 0.05
> #broj uzoraka
> K <- length(unique(db$segment.e))
> K
[1] 4
> #velicina uzoraka
> nj <- tapply(db$LGDR, db$segment.e, length)
> nj
segment1 segment2 segment3 segment4 
      41      233      197       29 
> #varijansa LGDR unutar uzoraka
> s2j <- tapply(db$LGDR, db$segment.e, var)
> nj
segment1 segment2 segment3 segment4 
      41      233      197       29 
> #ponder
> wj <- nj / s2j
> wj
 segment1  segment2  segment3  segment4 
 21740.97 104153.85  81096.60  12038.12 
> #suma pondera
> w <- sum(wj)
> w
[1] 219029.5
> #prosjecni LGDR unutar uzoraka
> xbar.g <- tapply(db$LGDR, db$segment.e, mean)
> xbar.g
 segment1  segment2  segment3  segment4 
0.3027159 0.4555191 0.6006755 0.7529565 
> #zajednicki prosjecni LGDR (ponderisan u odnosu na pojedinacne uzorke)
> xbar.t <- sum(wj * xbar.g / w)
> xbar.t
[1] 0.5104441
> #suma varijabiliteta izmedju i unutar uzoraka
> var.between <- sum(wj * (xbar.g - xbar.t)^2) / (K - 1)
> var.within <- 1 + (2 * (K - 2) / (K^2 - 1)) * sum((1 / (nj - 1)) * (1 - wj / w)^2)
> var.between 
[1] 873.5349
> var.within
[1] 1.01477
> #F statistika
> f.stat <- var.between / var.within
> f.stat
[1] 860.8203
> #stepeni slobode     
> df1 <- K - 1
> df2 <- (K^2 - 1) / (3 * sum((1 / (nj - 1)) * (1 - wj / w)^2))
> df1
[1] 3
> df2
[1] 90.27112
> #p vrijednost
> p.val <- pf(f.stat, df1, df2, lower.tail = FALSE)
> p.val
[1] 2.905626e-66
> p.val < alpha
[1] TRUE
> #r funkcija
> w.anova.r <- oneway.test(formula = LGDR ~ segment.e, data = db)
> w.anova.r

        One-way analysis of means (not assuming equal variances)

data:  LGDR and segment.e
F = 860.82, num df = 3.000, denom df = 90.271, p-value < 2.2e-16

> w.anova.r$p.val
[1] 2.905626e-66

Kako bi se utvrdili parovi prosječnih vrijednosti LGD-ija kod kojih postoji statistički značajna razlika i pod pretpostavkom da uslov homogenosti varijanse nije ispunjen (Welch ANOVA inicijalno primijenjena), primjenjujemo tzv. Games-Howell test, čija test statistika ima sljedeću formu:

$$GH = \frac{\bar{X_i} - \bar{X_j}}{\sqrt{\frac{s_i^{2}}{n_i} + \frac{s_j^{2}}{n_j}}}$$

gdje $\bar{X_i}$ i $\bar{X_j}$ predstavljaju par aritmetičkih sredina uzoraka koji se porede, dok $s_i^{2}$, $s_j^{2}$, $n_i$i $n_j$ predstavljaju varijansu i veličinu $i$-tog odnosno $j$-tog uzorka. Izračunata test statistika prati Tukeyev raspored sa standardnom devijacijom $\sigma$ i stepenima slobode $df1$ i $df2$ .

$$\sigma^{2} = \frac{\frac{s_i^{2}}{n_i} + \frac{s_j^{2}}{n_j}}{2}$$

$$df1 = K$$

$$df2 = \frac{(\frac{s_i^{2}}{n_i} + \frac{s_j^{2}}{n_j})^{2}}{\frac{(\frac{s_i^2}{n_i})^2}{n_i-1} + \frac{(\frac{s_j^2}{n_j})^2}{n_j-1}}$$

Primjer 48: Za podatke iz Primjera 47, ispitati kod kojih parova segmenata prosječnog LGD-ija postoji statistički značajna razlika. Zbog pretpostavke da uslov homogenosti varijanse nije ispunjen, primijeniti Games-Howell test.

> #korisnicka funkcija
> gh.test <- function(samples, values, sig.level) {
+ comb <- combn(sort(unique(samples)), 2)
+ ncomb <- ncol(comb)
+ K <- length(unique(samples))
+ n.g <- tapply(values, samples, length)
+ xbar.g <- tapply(values, samples, mean)
+ s2.g <- tapply(values, samples, var)
+ res <- vector("list", ncomb)
+ for(i in 1:ncomb) {
+ #par prosjecnih LGD-ijeva
+ pair <- comb[, i]
+ pair.d <- paste0(pair, collapse = " - ")
+ #razlika para prosjecnih LGD-ijeva
+ mean.d <- as.numeric(xbar.g[pair[1]] - xbar.g[pair[2]])
+ #standardna greska testa
+ se <- as.numeric(sqrt(s2.g[pair[1]] /  n.g[pair[1]] + s2.g[pair[2]] /  n.g[pair[2]]))
+ #GH test statistika
+ gh.stat <- abs(mean.d) / se
+ #stepeni slobode
+ df1 <- K
+ df2 <- se^4 / ((s2.g[pair[1]] / n.g[pair[1]])^2 / (n.g[pair[1]] - 1) + 
+ (s2.g[pair[2]] / n.g[pair[2]])^2 / (n.g[pair[2]] - 1))
+ #Tukey-ijeva studentizovana p vrijednosti
+ p.val <- ptukey(gh.stat * sqrt(2), df1, df2, lower.tail = FALSE)
+ #rezultat testa i poredjenje izracunare p vrijednosti i odabranog nivoa znacajnosti
+ significance <- p.val < alpha
+ #sumiranje rezultata
+ res.l <- data.frame(pair = pair.d, mean.d, se, GH = gh.stat, df1, df2, p.val, significance)
+ row.names(res.l) <- NULL
+ res[[i]] <- res.l
+ }
+ res <- do.call("rbind", res)
+ return(res)
+ } 
> gh.test(samples = db$segment.e, values = db$LGDR, sig.level = alpha)
                 pair     mean.d          se       GH df1       df2        p.val significance
1 segment1 - segment2 -0.1528031 0.007456360 20.49300   4  58.00623 1.490907e-11         TRUE
2 segment1 - segment3 -0.2979595 0.007637217 39.01415   4  63.39201 1.331280e-11         TRUE
3 segment1 - segment4 -0.4502406 0.011360703 39.63140   4  55.64907 5.987655e-12         TRUE
4 segment2 - segment3 -0.1451564 0.004683178 30.99528   4 410.03454 0.000000e+00         TRUE
5 segment2 - segment4 -0.2974374 0.009626559 30.89759   4  34.79044 0.000000e+00         TRUE
6 segment3 - segment4 -0.1522810 0.009767314 15.59088   4  36.81383 0.000000e+00         TRUE

Često su u praksi pretpostavke normalnosti i homogenosti varijanse narušene. U tim situacijama i posebno u kombinaciji sa uzorcima male veličine (30 ili manje opservacija) veću moć otkrivanja statistički značajnih razlika imaju neparametarski testovi. Kod većih uzoraka testovi jednakosti aritmetičkih sredina usljed djelovanja centralne granične teoreme prilično su stabilni, ali testovi disperzije su izuzetno osjetljivi. U narednim primjerima biće prikazani neparametarski testovi koji su analogni prethodnom predstavljenim parametarskim testovima.