NMST539 | Cvičenie 4

zdrojový Rmd súbor (kódovanie UTF8)

---

Outline štvrtého cvičenia:

• náhodné matice a niektoré mnohorozmerné rozdelenia – Wishartove, Hottellingovo a Wilkovo Lambda rozdelenie;
• mnohorozmerné testy pre vektor strednej hodnoty (mnohorozmerné analógie jednovýberových a dvojvýberových t-testov);

Štatisticky program R je k dispozícii (GNU public licence) na adrese https://www.r-project.org

RStudio – “user-friendly” interface (jeden z mnohých, ktoré na internete nájdete): RStudio.

Užitočné návody a manuály pre prácu s programom R:

• Bína, V., Komárek, A. a Komárková, L.: Jak na jazyk R. (PDF súbor)
• Komárek, A.: Základy práce s R. (PDF súbor)
• Kulich, M.: Velmi stručný úvod do R. (PDF súbor)
• De Vries, A. a Meys, J.: R for Dummies. (ISBN-13: 978-1119055808)

Odborná literatúra:

• Hardle, W. and Simar.L.: Applied Multivariate Statistical Analysis. Springer, 2015
• Mardia, K., Kent, J., and Bibby, J.:Multivariate Analysis, Academic Press, 1979.

---

1. Wishartovo pravdepodobnostné rozdelenie

Wishartovo rozdelenie je mnohorozmerným zobecnením jednorozmerného $\chi^2$ rozdelenia, ktoré sa štandardne používa pre inferenciu o neznámom parametri rozptylu na základe jednorozmerného náhodného výberu. Rozdelenie nesie názov po svojom autorovi, John-ovi Wishart-ovi, ktorý toto rozdelenie formuloval v roku 1928. Wishartovo rozdelenie je základným nástrojom pre analýzu/inferenciu ohľadom variančnej-kovariančnej matice na základe mnohorozmerného náhodného výberu $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n}$, pre $p$ rozmerné náhodné vektory $\boldsymbol{X}_{i} = (X_{i 1}, \dots, X_{i p})^\top \in \mathbb{R}^p$ a rozsah náhodného výberu $n \in \mathbb{N}$.

Wishartovo rozdelenie je mnohorozmerným zobecnením jednorozmerného $\chi^2$ rozdelenia v nasledujúcom zmysle: pre mnohorozmerný náhodný výber $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n}$ (stĺpcové náhodné vektory) z $p$ rozmerného normálneho rozdelenia $N_{p}(\boldsymbol{0}, \Sigma)$ s nulovým vektorom stredných hodnôt a variančnou-kovariančnou maticou $\Sigma$ (symetrická a pozitívne definitná matica) má príslušná kvadratická forma Wishartovo rozdelenie

$$\mathbb{X}^{\top}\mathbb{X} \sim W_{p}(\Sigma, n),$$ s parametrami $p \in \mathbb{N}$ (dimenzia), $n \in \mathbb{N}$ (prozsah náhodného výberu) a $\Sigma$ (variančná-kovariačná matica náhodného výberu $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n}$ kde $\mathbb{X} = (\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n})^{\top}$). Podobne ako v jednorozmernom prípade pre náhodný výber $X_{1}, \dots, X_{n}$ (náhodné veličiny) z normálelneho rozdelenia $N(0,1)$, má príslušná kvadratická forma $$\boldsymbol{X}^{\top}\boldsymbol{X} \sim \chi_{n}^2,$$ $\chi^2$ rozdelenie, pre $\boldsymbol{X} = (X_{1}, \dots, X_{n})^\top$. Je zrejmé, že pre náhodný výber s rozsahom $n \in \mathbb{N}$ z jednorozmerného normálneho rozdelenia $N(0,1)$ má náhodná veličina $\mathbb{X}^\top\mathbb{X}$ (jednorozmerné) Wishartovo rozdelenie $W_{1}(1, n)$, ktoré je ekvivalentné s $\chi_{n}^2$ rozdelením.

Wishartovo rozdelenie predstavuje rodinu náhodných rozdelení pre definovaných na symetrických, pozitívne semi-definitných .náhodných maticiach. Príslušná hustota Wishartovho rozdelenia má nasledujäci tvar:

$$f(\mathcal{X}) = \frac{1}{2^{np/2} |\Sigma|^{n/2} \Gamma_{p}(\frac{n}{2})} \cdot |\mathcal{X}|^{\frac{n - p - 1}{2}} e^{-(1/2) tr(\Sigma^{-1}\mathcal{X})},$$ kde $\mathcal{X}$ je náhodná matica typu $p\times p$ a $\Gamma_{p}(\cdot)$ predstavuje mnohorozmerné zobecnenie jednorozmernej Gamma funkcie $\Gamma(\cdot)$. V programe R sú k dispozícii rôzne možnosti (knižnice a príkazy) pre prácu s Wishartovým rozdelením. Niektoré z týchto nástrojov využijeme v následujúcej časti.

Samostatne (teoretické a praktické úlohy)

---

• Uvažujte mnohorozmerný náhodný výber $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n}$ z $p \in \mathbb{N}$ mnohorozmerného normálneho rozdelenia $N_{p}(\boldsymbol{\mu}, \Sigma)$, s vektorom stredných hodnôt $\boldsymbol{\mu} \in \mathbb{R}^p$ a variančnou-kovariančnou (pozitívne definitnou) maticou $\Sigma$. Ukážte, že náhodná matica definovaná predpisom $n\mathbb{X}^\top \mathcal{H}\mathbb{X}$ má Wishartovo rozdelenie $W_{p}(\Sigma, n - 1)$, kde $\mathcal{H}$ idenpotentná, tzv. centrovacia matica $\mathcal{H} = \mathbb{I}_{n} - \frac{1}{n}\boldsymbol{1}_{n}\boldsymbol{1}_{n}^\top$.




• V štandardnej inštalácii programu R je pre Wishartovo rozdelenie určená funkcia rWishart() – t.j. generátor náhodných matíc s Wishartovým rozdelením. Pomocou helpu v programe R a uvedených príkladov naštudujte implementáciu daného príkazu.




• Ďalšie užitočné príkazy a funkcie pre prácu s Wishartovým rozdelením sú napr. funkcia Wishart() v knižnici ‘MCMCpack’, funkcia Wishart() v knižnici ‘mixAK’, funkcia rwishart() z knižnice ‘dlm’, alebo dist.Wishart() z knižnice ‘LaplacesDemon’.

---

Pre jednoduchú ilustráciu použijeme jednorozmerný generátor: náhodný výber z Wishartovho rozdelenia $W_{1}(\Sigma = 1, n = 10)$ (resp. ekvivalentne z $\chi^2$ rozdelenia s $n = 10$ stupňami voľnosti) získame pomocou pomocu funkcie rWishart() ako

S <- as.matrix(1) ### variancna-kovariancna matica (i.e., rozptyl)
sample <- rWishart(10, df = 10, S)

Pre obecný rozmer $p \in \mathbb{N}$ použijeme analogický príkaz s vhodne definovanou variančnou-kovariačnou maticou $\Sigma$:

Sigma <- cbind(c(1,0,0,0), c(0,1,0,0), c(0,0,1,0), c(0,0,0,1)) ### matica typu 4x4
rWishart(2, df = 100, Sigma) ### sample of size two

## , , 1
## 
##             [,1]       [,2]      [,3]       [,4]
## [1,] 102.2004664 -0.5779487 -2.731336  28.946426
## [2,]  -0.5779487 88.5347475 10.996492  13.012143
## [3,]  -2.7313359 10.9964921 98.400893   2.851688
## [4,]  28.9464262 13.0121427  2.851688 106.993455
## 
## , , 2
## 
##            [,1]       [,2]       [,3]        [,4]
## [1,]  88.428345 -3.1209897   5.870712 -16.9115631
## [2,]  -3.120990 97.4225809   8.583250   0.6412193
## [3,]   5.870712  8.5832501 112.053916  12.1145055
## [4,] -16.911563  0.6412193  12.114506 101.8873937

V prípade, že $p = 1$, môžeme ekvivalentnosť Wishartovho rozdelenia $W_{1}(\Sigma = 1, n)$ a $\chi^2$ rozdelenia s $n$ stupňami voľnosti jednoducho overiť aj vizuálne, napr. pomocou histogramu, alebo príslušného neparametrického odhadu hustôt.

set.seed(1234)
n <- 10 ### pocet stupnov volnosti
sampleWishart <- rWishart(5000, df = n, S)
sampleChiSq <- rchisq(5000, df = n) 

par(mfrow = c(1,2))
hist(sampleWishart, col = "lightblue", main = expression(paste("Wishart Distribution", sep ="")), xlim = c(0, 40), breaks = 30, freq = F)
lines(density(sampleWishart), col = "red", lwd = 2)
hist(sampleChiSq, col = "lightblue", main = expression(paste(chi^2, "Distribution", sep = "")), xlim = c(0,40), breaks = 30, freq = F)
lines(density(sampleChiSq), col = "red", lwd = 2)

Dvojrozmerné Wishartovo rozdelenie (rozdelenie náhodných symetrických matíc typu $2 \times 2$) je už ale výrazne náročnejšie vizualizovať vhodným spôsobom. Jedná z možnosti je napr. následujúci graf, ktorý zobrazuje $2 \times 2$ maticu ako usporiadaný vektor - t.j., náhodná matica $\boldsymbol{X}_i$ je tvorená dvoma vektormi $(x_{11}, x_{21})^\top$ a $(x_{12}, x_{22})^\top$. Tieto vektory definujú v dvojrozmernej $xy$ rovine ($xy$ scatterplot) počiatočný a koncový bod úsečky. Úsečky su následne vyzualizované v grafe (počiatočný bod – t.j. bod určený vektorom $(x_{11}, x_{21})^\top$ je v grafe zvýraznený).

Graf možno nie je úplne intuitívny, ale umožňuje napríklad jednoznačnu a priamočiaru reprodukciu pôvodných hodnôt náhodného výberu.

sample2DWishart <- rWishart(100, df = 10, cbind(c(1,1), c(1,3)))

plot(0,0, pch = "", xlim = c(min(sample2DWishart[1,,]), max(sample2DWishart[1,,])), ylim = c(min(sample2DWishart[2,,]), max(sample2DWishart[2,,])),
     xlab = "First column vector", ylab = "Second column vector")
for(i in 1:100){
  lines(sample2DWishart[1,,i], sample2DWishart[2,,i], col = "red")
  points(sample2DWishart[1,1,i], sample2DWishart[2,1,i], pch = 21, bg = "red", cex = 0.8)
}

Iné možnosti grafickej reprezentácie náhodného výberu z Wishartovho rozdelenia by boli napríklad vizualizácie založené na vlastných čislach, vlastných vektoroch, alebo ďalšie. Väčšinou ale neumožňujú takto jednoduchú spätnú rekonštrukciu náhodných matíc na základe samotného grafu.

Samostatne

---

• V predcházajúcom príklade použijte odlišnú štruktúru pre variančnú-kovariančnú maticu $\Sigma$ alebo voľbu počtu stupňov voľnosti a pomocou obrázku sledujte vplyv týchto parametrov na celkovú štruktúru náhodného výberu reprezentovanú pomocou $xy$ scatterplotu. Zamerajte sa hlavne na vplyv jednotlivých parametrov variančnej-kovariančnej matice $\Sigma$.




• Uvažujte obecnú maticu $A \in \mathbb{R}^{p \times q}$ pre nejaké vhodné $q \in \mathbb{N}$, tak, aby $p \neq q$. Opäť pomocou grafickej vizualizácie oveřte (napríklad pomocou histogramu, alebo KDE), že $A^\top \mathbb{M} A \sim W_{q}(A^\top\Sigma A, n)$, kde $\mathbb{M} \sim W_{p}(\Sigma, n)$;

---

2. Hotellingovo $\boldsymbol{T^2}$ rozdelnie

Podobným spôsobom, ako je v jednorozmernom prípade definované $t$-rozdelenie s $n$ stupňami voľnosti, môžeme defnovať aj jeho mnohorozmerné zobecnenie. V jednorozmernom prípade sa jedná o náhodnú veličinu s normálnym $N(0,1)$ rozdelením, ktorá je podelená inou, na nej nezávislou náhodnou veličinou, ktorej kvadrát má $\chi^2$ s $n$ stupňami voľnosti (a navyše je štandardizovaná odmocninou z týchto stupńov voľnosti).

Mnohorozmerné zobecnenie (t.j., náhodná veličina s Hotellingovým $T^{2}$ rozdelením) je definované predpisom

$$n \boldsymbol{Y}^{\top} \mathbb{M}^{-1} \boldsymbol{Y} \sim T^{2}(n, p),$$ kde $p \in \mathbb{N}$ je dimenzia/rozmer náhodného vektoru $Y \sim N_{p}(0, \mathbb{I})$ a $n \in \mathbb{N}$ je parameter Wishartovho rozdelenia náhodnej matice $\mathbb{M} \sim W_{p}(\mathbb{I}, n)$. Analogicky ako v jednorozmernom prípade, prepokládame, že náhodný vektor $\boldsymbol{Y}$ je nezávsilý od náhodnej matice $\mathbb{M}$.

V špeciálnom prípade, pre $p = 1$, dostaneme jednorozmerné Fisherovo F rozdelenie s jedným a s $n$ stupňami voľnosti (čo je vlastne ekvivalentné s rozdelením kvadrátu náhodnej veličiny s $t$ rozdelením s $n$ stupňami voľnosti). Hotellingovo $T^2$ rozdelenie s prametrami $p, n \in \mathbb{N}$ preto možno považovať aj za mnohorozmerné zobecnenie Fisherovho F rozdelenia. Medzi oboma rozdeleniami je dokonca jednoznačná analytická súvislosť, ktorú je možné vyjadriť pomocou vzťahu

$$T^{2}(p, n) \equiv \frac{n p}{n - p + 1}F_{p, n - p + 1}.$$

Náhodný výber z jednorozmerného Fisherovho rozdelenia preto môže byť efektívne použitý aj k tomu, aby sme získali kritické hodnoty mnohorozmerného Hotellingovho $T^2$ rozdelenia. Príslušná transformácia medzi Hotellingovým $T^2$ rozdelením a Fisherovým F rozdelením závisí iba na parametroch $n, p \in \mathbb{N}$.

Úlohu týchto parametrov (resp. efekt týchto parametrov) v Hotellingovom $T^2$ rozdelení môžeme opäť jednoducho vizualizovať, napríklad pomocou náhodného generátora z Fisherovho F rozdelenia (v programe R príkaz rf()).

samples <- NULL
samples <- cbind(samples, rf(n = 1000, df1 = 1, df2 = 1))
samples <- cbind(samples, rf(n = 1000, df1 = 1, df2 = 10))

samples <- cbind(samples, rf(n = 1000, df1 = 10, df2 = 10))
samples <- cbind(samples, rf(n = 1000, df1 = 10, df2 = 100))

samples <- cbind(samples, rf(n = 1000, df1 = 100, df2 = 10))
samples <- cbind(samples, rf(n = 1000, df1 = 100, df2 = 100))

samples <- cbind(samples, rf(n = 1000, df1 = 1000, df2 = 100))
samples <- cbind(samples, rf(n = 1000, df1 = 1000, df2 = 1000))

par(mfrow = c(4,2))
for (i in 1:dim(samples)[2]){
  hist(samples[,i], xlab=paste("Sample no.", i, sep = ""), col = "yellow", main = "", freq = F, breaks = 20)
  lines(density(samples[,i]), col = "red")
}

Pre pripomenutie, náhodná veličina $X$ s Fisherovým F rozdelením $F_{df_1, df_2}$ má strednú hodnotu a rozptyl definované predpismi:

• $E X = \frac{df_{2}}{df_{2} - 2}$, pre $df_2 > 2$;
• $Var(X) = \frac{2 df_{2}^2 (df_1 + df_2 - 2)}{df_1 (df_{2} - 2)^2 (df_2 - 4)}$ pre $df_2> 4$;

Obecne platí, že pre mnohorozmerný náhodný výber $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n} \sim N_{p}(\boldsymbol{\mu}, \Sigma)$ s neznámym vektorom stredných hodnôt $\boldsymbol{\mu} \in \mathbb{R}^{p}$ a neznámou variančnou-kovariančnou maticou $\Sigma$, dostaneme pre výberový priemer $\overline{\boldsymbol{X}}$ rozdelenie

$$(n - 1)\Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}\Big)^\top \mathcal{S}^{-1} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}\Big) \sim T^2(p, n - 1),$$

čo lze ekvivalentne zapísať aj do tvaru

$$\frac{n - p}{p} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}\Big)^\top \mathcal{S}^{-1} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}\Big) \sim F_{p, n - p}.$$

Toto je následne možné aplikovať pri úlohach o inferencii pre neznámy vektor stredných hodnôt – napr. pri štatistických testoch, alebo pri konštrukcii konfidenčných oblastí pre vektor neznámych stredných hodnôt $\boldsymbol{\mu}$, kde $\boldsymbol{\mu} = (\mu_{1}, \dots, \mu_{p})^{\top}$.

V praxi sa často využíva namiesto confidenčného regiónu pre $\boldsymbol{\mu}$ (čo môže byť nepraktické hlavne pre vyšší rozmer dimenzie $p$) radšej súbor intervalov spoľahlivosti pre jednotlivé zložky $\boldsymbol{\mu}$, tak aby bola celková pravdepodobnosť pokrytia pod kontrolou – najčastejšie vyžadujeme pokrytie minimálne $(1 - \alpha)\times 100~\%$ pre vhodné a dostatočne malé hodnoty $\alpha \in (0,1)$.

Analogicky, v prípade štatistického testu definujeme dvojicu hypotéz

$$H_{0}: \boldsymbol{\mu} = \boldsymbol{\mu}_{0} \in \mathbb{R}^{p}$$

$$H_{1}: \boldsymbol{\mu} \neq \boldsymbol{\mu}_{0} \in \mathbb{R}^{p}$$

a využijeme testovú štatistiku

$$(n - 1)\Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}_{0}\Big)^\top \mathcal{S}^{-1} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}_{0}\Big),$$

ktorá má za platnosti nulovej hypotézy $T^2(p, n - 1)$ rozdelenie. Ekvivalentne, testová štatistika vyjadrená v tvare

$$\frac{n - p}{p} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}_{0}\Big)^\top \mathcal{S}^{-1} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}_{0}\Big)$$

má za platnosti nulovej hypotézy Fisherovo $F$-rozdelenie s $p$ a $n - p$ stupňami voľnosti.

V programe R je k dispozícii knižnica DescTools – nutné nainštalovať pomocou príkazu install.packages("DescTools")a následne inicializovať pomocou príkazu library(DescTools). Samotný test nulovej hypotézy vyššie je potom možné uskutočniť pomocou funkcie HotellingsT2Test() (pre konkrétnu implementáciu príkazu použijte help k danej funkcii).

Analogickým spôsobom je možné skonštruovať aj konfidenčnú oblasť - tzv. konfidenčný elipsoid $\boldsymbol{\mu} \in \mathbb{R}^p$. Keďže platí, že $$\frac{n - p}{p} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}_{0}\Big)^\top \mathcal{S}^{-1} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}_{0}\Big) \sim F_{p, n - p},$$ tak nasledujúca množina $$\left\{\boldsymbol{\mu} \in \mathbb{R}^p;~ \Big( \overline{\boldsymbol{X}} - \boldsymbol{\mu}\Big)^\top \mathcal{S}^{-1} \Big( \overline{\boldsymbol{X}} - \boldsymbol{\mu}\Big) \leq \frac{p}{n - p} F_{p, n- p}(1 - \alpha) \right\}$$ je konfidenčnou oblasťou pre neznámy vektor stredných hodnôt $\boldsymbol{\mu} \in \mathbb{R}^p$ s pravdepodobnosťou pokrytia $\alpha = 1 - \alpha$. Jedná sa o ellipsoid v $\mathbb{R}^p$.

Ilustračný príklad v 2D:

library(mvtnorm)
s=matrix(c(1,-0.5,-0.5,1),2);x=seq(-3,3,by=0.015) ### variance-covariance matrix

contour(x,x,outer(x,x, function(x,y){dmvnorm(cbind(x,y),sigma=s)}))
n <- 20
X <- rmvnorm(n,sigma=s)
m <- apply(X,2,mean)
S <- cov(X)
points(m[1],m[2],pch=8,col="red",cex=2)

S1=solve(S)
contour(x,x,outer(x,x,function(x,y){(n-2)*
                                      apply(t(t(cbind(x,y))-m),1,function(x){t(x)%*%S1%*%x})<
                                      2*qf(0.95,2,n-2)}),col="red",add=TRUE)
bodyx=m[1]+c(-1,1)*sqrt(S[1,1]*2*qf (0.95,2,n-2) /(n-2))
bodyy=m[2]+c(-1,1)*sqrt(S[2,2]*2*qf (0.95,2,n-2) /(n-2))
polygon(x=bodyx[c(1,1,2,2,1)],y=bodyy[c(1,2,2,1,1)],border="blue")

Poznámka

---

• Konfidenčný elipsoid je výhodnejší z geometrického hľadiska, keďže pri stejnej hladine spoľahlivosti definuje menšiu plochu;
• Súbor konfidenčných intervalov pre jednotlivé zložky – t.j. kartezský súčin intervalov je na druhu stranu výhodnejší z praktického hľadiska a poskytuje omnoho lepšiu intuíciu pre pochopenie a interpretáciu. Pri stejnej hladine spoľahlivosti ale pokrýva výrazne vačšiu plochu.

---

Príklad I: ozdiel dvoch vektorov stredných hodnôt / stejná variančná-kovariančná matica

Hotellingove $T^2$ rozdelenie použijeme najprv pre test o rovnosti dvoch neznámych vektorových parametrov stredných hodnôt.

K dispozícii máme náhodný výber $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n} \sim N_{p}(\boldsymbol{\mu}_{1}, \Sigma)$ a na ňom nezávislý náhodný výber $\boldsymbol{Y}_{1}, \dots, \boldsymbol{Y}_{M} \sim N_{p}(\boldsymbol{\mu}_{2}, \Sigma)$. Parametre stredných hodnôt – t.j., p rozmerné pevné, ale neznáme vektory – sú obecne rôzne $\boldsymbol{\mu}_{1} \neq \boldsymbol{\mu}_{2}$, ale variačna-kovariančná matica je pre oba náhodné výbery stejná (hoci neznáma). Rovnaká je samozrejme aj dimenzia náhodných výberov.

Zaujíma nás test nulovej hypotézy $$H_{0}: \boldsymbol{\mu}_1 = \boldsymbol{\mu}_{2}$$ oproti obecnej alternatíve, že nulová hypotéza neplatí. Je zrejmé, že pre výberové priemery $\overline{X}_{1}$ a $\overline{X}_{2}$ platí, že

$$(\overline{X}_{1} - \overline{X}_{2}) \sim N_{p}\Bigg(\boldsymbol{\Delta}, \frac{n + m}{n m} \Sigma\Bigg),$$

a tiež

$$n\mathcal{S}_{1} + m\mathcal{S}_{2} \sim W_{p}(\Sigma, n + m - 2),$$ kde $\mathcal{S}_{1}$ a $\mathcal{S}_{2}$ sú empriciké odhady variančnej-kovariančnej matice $\Sigma$ spočítané samostatne na základe prvého a druhého náhodného výberu. Kritický obor pre daný test je preto definovaný následovne:

$$\frac{nm(n + m - p - 1)}{p(n + m)^2}(\overline{\boldsymbol{X}}_{1} - \overline{\boldsymbol{X}}_{2})^{\top} \mathcal{S}^{-1} (\overline{\boldsymbol{X}}_{1} - \overline{\boldsymbol{X}}_{2}) \geq F_{p, n + m - p - 1}(1 - \alpha).$$

---

Príklad II: Rozdiel dvoch vektorov stredných hodnôt / rôzne variančné-kovariančné matice

Analogicky ako v predchádzajúcom príklade, opäť uvažujme náhodný výber $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n} \sim N_{p}(\boldsymbol{\mu}_{1}, \Sigma_{1})$ a na ňom nezávislý druhý náhodný výber $\boldsymbol{Y}_{1}, \dots, \boldsymbol{Y}_{M} \sim N_{p}(\boldsymbol{\mu}_{2}, \Sigma_{2})$, s obecne rôznymi parametrami stredných hodnôt $\boldsymbol{\mu}_{1} \neq \boldsymbol{\mu}_{2}$, ale tentokrát aj s potenciálne rôznymi (neznámymi) variančnými-kovariančnými maticami. Opäť nás zaujíma test nulovej hypotézy o rovnosti vektorových parametrov stredných hodnôt, t.j., nulová hypotéza v tvare

$$H_{0}: \boldsymbol{\mu}_1 = \boldsymbol{\mu}_{2}$$ oproti obecnej alternatíve, že mnulová hypotéza neplati. Je zrejmé, že platí $$(\overline{X}_{1} - \overline{X}_{2}) \sim N_{p}\Bigg(\boldsymbol{\Delta}, \frac{\Sigma_{1}}{n} + \frac{\Sigma_{2}}{m}\Bigg),$$

a taktiež $$(\overline{X}_{1} - \overline{X}_{2})^\top \Big(\frac{\Sigma_{1}}{n} + \frac{\Sigma_{2}}{m}\Big)^{-1} (\overline{X}_{1} - \overline{X}_{2}) \sim \chi_{p}^{2}.$$

S využitím vyššie uvedeného môžeme zostrojiť testovú štatistiku a analogickym spôsobom ako v predchádzajom prípade získame kritický obor na základe ktorého rozhodneme o zamietnuti, resp. nezamietnuti nulovej hypotézy.

Samostatne

---

• Nainštalujte si v programe R knižnice ‘Hotelling’ a ‘DescTools’;
• Pre jednovýberové a dvojvýberové Hotellingové testy využijte príkaz HotellingsT2Test() z knižnice ‘DescTools’;
• Pre dvojvýberové testy môžete využiť aj príkaz hotelling.test() z knižnice ‘Hotelling’;
• Pomocou helpu k príkazom hotelling.stat() a hotelling.test()zistite konkrétnu implementáciu jednotlivých funkcii;

---

Ilustračný príklad:

library("Hotelling")

Využijeme datový súbor container.df z knižnice ‘Hotelling’. Data reprezentujú konkrétne merania koncentrácie rôznych chemických prvkov (podrobnejšie v helpe ?container.df).

data(container.df)
summary(container.df)

##        gp            Ti                Al               Fe         
##  Min.   :1.0   Min.   :0.02800   Min.   :0.6590   Min.   :0.02700  
##  1st Qu.:1.0   1st Qu.:0.03175   1st Qu.:0.7405   1st Qu.:0.02975  
##  Median :1.5   Median :0.03300   Median :0.7850   Median :0.06950  
##  Mean   :1.5   Mean   :0.03350   Mean   :0.8047   Mean   :0.07665  
##  3rd Qu.:2.0   3rd Qu.:0.03700   3rd Qu.:0.8730   3rd Qu.:0.12200  
##  Max.   :2.0   Max.   :0.03900   Max.   :0.9410   Max.   :0.16800  
##        Mn               Mg               Ca              Ba         
##  Min.   :0.0020   Min.   :0.0900   Min.   :5.511   Min.   :0.00900  
##  1st Qu.:0.0020   1st Qu.:0.1130   1st Qu.:6.422   1st Qu.:0.01200  
##  Median :0.0030   Median :0.1595   Median :6.989   Median :0.01550  
##  Mean   :0.0034   Mean   :0.1759   Mean   :7.001   Mean   :0.01615  
##  3rd Qu.:0.0050   3rd Qu.:0.2362   3rd Qu.:7.578   3rd Qu.:0.02000  
##  Max.   :0.0050   Max.   :0.2700   Max.   :7.780   Max.   :0.02300  
##        Sr               Zr         
##  Min.   :0.0130   Min.   :0.00900  
##  1st Qu.:0.0150   1st Qu.:0.01000  
##  Median :0.0160   Median :0.01150  
##  Mean   :0.0162   Mean   :0.01125  
##  3rd Qu.:0.0170   3rd Qu.:0.01200  
##  Max.   :0.0190   Max.   :0.01400

Príslušné empirické odhady neznámych vektorov stredných hodnôt aj variačných-kovariačných matíc sú:

apply(container.df[,2:10], 2, mean) ### mean vector estimate

##      Ti      Al      Fe      Mn      Mg      Ca      Ba      Sr      Zr 
## 0.03350 0.80470 0.07665 0.00340 0.17590 7.00055 0.01615 0.01620 0.01125

cov(container.df[,2:10])            ### variance-covariance estimate

##               Ti           Al            Fe            Mn            Mg
## Ti  1.194737e-05 2.289474e-04  1.527632e-04  4.421053e-06  1.975789e-04
## Al  2.289474e-04 7.775589e-03  3.868363e-03  1.112842e-04  5.884863e-03
## Fe  1.527632e-04 3.868363e-03  2.465608e-03  7.072632e-05  3.294963e-03
## Mn  4.421053e-06 1.112842e-04  7.072632e-05  2.147368e-06  9.583158e-05
## Mg  1.975789e-04 5.884863e-03  3.294963e-03  9.583158e-05  4.784832e-03
## Ca  6.710526e-06 2.654805e-03  2.015711e-04  1.297895e-05 -4.767895e-05
## Ba  1.355263e-05 3.768895e-04  2.209500e-04  6.410526e-06  3.091211e-04
## Sr -2.105263e-07 1.484211e-06 -1.876842e-05 -6.631579e-07 -1.371579e-05
## Zr  1.921053e-06 6.907895e-05  2.451316e-05  6.842105e-07  4.423684e-05
##               Ca            Ba            Sr            Zr
## Ti  6.710526e-06  1.355263e-05 -2.105263e-07  1.921053e-06
## Al  2.654805e-03  3.768895e-04  1.484211e-06  6.907895e-05
## Fe  2.015711e-04  2.209500e-04 -1.876842e-05  2.451316e-05
## Mn  1.297895e-05  6.410526e-06 -6.631579e-07  6.842105e-07
## Mg -4.767895e-05  3.091211e-04 -1.371579e-05  4.423684e-05
## Ca  4.574089e-01  5.564921e-04 -3.695368e-04 -1.588158e-05
## Ba  5.564921e-04  2.139737e-05 -1.821053e-06  2.644737e-06
## Sr -3.695368e-04 -1.821053e-06  2.484211e-06  1.210526e-06
## Zr -1.588158e-05  2.644737e-06  1.210526e-06  1.671053e-06

Zaujíma nás test nulovej hypotézy, či (neznáme) stredné koncentrácie v prvom kontajneri (gp == 1) sú rovnaké, ako v druhom kontajneri (gp == 2).

split.data = split(container.df[,-1],container.df$gp)
x <- split.data[[1]]
y <- split.data[[2]]
statistic <- hotelling.stat(x, y)

Samotný test:

hotellingTest <- hotelling.test(.~gp, data = container.df)
hotellingTest

## Test stat:  126.11 
## Numerator df:  9 
## Denominator df:  10 
## P-value:  4.233e-09

Otázky

---

• Ako by ste na základe prezentovaných čísel interpretovali výsledok testu?
• Ako by ste výsledky vizuálne reprezentovali pomocou nejakého obrázku?
• Pre jednotlivé zložky aplikujte jednorozmerý dvojvýberový $t$-test(s). Ako sa líšia závery?
• Ako by ste zostrojili simultánny interval spoľahlivosti pre jednotlivé zložky neznámych vektorových parametrov stredných hodnôt?

---

Príslušné simultánne intervaly spoľahlivosti pre (všetky) lineárne kombinácie jednotlivých zložiek $\boldsymbol{\mu} \in \mathbb{R}^{p}$ (lineárne kombinácie definované parametrami $\boldsymbol{a}^\top$ získame pomocou vzťahu

$$P\Big(\forall \boldsymbol{a} \in \mathbb{R}^{p};~ \boldsymbol{a}^\top \boldsymbol{\mu} \in \big( \boldsymbol{a}^\top\overline{\boldsymbol{X}} - \sqrt{K_{\alpha} \boldsymbol{a}^\top \mathcal{S} \boldsymbol{a}}, \boldsymbol{a}^\top\overline{\boldsymbol{X}} + \sqrt{K_{\alpha} \boldsymbol{a}^\top \mathcal{S} \boldsymbol{a}} \big) \big)\Big) = 1 - \alpha,$$ kde $K_{\alpha}$ je príslušný transformácia kvantilu Fisherovho $F$ rozdelenia v tvare $K_{\alpha} = \frac{p}{n - p} F_{p, n - p}(1 - \alpha)$ a $\mathcal{S}$ výberová variančná-kovariančná matica.

Fquant <- (9/(20 - 9)) * qf(1 - 0.05, 9, 11)
meanVec <- apply(container.df[,2:10], 2, mean)
LB <- diag(9) %*% meanVec - sqrt(Fquant * diag(diag(9) %*% cov(container.df[,2:10])  %*% diag(9)))
UB <- diag(9) %*% meanVec + sqrt(Fquant * diag(diag(9) %*% cov(container.df[,2:10])  %*% diag(9)))

confInt <- data.frame(LB, UB, meanVec,  row.names = names(container.df)[2:10])
names(confInt) <- c("Lower Boundary", "Upper Boundary", "Mean Estimate")

confInt

##    Lower Boundary Upper Boundary Mean Estimate
## Ti   0.0281791989    0.038820801       0.03350
## Al   0.6689600906    0.940439909       0.80470
## Fe   0.0002131292    0.153086871       0.07665
## Mn   0.0011442333    0.005655767       0.00340
## Mg   0.0694184851    0.282381515       0.17590
## Ca   5.9594482060    8.041651794       7.00055
## Ba   0.0090293265    0.023270674       0.01615
## Sr   0.0137737525    0.018626247       0.01620
## Zr   0.0092600784    0.013239922       0.01125

3. Wilkovo Lambda rozdelnie

Toto rozdelenie je odvodené z dvoch nezávislých náhodných matíc s Wishartovým rozdelenim. Vpodstate sa jedná o mnohorozmerné zobecnenie jednorozmerného Fisherovho F rozdelenia a používa sa hlavne k inferencii dvoch variačných-kovariačných matíc (analogicky, ako sa jednorozmerné Fisherovo F rozdelenie používa k inferencii o dvoch parametroch rozptylu).

Pre dve náhodné a vzájomné nezávislé náhodné matice $$\mathbb{A} \sim W_{p}(\mathbb{I}, n) \quad \textrm{a} \quad \mathbb{B} \sim W_{p}(\mathbb{I}, m)$$

môžeme definovať náhodnú veličinu $\frac{|\mathbb{A}|}{|\mathbb{A} + \mathbb{B}|}$ ktorej rozdelenie je Wilkovo Lambda, t.j., $\Lambda(p, n, m)$. Toto rozdelenie sa tiež používa v súvislosti s testom pomerom vierohodnosti. Analogickym spôsobom, ako funguje metóda analýzy rozptylu (ANOVA) na základe Fisherovho F rozdelenia, v mnohorozmernom prípade používame mnohorozmerné zobecnenie (MANOVA), ktorá vychádza z Wilkovho Lambda rozdelenia. V programe R príkaz manova().

Ilustračný príklad

---

• K ilustrácii Wilkovho Lambda rozdelenia využijeme priamo definíciu náhodnej veličiny s týmto rozdelením. Generovať budeme náhodné matice z Wishartovho rozdelenia a vhodnou transformáciou získavame požadované Wilkovo Lambda rozdelenie.
  
  

  set.seed(1234)
  I <- cbind(c(1,0,0), c(0,1,0), c(0,0,1))
  A <- rWishart(100, df = 10, I)
  B <- rWishart(100, df = 4, I)
  
  WilksLambda <- apply(A, 3, det)/apply(A + B, 3, det)
  
  hist(WilksLambda, col = "lightblue", freq = F, xlab = "Wilk's Lambda", ylab = "Density", ylim = c(0, 4.5))
  lines(density(WilksLambda), col = "red", lwd = 2)

  

---

---

---