Tuloksia 25.5

Päivä           Smoltteja
22.5             0
23.5             2
24.5             0
25.5             0

Päiviteltiin taas hieman mallia viimeisimmällä, joskaan ei niin kovin tuoreella saatikka runsaalla datalla. Smoltittomia päiviä on sattunut tässä jo useampi, joko rupeaa kalat vedestä loppumaan??

Kokonaismäärän N keski- ja hajontalukuja

Kokonaismäärän N keski- ja hajontalukuja

Arvoi saatiin noin 120 000 iteraatiolla, joista poistettiin 24 000 ensimmäistä ja valittiin joka 15:sta iteraatio konvergoitumisen varmistamiseksi ja autokorrelaation pienentämiseksi. Kuvaajassa näkyvä piikikkyys johtuu arvojen pyöristämisestä mallin sisällä.

N255

Kokonaismäärän N posteriorijakauma

Uusin mallin antama arvio kokonaismäärän todennäköisyysjakauman keskiarvosta on 2272, keskihajonnalla 1181. 95 % todennäköisyysväli löytyy arvojen 582 ja 4975 väliltä. Arvio on siis edelleenkin hyvin varovainen, mutta arviot tarkentuvat kun tieto lisääntyy!

Rysän uuden pyyntipaikan huomioiminen mallin sisällä

Tänään tuli puheeksi rysän uusi sijainti ja sen huomioiminen mallin sisällä. Rysä on aikaisemmin sijainnut joen haarautumiskohdassa ja on nyttemmin (13.5.) siirretty joen länsihaaran puolelle, hieman haarautumiskohdan yläpuolelle.

Saimme myös uutta mandariinikoetietoa, jonka mukaan nyt kaikkien kolmen kokeen aikana “vapautetusta” 188 mandariineista 139 on päätynyt joen länsihaaraan. Olemme käyttäneet tätä estimaattina smolttien vaeltamisesta, joten voisimme myös olettaa rysän pyytävän nyt tehokkaammin sen ollessa voimakkaimmin virtaavassa paikassa.

Hahmottelin helpohkon pikaratkaisun ongelmalle, jossa uutta rysää lähdetään yksinkertaisesti arvoimaan kokonaan uudella priorijakaumalla. Tämä ratkaisu ei edellyttäisi suuria muutoksia malliin:

Fyke_q[d]<-q_F*Fyke_op[d]*step(32-d)+q_F2*Screw_op[d]*step(d-32)

q_F2~dbeta(1,3)

Step-funktiot palauttavat yksinkertaisesti nollan erotuksen ollessa <0, jolloin summattava saa arvon nolla. Siis päivillä 1-32 käytetään rysän aiempaa priorijakaumaa q_F ja 32-60 uutta q_F2:ta. q_F2:lle annettu arvio hieman suuremmasta pyyntitehosta.

Tämä ratkaisu ei välttämättä ole ongelmaton. Rysän pyydystävyys päivittyy mielestäni datalla samaan tapaan kuin aikasemminkin. Se tarkoittaa, että kaikki rysän aikasemmassa paikassa sijaitessaan keräämä data pyydystettävyydestä siirtyy nyt “uudelle rysälle”, eli uudelle sijainnille. Tämä on mielestäni ristiriitaista sen kanssa, että rysän uuteen sijaintiin perustuva priori antaa nimenomaan ymmärtää, ettei sen pyydystävyys ole vain sen ominaisuus vaan riippuvainen myös sen sijannista. Tähän on syytä paneutua lisää ja katsoa miten tarkalleen riippuvuudet kulkevat pyydystävyys-mallin sisällä.

 

Pohdintaa länsi- ja itähaaran mallintamisesta

Pohdin aikaisemmin mallin kehittämisen näkökulmasta mahdollisuutta mallintaa erikseen joen länsi- ja itähaaraan päätyvien smolttien määrää. Tämä olisi hyödyllistä mikäli halutaan jatkossa selvittää esimerkiksi kuinka suuri osa smolteista joutuu länsihaaraan ja kuinka paljon voimalan läpi kulkeminen vaikuttaa niiden kuolleisuuteen. Olen nyt yrittänyt hieman hahmotella mitä tämänkaltainen mallinnus edellyttäisi. Ongelma ei sinänsä varmastikkaan olet teknisesti mitenkään hankala, mutta omalle hahmottamiselleni löytyy siitäkin runsaasti haastetta ja on siksi hyvää pohdintaa pelkästään oppimista ajatellen.

Tällä hetkellä vedessä on taas kaksi pyydystä, rysä ja smolttiruuvi. Ruuvi pyyntää joen itähaaraa, jossa siis kulkevat kalaportaat. Ruuvin sijainti on haarautumiskohdasta alavirtaan. Rysä on nyt siirretty pyytämään länsihaaraa, eli voimalan puolta, mutta se sijaitsee haarautumiskohdasta hieman ylävirtaan. Sijainnista johtuen ei voida olettaa varmasti, etteivät kalat voisi vielä rysän ohitettuaan siirtyä itähaaraan ja otaksuisin että voi olla jopa mahdollista smolttien ruuvin ohitettuaan palata vielä ylävirtaan ja mahdollisesti toiseen haaraan. Tämä varmasti lisää tuntematonta vaihtelua, mutta pohdintani oletuksena onkin nyt tilanne, jossa molemmat pyydykset sijaitsisivat selkeästi haarautumiskohdasta alavirtaan.

Millä tavalla tätä ongelmaa voitaisiin sitten mielestäni mallintaa?

Päivittäin vaeltavan smolttimäärän muutujan n[i]:n suhteen ei muutoksia tarvittaisi. Smoltteja saapuisi siis haarautumiskohtaan n[i]:n suuruinen määrä, jonka jälkeen ne valitsivat joko itäisen, tai läntisen haaran. Näille voisimme valita priorijakaumat tutkimuksen alussa tehdyn mandariinkokeen perusteella. Jakaumana käyttäsin tuttua beta-jakaumaa:

prob_w~dbeta(7,3)       #prioriodennäköisyys länsihaaran päätymiselle
prob_e~dbeta(3,7)       #prioritodennäköisyys itähaaraan päätymiselle

Nämä valitisin parametreiksi eri haaroihin päätyvien smolttien määrää kuvaavan binomijakaumaan:

nw[i]~dbin(prob_w, n[i])     #jakauma länsihaaraan matkaaville smolteille
ne <- n[i]-nw[i]                  #itähaaraan matkaavat smoltit

Saalis mallinnettaisiin molemmille haaroille erikseen binomijakaumana:

FykeX[i]~dbin(Fyke_q[i], nw[i])         #jakauma länsihaaraan saaliismäärälle
ScrewX[i]~dbin(Screw_q[i], nw[i])    #jakauma itähaaran saalismäärälle

Pyydystävyyden prioriarvoja pitäisi hiukan muuttaa. Rysälle prioritietoa mandariinikokeesta ei ole, ruuvin suhteen on huomiotava että nyt pyydystävyys mallintaisi vain itähaaran osuutta, ei smolttien kokonaismäärän pyydystävyyttä. Mandariinkokeen mukaan itähaaraan päätyi kaksi mandariinia, joista kaksi jäi pyydykseen:

q_S~dbeta (2,1)       #prioritodennäköisyys ruuvin pyydystävyydelle
q_F~dbeta (5,3)       #prioritodennäköisyys rysän pyydystävyydelle, mielivaltainen arvio

Näiden päivittäminen voisi tapahtua mielestäni samaan tapaan kuin ennenkin.

Priorien miettiminen on tässä tapauksessa suhteellisen helppoa. Itselleni alkoi tässä vaiheessa muodostua ongelmaksi juurikin ongelman laskennallinen ja looginen puoli. Itselleni on hankalaa hahmottaa jakaumien välillä tapahtuvia matemaattisia operaatiota ja yhdistää mallin osia loogisesti toisiinsa. Nyt minusta tuntuu, ettei malli näin päivitä ollenkaan eri haaroihin päätymismahdollisuuksien todennäköisyysjakaumia. Joudun vielä pohtimaan tarkemmin miten tämä tapahtuu. Myös kyseisten jakaumien suhteellisuus tulisi ottaa huomioon, kävisikö tämä vain jakamalla kukin molempien summalla? Näinkö saataisiin myös itähaaraan pyydystettävyys siirtymään posteriorijakaumiin?

Neuvot ovat tervetulleita!

Bayesiläinen mallikeskiarvoistaminen

Bayesiläinen keskiarvoistaminen yleisesti

Bayesilaiseen mallivalintaan kuuluu käytännössä kaksi osaa: ”havainnot generoineen” uskottavuusfunktion sekä tämän funktion parametrien priorin määrittäminen. Ensin mainittuun kuuluu epävarmuutta siinä missä jälkimmäiseenkin. Klassisessa tilastotieteessä on tapana kritisoida priorivalinnan subjektiivista luonnetta antaen ymmärtää tämän valinnan olevan jossain määrin tutkijan omista epätieteellisistä mieltymyksistä riippuvaista. Kuitenkin myös uskottavuusfunktion valinta on tutkijan enemmän tai vähemmän perusteltu valinta siinä missä priorikin. Molemmat valinnat perustuvat (yhtä lailla) tutkijan henkilökohtaiseen ennakkokäsitykseen havainnot tuottaneesta (fysikaalisesta) mekanismista. Uskottavuusfunktion valinta onkin harvoin jos koskaan aukottomasti perusteltu, vaan myös siihen liittyy käytännössä aina huomionarvoista epävarmuutta.
Bayesilaisessa mallin keskiarvoistamisessa (bayesian model averaging BMA) voidaan huomioida mallivalinnan epävarmuus kokonaisuudessaan pitäen sisällään sekä uskottavuusfunktion että priorin tutkijan subjektiiviseen harkintaan perustuneen mallivalinnan pätevyyden. Bayesilainen keskiarvoistaminen perustuu mallikokonaisuuden johtamiselle keskiarvoistamalla malli yli osamalliensa painotten kunkin osamallin vaikutusvaltaa niille datan valossa päivitettyjen posterioritodennäköisyyksien suhteessa (bayesin teoreeman avulla). Posterioritodennäköisyydeltään suurempi malli saa enemmän, pienempi malli vähemmän painoarvoa mallikokonaisuutta johdettaessa.
Tällöin oletetaan (kokonaistodennäköisyyden kaavassa), että jokin esitetyistä malleista on tosi. Esimerkiksi, jos jompikumpi kahdesta mallista M1 tai M2 on tosi ja on havaittu aineisto d, tällöin mallin yksi posterioritodennäköisyys saadaan kaavasta: p(M=1|d)=[p(d|M=1)*p(M=1)]/[p(d|M=1)*p(M=1)+ p(d|M=2)*p(M=2)]. BMA painottaa esimerkiksi posterioriennustavaa jakaumaa johdettaessa kutakin mallia edellä laskettujen posteriorien suhteissa.
Kun bayesilainen päättely yhdistetään toimintaan valitsemalla se toimintamalli, joka tuottaa parhaimman odotusarvoisen hyödyn, keskiarvoistetuilla malleilla saadaan oletettavasti parempi päätöksiä aikaisesti kuin yksittäisillä mallikokonaisuuden osamalleilla. Tämän voisi pyrkiä hahmottamaan siten, että keskiarvoistaminen huomioi paitsi useaan malliin liittyvän epävarmuuden, myös niihin liittyvän informaation. Suuremmalla määrällä informaatiota voidaan oletettavasti tehdä perustellumpi päätöksiä.
Kuten informaatioteoriaan perustuvat mallivalinnan kriteeri DIC ja BIC, myös BMA huomioi mallien yksinkertaisuuden (parsimony) arvona sinänsä mallin dataan sopivuuden (goodness-of-fit) rinnalla. Jälleen – kuten informaatioteoreettiset mallivalinnan kriteerit – myöskään BMA ei kykene huomioimaan sitä, jos kaikki estimointiin käytetyt/vertaillut mallit ovat huonoja. Roskan keskiarvo on keskimääräinen roska.
Siksi mallivalinnan perusteltavuuteen tarvitaan ennustavia jakaumia, joiden avulla voidaan mallin toimivuutta vertailla aitoon kerättyyn aineistoon. Ei voida tyytyä olettamaan, että joku osamalleista on tosi ja että mallikokonaisuus siksi toimii. Tulee tarkastella, kuinka hyvin (keskiarvoistetun) mallin ennusteet vastaavat todellisuutta. Paras tilanne on, jos voidaan ennustaa tulevia todellisia havaintoja, kuten Vantaanjoen smoltteja arvioitaessa johtamalla seuraavan päivän havaintojen posterioriennustava jakauma. Ennusteiden ja havaintojen yhteensopivuuden tarkastelu voidaan suorittaa paitsi visuaalisella tarkastelulla, myös määrämittaistamalla arvio esimerkiksi laskemalla niin sanottuja ”bayesilaisia p:n arvoja” (bayesian p-values). Logiikka toimii analogisesti nollahypoteesin testauksen logiikan kanssa, ainoastaan, että nollahypoteesi koskee (keskiarvoistettua) mallivalintaa sisällöttömän nollahypoteesin sijaan. Voidaan siis arvioida todennäköisyys havaita havaitun kaltainen (tai harvinaisempi/oudompi) aineisto, jos (keskiarvoistettu) mallivalintamme olisi oikea. Jos posterioriennusteet ovat huonoja, on malli todennäköisesti huono.
Smolttien totaalia estimoitaessa on keskeisenä ongelmana ollut juuri mallivalinta liittyen smolttien lähtötodennäköisyyksien jakauman muodolle tarkasteluaikavälin päivien yli. Tällä mallivalinnalla on myös kohtalaisen paljon vaikutusta saatuihin totaaliestimaatteihin. Kuten aina, yhtäkään mallivalintaa ei voida tässäkään varmuudella pitää totena tai epätotena, vaan ne ovat kaikki enemmän tai vähemmän epävarmoja. Lisäksi keskiarvoistamisen toimivuuden kannalta hyvänä asiana voitaneen pitää sitä, että mallit huomioivat ikään kuin hieman erilaisia vaihtoehtoisia näkökantoja smolttien muuttokäyttäytymiseen. Näin ollen on oletettavaa, että BMA -menetelmää soveltamalla saadaan mallivalinnan epävarmuuden paremmin huomioiva ja oletettavasti luotettavampi estimaatti mereen Vantaanjoesta muuttavien smolttien kokonaismäärälle.

Malliparametrien vertailu BMA-menetelmällä

Osana päivän tehtävänantoa yritimme soveltaa oppimaamme Carling & Chib-menetelmää omaan malliimme. Muodostimme vertailevan mallin kolmesta vaihtoehtoisesta smolttien lähtemistodennäköisyyttä kuvaavasta p-jakaumastamme.

Jakaumat muotoiltiin seuraavasti:

#p for the linear distribution
p[i,1]<-1/60

#p for the normal distribution
pl[i] <- exp(- pow((i-myy_p)/sd_p,2)*0.5)
p[i,2] <- pl[i] / sum(pl[1:60])

#p for the log-normal distribution
pn[i] <- (1/i)*exp(- pow(log(i)-location, 2) / scale)
p[i,3] <- pn[i] / sum(pn[1:60])

# Number of leaving trouts at day i
n[i]~dbin(p[i,model],N)

Muodostettiin jakaumia vertaileva malli. Painot eri malleille annettiin omiin kokemuksiimme niiden toiminnasta perustuen. Log-normaali jakauma sai pienimmän painon, sillä sen toimivuus on ollut kyseenalaista.

#prior for the model BMA for p

model~dcat(z[1:3])
z[1]<-2/5 #weight for the uniform distribution
z[2]<-2/5 #weight for the normal distribution
z[3]<-1/5 #weight for the log-normal distribution
Z[1]<-equals(model, 1)
Z[2]<-equals(model, 2)
Z[3]<-equals(model, 3)

Lisättiin vielä alkuarvot ladattavaksi ketjuille ennen mallin ajoa.

#Inits: use different initial value for each chain
list(model=1)
list(model=2)
list(model=3)

Alkuarvojen antamisessa ennen ajoa oli pieni ongelma, joka poistui kuitenkin, kun ne generoitiin ensin ja ladattiin vasta sitten.

1

Ajon aikana todennäköisyyksien Z[1] ja Z[2] havaittiin selkeää ketjujen hyppimistä. Miksi Z[3] ei toiminut samoin?

2

Myös statistiikka Z[3]:n osalta näytti oudolta. Sen arvo oli tasan nolla.

Model-jakaumakin näytti sivuuttavan Z[3]:n kokonaan. Mistä tämä johtuu?

4

 

Opittavaa BMA:sta – teoria ja käytäntö

Tämän kappaleen alkuun on rehellisyyden nimissä todettava, että harva asia Bayeslaisen analyysin teoriassa tuskin tuntuu yhdellekään kurssilaiselle vielä täysin selkeältä. Bayesin kaavaan ja todennäköisyyksien päivittämiseen perustuva Bayeslaisen mallin perusajatus on toki suhteellisen suoraviivainen, mutta käytännössä on vaikeaa intuitiivisesti ymmärtää, miten esimerkiksi erilaiset priorivalinnat vaikuttavat posteriorijakaumiin. Lisäksi, koska käytännössä laskukaavoja ei koskaan voida soveltaa suoraan ja joudutaan turvautumaan laskennallisiin menetelmiin, tulee rakennetun mallin ymmärtämiseen mukaan uusi välivaihe – ohjelmisto, jolla itse laskenta suoritetaan, kuten projektissa käyttämämme BUGS.
Erilaisten mallin parametrien priorijakaumien vaikutus mallien posterioritodennäköisyyksiin BMA-analyysissa aiheutti ryhmässä pohdintaa. Lukemamme mukaan malleja vertailtaessa epäinformatiiviset priorijakaumat johtavat ’hyvin suuriin’ posterioritodennäköisyyksiin yksinkertaisten mallien tapauksessa. Miksi näin on? Millä perusteella eri mallien prioritodennäköisyydet täsmälleen päivittyvät?
Pohdimme myös erilaisten mallien prioritodennäköisyyksiä ja sitä, millä perusteella ne kannattaa määrittää. Kannattaako esimerkiksi aina käyttää tasaista prioria, vai kannattaako mieluummin käyttää omaa asiantuntija-arviota eri mallien sopivuudesta? Jos on jo havaittu tuloksia, voiko näitä huomioida arviossa vai pitääkö ne jättää huomioimatta?
Bayeslaisten p-arvojen käsite jäi myös toistaiseksi hieman hämäräksi.
Carlin & Chib – metodin käyttö BMA:n implikoinnissa BUGS – koodiin vaikutti mielestämme jokseenkin intuitiiviselta, mutta herätti silti useita kysymyksiä BUGSin muuttujiin ja funktioihin liittyvistä yksityiskohdista. Ymmärsimme, että eri mallien todennäköisyyksien päivittyminen varmasti tapahtuu samalla tapaa, kuin mallin muidenkin prioritodennäköisyyksien päivittyminen, mutta silti tähän liittyvä uusi syntaksi aiheutti pohdintaa siitä, miksi ja mitä kautta päivittyminen oikeastaan tapahtuu.
model –muuttujan määrittely tuntui esimerkkikoodeissa epäintuitiiviselta ja tuntui tapahtuvan hieman salaa allaolevan rakenteen mukaisesti
x ~ dnorm(mu[model],tau)
mu[1] < – a
mu[2] <- b

Tulkintamme mukaan ylläoleva määrittele model nimisen (1,2) vektorin.

Tämän lisäksi pohdimme BUGSin equals() –funktion toimintaa. Mitä täsmälleen tarkoittaa koodi P <- equals(model,1)

Luentokalvoilla myös sanotaan, että malli-indokaattoria ei yleensä tarvita sellaisille prioreille, jotka ovat malleille yhteisiä. Onko siis olemassa tilanteita, jossa tälle on kuitenkin tarve?

Pohdintaa mallista ja sen kehittämisestä

Mallissamme on vielä nykyisellään joitakin ratkottavia ongelmia mm. pyydystettävyyden suhteen. Toinen pyydys, hieman smolttiruuvista ylävirtaan sijoitettu rysä on mutkistanut mallinnusta huomattavasti. Emme vielä ole keksineet täysin toimivaa ratkaisua, jolla saisimme kerättyä informaation molempien pyydyksien suhteellisesta pyydystävyydestä. Tätä ratkotaan edelleen, mutta mielestäni on silti hedelmällistä pohtia myös mallin kehittämistä edelleen paremmaksi., vaikka teknisiä ratkaisuja en pystykkään oman osaamiseni puolesta juurikaan tarjoamaan.

Kalojen parveutuminen

Erilaisten biologisten prosessien mallinnus on hyvinkin haastavaa, sillä huomioon otettavia muuttujia on niin runsaasti. Esimerkiksi pyydystettävyyden arvioinnissa huomionarvoinen asia on taimenten vaelluspoikasten parveutumiskäyttäytyminen. Taimenen jokipoikasen reviirikäyttäytyminen muuttuu huomattavasti sen smolttiuduttua, eli sen käytyä läpi sitä mereen siirtymiseen valmistavan fysiologisen muutoksen. Tällöin taimen jättää reviirinsä ja yleensä parveutuu muiden vaelluspoikasten kanssa. Tämä muodostaa ongelman silloin, kun käytettävä malli olettaa kalojen olevan jakautuneen tasaisesti ja toimivan riippumattomasti yksilöinä. Silloin malli antaa liian optimistisen kuvan tarkkuudestaan. Kalojen parveutuminen tulisi ottaa huomioon, mikäli estimaattia smolttien määrästä halutaan tarkentaa.

Nykyinen mallimme käyttää pyydystettävyyden mallintamiseen beta-jakaumaa, jossa pyydyksen “onnistumiset” ja “epäonnistumiset” muodostavat jakauman tiheysfunktion, eli määräävät sen muodon. Tämä ei ota huomioon parveutumiskäyttäytymistä ja sen aiheuttamaa dispersiota. Ratkaisuja tähän ongelmaan on viisaampien toimesta kehitetty ja näitä menetelmiä voisi yrittää myös tämän projektin puitteissa soveltaa, mikäli resursseja siihen riittää. Näitä ovat mm. jonkinlaisen beta-binomiaalisen jakauman käyttäminen, joka sisältää erityisen dispersio-parametrin, tai negatiivisen-binomijakauman soveltaminen. Näiden aukeaminen ainakin allekirjoittaneelle tulee kuitenkin olemaan suuren työn takana.

Mitä tapahtuu joen länsihaarassa?

Mielestäni mielenkiintoinen ajatus mallin laajentamisesta olisi kerätä tarkemmin tietoa myös siitä, kuinka suuri osa kaloista matkaa mereen joen padotun länsihaaran kautta. Tieto olisi erityisen arvokasta nyt, kun padon kohtalosta kiistellään. Jatkossa olisi myös hedelmällistä selvittää kuinka suuri osuus taimenen vaelluspoikasista selviää hengissä voimalan turbiineista. Tällöin tarvittaisiin myös padon alajuoksulle jonkinlainen pyydys.

Tällä hetkellä voimma antaa vain jonkinlaisen estimaatin perustuen alun “mandariinikokeeseen”. Tällöin (hyvin epätarkka) arvio olisi, että länsihaaraan päätyy noin 80% mereen matkaavista taimenista. Tällä hetkellä estimaattia ei ole mahdollista tarkentaa, sillä länsihaaran puolella ei ole minkäänlaista havainnointivälinettä.

Mikäli tälläinen mahdollisuus joskus suotaisiin, olisi varmasti tarpeen mallintaa molempia joenhaaroja erikseen. Tällöin haaroille pitäisi muodostaa todennäköisyysjakaumat myös sen suhteen, kumpaan haaraan kala joutuu, pyydystettävyyden lisäksi. Tämä mutkistaisi mallia hieman, mutta tuottaisi varmasti arvokasta ja ajankohtaista lisätietoa.