UNIVERZA V LJUBLJANI ˇ BIOTEHNISKA FAKULTETA ODDELEK ZA ZOOTEHNIKO
Gregor GORJANC1
(Bio)informatika in Bayesovska statistika
Seminarska naloga pri metodoloˇskem predmetu “Bioinformatika” Podpiplomski ˇstudij genetike
Ljubljana, 2006 1
[email protected]
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
2
KAZALO VSEBINE str. 1
UVOD
3
2
ORIS BAYESOVSKE STATISTIKE
3
3
UPORABA BAYESOVSKE STATISTIKE V “BIOLOGIJI”
5
3.1
MATERIAL
5
3.2
METODE
5
3.2.1
Orodja
5
3.2.2
Prikaz uporabe
6
3.2.3
Iskanje
8
3.3
REZULTATI
8
4
IZBOR APRIORNE PORAZDELITVE
10
5
SKLEPI
16
6
VIRI
16
A
PRILOGA
17
entrezSearch
17
entrezCountByYear
23
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
1
3
UVOD
Od odkritja zapisa molekule DNA do danes se je koliˇcina informacij pri raziskavah v biologiji moˇcno poveˇcala. Tako kot so prepleteni pojavi, ki jih preuˇcujejo raziskovalci, so prepletena tudi podroˇcja raziskovanja. Prepletenost raziskav se v zadnjih letih vse bolj stopnjuje. Slednje je posledica do sedaj zbranega znanja in v veliki meri dostop do novih tehnologij, ki omogoˇcajo vse laˇzje pridobivanje podatkov v velikih koliˇcinah in na razliˇcnih podroˇcjih/nivojih. Na drugi strani se s tem poveˇcuje tudi kompleksnost analiz zbranih ˇ so bila vˇcasih posamezna podroˇcja raziskav “dobro” loˇcena, se te meje vse bolj briˇsejo. S podatkov. Ce pogledom na analizo podatkov lahko posamezna podroˇcja zdruˇzimo pod izrazom bioinformatika. Pri preuˇcevanju naravnih zakonitosti/pojavov je potrebno podatke tudi analizirati. Za statisitiˇcno analizo podatkov je na razpolago veˇc pristopov k t.i. statistiˇcnem sklepanju. Najbolj pogosta je uporaba klasiˇcne oz. frekvencistiˇcne statistike, ki je, kot ˇze ime pove, dobro poznana in uveljavljena. Na drugi strani lahko izberemo Bayesovski pristop, ki jo zagovorniki predstavljajo kot bolj naravni naˇcin pristopa k statistiˇcnem sklepanju. Namen tega prispevka je najprej orisati Bayesovsko statistiko in podati izhodiˇsˇce za osrednji temi in sicer: poizkus primerjalne ocene uporabe Bayesovske statistike v na podroˇcju bioloˇskih ved po letih s
pomoˇcjo novo razvitih orodij in preuˇcitev ter prikaz uporabe teorije informatike pri izbiri apriorne porazdelitve za analizo podatkov
z Bayesovskim pristopom.
2
ORIS BAYESOVSKE STATISTIKE
Raziskovalci ob prouˇcevanju naravnih zakonitosti spremljajo pojave/dogodke. Pojav lahko matematiˇcno zapiˇsemo kot skupno porazdelitev p(θ, y), kjer θ predstavlja vektor neznanih parametrov, ki so tema prouˇcevanja, in y vektor zbranih podatkov. Ker te skupne porazdelitve pravzaprav nikoli ne poznamo, skuˇsamo pri raziskavah oceniti neznane parametre na osnovi zbranih podatkov in znanja, ki je na voljo. Slednje matematiˇcno zapiˇsemo kot p(θ|y). Bayesovska statistika za izraˇcun te koliˇcine uporabi Bayesov izrek, ki ga v obliki za kontinuirane spremenljivke lahko zapiˇsemo kot:
p(θ|y) =
p(y|θ) × p(θ) , p(y)
(1)
kjer p(θ|y) predstavlja posteriorno porazdelitev neznanih parametrov, ki nam nudi vse informacije za sklepanje o parametrih glede na podatke in ostale predpostavke. Zbrani podatki za analizo so zastopani v p(y|θ). Ta ˇclen predstavlja model vzorˇcenja oz. nastajanja podatkov in je kot funkcija neznanih parametrov obiˇcajno zapisan tudi kot L(θ|y). V tej obliki ga imenujemo verjetje (ang. likelihood) in ˇ veˇc, princip verjetja predstavlja v statistiki pristop, ki predstavlja osrednji del Bayesovega izreka. Se ga lahko vmestimo med frekvencistiˇcno in Bayesovsko ˇsolo statistike. Metoda najveˇcjega verjetja (ang. maximum likelihood - ML) je najpogosteje uporabljena metoda za ocenjevanje parametrov. Predhodno znanje o neznanih parametrih je v Bayesov izrek vkljuˇceno preko apriorne porazdelitve, p(θ). Delitelj, p(y), predstavlja robno porazdelitev podatkov. Ta del ni odvisen od neznanih parametrov, saj predstavlja
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
4
robno porazdelitev podatkov preko vseh neznanih parametrov (Θ), tudi tistih, ki jih nismo vkljuˇcili v analizo.
p(y) =
Z
Θ
p(y|θ) × p(θ)dΘ
(2)
Tako lahko Bayesov izrek zapiˇsemo tudi kot:
p(θ|y) ∝ p(y|θ) × p(θ)
(3)
Predstavljeni izrek nudi celovit pristop k statistiˇcnemu sklepanju, a se kljub temu Bayesovska statistika v praktiˇcne namene v preteklosti ni veliko uporabljala ali pa je bila omejena na enostavne primere. Razlog za to je v zahtevnosti izraˇcuna zapisa [3]. Poleg tega skupna posteriorna porazdelitev vseh ˇ neznanih parametrov ni primerna za sklepanje. Stevilo parametrov je lahko zelo veliko s ˇcimer se poveˇcuje razseˇznost te porazdelitve. Predstavljivo informacijo lahko pridobimo iz robne porazdelitve posameznega parametra (θi ). To porazdelitev lahko teoretiˇcno izvrednotimo s pomoˇcjo veˇckratnega integrala skupne posteriorne porazdelitve preko preostalih parametrov (θ i− ):
p(θi |θ i− , y) =
Z
θi−
p(y|θ) × p(θ)dθ i− .
(4)
Obe naˇsteti operaciji (izvrednotenje izraza [3] in veˇckratni integral v [4]) sta v praktiˇcnih pogojih velikokrat analitiˇcno neizvedljivi. S spoznanjem, da lahko predstavljene probleme premostimo na stohatsiˇcen naˇcin z uporabo metode Monte Carlo z Markovskimi verigami (MCMC), je Bayesovska statistika postala moˇcno orodje za praktiˇcne in raziskovalne namene. Metode MCMC so zaradi svoje sploˇsnosti primerne za ˇsirok nabor problemov. Najbolj znan in najpogosteje uporabljen algoritem je nedvoumno algoritem Metropolis-Hastings (Metropolis et al., 1953; Hastings, 1970). Algoritem Gibbs (Geman and Geman, 1984; Gelfand et al., 1990; Casella and George, 1992) je poseben primer tega algoritma in zaradi enostavnosti tudi pogosto uporabljen. Green (1995) je s t.i. algoritmom “Reversible Jump MCMC” ˇse nadalje posploˇsil algoritem Metropolis-Hastings in s tem odprl moˇznost analize izbire modela za statistiˇcno sklepanje. Zaradi zmoˇznosti analize zelo razliˇcnih problemov/modelov, ki so lahko kompleksni in imajo veliko ˇstevilo (tudi slabo definiranih) parametrov, je Bayesovska statistika z metodami MCMC veˇc kot primerno orodje za analizo podatkov iz biologije. Za nadaljne branje o Bayesovski statistiki priporoˇcamo prispevek Bernardo (2003) in literaturo v njem, o uporabi Bayesovske statike in metod MCMC pa osrednji deli s tega podroˇcja: Gilks et al. (1998) in Gelman et al. (2004). Kljub prepodoru Bayesovsko statistiko ˇse vedno nekoliko omejuje uporaba apriorne porazdelitve. Ravno v tem je jedro nestrinjanja med Bayesovskimi in frekvencistiˇcnimi statistiki. Apriorna porazdelitev opiˇse (naˇse) predhodno znanje o neznanih parametrih ˇse predno zberemo podatke. Slednje je logiˇcno gledano nesmiselno ali nemogoˇce in s staliˇsˇca objektivnosti raziskovanja velikokrat sporno. Z izborom apriorne porazdelitve lahko tudi vplivamo na rezultat statistiˇcne analize. Da bi se temu izognili, smo pri izboru te porazdelitve kar se da pazljivi. V kolikor obstaja dovolj zanesljivega znanja s katerim se strinja veˇcina je
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
5
to znanje vsekakor koristno vkljuˇciti v analizo. Veˇcji problem nastane, ko takˇsnega znanja ali pogosteje soglasja ni. Takrat skuˇsamo v analizo vkljuˇciti apriorno porazdelitev, ki je ˇcimmanj informativna. Z intuitivnega vidika se kot idealna ponuja enakomerna (ang. uniform) porazdelitev, ki pravi, da so vse vrednosti parametra enako verjetne. Kasneje bomo pokazali, da temu ni vedno tako, in predstavili pristope za izbor “neinformativne” apriorne porazdelitve, ki temeljijo tudi na teoriji informatike. 3
UPORABA BAYESOVSKE STATISTIKE V “BIOLOGIJI”
3.1 MATERIAL Primerjalno oceno uporabe Bayesovske statistike v “biologiji”po letih smo opravili na podlagi zastopanosti izbranih kljuˇcnih besed v objavljenih delih. Za vir podatkov smo izbrali podatkovno bazo objav (ˇclankov) PubMed pri National Center for Biotechnology Information - (NCBI., 2005). Ta vir smo izbrali, ker po ˇ nam znanih podatkih zajema objave s podroˇcja bioloˇskih ved v najveˇcjem obsegu. Stevilo objav v tej bazi znaˇsa veˇc kot 15 milijonov in sicer iz veˇc kot 30 000 revij. Kot dodatni in zelo pomemben kriterij pri izboru tega vira sta bili ˇse dejstvi, da iskalni portal “Entrez” omogoˇca enostaven in prost dostop do podatkov ter nenazadnje tudi “programirano” iskanje. 3.2 METODE 3.2.1
Orodja
Poizvedovanje po podatkih bi lahko opravili s pomoˇcjo iskalnika na spletnih straneh portala Entrez, a bi bilo takˇsno delo zamudno in hkrati teˇzavno za ponovitev ali posodobitev. Za bolj avtomatizirano oz. programirano iskanje po vseh bazah pri NCBI lahko uporabimo t.i. Entrez programska orodja eUtils (Sayers and Wheeler, 2005). Ta orodja nam omogoˇcajo dostop do zapisov v podatkovnih bazah pri NCBI preko spleta. Dostop lahko izvedemo kar s spletnim brskalnikom in ustrezno sestavo URL-ja, ki zajema osnovni naslov, orodje in dodatne argumente, s katerimi izvedemo iskanje. Uporabo orodja EInfo nam prikazuje naslednji primer, ki nam izpiˇse osnovne podatke o podatkovni bazi (argument db) PubMed: http://eutils.ncbi.nlm.nih.gov/entrez/eutils/einfo.fcgi?db=PubMed. Prednost pri uporabi takˇsnega naˇcina iskanja je v enostavnem shranjevanju uporabljenih kljuˇcnih besed, hitra in enostavna ponovljivost in posodobitev ter navsezadnje moˇznost shranjevanja rezultatov v razliˇcnih formatih (TXT, XML, . . . ). Vse to je v veliki meri pogojeno z moˇznostmi, ki jih omogoˇca izbrano orodje eUtils. S poveˇcevanjem ˇstevila argumentov in njihovih vrednosti je takˇsen pristop zamuden in v veliki meri podvrˇzen napakam. Slednje lahko premostimo, ˇce navedena orodja uporabljamo v okviru nam priljubljenega programskega jezika. Za takˇsne naloge se zelo gosto uporablja programski jezik Perl. Mi smo si za to delo izbrali odprtokodni programski jezik in okolje za analizo podatkov R (Ihaka and Gentleman, 1996; R Development Core Team, 2005). Izbira tega jezika je bila pogojena z avtorjevim znanjem in v veliki meri z odprtokodnim projektom Bioconductor (2004), ki deloma ˇze nudi orodja za poizvedovanje po podatkovnih bazah pri NCBI v paketu annotate (Gentleman and Gentry, 2002). Projekt Bioconductor je v osnovi nastal zaradi potrebe po analizi podatkov s podroˇcja mikroˇcipov in danes predstavlja eno od osrednjih orodij za analizo takˇsnih podatkov. Funkcije v paketu annotate omogoˇcajo poizvedovanje in pridobivanje informacij na podlagi t.i. UID oznak, ki skrbijo za enoliˇcno oznaˇcitev znotraj sistema Entrez. Ker pri iskanju literature teh oznak nimamo, smo za naˇs namen dodatno razvili paket entrez , ki vsebuje funkciji entrezSearch in entrezCountByYear.
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
6
Funkcija entrezSearch je osrednja in predstavlja vmesnik med okoljem R in eUtils orodjem eSearch. Slednje omogoˇca iskanje po podatkovnih bazah pri NCBI s poplnoma istimi moˇznostmi, kot jih imamo pri bolj znani razliˇcici na spletnih straneh Entrez portala. Rezultat te funkcije so UID oznake, ki jih potem lahko uporabimo pri drugih funkcijah paketa annotate. Funkcija entrezCountByYear je nadgradnja funkcije entrezSearch in omogoˇca iskanje po letih in predstavi rezultat v obliki, ki je enostavna za tabelarni ali grafiˇcni prikaz. Za uporabo funkcij potrebujemo nameˇsˇcen program R, osnovo projekta Bioconductor in paket entrez , ki je odvisen od paketa XML. Za nadalnje delo je potreben ˇse paket annotate. Dodaten opis in pregled funkcij je v prilogi. Informacije o samem paketu entrez so izpisane spodaj. Package: entrez Type: Package Title: Entrez search in R Version: 0.1 Date: Check NEWS file for changes Author: Gregor Gorjanc Maintainer: Gregor Gorjanc
URL: http://www.bfro.uni-lj.si/MR/ggorjan Description: Proof of concept for Entrez search (via esearch) in R Depends: XML Suggests: annotate License: GPL version 2 or later Built: R 2.3.1; ; 2006-09-13 12:41:38; unix -- File: /usr/local/lib/R/site-library/entrez/DESCRIPTION
3.2.2
Prikaz uporabe
Ko uspeˇsno namestimo potrebno programsko opremo (glej konec predhodnega odstavka) in ˇce imamo dostop do spleta, lahko uporabimo razvite funkcije. V tem delu bomo prikazali njihovo uporabo. Po zagonu programa R moramo najprej naloˇziti paket entrez . R> library(package = "entrez") Loading required package: XML Osnovno iskanje lahko opravimo s pomoˇcjo funkcije entrezSearch npr. po kljuˇcni besedi “Bayes”. Iskanje smo za prikaz omejili le na prvih 5 zadetkov (retMax=5). R> entrezSearch(term = "Bayes", retMax = 5) [1] "Number of matches is greater than 'retMax'" $term [1] "Bayes" $URL [1] "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?term=Bayes&retstart=0&retmax=5"
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
7
$retStart [1] 0 $retMax [1] 5 $retType [1] "uilist" $count [1] 7009 $PMID [1] "16984315" "16979071" "16970551" "16970511" "16969519" V izpisu so podani: kljuˇcne besede (term), uporabljen URL, ki ga lahko uporabimo tudi kasneje npr. v brskalniku, ˇstevilka prvega zadetka (retStart), ˇstevilo zadetkov (retMax), format rezultata (retType), ˇstevilo vseh zadetkov (count) v bazi, ki smo jo uporabili - privzeta je baza PubMed, in na koncu rezultat. Podrobnosti so opisane v prilogi. V primeru uporabe te funkcije na bazi PubMed dobimo kot rezultat PMID oznake, katere lahko uporabimo za razliˇcne namene. Poleg tega se na koncu izpiˇse ˇse opozorilo, da je ˇstevilo vseh zadetkov veˇcje kot argument retMax in da bi bilo smiselno poveˇcati vrednost za argument retMax. Za nadaljnje delo je priroˇcno rezultat shraniti v objekt npr. prikaz. Potem lahko s pomoˇcju funkcije pubmed (v paketu annotate) pretoˇcimo izvleˇcke zadetkov v formatu XML na naˇs raˇcunalnik in jih obdelamo s funkcijami iz paketa annotate: xmlRoot, xmlSApply in buildPubMedAbst. Za prikaz smo izpisali izvleˇcek prvega zadetka. Uporaba funkcije entrezCountByYear je prikazana niˇzje v okviru zadane naloge. R> prikaz <- entrezSearch(term = "Bayes", retMax = 5) [1] "Number of matches is greater than 'retMax'" R> prikaz$PMID [1] "16984315" "16979071" "16970551" "16970511" "16969519" R> library(package = "annotate") Loading required package: Biobase Loading required package: tools R> prikazAbst <- pubmed(prikaz$PMID) R> tmp <- xmlRoot(prikazAbst) R> abst <- xmlSApply(tmp, buildPubMedAbst) R> length(abst)
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
8
[1] 5 R> abst[[1]] An object of class 'pubMedAbst': Title: Parametric and Nonparametric FDR Estimation Revisited. PMID: 16984315 Authors: B Wu, Z Guan, H Zhao Journal: Biometrics Date: Sep 2006 R> writeLines(strwrap(abst[[1]]@abstText)) Nonparametric and parametric approaches have been proposed to estimate false discovery rate under the independent hypothesis testing assumption. The parametric approach has been shown to have better performance than the nonparametric approaches. In this article, we study the nonparametric approaches and quantify the underlying relations between parametric and nonparametric approaches. Our study reveals the conservative nature of the nonparametric approaches, and establishes the connections between the empirical Bayes method and p-value-based nonparametric methods. Based on our results, we advocate using the parametric approach, or directly modeling the test statistics using the empirical Bayes method.
3.2.3
Iskanje
Za pridobitev primerjalne ocene uporabe Bayesovske statistike v “biologiji” po letih smo v podatkovni bazi objav PubMed izvedli iskanje po prispevkih iz obdobja od 1980 do vkljuˇcno 2006. Za primerjavo smo iskanje izvedli tudi za klasiˇcni oz. frekvencistiˇcni pristop k statistiˇcnem sklepanju, kakor tudi za pristop na osnovi verjetja. Kljuˇcne besede za iskanje (tabela 1) smo pripravili s pomoˇcjo tezavra MeSH v okviru sistema Entrez. Takˇsno iskanje je enostavnejˇse in omogoˇca hkratno iskanje po veˇcjem ˇstevilu podobnih kljuˇcnih besed (tabela 1), ki jih doloˇca MeSH sistem klasifikacije prispevkov. Predhodna t.i. “enostavna”iskanja po naslovih in izvleˇckih prispevkov s kljuˇcnimi besedami“bayes*[tiab]”, “likelihood*[tiab]” in “ANOVA[tiab]” so pokazala, da izbrani MeSH izrazi niso vedno najboljˇsi izbor (tabela 2). Najbolj oˇcitna razlika (petkratna) je bila pri izrazu ”likelihood*[tiab]”. Da bi poveˇcali obˇcutljivost oz. iskanja, smo kljuˇcne besede iskali v naslovih in izvleˇckih prispevkov skupaj z izrazi po MeSH klasifikaciji (tabela 2). S tem smo poveˇcali ˇstevilo zadetkov in se najverjetneje bolje ocenili uporabo posameznih statistiˇcnih pristopov.
3.3 REZULTATI Uporaba vseh treh statistiˇcnih pristopov s ˇcasom naraˇsˇca (slika 1), razen padca v zadnjem letu (2006), ker to leto ˇse ni zakljuˇceno. Pri interpretaciji grafa je potrebno biti nekoliko pozoren zaradi uporabe logaritemske skale. Priˇcakovano smo dobili najveˇc zadetkov za frekvencistiˇcno statistiko, ki tako upraviˇcuje svoje drugo ime, klasiˇcna statistika. Za frekvencistiˇcno statistiko sledi pristop z verjetjem in nato
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
Tabela 1: MeSH iskalni izrazi in vnosi po statistiˇcnih pristopih Statistiˇcni pristop MeSH izraz Bayesovski Bayes Theorem
Verjetje
Likelihood Functions
Frekvencistiˇcen
Analysis of Variance
9
MeSH vnosi Theorem, Bayes Bayesian Analysis Analysis, Bayesian Bayesian Forecast Forecast, Bayesian Bayesian Method Method, Bayesian Bayesian Prediction Prediction, Bayesian Function, Likelihood Functions, Likelihood Likelihood Function Maximum Likelihood Estimates Estimate, Maximum Likelihood Estimates, Maximum Likelihood Maximum Likelihood Estimate Analysis, Variance Analyses, Variance Variance Analyses ANOVA Variance Analysis
Tabela 2: Skupno ˇstevilo zadetkov v bazi PubMed glede na uporabljene iskalne izraze ˇ Iskalni izraz Stevilo zadetkov Bayesovski pristop ”Bayes Theorem”[MeSH] 6110 bayes*[tiab] 7562 ”Bayes Theorem”[MeSH] OR bayes*[tiab] 9290 Pristop z verjetjem ”Likelihood Functions”[MeSH] 8419 likelihood*[tiab] 40366 ”Likelihood Functions”[MeSH] OR 44734 likelihood*[tiab] Frekvencistiˇcen pristop ”Analysis of Variance”[MeSH] 149741 ANOVA[tiab] 16221 ”Analysis of Variance”[MeSH] OR ANOVA[tiab] 160441 Bayesovska statistika. Rast uporabe ni enaka med pristopi. Po letu 1988 opazimo znaten porast uporabe frekvencistiˇcne statistike. Druga dva pristopa po uporabi prav tako naraˇsˇcata, a ne tako moˇcno. Uporaba pristopa z verjetjem naraˇsˇca poˇcasneje kot uporaba Bayesovskega pristopa, ˇse posebej po letu 2000. Razloge za to smo nakazali ˇze pri orisu tega pristopa. V tem obdobju je rast uporabe Bayesovske statistike celo veˇcja kor rast uporabe frekvencistiˇcne statistike. Na podlagi pridobljenega ˇstevila zadetkov lahko sklepamo, da uporaba Bayesovske statistike naraˇsˇca in da bo v prihodnjih letih uporaba ˇse narasla. Kljub temu, je absolutno gledano zastopanost Bayesovske statistike velikokrat manjˇsa kot pri frekvencistiˇcni statistiki. Uporabljena metoda za primerjalno oceno uporabe Bayesovske statistike vsekakor ni najboljˇsa, nam pa daje grobo oceno. Najveˇcji problem so iskalni izrazi. To velja ˇse posebej za frekvencistiˇcno statistiko.
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
10
50
100
200
Stevilo zadetkov v PubMed 500 1000 2000
5000 10000
Bayes Verjetje ANOVA
1980
1985
1990
1995
2000
2005
Leto
ˇ Slika 1: Stevilo zadetkov v bazi PubMed glede na statistiˇcni pristop in leto
Za slednjo namreˇc kluˇcna beseda ANOVA ni dovolj reprezentativen, a je teˇzko poiskati kljuˇcne besede, ki bi zajele vse moˇzne analize in ne bi hkrati zajemal tudi druga dva pristopa. Tako lahko upraviˇceno sklepamo, da je ocena za frekvencistiˇcen pristop najbolj podcenjena.
4
IZBOR APRIORNE PORAZDELITVE
Uporaba apriorne porazdelitve (p(θ), 1) je osrednji problem pri Bayesovski statistiki, saj z njo podamo informacijo o neznanih parametrih ˇse preden smo zbrali podatke. S staliˇsˇca Bayesovske statistike je lahko ta informacija “objektivna” kakor tudi “subjektivna”. Vse to je s staliˇsˇca raziskovanja velikokrat sporno. V primeru, da nimamo zanesljive in sploˇsno sprejete predhodne informacije, si ˇzelimo uporabiti apriorno porazdelitev, ki ne bi vsebovala nobene informacije oziroma bi bila “neinformativna”. Doloˇciti takˇsno porazdelitev je vse prej kot enostavno. Bernardo (1997) poroˇca, da “neinformativnih” porazdelitev ni! Za primer vzemimo parameter θ, ki ima zalogo vrednosti na intervalu [a, b]. Z intuitivnega vidika bi bila za apriorno porazdelitev primerna uniformna porazdelitev (p(θ) ∼ U nif orm(a, b)). Zaradi nezadostnih
informacij namreˇc sklepamo, da so vse vrednosti parametra enako verjetne. To ni porazdelitev “brez informacije”, saj porazdelitev pravi, da so vse vrednosti enako verjetne in to je informacija. S staliˇsˇca
informatike je to porazdelitev, ki ima najveˇcjo entropijo2 . Temu pristopu izbora apriorne porazdelitve pravimo naˇcelo pomankljivosti (ang. principle of insufficient reason) in uporabil, ga je ˇze sam Bayes (1763) in kasneje tudi drugi raziskovalci npr. Laplace. Slednji je opazil, da v primeru da parameter θ 2 Entropija
predstavlja mero povpreˇ cne nezanesljivosti sluˇ cajne spremenljivke.
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
11
transformirano npr. ρ = f (θ), uniformna porazdelitev za θ ni nujno “neinformativna” za ρ. Slednje nam vzbuja nekoliko dvoma o izbiri te porazdelitve. V kolikor je zaloga vrednosti za parameter θ cel realni interval (−∞, ∞) ali le del tega intervala npr. (0, ∞) se prav tako pojavljajo problemi z uniformno porazdelitvijo. Porazdelitev na takˇsnem intervalu
enostavno ne more biti primerna (ang. proper) porazdelitev, saj se gostota verjetnosti ne seˇsteje v 1. Pri uporabi Bayesovega izreka sicer ni treba biti pozoren na lastnosti apriorne temveˇc posteriorne porazdelitve. Tako lahko uporabimo tudi neprimerno apriorno porazdelitev, ˇce le ta privede do primerne posteriorne porazdelitve (Kass et al., 1998). Pri izboru apriornih porazdelitev je poleg naˇstetih teˇzav ˇse cela vrsta drugih (Kass and Wasserman, 1996). Tako je tudi izbor “neinformativne” porazdelitve zelo ˇsiroko podroˇcje. Obstaja veˇc pristopov/naˇcinov za izbor apriornih porazdelitev. Pristopi in tudi rezultati le teh se razlikujejo. Naˇstejmo nekaj pristopov (Kass and Wasserman, 1996): Jeffreysova metoda: apriorna porazdelitev mora zadostiti naˇcelu nespremenljivosti (ang. invariance
principle) zaradi transformacij; naˇcelo izbira t.i. privzete/dogovorjene apriorne porazdelitve in predstavlja izhodiˇsˇce za ostale pristope referenˇcne apriorne porazdelitve: Bernardo (1979) je razvil metodo, ki temelji na maksimiziranju
“informacije” med apriorno in posteriorno porazdelitvijo - apriorna porazdelitev naj bi podala ˇcim manj informacije v primerjavi s podatki; ta metoda je najbolj sploˇsna a tudi zelo zahtevna maksimalna entropija: apriorne porazdelitve z veˇcjo entropijo so manj informativne.
Od zgoraj naˇstetih pristopov je Jeffreysova metoda najbolj osnovna. Hkrati je ta metoda tudi izhodiˇsˇce za druge, ki v enostavnih primerih dajejo enake rezultate kot Jeffreysova metoda. Zato si bomo v nadaljevanju ogledali le to metodo. Za primer bomo izbrali meritve (y), ki naj bi bile porazdeljene po normalni porazdelitvi [5]. Ta porazdelitev/model ima dva parametra, ki ju ˇzelimo oceniti. Prvi parameter µ je hkrati tudi priˇcakovana vrednost in lokacijski parameter in ima za zalogo vrednosti celo realno mnoˇzico (−∞, ∞). Drugi parameter σ predstavlja variabilnost porazdelitve in ima za zalogo vrednosti pozitiven del realne mnoˇzice (0, ∞).
p y|µ, σ
2
(y − µ)2 exp − =√ 2σ 2 2πσ 2 1
!
(5)
Za zanimivost dodatno predpostavljajmo, da lahko z neko drugo metodo namesto y izmerimo vrednost z, ki pa lahko ima le dve vrednosti in sicer z = [0, 1], kjer je 0 = y < µ in 1 = y ≥ µ. Sluˇcajna spremenljivka
z ima tako Bernoullijevo porazdelitev [6] s parametrom p = 12 .
p(z|p) = pz (1 − p)1−z
(6)
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
12
Za oceno izgube informacije pri z v primerjavi z y izraˇcunajmo entropijo H(X) ene in druge sluˇcajne spremenljivke. Za z, ki je diskretne narave, znaˇsa entropija:
H(z) = −
X
p(z)logp(z)
z
H(z) = − p0 (1 − p)1−0 log p0 (1 − p)1−0 + p1 (1 − p)1−1 log p1 (1 − p)1−1 p 1−p + log((p) ) = −((1 − p)log(1 − p) + p × log(p)) = − log (1 − p) 1− 21 21 !! 1 1 1−p p ∗ = − log 1− = − log (1 − p) ∗ (p) 2 2 !! !! r r 12 12 1 1 1 1 = − log = − log ∗ ∗ 2 2 2 2 1 = 1bit = − log2 2 1 = − loge = 0.7nat. 2
(7)
Za y, ki je zvezne narave, sta Cover and Thomas (1991) pokazala, da (diferencialna) entropija standardne normalne porazdelitve znaˇsa:
h(y) = −
Z
1 1 loge 2πeσ 2 = loge 2πe12 = 1.4nat. 2 2
p(y)logp(y)dy =
y
(8)
Tako lahko zakljuˇcimo, da je v sluˇcajni spremenljivki z shranjene za 50 % manj informacije kot v standardni normalni porazdelitvi. Vrnimo se na izbor apriorne porazdelitve za parametre normalne porazdelitve [5]. Bayesov izrek [3] za naˇs primer pravi:
p µ, σ 2 |y ∝ p y|µ, σ 2 × p µ, σ 2 , pri ˇcemer verjetje p y|µ, σ 2
(9)
za n podatkov znaˇsa (Sorensen and Gianola, 2002) y¯ =
n Y p yi |µ, σ 2 p yi=1,2,...,n |µ, σ 2 =
Pn
i=1
n
yi
:
i=1
∝ σ
n 2 −2
∝ σ
n 2 −2
exp −
Pn
!
exp −
Pn
(yi − y¯) + n(¯ y − µ) 2 2σ
2 i=1 (yi − µ) 2σ 2 i=1
(10)
2
Za izraˇcun (skupne) posteriorne porazdelitve p µ, σ 2 |y
2
!
.
potrebujemo sedaj ˇse apriorno porazdelitev
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
13
p µ, σ 2 . Kot smo ˇze omenili, je zaloga vrednosti za µ enaka (−∞, ∞), medtem, ko je za σ enaka (0, ∞). Jeffreys je za parametre normalne porazdelitev izbral (Kass and Wasserman, 1996):
π(µ, σ) ∝ σ −1 .
(11)
Najprej opozorilo za spremembo oznake p v π. Jeffreysove apriorne porazdelitve veˇcinoma niso primerne (ang. proper) porazdelitve ampak so zgolj funkcije, ki omogoˇcajo izraˇcun “objektivnih” posteriornih porazdelitev. Zanimivo je to, da nam pri tej funkciji ni potrebno vstavljati nobenih vrednosti kot je to obiˇcajno potrebno pri “pravih” apriornih porazdelitvah npr. srednjo vrednost in varianco za normalno porazdelitev ali pa spodnjo in zgornjo mejo za uniformno porazdelitev. Poleg tega je zgoraj podana apriorna porazdelitev za σ in ne σ 2 . Tako moramo opraviti transformacijo s ˇcimer bomo pokazali, da so Jeffreysove apriorne porazdelitve nespremenljive glede na transformacije. Uporabili bomo pravilo transformacije sluˇcajne spremenljivke:
∂x ∂f −1 (y) −1 py (y) = px (x) = px f (y) ∂y ∂y ∂f −1 σ 2 ∂σ π µ, σ 2 = π(µ, σ) 2 = π µ, f −1 σ 2 2 ∂σ ∂σ 1 2 2 1 1 − 1 1 1 2 −1 2 ∂ σ −2 2 1 = σ σ 2 2 = σ −1 σ −1 = σ −2 ∂σ 2 = σ 2 2 2
(12)
∝ σ −2 .
Dobili smo enako “obliko” porazdelitve kot v [11]. Ta Jeffreysova apriorna porazdelitev je pravzaprav uniformna porazdelitev za µ in uniformna za log(σ) (Kass and Wasserman, 1996). Prvo trditev lahko potrdimo z spodnjim izrazom, drugo pa na enak naˇcin kot smo pokazali transformacijo π(µ, σ) v pi µ, σ 2 . π µ, σ 2 = π µ|σ 2 π σ 2 = konst. × π σ 2 ∝ π σ 2 ∝ σ −2
(13)
Posteriorna porazdelitev je sedaj:
! (yi − y¯)2 + n ∗ (¯ y − µ)2 × σ −2 p µ, σ |y ∝ σ exp − 2σ 2 ! Pn 2 2 (y − y ¯ ) + n ∗ (¯ y − µ) i ∝ σ −n−2 exp − i=1 2σ 2 Z p µ, σ 2 |y p(µ|y) ∝ 2 Zσ p σ 2 |y ∝ p µ, σ 2 |y 2
n 2 −2
Pn
i=1
(14)
µ
Jeffreys je zgoraj omenjeno apriorno porazdelitev doloˇcil s pomoˇcjo Fisherjeve informacije [15].
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
14
Cover and Thomas (1991) sta pokazala, da medtem, ko je entropija mera za nezanesljivost sluˇcajne spremenljivke oziroma “volumen” seta podatkov, je Fisherjeva informacija povezana z minimalno napako ocene parametrov oziroma s “povrˇsino” seta podatkov. Fisherjevo infromacijo lahko zapiˇsemo kot:
I(θ) = Ey
∂ 2 log p y|µ, σ 2 − ∂θi ∂θj
!
.
(15)
Jeffreys je definiral apriorno porazdelitev kot kvadratni koren iz determinante Fisherjeve informacijske matrike (Kass and Wasserman, 1996):
π(θ) ∝
p |I(θ)|.
(16)
Za naˇs primer bomo najprej izpeljali Fisherjevo informacijska matriko:
I(θ) = Ey −
∂ 2 log(p(y|µ,σ2 )) ∂µ∂µ ∂ 2 log(p(y|µ,σ2 )) ∂σ2 ∂µ
∂ 2 log(p(y|µ,σ2 )) ∂µ∂σ2 ∂ 2 log(p(y|µ,σ2 )) . ∂σ2 ∂σ2
(17)
Logaritem verjetja 10 in posamezni odvodi znaˇsajo:
l = log p y|µ, σ
2
n = − log σ 2 − 2
Pn
i=1
(yi − y¯) + n(¯ y − µ) 2 2σ
n¯ y − nµ ∂l = ∂µ σ2 2 ∂ l n =− 2 ∂µ∂µ σ ∂2l = −σ −4 (n¯ y − nµ) ∂µ∂σ 2 ! n ∂l n 1 X (yi − y¯) + n(¯ y − µ) =− 2 + 4 ∂σ 2 2σ 2σ i=1 ! n X ∂ 2l n 6 (yi − y¯) + n(¯ y − µ) . =− 4 −σ ∂σ 2 ∂σ 2 2σ i=1
(18)
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
15
Od vsakega drugega odvoda moramo izraˇcunati ˇse Ey (−x), to je minus priˇcakovane vrednosti glede na y:
n n −Ey − 2 = 2 σ σ X Ey (y)/n − nµ y − nµ) = σ −4 (Ey (n¯ y) − nµ) = σ −4 n −Ey −σ −4 (n¯ = σ −4 (n(nµ)/n − nµ)
−Ey −
!!
n X n 6 (yi − y¯) + n(¯ y − µ) − σ 2σ 4 i=1
= σ −4 (nµ − nµ) = σ −4 (nµ − nµ) × 0 = 0 ! n X n 6 = 4 + σ Ey (yi − µ) 2σ i=1
n n X X n n 6 6 (E (y ) − µ) = (µ − µ) + σ + σ y i 2σ 4 2σ 4 i=1 i=1 n n = 4 + σ6 × 0 = 4 . 2σ 2σ (19)
=
Matrika Fisherjeve informacije je tako diagonalna:
I(θ) =
n σ2
0
0
n 2σ4 .
!
(20)
Sedaj pa ˇse kvadratni koren determinante te matrike in dobimo Jeffreysovo apriorno porazdelitev za naˇs primer:
π µ, σ
2
p = |I(θ)| =
= konst.σ −3
r
n n × 4 −0×0= σ2 2σ
r
√ n2 n2 2−1 σ −6 = n2×0.5 2−1×0.5 σ −6×0.5 = 2σ 6 (21)
π µ, σ 2 ∝ σ −3 . To pa ni enaka apriorna porazdelitev kot v [12]. Tukaj smo naleteli na nekonzistentnost Jeffreysove metode. Metoda namreˇc v primeru veˇcjega ˇstevila parametrov daje drugaˇcne apriorne porazdelitve kot pri enem samem parametru na enkrat. Ravno zaradi tega so se razvile tudi druge metode, ki smo jih ˇze naˇsteli. Zgornji problem je Jeffreys reˇseval z ad-hoc pristopi. Razˇclenimo skupno apriorno porazdelitev:
p µ, σ 2 = p µ|σ 2 p σ 2
(22)
in izraˇcunajmo apriorne porazdelitve za µ in σ 2 posebej. Pri p µ|σ 2 predpostavljamo, da poznamo σ 2 in tako je Jeffreysova apriorna porazdelitev za µ:
p p p µ|σ 2 = I(µ) = n/σ 2 ∝ konst.
(23)
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
16
kar predstavlja uniformno porazdelitev na celem realnem intervalu, kakor smo omenili ˇze prej [13]. Sedaj pa ˇse za σ 2 :
p p p σ 2 = I(σ 2 ) = n/σ 4 = n0.5 σ −4×0.5 ∝ σ −2 ,
(24)
kar je rezultat, ki smo ga iskali in predstavlja uniformno porazdelitev za log(σ) (Kass and Wasserman, 1996).
5
SKLEPI
Povzamemo lahko, da se Bayesovska statistika pogosteje uporablja za analizo podatkov (s podroˇcja biˇ vedno pa je njen deleˇz manjˇsi, kar je glede na zahtevnost izraˇcunov, predvsem pa izbora ologije). Se apriornih porazdelitev tudi razumljivo. Obstaja veˇc pristopov k izboru t.i. “neinformativnih” porazdelitev. Jeffreysova metoda je ena najbolj osnovnih in temelji na Fisherjevi informaciji. Medtem ko, entropija predstavlja mero nezanesljivosti sluˇcajne spremenljivke oziroma “volumen” seta podatkov, je Fisherjeva informcaija povezana z minimalno napako ocene parametrov oziroma s “povrˇsino” seta po datkov. Jeffreysovo metodo smo prikazali na primeru izbora apriornih porazdelitev parametrov µ, σ 2 normalne porazdelitve. Za p(µ) smo dobili uniformno porazdelitev, za p σ 2 pa “porazdelitev”, ki je ˇ Jeffreysova proporcionalna izrazu σ −2 . Slednji izraz je pravzaprav uniformna porazdelitev za log(σ). Zal metoda ni konzistentna pri hkratnem “izraˇcunu” skupne apriorne porazdelitve veˇcih parametrov.
6
VIRI
Bayes, T. (1763). An essay towards solving a problem in the doctrine of chances. Phil. Trans. Roy. Soc. London, 53:370–418. Reprinted in Biometrika, 45:293-315 (1958) and in Press (1989). Bernardo, J. M. (1979). Reference posterior distributions for bayesian inference (with discussion). J. R. Stat. Soc. B Stat. Methodol., 41:113–147. Bernardo, J. M. (1997). Noninformative priors do not exist: A discussion. J. Stat. Plan. Inference, 65:159–189. Bernardo, J. M. (2003). Probability and statistics, chapter Bayesian statistics. Encyclopedia of life support systems (EOLSS). Eolss Publishers, Oxford ,UK. Casella, G. and George, E. I. (1992). Explaining the Gibbs sampler. Am. Stat., 46(3):167–174. Cover, T. M. and Thomas, J. A. (1991). Elements of information theory. John Wiley & Sons. Gelfand, A. E., Hills, S. E., Racine-Poon, A., and Smith, A. F. M. (1990). Illustration of Bayesian inference in normal data using Gibbs sampling. J. Am. Stat. Assoc., 85(412):972–985. Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. (2004). Bayesian data analysis. Texts in statistical science. Chapman & Hall / CRC, 2nd edition. Geman, S. and Geman, D. (1984). Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. IEEE Trans. Pattern Anal. Mach. Intell., 6(6):721–741.
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
17
Gentleman, R. and Gentry, J. (2002). Querying pubmed. R News, 2(2):28–31. Gentleman, R. C., Carey, V. J., Bates, D. J., Bolstad, B. M., Dettling, M., Dudoit, S., Ellis, B., Gautier, L., Ge, Y., Gentry, J., Hornik, K., Hothorn, T., Huber, W., Iacus, S., Irizarry, R., Leisch, F., Li, C., Maechler, M., Rossini, A. J., Sawitzki, G., Smith, C., Smyth, G. K., Tierney, L., Yang, Y. H., and Zhang, J. (2004). Bioconductor: Open software development for computational biology and bioinformatics. Genome Biology, 5(10):R80. Gilks, W. R., Richardson, S., and Spiegelhalter, D. J., editors (1998). Markov Chain Monte Carlo in practice. Chapman & Hall / CRC. Green, P. J. (1995). Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika, 82:711–732. Hastings, W. K. (1970). Monte Carlo sampling-based methods using Markov chains and their applications. Biometrika, 57(1):97–109. Ihaka, R. and Gentleman, R. (1996). R: A language for data analysis and graphics. J. Comput. Graph. Stat., 5(3):299–314. Kass, R. E., Carlin, B. P., Gelman, A., and Neal, R. M. (1998). Markov chain Monte Carlo in practice: a roundtble discussion. Am. Stat., 52(2):93–100. Kass, R. E. and Wasserman, L. (1996). The selection of prior distributions by formal rules. J. Am. Stat. Assoc., 91(435):1343–1370. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. (1953). Equations of state calculations by fast computing machines. J. Chem. Phys., 21(6):1087–1091. NCBI. (2005). National center for biotechnology information. R Development Core Team (2005). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0. Sayers, E. and Wheeler, D. (2005). NCBI short courses, chapter Building Customized Data Pipelines Using the Entrez Programming Utilities (eUtils). The National Library of Medicine. Sorensen, D. and Gianola, D. (2002). Likelihood, Bayesian, and MCMC Methods in Quantitative Genetics. Statistics for Biology and Health. Springer. New York, Springer: 732 str.
A
PRILOGA
entrezSearch
Entrez search
Description entrezSearch performs search for primary IDs in PubMed via Entrez utility esearch at NCBI.
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
18
Usage
entrezSearch(term, db=NULL, WebEnv=NULL, useHistory=FALSE, queryKey=NULL, tool=NULL, email=NULL, field=NULL, relDate=NULL, minDate=NULL, maxDate=NULL, dateType=NULL, retMax=100, retStart=0, retType="uilist", xmlKeep=FALSE, entrezURL="http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?", warn=TRUE)
Arguments Arguments to this function are simply arguments for Entrez utility esearch and you should look under details and references therein for detailed explanation: term
character, search string
db
character, database name
WebEnv
character, web environment
useHistory
logical, maintain results
queryKey
integer, history search number
tool
character, resource which is using Entrez links
email
character, your contact
field
character, PubMed field
relDate
integer, number of days preceding today
minDate
character, limit results by date
maxDate
character, limit results by date
dateType
character, date type
retMax
integer, number of matches to retrieve
retStart
integer, sequential number of the first record retrieved
retType
character, return type (uilist for PMIDs or count for number of matches, which should be much faster with many matches)
xmlKeep
logical, keep retrieved XML file
entrezURL
character, URL for NCBI utilities
warn
logical, warn if there are more hits than retMax
Details This function is a wrapper around Entrez utility esearch http://eutils.ncbi.nlm.nih.gov/entrez/query/static/esearch_help.html
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
19
at NCBI http://www.ncbi.nlm.nih.gov/. For information on search string i.e. term see the PubMed or Entrez help at: http://eutils.ncbi.nlm.nih.gov/entrez/query/static/help/pmhelp.html and http://www.ncbi.nlm.nih.gov/entrez/query/static/help/helpdoc.html. For other Entrez utilities, user requirements and any additional information look at: http://eutils.ncbi.nlm.nih.gov/entrez/query/static/eutils_help.html, http://www.ncbi.nlm.nih.gov/books/bv.fcgi?rid=coursework.chapter.eutils, http://www.ncbi.nlm.nih.gov/Class/PowerTools/eutils/course.html. User requirements state that one should not send more than one request per 3 seconds! Term string term can use some characters that raise problems in R as well as in some browsers. Therefore some characters must be escaped, while for others the whole string in term is processed with URLencode to circumvent problems.
Value A list with used search options and results
Author(s) Gregor Gorjanc
See Also pubmed in Bioconductor package annotation and URLencode
Examples ## You need internet access for example to work if(capabilities()["http/ftp"] && .Platform$OS.type == "unix" && !is.null(nsl("cran.r-project.org"))) { ## Are there any articles on Bioconductor in PubMed? (search <- entrezSearch("Bioconductor[tiab]")) ## Fetch other data based on retrieved UID's library(annotate) x <- pubmed(search$PMID) a <- xmlRoot(x) BioCAbst <- xmlSApply(a, buildPubMedAbst) length(BioCAbst) }
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
# ## entrezSearch. R # ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # ## What : Entrez search # ## $ Id : entrezSearch.R , v 1 .8 2006 / 09 / 11 20:48:42 ggorjan Exp $ # ## Time - stamp : <2006 -09 -11 22:41:13 ggorjan > # ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - entrezSearch <- function ( term , db = NULL , WebEnv = NULL , useHistory = FALSE , queryKey = NULL , tool = NULL , email = NULL , field = NULL , relDate = NULL , minDate = NULL , maxDate = NULL , dateType = NULL , retMax =100 , retStart =0 , retType = " uilist " , xmlKeep = FALSE , entrezURL = " http : / / e u t i l s . n c b i . n l m . n i h . go v/ entrez / eutils / esearch.fcgi ? " , warn = TRUE ) { if ( missing ( term )) stop ( " ' term ' is missing " ) tmp <- list () tmp $ term <- term # # Build URL for search in PubMed for number of matched publications and # # fix special characters in it. tmp $ URL <- paste ( entrezURL , " term = " , URLencode ( term ) , sep = " " ) if ( ! is.null ( db )) { tmp $ URL <- paste ( tmp $ URL , " & db = " , db , sep = " " ) tmp $ db <- db } if ( ! is.null ( WebEnv )) { tmp $ URL <- paste ( tmp $ URL , " & WebEnv = " , WebEnv , sep = " " ) tmp $ UsedWebEnv <- WebEnv } if ( ! is.null ( queryKey )) { tmp $ URL <- paste ( tmp $ URL , " & query _ key = " , queryKey , sep = " " ) tmp $ UsedQueryKey <- queryKey } if ( useHistory ) { tmp $ URL <- paste ( tmp $ URL , " & usehistory = y " , sep = " " ) tmp $ useHistory <- useHistory } if ( ! is.null ( tool )) { tmp $ URL <- paste ( tmp $ URL , " & tool = " , tool , sep = " " ) tmp $ tool <- tool } if ( ! is.null ( email )) { tmp $ URL <- paste ( tmp $ URL , " & email = " , email , sep = " " ) tmp $ email <- email } if ( ! is.null ( field )) { tmp $ URL <- paste ( tmp $ URL , " & field = " , field , sep = " " ) tmp $ field <- field
20
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
21
} if ( ! is.null ( relDate )) { tmp $ URL <- paste ( tmp $ URL , " & reldate = " , relDate , sep = " " ) tmp $ relDate <- relDate } if ( ! is.null ( minDate )) { tmp $ URL <- paste ( tmp $ URL , " & mindate = " , minDate , sep = " " ) tmp $ minDate <- minDate } if ( ! is.null ( maxDate )) { tmp $ URL <- paste ( tmp $ URL , " & maxdate = " , maxDate , sep = " " ) tmp $ maxDate <- maxDate } if ( ! is.null ( dateType )) { tmp $ URL <- paste ( tmp $ URL , " & datetype = " , dateType , sep = " " ) tmp $ dateType <- dateType } if ( ! is.null ( retStart )) { tmp $ URL <- paste ( tmp $ URL , " & retstart = " , retStart , sep = " " ) tmp $ retStart <- retStart } if ( ! is.null ( retMax )) { tmp $ URL <- paste ( tmp $ URL , " & retmax = " , retMax , sep = " " ) tmp $ retMax <- retMax } if ( retType == " count " ) { tmp $ URL <- paste ( tmp $ URL , " & rettype = count " , sep = " " ) tmp $ retType <- retType } else { tmp $ retType <- retType } # # Query PubMed tmp $ xml <- xmlTreeParse ( tmp $ URL , isURL = TRUE ) # # Get various parameters tmp $ count <- as.integer ( xmlValue ( xmlChildren ( tmp $ xml $ doc $ children $ eSearchResult ) $ Count )) if ( retType ! = " count " ) {
tmp $ retMax <- as.integer ( xmlValue ( xmlChildren ( tmp $ xml $ doc $ children $ eSearchResult ) $ RetMax ) } else { tmp $ retMax <- NA } if ( warn & & retType ! = " count " & & tmp $ count > tmp $ retMax ) print ( " Number of matches is greater than ' retMax ' " ) if ( useHistory ) { tmp $ WebEnv <- xmlValue ( xmlChildren ( tmp $ xml $ doc $ children $ eSearchResult ) $ WebEnv ) tmp $ queryKey <- xmlValue ( xmlChildren ( tmp $ xml $ doc $ children $ eSearchResult ) $ QueryKey )
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
} if ( tmp $ count ! = 0) { # Proceed only if there are any matches # # Extract PMIDs if ( retType ! = " count " ) { tmp $ PMID <- xmlChildren ( tmp $ xml $ doc $ children $ eSearchResult ) $ IdList tmp $ PMID <- xmlSApply ( tmp $ PMID , xmlValue ) names ( tmp $ PMID ) <- NULL } # # Keep the whole XML if wanted if ( ! xmlKeep ) tmp $ xml <- NULL } tmp } # ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # ## entrezSearch. R ends here
22
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
entrezCountByYear
23
Number of publications in NCBI Entrez by year
Description entrezCountByYear counts number of publications in NCBI Entrez database by year for given search term and date of publication limits. Usage
entrezCountByYear(term, dp, sleep=3, db="pubmed", echo=TRUE, retType="count", ...)
Arguments term
character, search term (look under details)
dp
character or numeric, publication date limits in PubMed format (look under details)
sleep
integer, number of seconds between requests with NCBI (look under details)
db
character, database name
echo
logical, print progress
retType
character, return type as described in entrezSearch - defaults to ˇ count" due to speed
...
further arguments for entrezSearch
Details This function is a wrapper around entrezSearch function and all the details described in its help page apply also here! Value Table with number of publications per year and attribute term with used search term. Author(s) Gregor Gorjanc See Also pubmed in Bioconductor package annotation and entrezSearch
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
Examples ## You need internet access for example to work if(capabilities()["http/ftp"] && .Platform$OS.type == "unix" && !is.null(nsl("cran.r-project.org"))) { ## Number or "prion" hits by year prion <- entrezCountByYear(term="prion[tiab]", dp=1985:1990) ## Plot hits plot(prion, typ="l", xlab="Year", ylab="Number of publications in PubMed") }
24
Gorjanc G. (Bio)informatika in Bayesovska statistika. Semin. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
# ## entrezCount B yY e ar . R # ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # ## What : Count number of publications in Entrez by year for given search # ##
term
# ## $ Id : entrezCountByYear.R , v 1 .9 2006 / 09 / 11 20:48:06 ggorjan Exp ggorjan $ # ## Time - stamp : <2006 -09 -13 01:14:35 ggorjan > # ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - entrezCountB y Ye ar <- function ( term , dp , sleep =3 , db = " pubmed " , echo = TRUE , retType = " count " , ... ) { # # Checks if ( missing ( term )) stop ( " missing search term " ) if ( missing ( dp )) stop ( " missing date limits " ) # # The machine result <- vector ( mode = " integer " , length = length ( dp )) names ( result ) <- dp for ( i in seq ( along = dp )) { # # Hmm , eUtils documentation says that one needs to use + AND + , but it does # # not work this way string <- paste ( " ( " , term , " ) " , " AND " , dp [ i ] , " [ dp ] " , sep = " " ) result [ i ] <- entrezSearch ( term = string , db = db , retType = retType , ... ) $ count if ( interactive () & echo ) cat ( dp [ i ] , " \ b " ) Sys.sleep ( sleep ) } if ( interactive () & & echo ) cat ( " \ n " ) attr ( result , " searchTerm " ) <- term as.table ( result ) } # ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # ## entrezCount B yY e ar . R ends here
25