Onko mandariini kala?

Ruuvin pyydystettävyyden prioria mallinnamme binomijakauman konjugaatilla beta -jakaumalla. Valinta on kätevä laskettavuuden lisäksi sisällöllisesti siksi, että täten syntyvässä betabinomimallissa beta-(priori)jakauman parametrit voidaan ymmärtää niin sanottuina ”pseudohavaintoina”. Tämä tarkoittaa, että esimerkiksi alun perin käyttämämme priori Beta (2,8) muuntuu ajatusleikiksi, jossa olisimme saaneet kiinni ennen taimensmolttien muuttoon perustuvan ”varsinaisen” aineiston keruuta 2 kalaa (onnistumista) 10 kokeesta (eli 8 epäonnistumista).

 
Pseudohavaintopriorimme on perustunut mandariinikokeeseen, jossa on tarkasteltu, kuinka moni eri paikoista jokeen päästetyistä mandariineista on jäänyt ruuviin. Pseudohavaintojen painoarvoa olemme häivyttäneet suhteessa vapautettujen mandariinien todelliseen määrään siten, että lähes kahta sataa mandariinia koskeva aineisto on muunnettu 9 pseudohavainnoksi (1, 8). Tämä kielii mandariinien ja kalojen vertailtavuutta kohtaan tuntemastamme skeptisismistä.

 
Tämä skeptisyys on näkynyt pseudohavaintojen tarkkuudessa, ei oletuksena mandariinikokeen harhaisuudesta kalojen liikkeiden mallina. Olemme nimittäin määrittäneet pseudohavaintopriorimme siten, että onnistumisten (ruuviin jääneet kalat) ja epäonnistumisten suhde vastaa suoraan mandariinihavaintojen suhteita.

 
Pseudohavaintopriorimme sisältää siis piilo-oletuksen, jonka perusteella kalat ovat enemmän tai vähemmän kuin mandariinit. Eritoten oletamme, etteivät kalat kykene tietoisesti väistämään rysää tai ruuvia. Nimittäin jos kykenevät, tulisi yleistettäessä mandariinit kaloihin mandariinikokeen todellisten onnistumisten määrää pienentää suhteessa havaittuihin epäonnistumisten määrään.

 

Kuinka paljon priori sitten vaikuttaa tekemiimme päätelmiin?
Tarkastellaan asiaa kahdella heuristisella tavalla (tähänhän olisi käytössä mm. informaatioteoreettisia tarkkoja mittoja): vertailemalla ruuvin pyydystettävyyden posterioritodennäköisyyden moodia suurimman uskottavuuden estimaattiin sekä posterioritodennäköisyyteen, joka saataisiin, jos mallin rakentamiseen käytettäisiin informatiivisen sijaan maksimaalisen epäinformatiivista Jeffreysin prioria.

 
Suurimman uskottavuuden estimaatti määräytyy aineiston perusteella. SU -estimoinnissa selvitetään uskottavuusfunktion logiikan mukaisesti, millä parametrin arvolla havaittu aineisto olisi todennäköisin. Jos tässä ”aineisto” ymmärrettäisiin kalojen uudelleenpyyntidatana (jotta vertailu mandariineihin olisi mahdollisimman puhdas), olisi binomikokeen onnistumistodennäköisyyden suurimman uskottavuuden estimaatti yksinkertaisesti kiinniotettujen kalojen suhde kaikista merkatuista kaloista. Uusimmalla aineistolla tämä suhde on ruuvin kohdalla 2/72=1/36=0.0278.

 
Tällä hetkellä ruuvin pyydystystehon posteriorimme näyttää seuraavalta:

Pyydystettävyyden posteriori
Kuvion punainen viiva näyttää suurimman uskottavuuden estimaatin, vihreä viiva posteriorimoodin ja oranssi  priorimoodin paikan. Posteriorimoodimme on ymmärrettävästi suurimman uskottavuuden estimaattia korkeampi, onhan tätä myös priorimme (1/9). Posterioriestimaattimme asettuu SU -estimaatin ja priorimoodin väliin, jonkin verran lähemmäksi ensin kuin viimeksi mainittua.

 

Käytettäessä Jeffreysin epäinformatiivista prioria Beta (0.5;0.5) sekä ruuville että rysälle, saadaan ruuvin pyydystettävyydelle vain hieman vasemmalle varsinaisesta mallistamme siirtynyt posteriorijakauma.

 

Käyttämämme pyydystettävyysmallin posteriorin keskiarvo on 0.064 (noin 1/16) ja 95 % luottamusväli (0.022; 0.136). Pyydystettävyyden posteriorin keskiarvo Jeffreysin priorilla on puolestaan 0.055 ja 95 % luottamusväli (0.018; 0.129). Jeffreysin priori nostaa aineiston vaikutusvaltaa ja hilaa täten posterioria kohti suurimman uskottavuuden estimaattia (1/36) ja pois priorimoodista (1/9).

 

Onko mandariini siis kala? Vaikuttaisi siltä, että informatiivinen beta -priorimme on ”ristiriidassa” havaitun aineiston kanssa ja mandariini ei tämän perusteella vaikuttaisi olevan kala. Todennäköisyys havaita 2 tai vähemmän kalaa 72 kokeella, jos mandariinipseudohavaintomme 1/9 vastaisi todellista pyydystystehoa, on noin 1 % (nollahypoteesin testauksen hirviölogiikalla).
Huomionarvoista on, ettei ero pyydystettävyyden posteriorissa ole kuitenkaan kovin suuri, vaikka siirryttäisiinkin käyttämään epäinformatiivista Jeffreysin prioria (posteriorin keskiarvo siirtyy 1/16 -> 1/18).

 

Yksi tapa sisällyttää mandariinihavaintojen kalojen priorina käyttämiseen liittyvät epäilyksemme olisi rakentaa hierarkinen malli, jossa myös betapriorin parametrimme noudattaisivat omaa priorijakaumaansa. Tämä olisi oletettavasti perustellumpi ratkaisu haasteeseen kuin Jeffreysin priori. Onhan niin, että mandariinit välittävät jotain tietoa kalojen käyttäytymisestä, vaikka olisikin luonteeltaan harhaista ja epätarkkaa.

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!

Päätöksen tekoa

Kurssilla ollaan nyt tilanteessa, jossa mallin rakentamisesta on siirrytty päätöksentekoanalyyseihin. Suurimmalta osin ainakin. Rysän ollessa padon edessä, saadaan kerättyä tietoa, kuinka moni taimen joutuu mahdollisesti turbiinien rattaisiin ja sushina padon alla odottavien kalojen suihin. Tai näin ainakin oli ajatus kurssilaisten mielissä.

 

Rysä otettiin pois käytöstä jo muutamia päiviä sitten, eikä dataa saatu kerättyä juuri yhtään. Rysä kerkesi olla padon edessä vain kahdeksan päivän ajan. Tänä aikana n. 25 % pyydetyistä kaloista päätyi padon edessä olevaan rysään. Mikäli luku olisi sama myös todellisuudessa, olisi merelle vaeltavilla taimenilla suuri mahdollisuus päättää vaelluksensa tähän esteeseen.

 

Padon purkamisen päätöksentekoanalyysia rakentaessamme törmäsin itse monesti ongelmaan, jossa osa osioista tuli ehkä ylimietittyäkin, tämän johtaen hyvin yksityiskohtaiseen ongelman tarkasteluun, osan taas jäädessä hyvin yleiselle tasolle. Jotta päätöksenteko ”kartta” olisi toimiva, pitäisi kaikki kohteet saada suunnilleen samalle tarkastelutasolle.?!

 

Toinen ongelma muodostui hyödyn (utility) ajattelusta. Analyysissä pitäisi pystyä arvottamaan kartan eri osiot erilaisilla hyödyillä tai vaihtoehtoisesti haitoilla. Jotta tämä onnistuisi, pitäisi mallissa pystyä jotenkin valitsemaan vaalikonemaisesti näkeekö asiat hyötyinä vai haittoina ja antamaan näille, erilaisia painoarvoja.

Tulokset 21.5.2015

Päivä: 42

Pyydettyjä taimenen smoltteja kaiken kaikkiaan: 70

Merkattuja taimenen smoltteja kaiken kaikkiaan: 70

Uudelleen pyydettyjä taimenen smoltteja kaiken kaikkiaan: 2

Vähän pienempien ja suurempien teknisten ongelmien jälkeen valmistui päivän 42 (21.5.) tulokset. Tässä mallissa olemme päivittäneen pyydystettävyyden prioria, koska pyydystettävyyttä testaava mandariinikoe uusittiin. Kuten myös edellisessä tulospostauksessa, kaikki seuraavat tulokset perustuvat yli neljän mallimme keskiarvostettuihin posterioriestimaatteihin.

Malli tuottaa lähtemisintensiteetille totaaliestimaatin (Kuva 1), jossa 95 % todennäköisyydellä vaeltavien taimensmolttien kokonaismäärä on välillä 483–5035 (Taulukko 1). Todennäköisin yksittäinen arvio kokonaismäärästä on 500 ja 1500 välissä.

N26052015

Kuva 1. Posteriorijakauma vaeltavien taimensmolttien kokonaismäärästä 2015.

Perjantaina 22.5. vaeltavien taimensmolttien lukumäärä (n[43]) on 95% todennäköisyydellä välillä 3–81 (Taulukko 1). Todennäköisin yksittäinen arvio vaeltavien smolttien lukumäärästä on 9 (Kuva 2).

pikkun26052015

Kuva 2. Posteriorijakauma vaeltavien taimensmolttien määrästä perjantaina 22.5.2015.

Perjantaina 22.5. pyydyksiin jäävien taimensmolttien lukumäärä (x[26]) on 95% todennäköisyydellä välillä 0–5 (Taulukko 1). Jälleen kerran todennäköisin yksittäinen arvio on, että pyydyksiin ei jää yhtään taimensmolttia (Kuva 3).

Edellisellä kurssikerralla puhuimme toisen pyydyksen (rysän) uuden sijainnin huomioimisesta mallissa. Aihetta on käsitelty myös aiemmin blogissa, mutta ehdotusta ei ole huomioitu vielä tässä mallissa.

TotalX26052015

 

Kuva 3. Posteriorijakauma pyydyksiin jäävien taimensmolttien määrästä perjantaina 22.5.2015.

 

Taulukko 1. Parametrien N, n[26] ja TotalX[26] posteriorijakaumien tunnusluvut.

mean sd MC_error val2.5pc median val97.5pc start sample
N 1985.0 1212.0 35.16 483.0 1694.0 5035.0 51 164900
n[43] 25.76 20.81 0.5681 3.0 20.0 81.0 51 164900
TotalX[43] 1.363 1.66 0.02455 0.0 1.0 5.0 51 164900

 

Miksei tasajakauma?

Minkä muotoisen käppyrän, suoran tai niin edespäin sovittaisit seuraavaan aineistoon?

Smoltit_kiinni

Pisteet kuvaavat päivittäisiä kiinnijääneitä taimenen smoltteja päivään 38 asti.

Mielestäni prima facie ja äkkiseltään kohtalaiselta sovitteelta vaikuttaisi horisontaalinen suora, joka kulkisi noin yhdessä päivittäisessä havainnossa yli kaikkien seuranta-ajanjakson päivien.

Mitään trendiä on nimittäin ainakin minun mielestäni kuviosta kovin vaikea havaita. Herkällä tulkinnalla voitaisiin kenties todeta havaintojen hieman lisääntyneen päivien kuluessa. Tällöin horisontaalisen sijaan varovaisesti nouseva suora saattaisi toimia.

Entä onko kuviossa useamman havainnon muodostamia havaintohuippuja (useamman siksi, että näissä olisi tulkittavissa jotain systematiikkaa)? Edes tällaisia on mielestäni vaikea päätellä kuviossa esiintyvän: ehkä korkeintaan päivinä 20, 21 ja 22, kenties jälleen päivinä 36,37 ja 38.

Onko näissä päivissä jotain muuta erityistä, kuin mitä mallien teoriapohjalta voitaisiin ennustaa? No, rysä on ollut yhdessä ruuvin kanssa toiminnassa kaiken kaikkiaan päivinä 16-21 ja 35-38. Kun tämän huomioi, vaikuttaisi hajontakuvioon sopivan erityisen hyvin mainittu horisontaalinen suora, josta muutamat rysän vaikutuksen poistamisen jälkeen jäljelle jäävät poikkeavat havainnot (outlierit), ovat suurusluokkaa 3-4, eivätkä siksi kovin epätodennäköisiä ehdolla horistontaalinen sovite.

Tällainen “horisontaalinen sovite” sopisi yhteen muuttointensiteetin mallin tasajakauma kanssa. Havaintommehan oletetaan tulevan (muuten paitsi pyydysten määrältä) yli päivien vakiolla pyydystettävyysteholla, tietyn mutta  tuntemattoman kokoisesta smolttipopulaatiosta, yli päivien mallinnetuilla muuttointensiteeteillä. Siispä ennustettu muuttointensiteetti määrittää yhdessä pyydysten määrän kanssa suurimman osan siitä, minkä “muodon” oletamme hajontakuviossa esiintyvän.

Mielenkiintoisen poikkeuksen tähän tekee totaaliestimaattimme. Jos havaintomäärät “raahaavat perässä” totaaliprioria, tulee mallin kyetä ennustamaan havaintopiikkejä tai systemaattisesti lisääntyviä havaintoja jäljellä oleville seuranta-ajanjakson päiville. Tasajakauma ei tähän juuri kykene.

Mutta selittääkö tämä sen, että tasajakauma saa mallivertailussamme (https://blogs.helsinki.fi/taimenlaskenta/?p=221) posterioritodennäköisyydekseen 0? Lähes varmoina tulosten perusteella pidettävät (joko tai) normaalijakauma ja lognormaalijakauma eivät vaikuttaisi sopivan aineistoon yhtään tasajakaumaa paremmin. Tekisi mieli sanoa päinvastoin.

No, asia ei ole näin yksinkertainen. Mallia ei tule suosia vain datan perusteella, minkä vuoksi postauksen aloittama kysymys on ikään kuin nurinkurinen. Pitäisi pohtia, mikä priorikäsityksen mukaan perusteltu malli on myös perusteltu ehdolla aineisto. Mallien posterioreihin vaikuttavatkin niiden uskottavuuden lisäksi sekä mallien parametrien määrä, parametrien priorit että mallien priorit.

Koska viimeksi mainitut ovat meillä kategorisesti tasan (pun intended) eli 1/4, eivät erot mallien posterioreissa voi tästä johtua. Siispä ne vaikuttaisivat johtuvan ennen kaikkea mallien parametreista ja näiden prioreista.

Muissa malleissa on enemmän parametreja kuin tasajakaumassa, jossa muuttointensiteetin määrittää 1/60 jokaisena päivänä. Useammalla parametrilla voidaan paremmin myötäillä saatua dataa. Sopivilla parametrien prioreilla voitaneen mallista saada tasajakaumaa selvästi joustavampi ja vielä kun käy niin, että havaintopiikki osuu ennalta määrättyyn priorikohtaan normaalijakaumaa, voidaan juuri kyseinen piikki selittää erinomaisesti jakauman joustavalla muodolla. Mutta mitä tapahtuu, jos ja kun “piikkejä” tulee useampia ja ne ovat vain yhden päivän mittaisia?

Ja eikö parametrien määrästä pitäisi rankaista? Eikö hyvä malli ole SEKÄ sopiva aineistoon (goodness of fit) että yksinkertainen (parsimonious)? Tällöin tasajakauma tuntuisi edellisen kuvion ja yksinkertaisuudensa vuoksi vähintään intuitiivisesti hyvältä valinnalta.

Mutta. Teoriapohjalta varmasti tiedetään, ettei tasajakauma voi olla oikea malli. Smoltit eivät muuta tasaista virtaa. Jos ja kun näin on, tulisi tämän näkyä mallivalintojen prioreissa ja vaikuttaa pääasiassa sitä eikä muuta kautta malliemme posterioritodennäköisyyksiin.

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ä.

 

Mallin päivitys 20.5.2015

Saimme uutta prioridataa!:) Tämä kummajainen pääsi tapahtumaan, kun ruuvin prioripyydystettävyyden määräämä mandariinikoe toistettiin. Alkuperäisessä kokeessa olimme saaneet tulokset, joiden perusteella 70 mandariinia päätyi pato- ja 30 ruuviuomaan. Ruuviuomaan päätyneistä 20 päätyi ruuviin.

Nyt saimme aineistoa sekä yläuomasta että ainoastaan ruuviuomaan päästetyistä mandariineista. Kuinka yhdistää saatu aineisto?

Olkoon P(ru) mandariinin todennäköisyys päätyä ruuviuomaan, P(pu) todennäköisyys päätyä patouomaan ja P(r) todennäköisyys päätyä ruuviin. Tällöin kokonaistodennäköisyyden kaavasta P(r)=P(r|ru)*P(ru)+P(r|pu)*P(pu).

(Viattomalla apuoletuksella, että mandariini tosiaan päätyy jompaankumpaan uomaan, eikä jää matkan varrelle jonnekin jumiin tai tule syödyksi. :))

Koska todennäköisyys, että mandariini päätyy ruuviuomassa sijaitsevaan ruuviin ehdolla, että mandariini päätyy patouomaan, on toki 0, eli P(r|pu)=0, niin tästä seuraa, että P(r)=P(r|ru)*P(ru). Mandariinin todennäköisyys päätyä ruuviin on siis sama kuin mandariinin todennäköisyys päätyä ruuviuomaan kertaa tämän todennäköisyys päätyä ruuviin ehdolla, että on päätynyt ensin ruuviuomaan. Tämä on tietenkin sama kuin todennäköisyys päätyä ruuviuomaan ja ruuviin, P(r ja ru).

Kun huomioidaan kaikista mandariinikokeista saatu data, ovat parhaat yksittäiset pistearviomme P(ru)=49/188, ja P(r|ru)=30/66. Tällöin P(r)=0,118.

Tässä ratkaisussa arveluttaa, ovatko eri kohtia yläjuoksua päästettyjen mandariinikokeiden tulokset vertailukelpoisia? Tämän edellä koetuloksia yhdistettäessä implisiittisesti tehdyn oletuksen paikkansapitävyys vaikuttaisi epätodennäköiseltä siksi, että alkuperäisessä kokeessa ruuviuomaan päätyneiden osuus oli 30/100, kun taas 1.5.2015 neljännen tien sillalta päästetyistä 50 mandariinista KAIKKI päätyivät patouomaan. Todennäköisyys, että 1.5. tehdyssä binomikokeessa saataisiin 50 toistossa yhteensä 50 epäonnistumista ehdolla, että todennäköisyys yksittäisen toiston onnistumiselle olisi ensimmäisen kokeen 0.3 on häviävän pieni (0.7^50). Oletettavasti kokeet eivät siis ole vertailukelpoisia. Koska emme kuitenkaan tiedä, mikä kokeista vastaa parhaiten kalojen käyttäytymistä, oletamme kaikki koetulokset kokeiden erilaisesta luonteesta huolimatta yhtä tärkeiksi emmekä painota tuloksia yhdistettäessä yhtä koetta toisen yli. (Jos 1.5.2015 tehty koe vastaisi parhaiten kalojen liikkeitä, olisi pyydystettävyyden arviomme selvästi esitettyä alhaisempi ja totaaliestimaattimme vastaavasti korkeampi. Tämä selittäisi osin myös Miken blogissa aprikoimat pienet havaintomäärät.)

Tällöin saadut havainnot vastaavat pseudohavaintoina noin suurin piirtein beta-jakaumaa, jossa onnistumisten ja epäonnistumisten suhde on 1/8.

Emme kuitenkaan missään nimessä halua käyttää mandariinikoeaineistoa siten, että yhteenlasketut havainnot vastaisivat todellisten laskettujen mandariinien määrää. Mandariinit eivät ole kaloja, ja kaloilla kertynyt tieto pyydystettävyydestä on mandariineja tärkeämpää. Jos käyttäisimme kaikkia noin 200 mandariinipseudohavaintoa, olisi (pseudohavainto)priori aivan liian vaikutusvaltainen (varsinaisiin kala)havaintoihin nähden.

Siispä päätimme käyttää samankaltaista, joskin hieman vahvempaa mandariiniaineiston merkitystä suhteessa kalahavaintoihin häivyttävää painoa kuin aikaisemmin (noin 1/20). Täten päivitimme beta-priorin pyydystettävyydelle vastaamaan 1 onnistumista ja 8 epäonnistumista, eli Beta(1,8)

Missä kalat on?

Pyydettyjä smoltteja on tällä hetkellä vain 60 kappaletta vaikka oletettu vaellus mereen piti asiantuntija tiedon mukaan käynnistyä vedenlämpötilan saavuttaessa 8 astetta. Ympäristö.fi sivuilta katsottaessa vesi on tällä hetkellä Pirttinrannan mittauspisteessä jo 11 asteessa. Lämpötilan nousu olisi siis pitänyt jo näkyä pyydyksissä, suurimpina saalismäärinä. Tällä hetkellä mitään piikkiä ei ole havaittavissa, vaan pyydettyjen smolttien määrä vaihtelee satunnaisesti 1 – 6 smoltin välillä, vaikka pyydysten määrä on tuplattu. Alkujaanhan pyydyksenä toimi vain joensivuhaarassa oleva ruuvi. Nykyään käytössä on myös rysä. Toisen pyydyksen lisääminen luulisi jotenkin näkyvän pyydettyjen smolttien määrässä.

Vanhassa mallissamme päivitimme pyydysten ”tehoa” sitä mukaan kun kaloja saatiin pyydyksiin. Pyydetty kala oli ensin failure ja sitten succes mikäli merkattu kala saatiin pyydettyä uudestaan. Priorit saatiin malliin Vesi- ja Kalatutkimus Oy:n mandariinikokeesta. Mandariini kokeessa mandariineista pyydykseen päätyi 20 %. Jos mandariineja ajateltaisiin merkattuina smoltteina, olisi uudelleenpyydettävyys riittävällä tasolla, jotta voitaisiin arvioida smolttien kokonaismäärää luotettavasti. Tällä hetkellä uudelleen pyydettyjä kaloja on vasta 1, joten kahdella pyydyksellä päästään vain 1,7 % uudelleenpyydettävyyteen. Kurssilla oli tosin puhetta, että smoltit saattavat kuitenkin tulla pyydykseen uudestaan vasta 1-2 viikon päästä siitä, kun ne on edellisen kerran saatu pyydettyä. Pyydykset ovat nyt kuitenkin pyytäneet jo 40 päivän ajan, joten viiveen ei enää pitäisi vaikuttaa siihen, että merkattuja kaloja ei onnistuta pyytämään takaisin. Tällä hetkellä vanhan mallin pyydysteho olisi niin alhaalla, että en tiedä olisiko sillä voinut antaa luotettavaa arvioita vaelluskannan koosta.

Itseäni on myös ihmetyttänyt pyydysten sijoittaminen. Jo kokeen alkaessa puhuttiin, että ruuvi ei toimi kunnolla, koska virtaus on niin pieni, että ruuvi ei pyöri kunnolla. Sitten rysä asetettiin aluksi joen haaraumaan. Ymmärtääkseni pyydykset toimivat sitä paremmin mitä kovempi virtaus on. Tällöin smolttien on vaikeampi välttää ”omituista kapistusta”. Tästä syystä pyydys tulisi sijoittaa joen ”päävirran” kohdalle. Toki pyydysten sijoittamiseen on vaikuttanut moni asia, ja juurikin kovasta virrasta johtuen sitä, ei ehkä ole voitu sijoittaa juuri sinne minne se olisi haluttu.

Uudessa mallissamme otimme mallia kurssinvetäjämme pyydysmallista ja siinäkin pyydysten teho päivittyy ymmärtääkseni samaan tapaan, kuin vanhassakin mallissa, näin ainakin teoriassa. Koodia katsoessani, ei jotkut kohdat tästä kuitenkaan aukea minulle kuten ei myöskään ehkä muillekaan ja toivoisinkin, että seuraavalla kerralla palattaisiin vielä tähän, vaikka olemmekin siirtyneet jo päätöksentekoanalyysin maailmaan. Tässä kohta, mikä on vielä epäselvä tagged[d]~dbin(1,1000)}. Mistä nämä numerot tulevat?

Tulokset 15.5.2015

Päivä: 36

Pyydettyjä taimenen smoltteja kaiken kaikkiaan: 60

Merkattuja taimenen smoltteja kaiken kaikkiaan: 60

Uudelleen pyydettyjä taimenen smoltteja kaiken kaikkiaan: 1

Viimeisen tulospostauksen jälkeen olemme lisänneet neljännen mallin smolttien muuttojakaumaa estimoimaan. Tässä mallissa oletetaan muuton noudattavan seuranta-ajanjakson päivien yli Beta –jakaumaa. Ajatuksena on, että koska Beta -jakauma on ”joustava”, kykenee se riippuen parametriensa a ja b arvoista mallintamaan hyvin erilaisia muuttointensiteettien muotoja aina tasajakaumasta normaalijakauman  kautta vinoon jakaumaan. Näin ollen toivomme tämän huomioivan muuttointensiteetin muotoon liittyvän epävarmuuden entistä paremmin. Beta –jakauman määrittelyväli on [0,1], joten olemme tehneet yksinkertaisen muuttujamuunnoksen sovittaaksemme aineiston väliltä [1,60] sopimaan jakauman määrittelyvälille.

Nyt käytössämme on siis neljä vaihtoehtoista muuttointensiteettiä mallintavaa jakaumaa: tasajakauma, normaalijakauma, log-normaali- sekä beta -jakauma. Vielä merkittävämpi muutos on kuitenkin usean mallin keskiarvoistaminen. Tällä muutoksella kykenemme kontrolloimaan muuttointensiteetin mallivalintaan liittyvää epävarmuutta. Päivitetystä mallista sekä mallikeskiarvoistamisen periaatteista löydät lisätietoa täältä: https://blogs.helsinki.fi/taimenlaskenta/?cat=75422

Hieman diagnostiikkaa                                                                                                

Malli on melko raskas. 100 000 Gibbsin otantamenetelmällä kerätyn otoksen kokoamiseen meni lähes 3 tuntia. Poimittua suurempi otos saattaisi jatkossa olla perusteltu, sillä vielä 50 000 otoksen kohdalla vaihtoehtoisten MCMC -ketjujen vertailuun perustuva bgr -ketju ei ole täysin konvergoitunut tavoitearvoonsa yhteen:

 

Bgr

 

Ketjut poukkoilevat historiansa perusteella ainoa lähes 20 000 havaintoon asti. Täten on syytä polttaa kyseiset ketjun alun ensimmäistä 20 000 otosta.

Lisäksi havaitaan, että autokorrelaatio ketjujen välillä on korkea:

autokorrelaatio

MCMC -ketjun peräkkäiset otokset ovat toisistaan riippuvaisia. Liiallisesta autokorrelaatiosta pääsee useimmiten eroon, kun valitsee poimituista otoksista vain esimerkiksi joka kymmenennen. Näin saadaan autokorrelaatio seuraavanlaiseksi:

Aukorrelaatio2

Polttamalla alusta 20 000 ensimmäistä havaintoa ja valikoimalla vain joka kymmenen otoksen saadaan numeerisesti estimoidut posteriorijakauman arvot paremmin vastaamaan todellisen posteriorijakauman arvoja.

Tulokset

Kaikki seuraavat tulokset perustuvat yli neljän mallimme keskarvoistettuihin posterioriestimaatteihin.

Saadaan seuraava totaaliestimaatti muuttavien meritaimenten kokonaismäärälle seuranta-ajanjakson aikana.

Totaalin posterioriestimaatti

Jakauma on muodoltaan epämääräisen rosoinen. Epävarmuus on myös suuri. 95% prosentin todennäköisyydellä muuttavien taimenten smolttien kokonaismäärä on välillä [550, 4700] (Valinta vastaa bayesilaista luottamusväliä,  joka on  poimittu posteriorijakauman kvantiilien väliltä [0.025, 0.975].) Totaalin posteriorijakauman odotusarvo on 1960 ja mediaani 1695.

Seuraavan päivän (37) pyydettyjen smolttien kokonaismäärän posterioriennustava jakauma näyttää puolestaan tältä:

Ennustava

Tämä ennustava jakauma EI mielestäni huomioi riittävällä tavalla viimeisen parin päivän havainto”trendiä”, jonka perusteella 5 tai yli havaintoa ei tulisi päivänä 37 pitää aivan niin epätodennäköisenä, kuin miltä se koko seuranta-ajanjakson yli arvioitaessa vaikuttaisi.

Käytännössä kaikki havaintoja ennustavan mallin todennäköisyysmassa on välillä nollasta kymmeneen. Posterioriennustavan jakauman keskiarvo on 2,5 ja odotusarvo 2.

Jos oletetaan, että joku malleista 1,2,3 tai 4 on tosi ja että kaikkien mallien prioritodennäköisyys on ¼, muuttointensiteettiä mallintavien vaihtoehtoisten mallien posterioritodennäköisyydet ovat seuraavat:

Malli 1 tasajakauma: 0.0%

Malli 2 normaalijakauma: 73.5%

Malli 3 log-normaalijakauma: 25%

Malli 4 Beta -jakauma: 1.5%

Tämän perusteella vaikuttaisi siltä, että mallit 1 ja 4 eivät vastaa havaintoja, kun taas malli 2 saa paljon tukea aineistolta ja malli 4 jonkun verran.

Seuraavassa tulospostauksessa lisää vaihtoehtoisten mallien arviointia mm. näiden ennustamien muuttointensiteettien jakaumia tarkastelemalla.

 

Mallin hiomista: ympäristötekijöiden huomioiminen pyydystettävyydessä

Laskentamallimme Vantaanjoen taimensmolttien kokonaismäärälle on saavuttanut pisteen, jossa mallin kehittäminen on päätetty jättää vähemmälle huomioille. Nykyisessä mallissa yhdistyy neljä mallia, joissa smolttien muuttointensiteettiä kuvataan neljällä erilaisella todennäköisyysjakaumalla. Malli on hyvässä tilassa ja uskomme sen tarjoavan relevanttia informaatiota Vantaanjoen smolttien kokonaismäärästä. Parannettavaa ja viilailtavaa kuitenkin tottakai riittäisi useallakin mallin osa-alueella. Tässä postauksessa keskityn pyydystettävyyteen.

Mallin pyydystettävyyttä vastaava BUGS-koodi on alla. Pyydyksiä on ollut käytössä yhteensä kaksi; rysä ja ruuvi, joista ruuvi on ollut toiminnassa koko ajan ja rysä muutaman päivän.

# Catch
for(d in 1:days+1){

# Daily catch and total catchability
TotalX[d]~dbin(Total_q[d],n[d])
Total_q[d]<-Screw_q[d]+Fyke_q[d]-Screw_q[d]*Fyke_q[d]

Yllä määritellään kokonaispyydystettävyys, joka vaikuttaa päivittäisen saaliin todennäköisyysjakaumaan. Alla määritellään ruuvin suhteellinen pyydystettävyys, jonka avulla malliin välittyy tietoa pyydysten välisistä eroista.

ScrewX[d]~dpois(screw_muX[d])
screw_muX[d] <- Screwprop[d]*TotalX[d]
Screwprop[d]<-Screw_q[d]/Total_q[d]

Screw_q[d]<-q_S*Screw_op[d]
Fyke_q[d]<-q_F*Fyke_op[d]

# prior for the number of tagged fish
tagged[d]~dbin(1,1000)
}

Uudelleenpyydystettyjen lukumäärää mallinnetaan erikseen, mikä päivittää pyydysten pyydystettävyyden prioreja.

# Catchability update
for(d in 2:days+1){

FykeR[d]~dbin(Fyke_q[d],tagged[d-1])
Screw_nR[d]<-tagged[d-1]-FykeR[d]
ScrewR[d]~dbin(Screw_q[d],Screw_nR[d])

}

# prior for catchability
q_F~dbeta(1,4)
q_S~dbeta(2,8)

 

Huomioitavaa pyydystettävyyden määrittelemisessä on se, että malli olettaa pyydystettävyyden pysyvän vakiona yli vaelluspäivien. Kuitenkin käsityksemme mukaan ruuvin pyydystettävyys kyseenalaistettiin joinakin päivinä sääolosuhteiden johdosta. Tarkemmin: joen pinnan laskeminen esti ruuvin pyörimisen, jonka arvioitiin vaikuttavan pyydystettävyyteen. Blogissaan Henri pohti jo aikaisemmin hieman sitä, miten pyydystettävyys voitaisiin määritellä päiväkohtaisesti sen sijaan, että käytettäisiin ylläolevaa päivistä riippumatonta määrittelyä.

Yksi tapa lähestyä asiaa olisi muodostaa pyydystystodennäköisyyttä selittävä regressiomalli, jossa pyydystettävyyttä selitettäisiin ympäristötekijöillä. Regressiomallia voidaan yleisesti pitää mallina havaintojen odotusarvoille. Pyydystettävyyden tapauksessa data on bernoulli-jakautunutta – kala joko jäi pyydykseen, tai ei jäänyt. Nämä havainnot koodataan yleisesti 1: onnistuminen, 0: epäonnistuminen, jolloin havaintoja vastaavat odotusarvot ovat yhtä kuin onnistumistodennäköisyys.

Koska todennäköisyysmitta on määritelty välillä [0,1], ei lineaarista regressiomallia voida suoraan soveltaa binomimuotoisen datan odotusarvon mallintamiseen, sillä tavallisessa lineaarisessa regressiomallissa kohdemuuttujan arvojoukoksi oletetaan (-ääretön, ääretön). Tämä johtuu siitä, että kohdemuuttujan odotusarvo oletetaan normaalijakaantuneeksi.

Asia voidaan kuitenkin hoitaa sopivalla linkkifunktiolla. Yleisesti käytetty linkkifunktio on logit-funktio, jossa todennäköisyys q muunnetaan sitä vastaavan odds-muuttujan logaritmiksi. Odds tarkoittaa onnistumistodennäköisyyden suhdetta epäonnistumistodennäköisyyteen eli odds = q / 1-q. Tämän logaritmilla on se toivottu ominaisuus, että se saa arvoja negatiivisesta äärettömästä positiiviseen äärettömään, kun q lähestyy nollaa ja yhtä. Tällöin kohdemuuttuja (tai regressiomallin virhetermi) voidaan olettaa normaalijakaantuneeksi.

Muotoillen Mäntyniemen ja Romakkaniemen (2002) merkintöjä, voitaisiin pävittäistä pyydystystodennäköisyyttä (q_j) päivänä j mallintaa seuraavasti:

log( q_j / 1 – q_j ) | μ_j, ξ ~N( μ_j , ξ )

μ_j = α0 + α1WT_j + α2WL_j

WT = veden lämpötila

WL = veden korkeus

Mallin mukaan logit muunnos pyydystettävyyksistä q_j, noudattaa normaalijakaumaa, jonka odotusarvo on ympäristötekijöiden lineaarikombinaatio μj ja varianssi ξ. Ylläolevassa mallissa parametreille α0, α1, α2 ja ξ tulee asettaa priorijakaumat, minkä jälkeen malli on määritelty ja selittää päivittäistä pyydystystodennäköisyyttä veden lämpötilalla ja korkeudella.