Tradiční míry diverzity a citlivost mocninných entropií

cs17

Původní práce

Tradiční míry diverzity a citlivost mocninných entropií Martin Horáček1,2 , Jana Zvárová1,2 1 2

Centrum biomedicínské informatiky, Ústav informatiky AV ČR, v.v.i., Praha, Česká republika Ústav hygieny a epidemiologie 1. lékařské fakulty Univerzity Karlovy v Praze, Česká republika

Souhrn Cíle: Zabývali jsme se tradičními mírami diverzity a jejich odhady. Zkoumali jsme způsoby, jak porovnávat citlivost různých měr diverzity na změny. Metody: Navrhli jsme nový typ odhadu pro míry diverzity. V simulační studii jsme srovnali jeho vlastnosti, zejména vychýlení a rozptyl, se třemi zavedenými odhady. Zavedli jsme funkci, kterou jsme nazvali citlivost na změny míry diverzity H a studovali jsme její vlastnosti. Výsledky: Navržený odhad je v mnoha situacích srovnatelný či lepší než zavedené typy odhadů. Citlivost na změny má srozumitelnou interpretaci a jednoduchý tvar. Závěry: Citlivost míry diverzity na změny může být využita ke srovnání chování různých měr diverzity a k vybrání těch měr, které jsou nejvhodnější pro daný problém.

Mgr. Martin Horáček Klíčová slova diverzita, entropie, odhady diverzity, citlivost

Kontakt: Mgr. Martin Horáček Centrum biomedicínské informatiky, Ústav informatiky AV ČR, v.v.i. Adresa: Pod Vodárenskou věží 2, 182 07 Praha E–mail: [email protected]

1

Úvod

EJBI 2011; 7(1):17–21 zasláno: 29. září 2011 přijato: 24. října 2011 publikováno: 20. listopadu 2011

která splňuje • je nezáporná,

V tomto článku se zabýváme funkcemi, jejichž cílem • je symetrická vzhledem k permutacím, je vhodně a jednoduše zachytit míru diverzity zvolené po• je minimální, když pro nějaký index i platí pi = 1 pulace. Může nás zajímat například genetická diverzita (v populaci se vyskytuje jen jeden znak) diverzita alel zvoleného genu, druhová diverzita ve zvolené lokaci, ale také jazyková či ekonomická diverzita. Zejména • je maximální, když pi ≡ 1/r (znaky jsou rovnoměrně se budeme zabývat situací, kdy míra diverzity populace zastoupeny), závisí pouze na pravděpodobnostech pi , že náhodné vybraný člen populace má právě znak i z r možných na• když se zvetší jedna pravděpodobnost pi na úkor vzájem se vylučujících znaků. Při splnění několika dalších jiné menší pravděpodobnosti pj , hodnota H(p) se předpokladů (popsaných v následujícím odstavci), které nezvětší. lze intuitivně očekávat od funkce, která měří populační diverzitu, se tyto funkce často souhrně nazývají tradiční Pro ověřování těchto požadavků může být užitečné si všimnout, že když je funkce H : ∆r → R+ nezáporná míry diverzity. a Schur-konkávní, splňuje všech pět požadavků (viz [3]). Formálně zapsáno, tradiční míra diverzity je reálná Existuje několik často používaných tradičních měr difunkce H definovaná na verzity. Nejznámější z nich uvádíme v následující sekci. ( ) Většinu z nich lze zahrnout pod takzvané f -entropie, nar X vržené Zvárovou [1], nebo jsou s touto rodinou zobecněr ∆ = p = (p1 , . . . , pr ) : pi = 1, pi ≥ 0 ∀ i , ných entropií úzce spjaty. Tyto f-entropie, jejich vlastnosti i=1 c

2011 EuroMISE s.r.o.

EJBI – Ročník 7 (2011), číslo 1

cs18

Horáček, Zvárová – Tradiční míry diverzity a citlivost mocninných entropií

a způsob, jakým je lze využít k měření diverzity, je stu- a Rényiho entropie řádu α dován a popsán například v pracích Zvárová, Vajda [2] ! r X a Horáček [3]. −1 α HR,α (p) = (1 − α) ln pi ,

pro α > 0, α 6= 1.

i=1

2

Tradiční míry diverzity a jejich odhady

Limitním dodefinováním všech těchto indexů pro α = 1 dostaneme Shannonovu entropii. Tyto a další zobecňující parametrické jsou užitečné různými způsoby. Lze z nich vybírat vhodný index pro měření zvoleného druhu diverzity, lze je však také použít ke zlepšení procedur založených na tradiční Shannonově entropii. Nahradíme-li tuto entropii vhodným parametrických indexem, můžeme volbou parametru docílit zlepšení dotyčné procedury. Tento postup je použit například v práci Andrade, Wang [5].

V této sekci zavedeme nejznámější tradiční míry diverzity, jako například Simpsonův index, Shannonovu entropii, Rényiho entropii řádu α, Hillův index a další. Uvedeme problematiku odhadů tradičních měr diverzity a odvodíme nový druh odhadu. Část sekcí 2.2 a 3.1 byla publikována ve sborníku k 7. Letní škole výpočetní biologie 2.2 [4].

Odhady tradičních měr diverzity

Nechť p = {p1 , . . . , pr } ∈ ∆r je vektor neznámých pravděpodobností pi , že jedinec náhodně vybraný z po2.1 Příklady tradičních měr diverzity pulace má znak Ai z r možných znaků. V této situaci se odhad míry diverzity H(p) často provádí na základě repi , . . . , pˆr ), Mezi nejčastěji zmiňované a používané tradiční míry lativních frekvencí pˆn = (X1 /n, . . . , Xr /n) = (ˆ napozorovaných v náhodné výběru s vracením o rozsahu diverzity patří počet znaků (napříkal alel nebo druhů) n. Rozdělení vektoru X = (X1 , . . . , Xr ) pak je multir nomické M(n, p). Míru diverzity lze odhadnout různými X I(0,1] (pi ) − 1, H0 (p) = způsoby a vlastnosti zvoleného odhadu, zejména jeho vyi=1 chýlení a rozptyl, případně jejich střední čtvercová chyba, závisí jak na použité míře diverzity, tak na dané populaci. kde I značí funkci identity, Simpsonův index Zřejmě nejčastěji používaný je takzvaný plug-in odr had. Jeho princip spočívá v jednoduchém nahrazení neX p2i H2 (p) = 1 − známých pravděpodobností pi napozorovanými frekveni=1 cemi pˆi a jejich dosazením H(ˆ p). Ačkoli jsou frekvence pˆi nestranným odhadem pravděpodobností pi , plud-in odhad a Shannonova entropie není obecně nestranný. r X V některých případech lze jeho vychýlení opravit. Je to H1 (p) = − pi ln pi . možné například u odhadu Simpsonova indexu, pro jehož i=1 střední hodnotu platí Tyto tři míry jsou zobecněny parametrickou rodinou mocninných entropií ! r X −1 α Hα (p) = (α − 1) 1− pi , pro α > 0, α 6= 1, i=1

limitně dodefinovanou pro α = 0 (pak dostaneme počet znaků) a pro α = 1 (dostaneme Shannonovu entropii). Když je α = 2, dostaneme Simpsonův index.

EH2 (ˆ pn )

=

1 − n−2

r X

EXi2

i=1

=

1 − n−2

r X

varXi + (EXi )2

i=1

=

1 − n−2

r X

npi (1 − pi ) + n2 p2i

i=1 −1

= 1−n H2 (p). Mezi další často používané indexy patří γ-entropická Nevychýlený odhad Simsonova indexu je tedy funkce " !γ # ˆ 2 (ˆ r H pn ) = n(n − 1)−1 H2 (ˆ pn ). X 1/γ −1 HA,γ (p) = (1 − γ) 1− pi , γ > 0, γ 6= 1, U jiných měr diverzity je ale většinou obtížné či dokonce i=1 nemožné nalézt vhodnou korekci plug-in odhadu. Lze naHillův index příklad ukázat, že pro Shannonovu entropii nevychýlený 1 ! odhad neexistuje (Blyth [6]). Několik autorů se zabývalo 1−α r X tímto problémem a navrhli sofistikovanější typy odhadů. α , pro α > 0, α 6= 1 HH,α (p) = pi Uvádíme zde odhad navržený Bonachelou a kol. [7] zvaný i=1 EJBI – Ročník 7 (2011), číslo 1

c


cs19


Obrázek 1: Výběrový průměr a výběrový rozptyl odhadů mocninné entropie H3/2 . n i X balanced odhad. Navrhujeme modifikaci tohoto odhadu, + P (Xi = j)ζ 2 (j) w(pi )dpi = 0 ktérá bere v úvahu rozložení hodnot pi ve vektoru p. j=0 Předpokládejme, že zvolenou míru diverzity lze zapsat ve tvaru což se dá zjednodušit do tvaru ! r X Z 1h i H(p) = F h(pi ) , ζ(k)P (Xi = k) − h(pi )P (Xi = k) w(pi )dpi = 0. i=1 0

kde F a h jsou libovolné reálné spojité funkce. Tento tvar zahrnuje mimo jiné všechny indexy uvedené v předchozí Protože části kromě počtu znaků. Bonachela a kol. navrhli svůj odhad ve tvaru ! máme r X ˆ H(X) = F ζ(Xi ) , i=1

P (Xi = k) =

ζ(k)

n k p (1 − pi )n−k , k i

R1 k n−k h(pi )w(pi )(n dpi k )pi (1−pi ) . = 0R 1 n k w(pi )( k )pi (1−pi )n−k dpi 0

(2)

kde funkce ζ je zvolena tak, aby se minimalizovala hodBonachela a kol. [7] odvodili podobu balanced odhadu pro nota výrazu mocninné entropie s váhovou funkcí rovnou jedné na celém intervalu [0, 1]. Φ2ζ (pi ) = [E(ζ(Xi ) − h(pi ))]2 + var(ζ(Xi )) Nemáme-li žádnou apriorní znalost o hodnotách slos váhovou funkcí w(pi ) umožnující při minimalizaci měnit žek vektoru p, je přirozené pro tento vektor nepreferovat důraz kladený na různé hodnoty pi z intervalu [0, 1]. Po- žádný bod z ∆r . Na vektor p tak můžeme nahlížet jako kud zanedbáme vliv funkce F a korelací, můžeme zároveň na realizaci náhodného vektoru Y = (Y1 , . . . , Yr ), která redukovat rozptyl a druhou mocninu vychýlení odhadu. má rovnoměrné rozdělení na ∆r . Váhu w(pi ) je vhodné Vážená průměrná chyba je pak dána vztahem zvolit proporcionálně k očekávaným hodnotám složek pi , v tomto případě tedy jako marginální hustotu náhodné R ¯ 2ζ (pi ) = 1 Φ2 (pi )w(pi )dpi . Φ (1) veličiny Yi . ζ 0 Tu lze spočítat jako Nutná podmínka pro minimalitu chyby je nulová hodnota Z 1−y1 Z 1−y1 −...−yr−2 derivací f (y1 ) ∝ ... dyr−1 . . . dy2 δ ¯2 0 0 Φ (pi ) = 0, k ∈ {0, . . . , n} . (1 − y1 )r−2 δζ(k) ζ = , (r − 2)! Zvolíme proto ζ tak, aby což je až na konstantu hustota Beta rozdělení B(1, r − 1). Z 1h n X δ 2 Zvolili jsme proto váhovou funkci jako w(pi ) = (1 − h (pi ) − 2h(pi ) P (Xi = j)ζ(j) + r−2 δζ(k) 0 p ) a odvodili odpovídající funkci ζ. Odhad založený i j=0 c



cs20


na této funkci jsme nazvali β-odhad. Popíšeme zde odvo- pií H1/2 a H5/2 v populacích s p1 = (13, 9, 2, 2, 1)/27 zení β-odhadu pro Shannonovu entropii, tedy pro situaci, a p2 = (12, 8, 5, 3, 3, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1)/46, na zákdy h(pi ) = −pi ln pi a F (x) = x. kladě 1000 výběrů o velikosti n = 50. Symboly Γ, B a Ψ označují gama, beta and digamma funkce. Nejprve jsme nahradili h(pi ) a w(pi ) odpovídajícími výrazy a spočítali integrál ve jmenovateli výrazu zeta1. R1 X h(pi )p i (1−pi )n−Xi +r−2 dpi . ζ(Xi ) = 0 B(Xii +1,n−Xi +r−1)

Tabulka 1: Absolute sample bias and sample variance of estimates.

p1

Pro parciální derivace beta funkce platí δ B(x, y) = B(x, y)[Ψ(x) − Ψ(x + y)], δx a čitatel lze tedy vyjádřit jako Z 1 n−Xi +r−2 i dpi h(pi )pX i (1 − pi ) 0

Z =

−

n−Xi +r−2 i pi ln(pi )pX dpi i (1 − pi ) 1

Z

= =

− lim

α→0

pα i

0

− 1 Xi +1 pi (1 − pi )n−Xi +r−2 dpi α

1 lim [B(Xi + 2, n − Xi + r − 1) − α→0 α − B(Xi + α + 2, n − Xi + r − 1)] B(Xi + 2, n − Xi + r − 1)[Ψ(n + r + 1) − − Ψ(Xi + 2)].

Funkce ζ tedy splňuje ζ(Xi )

3

H5/2 |bias| var 0,0094 0,0010 0,0169 0,0008 0,0009 0,0006 0,0056 0,0010 0,0075 0.0002 0,1673 0,0002 0,0038 0,0001 0,0052 0,0002

Citlivost na změny

1

0

=

p2

plug-in balanced beta shrink plug-in balanced beta shrink

H1/2 |bias| var 0,0843 0,0339 0,0059 0,0143 0,0052 0,0135 0,0354 0,0211 0,7476 0,1486 0,2867 0,0281 0,1288 0,0218 0,0651 0,0817

= =

Xi + 1 [Ψ(n + r + 1) − Ψ(Xi + 2)] n+r n+r Xi + 1 X 1 n+r

k=Xi +2

k

Vlastnosti indexů používaných k měření diverzity se více či méně liší. Často by nás přitom zajímalo, jak by daná míra diverzity reagovala na změny frekvencí jednotlivých znaků v populaci. Tímto problémem se zabývalo několik autorů, zejména Boyle a kol. [10], kteří se zaměřili zejména na empirickou stránku věci, a Izsak [11], který navrhnul míru citlivosti na teoretickém základě. Na základě návrhu I. Vajdy zde navrhujeme míru citlivosti měr diverzity na změny ve frekvenci znaků, kterou je snadné spočítat a má intuitivní interpretaci. Míru citlivosti míry diverzity H na změny ve frekvenci j-tého znaku definujeme jako H(pj, ) − H(p) →0 H(p)

SH (p|j) = lim kde

a β-odhad Shannonovy entropie má tvar ˆ 1 (X) H

=

r P i=1

Xi +1 n+r

n+r P k=Xi +2

1 k.

pj, =

(p1 , . . . , pj−1 , pj + pj , pj+1 , . . . , pr ) . 1 + pj

Citlivost je tedy definována jako limita výrazu

relativní změna H β-odhad mocninných entropií, pro které platí F (x) = x relativní změna pj a h(pi ) = (α − 1)−1 (pi − pα i ), lze odvodit podobným způsobem. U mocninných entropií nabývá tvaru a odráží poměr relativních změn H(p) a pj , když se pj změní o malou hodnotu. r P B(n+r,α) ˆ α (X) = (α − 1)−1 1 − H B(Xi +1,α) . i=1

3.1

Citlivost mocninných entropií

Na obrázku 1 můžeme vidět srovnání β-odhadu, Bonachelova balanced odhadu, plug-in odhadu a JamesovaCitlivost mocninných entropií byla spočítána v práci Steinova shrinkage odhadu [8] v populaci se šesti znaky, Horáček [3]. Pokud platí pi > 0 pro všechna i ∈ {1, . . . , r}, rozdělenými v poměru 24 : 11 : 9 : 3 : 2 : 1. Zob- pro tuto citlivost platí vztah razeny jsou výběrové průměry a výběrové rozptyly odr P hadů spočítané z 300 pokusů. Obrázky byly vytvořeny piα−1 (pi − δij ) v programu R [9]. Další srovnání vývěrových rozptylů SHα (p|j) = α i=1 , r P a absolutních hodnot výběrového vychýlení je v tabulce α 1 − p i 1. Tentokrát jsme srovnávali odhady mocninných entroi=1 EJBI – Ročník 7 (2011), číslo 1

c



pokud α 6= 1, a

Poděkování r P

SH1 (p|j) =

cs21

(pi − δij ) ln pi

i=1 r P

. pi ln pi

Tato práce byla podporována projektem 1M06014 Ministerstva školství, mládeže a sportu ČR a projektem SVV-2011-262514 Univerzity Karlovy v Praze.

i=1

Srovnání citlivostí mocninných entropií jsme provedli na populaci s p = (24/50, 11/50, 9/50, 3/50, 2/50, 1/50) a je zobrazeno na obrázku 2. Můžeme vidět, že s klesajícím α jsou mocninné entropie více citlivé na změny ve znacích, které jsou v populaci zastoupeny vzácně. Podíváme-li se například na Shannonovu entropii, řekněme desetiprocentní nárůst v pravděpodobnosti p5 = 1/25 by vedl ke zhruba desetiprocentnímu nárůstu H1 (p), zatímco desetiprocentní nárůst v pravděpodobnosti p1 = 24/50 by vyústil v asi tříprocentní snížení hodnoty H1 (p). Malé změny v pravděpodobnosti p2 = 11/50 by pak na hodnotu H1 (p) měly jen zanedbatelný vliv.

Literatura [1] Zvárová J.: On Measures of Statistical Dependence. Časopis pro pěstování matematiky 1974; 99: 15–29 [2] Zvárová J., Vajda I.: On Genetic Information, Diversity and Distance. Methods of Inform. in Medicine 2006; 2: 173–179 [3] Horáček, M.: Measures of biodiversity and their applications. Master thesis, Charles university, Prague, supervisor J. Zvárová 2009 [4] Horáček, M., Zvárová J.: Traditional Measures of Diversity, Their Estimates and Sensitivity to Changes. Proceedings of the 7th Summer School on Computational Biology. 2011; 73– 81. [5] Andrade, M. de, Wang, X.: Entropy Based Genetic Association Tests and Gene-Gene Interaction Tests. Statistical Applications in Genetics and Molecular Biology 2011; 10: Iss. 1, Article 38 [6] Blyth, C. R.: Note on estimating information. Annals of Math. Stat. 1959; 30: 71–79 [7] Bonachela, J. A., Hinrichsen, H., Muñoz, M. A.: Entropy estimates of small data sets. J. of Phys. A: Math. and Theor. 2008; 41: 1–9 [8] Hausser, J., Strimmer, K.: Entropy Inference and the JamesStein Estimator, with Application to Nonlinear Gene Association Networks. Journal of Machine Learning Research 2009; 10: 1469-1484. [9] R Development Core Team: R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org 2011

Obrázek 2: Srovnání citlivostí na změny u mocninných entropií.

[10] Boyle, T. P., Smillie, G. M., Anderson, J. C., and Beeson, D. R.: A sensitivity analysis of nine diversity and seven similarity indices. Research Journal Water Pollution Control Federation 1990; 62: 749–762 [11] Izsak, J.: Sensitivity Profiles of Diversity Indices. Biom. J. 1996; 38: 921–930

c



Tradiční míry diverzity a citlivost mocninných entropií

Recommend Documents