Laskentamalli goes live?

BUGS ei ole ainoa tapaa harjoittaa markov chain monte carlo simulointia, johon bugsin numeerinen laskenta perustuu. R-ohjelmistosta lötyy esimerkiksi paketti rcppbugs, jonka avulla voidaan mcmc simulointi suorittaa pelkästään R-ympäristössä. R on myös sen verran monipuolinen ohjelmointikieli, että simulointi voitaisiin myös koodata sen avulla kokonaan alusta alkaen.

Pelkän R-ohjelmiston käytöstä on tiettyjä etuja. R-ohjelmistoon on saatavilla esimerkiksi shiny-lisäosa, joka mahdollistaa helpon web-appien kehityksen. Kalamallista löytyykin nyt toistaiseksi julkaisematon web-appina toteutettu versio, jonka voi nyt ajaa R-ohjelmiston kautta.

Tämän version laskentatulokset ovat toistaiseksi epäluotettavia ja simuloinnin koodaus saatetaan vielä uudistaa myöhemmin. Web-appi kuitenkin mahdollistaa esimerkiksi helpon priorijakaumien kuvaajien kanssa leikkimisen.

Projektin koodit löytyvät GitHubista ja “Smolt” nimisen appin voi käynnistää R-ohjelmiston kautta seuraavasti:

if(!require(“shiny”)) {
install.packages(“shiny”)
library(shiny)
}

runGitHub(“Smolt”,”hasiOne”)

 

 

Logistinen regressiomalli pyydystettävyydelle

Tämä postaus on jatkoa muutamalle edelliselle postaukselleni, joissa käsittelin pyydystettävyyden mallintamista logistisella regressiomallilla. Ideana ja motivaationa tälle oli se, että käytetyn pyydyksen tehon epäiltiin riippuvan erityisesti veden korkeudesta.

Käsittelin logistista regressiomalli teoreettisesti ja kirjoitin hahmotelmaa siitä, miten malli voitaisiin määritellä BUGS-kielellä. Kuitenkin epäselväksi jäi se, miten prioritieto pyydystettävyydestä saataisiin mukaan malliin. Tämä postaus käsittelee tätä ongelmaa.

Logistinen regressiomalli on muotoa

logit(q) = a + beta[1]*VK +beta[2]*VL + E

VK=veden korkeus, VL=veden lämpötila, E=virhetermi , a= vakiotermi, beta[1] = veden korkeutta vastaava regressiokerroin, beta[2] = veden lämpötilaa vastaava regressiokerroin.

Haluaisimme siis ‘valita’ tämän mallin parametrit niin, että pyydystettävyyden priorijakauma olisi edelleen suurinpiirtein tietty beta-jakauma, kuten aiemminkin. Nyt, koska mallissa pyydystettävyys vaihtelee päivittäin, tarkoittaa tämä oikeastaan sitä, että arvioituna kaikkien päivien yli (tai oikeammin, arvioituna mallin selittäjien jakauman yli), noudattaisi pyydystettävyys haluttua beta-jakaumaa.

Aivan täsmälleen tätä haluttua tulosta emme voi saavuttaa, mutta niin sanotun informatiivisen g-priorin avulla pääsemme hyvin lähelle.

Ratkaistessa todennäköisyys logit-muunnoksesta, saadaan yhtälö

q = exp(x’B) / (1 + exp(x’B) )

jossa X on mallin selittäjät sisältävä matriisi (1, VK, VL) ja B näitä vastaavat kertoimet (a, beta[1], beta[2]).

On näytetty, että jos

B ~ N(b*e, g*n*(X’X)^-1)

ja

x ~ N(u, Z)

missä

b = digamma(a) – digamma(b)
g = [trigamma(a) + trigamma(b)] / p
e = c(1,0,..0) <- p-pituinen vektori
p = parametrien lukumäärä
n = havaintojen lukumäärä

N() = normaalijakauma

x = selittävät muuttujat

niin silloin

exp(x’B) / (1 + exp(x’B) )

vastaa suunnilleen beta-jakaumaa parametrein a ja b.

Tämä siis tarkoittaa, että on mahdollista rakentaa logistinen regressiomalli, jonka priorijakauma vastaa haluttua beta-jakaumaa.

Tämän voi todeta esimerkiksi R-ohjelmiston avulla seuraavasti. Alla olevassa koodissa X on niin kutsuttu design-matrix määriteltynä R-objektina (eli esim matriisi tai data frame, jonka ensimmäinen kolumni on 1, ja toinen ja kolmas esimerkiksi päivittäiset veden korkeudet ja lämpötilat). Koodin ajamiseksi tämä matriisi tulee siis ensin määritellä.

library(MASS)

X <- …

get.bg <- function(alpha, beta, parameters = 3) {
b <- digamma(alpha) – digamma(beta)
g <- trigamma(alpha) + trigamma(beta)
g <- g/parameters

return(c(b = b, g = g))
}

bg <- get.bg(alpha = 3, beta = 50)
Sigma <- bg[[‘g’]]*nrow(X)*solve(t(X) %*% X)
mu <- c(1,0,0)*bg[[‘b’]]

# B noudattaa multinormaalijaumaa, tässä otetaan siitä satunnaisotos
B <- mvrnorm(1000, mu = mu, Sigma = Sigma)

#  vektori, jossa todennäköisyyksiä logit-muunnoksesta
p <- apply(B,1,function(bi) {
apply(X,1,function(xi) {
exp(xi%*%bi)/(1+exp(xi%*%bi))
})
})
p <- c(p)

plot(density(p))

 

Kuvasta nähdään, että tuloksena on haluttu beta(a=3, b=50) jakauma

Pyydystettävyyden mallintaminen ja BUGS

Viimeisessä henkilökohtaisessa postauksessani kirjoitin pyydystettävyydestä ja sen mallintamisesta. Taimenruuvin pyydystettävyys on tämänhetkisessä mallissa mallinnettu niin, että sen jakauma on riippumaton päivistä. Toisinsanoen ruuvi pyydystää samalla teholla joka päivä.

Esitin viimeksi tälle toisen teoreettisen mahdollisuden; logistinen regressiomalli, jossa päivittäistä pyydystettävyyttä selitetään veden korkeudella ja lämpötilalla.

logit(q) = a + beta[1]*VK +beta[2]*VL + E

VK=veden korkeus, VL=veden lämpötila, E=virhetermi , a= vakiotermi, beta[1] = veden korkeutta vastaava regressiokerroin, beta[2] = veden lämpötilaa vastaava regressiokerroin.

Korkeuden mahdollinen vaikutus on siinä mielessä selkeä, että veden korkeus vaikuttaa havaittavasti pyydyksen toimintaan. Lämpötila taas vaikuttaa taimenten vireyteen, joka saattaa vaikuttaa niiden kykyyn väistää pyydys.

Nyt esimerkiksi veden korkeutta vastaavan regressiokertoimen (b) tulkita on seuraava: Veden korkeuden noustessa yhden yksikön, muuttuu pydystettävyyden logit-muunnos keskimäärin (b) yksikköä silloin, kun veden lämpötila pysyy vakiona.

On eri asia määrittää tilastollinen malli matemaattisesti ja soveltaa kyseistä mallia. Soveltaminen tapahtuu tilastollisella ohjelmistolla, jolla on oma kielensä. Miten logistinen regressiomalli voitaisiin määritellä BUGS-kielellä? Alla on hahmotelma mallista vaihe vaiheelta.

# (1) Oletamme, että logit-muunnos pyydystystodennäköisyydestä (q), (eli log(q / 1- q) ), noudattaa normaalijakaumaa parametrein nyy ja tao. Indeksi d vittaa päivään.

logit_q[d]~dnorm(nyy[d],tao)

# (2) … jakauman odotusarvo on lineaarikombinaatio nyy = a + beta1*VK + beta2*VL

nyy[d] <- alpha + beta1*VK[d] + beta2*VL[d]

# (3) BUGS:ssa normaalijakauman toisena parametrina on varianssin käänteisluku. Varianssi vastaa mallin virhetermiä E. Keskihajonta on yleensä hieman intuitiivisempi tunnusluku hahmottaa, joten priori voidaan määrittää sille ja muunta se varianssin käänteisluvuksi. Voisimme esim olettaa, että emme tiedä virheen suuruudesta juuri mitään.

tao<-pow(E,-2)
E~dunif(0.01,10000)

# (4) mallinnettu logit-muunnos täytyy vielä muuntaa takaisin todennäköisyydeksi. Ratkaistaan siis q yhtälöstä log(q/1-q) = nyy ja käytetään saatua kaavaa. Alla ensin BUGSissa käytettävä tulos ja sen alla myös itse lasku, joka on toistettavissa lukiomatematiikan avuin.

q[d] <- exp(nyy[d]) / (1+exp(nyy[d]))

log(q/1-q) = nyy | exp()
q/1-q = exp(nyy) | *(1-q)
q = exp(nyy) – q*exp(nyy) | + q*e^nyy
q*(1+exp(nyy)) = exp(nyy) | / (1+e^nyy)
q = exp(nyy)/(1+exp(nyy))

# (5) mallin parametreille a, beta1, beta2 tarvitaan vielä priorijakaumat. Alla mallin hahmotelma kokonaisuudessaan ja sen jälkeen muutama huomio ja kehitysidea.

for(d in 1:days) {

logit_q[d]~dnorm(nyy[d],tao)
nyy[d] <- alpha + beta1*WT[d] + beta2*WL[d]

}

tao<-pow(sd,-2)
sd~ dunif(0.01,10000)
a~?
beta1~?
beta2~?

Parametrin a priorijakauma tulisi muodostaa niin, että se huomioi prioritiedon pyydystettävyydestä. Tässä on vielä pohdittavaa.

Lisäksi pitäisi huomioida se, että vaikkakin veden korkeuden nouseminen parantaa pyydystettävyyttä johonkin korkeuteen asti, on oltava sellainen piste, jonka jälkeen pyydystettävyys alkaakin huononemaan. Nimittäin kun veden tilavuus kasvaa niin kalat menevät todennäköisemmin pyydyksestä ohi.

Mallista pitäisikin luultavasti muodostaa polynominen versio, johon lisättäisiin selittäjäksi veden korkeuden neliö. Tällöin malli pystyy muotoutumaan paraabeliksi ja huomioimaan yllämainitun.

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.

Huomioita tulosten raportoinnista

Tilastollista päättelyä voidaan harrastaa kahdesta hieman filosofisesti toisistaan poikkeavasta näkökulmasta. Tämä blogi keskittyy ns. Bayeslaisen koulukunnan tarjoamiin menetelmiin. Toinen päättelyn koulukunta on ’perinteisempi’ ja suositumpi frekventistinen päättely. Teoreettisesti ero näiden kahden välillä muodostuu pienestä yksityiskohdasta: Bayes-päättelyssä kiinnostuksen kohteena oleviin mallin parametreihin kiinnitetään todennäköisyysjakauma, kun taas frekventisessä päättelyssä ne ajatellaan tuntemattomiksi, mutta kiinteiksi.

 

Tämä yksityiskohta osoittautuu hyvin merkitykselliseksi ja määrittelee oikeastaan kaksi hyvin toisistaan poikkeavaa näkökulmaa tilastolliseen päättelyyn. Frekventistisessä näkökulmassa aineiston tuottaneen tilastollisen mallin parametrit ajatelleen siis kiinteiksi, jolloin kaikki käsiteltävä satunnaisuus liittyy itse aineistoon. Aineisto ajatellaan satunnaiseksi leikkimällä sillä mahdollisuudella, että oltaisiin voitu kerätä/havaita myös toisenlaisia aineistoja, mikäli sama koe oltaisiin toistettu samoissa olosuhteissa.

 

Bayeslaisessa näkökulmassa tilastollisen mallin parametreille kiinnitetään todennäköisyysjakauma, jota päivitetään havaitun aineiston perusteella. Tämä ns. priorijakauma kvantifioi ilmiöön liittyvän enakkotiedon ja siihen liittyvän epävarmuuden. Asetelma on siis päinvastainen kuin frekventistisessä päättelyssä; aineisto on kiinteä ja sen tuottaneen tilastollisen mallin parametrit satunnaisia (epävarmoja).

 

Tyypillisesti tilastotieteilijöille opetetaan lähinnä pelkästään frekventististä päättelyä. Sellainen asetelma, jossa kaikki onkin toisinpäin aiheuttaa hieman totuttelemista. Olemme esimerkiksi tässä blogissa raportoineet laskentamallimme tuloksia osittain epäselvästi johtuen omalta osaltani kokemattomuudesta analysoida bayslaisen laskentamallin tuloksia. Puutun nyt muutamaan virheelliseen ilmaisuun ja selvennän, mitä yritettiin sanoa ja mitä olisi pitänyt sanoa.

 

Kokonaismäärän odotusarvon piste-estimaatti

 

Laskentamallimme tarkoituksena on tuottaa taimenten kokonaismäärän todennäköisyysjakauma. Kun olemme onnistuneet tämän jakauman tuottamaan, on sen odotusarvo meille tunnettu. Käytämme kuitenkin laskennallisia menetelmiä, joten emme aivan täsmälleen pysty tuottamaan haluamaamme jakaumaa, vaan BUGS-ohjelmiston tuottaman estimaatin siitä. On aivan oikein puhua piste-estimaatista, mutta on syytä ymmärtää, että kysymyksessä on estimaatti laskennallisista syistä. Tämä liittyy alempaan kohtaan.

 

Odotusarvon piste-estimaatin 95% luottamusväli

 

Luottamusväli on frekventistisen päättelyn termi. Frekventistisessä päättelyssä 95% luottamusväli tarkoittaa, että muodostettu väli sisältää todellisen parametrin arvon 95 kertaa sadasta, perustuen sellaiseen mielikuvitusleikkiin, jossa kerätään vastaavia aineistoja ja muodostetaan uusia välejä.  Tällä nimikkeellä olemme raportoineet kokonaismäärään liittyviä todennäköisyysvälejä (Näistä käytetään Bayeslaisessa analyysissa joskus termiä uskottavuusväli, engl. credible interval), eli sellaisia välejä, joiden sisällä kokonaismäärä on 95% todennäköisyydellä.

 

Vaikkakin oletettavasti lukijat ovat tienneet mitä näillä luvuilla on tarkoitettu, on niiden kohdalla tarkalleen ottaen tehty paha virhe liittyen edellisen kohdan huomioihin. Kuten todettu, liittyy odotusarvon piste-estimaattiin laskennallista epävarmuutta, joka on tunnettua, sillä laskentaohjelmistot kuten BUGS raportoivat samplaykseen liittyvän epävarmuuden (MCMC error), joka kunkin tunnusluvun kohdalla kuvaa siihen liittyvää laskennallista epävarmuutta. Tämä epävarmuus voitaisiin (ja kuuluisikin) raportoida, jolloin voitaisiin esimerkiksi kertoa piste-estimaattiin liittyvä 95% todennäköisyysväli (uskottavuusväli). Tämä väli on kuitenkin aivan eri asia, kuin mallin antama todennäköisyysväli kokonaismäärälle, sillä piste-estimaatin todennäköisyysvälissä kysymyksessä on parametriin liittyvä laskennallinen epävarmuus, ei mallin epävarmuus.

 

Raportoimamme luvut ovat siis olleet kokonaismäärän 95% todennäköisyysvälejä – välejä, joiden sisällä kokonaismäärä on mallimme mukaan 95% todennäköisyydellä.

 

Bayes-oppia kalojen kanssa

Tässä blogipostauksessa pohdin laskentaprojektiamme yleisesti, siihen liittyvää omaa oppimistani, sekä muutamaa projektin yksityiskohtaa, jotka tässä vaiheessa mietityttävät minua.

Olemme nyt puuhastelleet Vantaanjoen meritaimenkannan arvioinnin kanssa muutaman viikon, ja on selvää, että mielekkään projektin parissa työskentely tarjoaa loistavan tilaisuuden oppia. Olen todella iloinen siitä, että projektia suoritetaan suurilta osin ryhmätyönä ja siinä on mukana paljon eritaustaisia ihmisiä. Tämä tarjoaa jokaiselle hyvän mahdollisuuden oppia toisiltaan. Itselläni on tausta tilastotieteessä ja koen aina hyödylliseksi jos saan tilaisuuden selittää muille sellaisia teoreettisia peruskäsitteitä, jotka minulla pitäisi olla jo hyvin hallussa. Tämä on hyvä tapa varmistaa, että olen todella itse ymmärtänyt kyseiset asiat ja samalla tilaisuus syventää omaa ymmärrystäni niihin. Tieto on paljon hyödyllisempää, jos sen osaa selkeästi välittää eteenpäin, joten kaikki kokemus oman tiedon jakamisesta on minusta arvokasta.

Käytännöllisen projektin kanssa työskentely on myös tarjonnut mahdollisuuden peilata omaa osaamistaan ‘tosielämän tilanteiden’ valossa. Toisaalta keksiessään ratkaisuja projektin osaongelmiin voi tyytyväisenä huomata oppineensa jotain viimeisen kolmen vuoden aikana yliopistolla, mutta toisaalta taas joidenkin asioiden kanssa ymmärtää, että teorian soveltaminen käytäntöön ei ole helppo eikä itsestäänselvä harppaus. Välillä on vaikeaa tunnistaa miten käytännön ongelma muotoillaan matemaattisesti, vaikka ongelma tuntuisikin yksinkertaiselta. Käytännössä tutuilla ongelmilla on aina päällä uudenlaiset vaatteet.

 

eDefined

The definition of e on a T-shirt

 

Koen saaneeni paljon apua projektin biologisen ja tietoteknisen puolen asioihin muilta projektilaisilta. Käyttämämme BUGS -ohjelmisto oli minulle jokseenkin tuttu ennen projektia, mutta nyt koen, että muiden avulla olen saanut hyvän otteen ohjelmasta ja olen kiinnostunut oppimaan lisää mallidiagnostiikasta ja muista ohjelmiston tarjoamista ominaisuuksista. Erityisesti minua mietityttää esimerkiksi se, minkälaisia työkaluja on olemassa eri mallien vertailuun ja miten parhaiten vertailla erilaisten priorien vaikutusta posteriorijakaumiin. Tähän liittyen myös teoria priorijakaumien valinnasta kiinnostaa.

 

JeffreysPrior

Jeffreys prior

 

Yksi haastavimmista osista mallia on ollut kalojen muuttointensiteetin mallintaminen, jota kuvataan mallissa kalojen todennäköisyytenä lähteä liikkeellee tiettynä päivänä. Tämän mallinosan suunnittelu alkoi yksinkertaisesta diskreetistä tasajakaumapriorista, eteni siitä ’normaalin’ näköiseen pistetodennäköisyysfunktioon ja tästä edelleen ’log-normaaliin’ pistetodennäköisyysfunktioon. Pidimme alunperin jälkimmäistä ’parhaana’ ja ’realistisimpana’ priorivalintana, sillä se vaikutti sopivan yhteen vaelluksen dynamiikasta kertovat asiantuntijatiedon kanssa.

 

640px-Lognormal_distribution_PDF.svg

log-Normal distribution

 

Vaikutti kuitenkin siltä, ettei kyseinen priorivalinta saanut mallia käyttäytymään haluamallamme tavalla; lähtöintensiteetin huippu ajoittui ei-haluttuun kohtaan. Kyseinen priorivalinta vaikutti myös huomattavasti tuloksiin; ajattelen tämän johtuvan ainakin osittain siitä, että lähtöintensiteetin jakauma ei päivity mallissa kovinkaan nopeasti, sillä se kuvaa vaelluksen kulkua yli kaikkien vaelluspäivien ja alussa on tietysti dataa vain ensimmäisistä päivistä. Yritän hakea tästä epäonnistumisesta opetusta. Meillä ei ollut kovinkaan paljoa tietoa ilmiöstä ja silti valitsimme melko informatiivisen priorin. Lisäksi tiesimme, ettei data tulisi päivittämään tätä priorijakaumaa kovinkaan tehokkaasti. Ehkäpä voidaan yleisesti sanoa, että tällaisessa tilanteessa tulisi valita epäinformatiivisempi priori?

 

Tuloksia tuplateholla

  • Päivä: 20
  • Smoltteja saatu yhteensä: 21

 

Tässä postauksessa esittelemme lyhyesti uuden mallin antamia tuloksi. Mallin nykytilasta voit lukea tarkemmin täältä: https://blogs.helsinki.fi/taimenlaskenta/?p=32

Tärkeät rakenteelliset muutokset malliin ovat: 1) toisen pyydyksen (taimenrysä) huomioiminen, sekä 2) muuttointensiteetin jakauman uudelleenarviointi. Muuttointensiteettiä, eli todennäköisyyttä muuttaa tiettynä päivänä, mallinnetaan nyt log-normaalin jakauman sijaan tasajakaumalla ja normaalijakaumalla. Syyt muutoksiin löytyvät ylläviitatusta postauksesta.

Lisäksi päivitimme asiantuntijaprioriamme muuttajien kokonaismäärän odotusarvosta. Meillä on nyt käytössä kolme asiantuntija-arviota, jotka ovat ’5000’, ’1940’ ja ’alle 1000’. Optimistisesti tulkitsimme alle 1000 tarkoittamaan 999 ja saimme näin kokonaislukumäärän priorijakauman odotusarvoksi (5000+1980+999) / 3 = 2600 ja keskihajonnaksi sd(5000,1980,999) = 2000. Mallina tälle priorille on edelleen kokonaisluvuiksi pyöristetty normaalijakauma, joka on katkaistu tuottamaan ainoastaan positiivisia arvoja.

Sitten tuloksiin!

 

Malli 1. normaalimalli muuttointensiteetille

 

 

Taulukko 1: mallin 1 antamat estimaatit taimenten kokonaismäärän jakauman tunnusluvuille

mean sd MC error 95% conf
N 1086.0 814.3 18.81 219 – 3307

 

Kuvio 1: kokonaismäärän posteriori-jakauma mallissa 1

  normal_N

 

 

Normaalimalli lähtemisintensiteetille tuottaa taimenten kokonaismäärän odotusarvon piste-estimaatiksi 1086 kappaletta, jolle 95% luottamusväli on 219 – 3307. Luonnollisesti uusi pienempi asiantuntijapriori vetää odotusarvoa alaspäin. Myöskään data ei toistaiseksi tue normaalioletusta lähtemisintensiteetille kovinkaan hyvin, sillä taimenia on pyydystetty tasaisen pieniä määriä, minkä johdosta tasajakaumapriori lähtemistodennäköisyydelle tuottaa (vielä toistaiseksi) suuremman estimaatin kokonaislukumäärälle, kuin normaalipriori. Tästä lisää alempana.

 

Taulukko 2: Kokonaispyydystettävyys mallissa 1.

mean sd MC error 95% conf
q 0.139 0.074 0.001 0.034 – 0.318

 

Kokonaispyydystettävyyden  estimaatti on mallissa 0.139, eli odotamme, että 13.9% kaloista joutuu jompaankumpaan pyydykseen. Toistaiseksi pyydykset eivät tosin ole onnistuneet uudelleenpyydystämään ainuttakaan 17 merkatusta kalasta, joten pyydysten pyydystettävyys päivittyy priorioletuksista jatkuvasti huonompaan suuntaan. Tämän(kin) johdosta kokonaismäärän odotusarvon estimaatti tulee todennäköisesti päivittymään suuremmaksi tulevina päivinä, kun dataa saadaan lisää.

 

Malli 2. Vantaanjoesta muuttavien smolttien totaalin posterioriestimaatti, kun muuttointensiteettien oletetaan noudattavan tasajakaumaa päivien funktiona.

 

Olettaen muuttointensiteetin seuraavan tasajakaumaa saadaan seuraava posteriorijakauma seurantajakson aikana muuttavien smolttien kokonaismääräksi N.

Taulukko 3. Vantaanjoesta muuttavien smolttien kokonaismäärän posteriorijakauma.

Mean Sd MC_error val2.5pc median val97.5pc start sample
N 1346.0 873.5 9.441 364.0 1105.0 3694.0 66301 2827470

 

Kuvio 2. Muuttavien smolttien kokonaismäärän posteriorijakauma

Unif_N

 

Paras yksittäinen posterioriestimaatti liikkuu jossain reilun tuhannen paikkeilla. Mielestäni tämä vaikuttaa uskottavalla yksittäiseltä arviolta. Voidaan kuitenkin todeta, että arvio smolttien kokonaismäärästä on vielä kovin epävarma.  95% todennäköisyydellä totaali olisi mallin perusteella välillä (360; 3700).

Samaan hengenvetoon voidaan kyseenalaistaa, onko posterioriarvio riittävän epävarma? Aineistoa on nimittäin saatu kovin vähän, eikä sen ja prioritiedon valossa uskaltaisi välttämättä mennä väittämään, ettei smoltteja varmasti muuta esimerkiksi yli 6000. Itse en menisi tästä takuuseen: hyvin paljon nimittäin riippuu tehdyistä mallioletuksista. Kenties luotettavin arvio saataisiinkin keskiarvoistamalla posterioriestimaatit yli rakennettujen mallien. Palataan asiaan myöhemmässä blogipostauksessa.

Näkisin mallivalinnassa esiintyvän kaksi erityisen suurta epävarmuustekijää: 1) mitä jakaumaa smolttien muuttointensiteetti (ajan funktiona) noudattaa (ks. viimeisin mallit -osion blogipostaus)? sekä 2) mikä on ruuvin ja rysän pyydystystodennäköisyys?

Molempien suhteen priorioletuksemme päivittyy koko ajan melko vinhaa vauhtia. Ensimmäistä koskien olemme joutuneet jo luopumaan epärealistiseksi todetusta lognormaalijakaumasta: kyseinen priorioletus asetti liian paljon todennäköisyysmassaa ensimmäisten havaintopäivien kohdalle. Jälkimmäistä koskien oletuksemme ovat siirtyneet vahvasti kohti heikompaa pyydystettävyyttä, mikä on varsin perusteltua, sillä vielä emme nimittäin ole saaneet yhtään smolttia uudelleenpyydettyä. Lähtökohtainen joka viidennen smoltin jääminen ruuviin on muuttunut noin 7% odotusarvoksi kunkin pyydyksen pyydystettävyydelle. Kenties tämän suhteen beta-jakauma priorimme oli turhan informatiivinen vastatessaan kymmentä pseudosmolttihavaintoa.

 

Taulukko 4. Posterioriestimaatti yhteenlasketulle (ts. smoltti jää rysään TAI ruuviin) pyydystettävyydelle

Mean Sd MC_error val2.5pc median val97.5pc start sample
q 0.1248 0.0666 5.534E-4 0.03345 0.1122 0.2854 152121 2570010

 

Pyydystettävyyden estimointia tulisi parantaa varsinkin suhteessa rysän ja ruuvin suhteelliseen pyydystehoon. Nythän oletamme rysän ja ruuvin saalistavan samalla teholla: oletus, jonka paikkansapitävyydelle ei ole juuri perusteita. Samanaikaisesti pyydystettävyystehot päivittyvät ainoastaan uudelleenpyydettyjen smolttien aineistolla. Voidaan kuitenkin hyvin olettaa, että kaikki (muutkin kuin merkatut smoltit) kalat tuovat mukanaan informaatiota pyydystettävyydestä jäädessään kuhunkin pyydykseen. Jos esimerkiksi alkaisi vaikuttaa siltä, että yhteen pyydykseen jää kaksi kertaa niin paljon kaloja kuin toiseen, tulisi mallin kyetä päivittämään pyydysten suhteellisia tehoja havainnon mukaisesti. Tämä voitaisiin kenties tehdä laskemalla aineiston mukaan päivitettävät painot pyydyskohtaisille todennäköisyyksille. Näiden painojen keskiarvon tulisi olla 1. Jatkossa tämän näkökulman huomioimista on syytä edelleen kehittää.

 

 

 

Lisäpyydyksen huomioiminen ja lähtemisintensiteetin uudelleenmallintaminen

Muutoksia malliin: uusi pyydystystapa sekä priorioletusten päivittäminen

Pyydystystavan muutoksen huomioiminen mallissa

Päivän 16 kohdalla tapahtui tutkimusasetelman kannalta mielenkiintoinen muutos, kun datan lähteeksi joelle lisättiin Taimenruuvin lisäksi Taimenrysä, joka sijaitsee joen haaraumassa suunnilleen kuvan osoittamalla tavalla.

Kuva 1. Pyydysten sijainnit

Position of the samplers

 

Tämä tietysti vaikuttaa oleellisesti päivittäisen saalismäärän jakaumaan. Koko mallimme rakenne on karkeasti kuvattuna suunnilleen seuraavanlainen:

N & p -> n , n & q -> x

N = lähtevien taimenten kokonaismäärä

p = päivittäinen lähtötodennäköisyys

n = päivittäisten lähtijöiden lukumäärä

q = pyydystystodennäköisyys

x = pyydystettyjen taimenten lukumäärä

 

Eli kokonaislukumäärä ja päivittäinen lähtötodennäköisyys määräävät päivittäisten lähtijöiden lukumäärän jakauman, joka taas yhdessä pyydystystodennäköisyyden kanssa määrittää kiinnijääneiden lukumäärän jakaumaan. Tämä hierarkisen rakenteen ansiosta pyydykseen jääneiden kalojen lukumäärä päivittää esimerkiksi arviotamme kokonaislukumäärästä (N).

Taimenrysän tullessa peliin mukaan, täytyy pyydystystodennäköisyys määrittää uudestaan. Rysä ja ruuvi eivät vaikuta toistensa toimintaan, mutta sama kala ei voi jäädä kiinni molempiin. Siis:

P(kala jää ruuviin TAI rysään) = P(kala jää rysään) + P(kala jää ruuviin) – P(kala jää ruuviin JA rysään)

BUGS-koodilla tämä näyttää seuraavalta. (En käy läpi koodin yksityiskohtia; jos BUGS ei sinua kiinnnosta, siirry alas muuttointensiteettiä käsittelevään kappaleeseen!)

# Total cathability

q <- q_f + q_s – q_f * q_s

 

#separate catchabilities for both catching methods

# a=number of recaptures

# b=number of marked smolts not recaptured on the next day

# q_s=catchability of the smolt screw

q_s~dbeta(alpha,beta)

alpha<-(2 + sum(a[1:m[2]]))

beta<-(8 + sum(b[1:m[2]]))

# q_f=catchability of the fyke, that has been operational from day 15

q_f~dbeta(alpha2,beta2)

alpha2<-(2+(sum(a2[15:m[2]])))

beta2<-(8+(sum(b[15:m[2]])))

 

Nyt kun pyyntitodennäköisyys on uudelleenmääritelty, täytyy vielä mallissa huomioida, että muutos tapahtui tiettynä päivänä. Toteutimme tämän step-funktion avulla seuraavasti:

for(i in 1:m[1]) {

.

.

.

# catchability changed from day 16

lambda[i] <- q_s*n[i]*step(15-i) + q*n[i]*step(i-16)

 

# x[i] = the daily catch (num. of fish caught at day i)

x[i]~dpois(lambda[i])

}

Data näyttää mallissa nyt seuraavalta:

list(m=c(60,20))

x[]  a[]  a2[]      b[]

0    0    0    0

0    0    0    0

0    0    0    0

1     0    0    0

2     0    0    0

2     0    0    3

0    0    0    2

0    0    0    0

0    0    0    0

1     0    0    0

1     0    0    1

1     0    0    1

0    0    0    1

2     0    0    0

1     0    0    2

1     0    0    1

1     0    0    1

1     0    0    1

2     0    0    1

5     0    0    3

 

Muuttointensiteettiä koskevan priorioletuksen päivittäminen

Yhtenä oleellisena haasteena – sekä tarkkailuaikana kokonaisuudessaan Vantaanjoesta mereen siirtyvien smolttien totaalin posterioriestimoinnissa, että seuraavan päivän muuttomäärän estimoinnissa – on määrittää, kuinka muuttointensiteetti vaihtuu päivien kuluessa. Tämän suhteen asiantuntijatietoon perustuvana priorina käytettiin arviota, jonka mukaan a) muutto lähtee (luonnon olosuhteissa) toden teolla käyntiin noin 7 celsius -asteessa sekä b) muutto saavuttaa huippunsa 8 asteen tienoilla. Arvelimme intensiteetin huipun tämän perusteella saavutettavan siis jo seurannan alkuvaiheilla hahmotellen, että päivittäiset muuttomäärät laskevat huipun saavutettuaan hitaammin kuin olivat nousseet, hiipuen pikkuhiljaa lähtötasolle seurannan loppuessa noin 60 päivän jälkeen. Oletamme, että kaikki joessa majailevat smoltit lähtevät liikkeelle seuranta-ajan kuluessa. Tämä olettamus tuskin kirjaimellisesti pitää paikkansa, sillä jo ensimmäisinä päivinä saimme haaviin muutamia havaintoja, mikä antaa ymmärtää, että oletettavasti smoltteja on kulkenut vanhankaupunginkosken ohi jo ennen seurannan aloittamista.

Ennen kaikkea yllä tehdyn priorioletuksen b) perusteella päätimme arvioida yksittäisen smoltin todennäköisyyttä lähteä liikkeelle päivänä i vahvasti oikealle vinolla funktiolla. Tämän kaltaista jakaumamuotoa voitaisiin lähteä approksimoimaan esim. jonkinlaisella muunnelmalla Pareto-, Gamma- tai lognormaalijakaumasta. Päädyimme lognormaalin kaltaiseen funktioon, jonka varmistimme täyttävän todennäköisyysjakauman määritelmän määrittelyjoukossaan jakamalla funktion antamat arvot kaikkien arvojen summalla.

Havaitsimme kuitenkin pian, että lognormaalijakaumalla tullaan olettaneeksi aivan liian suuren osuuden smolteista lähteneen liikkeelle seurannan alkupäivinä. Tämän oletuksen seurauksena yksittäisten päivien muuttoestimaatit vaikuttivat liian korkeilta jo havaituille päiville. Samaisesta syystä totaaliestimaatti vaikutti arvioivan todellista muuttomäärää alakanttiin. Uumoilimme, että muuttohuippu ollaan vasta saavuttamassa ja joko että 1) kasvu huipulle tultaessa ei ole niin jyrkkään – ainakaan seurannan ensimmäisinä päivinä – kuin lognormaalijakauma olettaa, tai 2) että kasvu on vielä jyrkempää, mutta täysin äkillistä. Ehkäpä tiettyjen biologisten ennakkoehtojen kuten sopivan lämpötilan ja valonmäärän osuessa kohdalle, lukemattomat smoltit lähtevät liikkeelle samanaikaisesti.

Kuvio 1. Smolttien ennustetut odotusarvoiset muuttomäärät päivien funktiona mallivalinnalla lognormaali

log_n

Kuten kuviosta 1 havaitsee, mallioletuksella lognormaali smolttien muuton huipun pitäisi jo olla takanapäin, sillä tänään 29.5. olemme päivässä nro 21. Tämä on ristiriidassa asiantuntijatiedon kanssa, jonka mukaan muuttaminen saavuttaa huippunsa noin kahdeksassa asteessa, mikä lämpötila saavutettiin eilen. Samaisesta syystä mallin estimoima suhteellinen muuttaneiden osuus kaikista joen smolteista vaikuttaa liian korkealta.

Korjausliikkeenä em. haasteeseen 1) päätimme muuttaa oletustamme muuttointensiteetistä. Lähdimme mallintamaan asiaa kahdella vaihtoehtoisella tavalla. Ensinnäkin datavetoisesti olettaen muuttotodennäköisyyden vakioksi yli päivien (ts P[i] = 1/60). Toiseksi olettaen jakauman normaaliseksi saavuttaen huippunsa odotusarvoisesti 25 päivän kohdalla. Saimme seuraavat estimaatit päiväkohtaisille muuttomäärien odotusarvoille:

Kuvio 2. Smolttien ennustetut odotusarvoiset muuttomäärät päivien funktiona mallivalinnalla tasajakauma ja normaalijakauma

 Unif_nNormal_n

Nähdäksemme kyseisillä mallivalinnoilla on realistisemmat oletukset jo mereen muuttaneiden smolttien suhteellisesta osuudesta. Datavetoisella tasajakaumapriorilla saamme suhteessa havaittujen smolttien aineistossa esiintyneeseen tasaiseen trendiin perustellun oloisia estimaatteja. Normaalijakaumaoletuksella muuttohuippu on vasta edessä, minkä lisäksi pienempi osuus smolteista oletetaan jo muuttaneen kuin lognomaalioletuksella.

Edellä arveltiin, että menneiden tutkimusten perusteella saattaisi käydä niin, että muutto tapahtuu muutamassa harppauksessa sitten, kun se lähtee käyntiin. Näin ollen päivien yli oletetun muuttointensiteetin jakaumaan saattaa olla perusteltua tehdä yhä joitain kyseiset arviot paremmin huomioonottavia muokkauksia. Esimerkiksi voitaisiin olettaa, että muutto on tasaista (sitä kuvaa tasajakauma) tiettyyn päivään asti, jolloin jakauman muoto muuttuu ikään kuin kategorisena muuttointensiteetin tilan muutoksena. Tämä voidaan sisällyttää malliin olettamalla eri muuttoajanjaksoille (korkea / alhainen intensiteetti) omat niitä parhaiten kuvaavat todennäköisyysjakaumat esimerkiksi edellä jo esiintyneen step -funktion logiikkaa hyödyntäen.