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”)

 

 

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.

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)

Mallipäivitys 13.5.

Viimeisimmät muutokset käytössä olevaan malliin

  • Vaihtoehtoisen smolttien lähtemistodennäköisyyden luonnetta kuvaavan Beta-jakaumamallin lisääminen
  • Mallitekniset muutokset pyydystettävyyden laskemisessa

Käytössä oleva mallimme perustuu pääpiirteittäin edelleen alkuperäiseen rakenteeseen ja eri malliversioiden yhdistämiseen BMA:ta hyödyntäen. Seuraavassa esitellään mallin eri osat ja tehdyt muutokset (kiitos Leolle selkeistä kommenteista BUGS-skriptissä!).

Vuosittainen kokonaisvaellusmäärä

Muuttavien smolttien kokonaismäärän (N) odotusarvon priorina käytössä on edelleen kolmen asiantuntijaarvion keskiarvo (2600), jonka on odotettu tulevan ainoastaan positiivisia arvoja tuottavasta katkaistusta, kokonaisluvuiksi pyöristetystä normaalijakaumasta.

N<-round(cN)

cN~dnorm(2600, prec_N)I(1,)

prec_N<-pow(sd_N, -2)

sd_N<-2000

Kokonaisvaellusmäärän priorin määrittelemiseksi voisi myös asiantuntijaarvioiden sijaan tai lisäksi hyödyntää tietoa edellisvuosina jokeen nousseiden meritaimenten määristä, joko laskemalla nousukkaiden smolttituotantopotentiaali mortaliteettitekijät huomioiden, tai takautuvasti merikuolleisuuden kautta arvioiden aikaisempien vuosien mereen päätyneiden smolttien määriä.

Lähtemistodennäköisyys

Muuttointensiteetin tai päivittäisen lähtemistodennäköisyyden (p) kuvaamiseksi on tähän mennessä käytetty kolmea eri todennäköisyysjakaumaa; tasajakauma, normaali ja log-normaali. Vaihtoehtoisia malleja on ajettu erikseen rinnakkain ja myöhemmin yhdistettynä mallikeskiarvoistamista (BMA) käyttäen. Nykyisessä mallissa näiden kolmen jakauman lisäksi neljäntenä vaihtoehtona on käytetty beta-jakaumaa. Beta-jakauma (Beta distribution – Wikipedia) on intervallille [0,1] määritelty, melko mukautuva todennäköisyysjakauma, jonka muodon määrittää tiheysfunktion alfa ja beta parametrit. Beta-jakauma sopii usein osuuksien sattumanvaraisen käytöksen mallintamiseen.

d1c8bb0654c111cd0a16d1aafd8b970a

Mallissa on tarkastelujakson jokaisen päivän (1,2,3,…60) huomioivan silmukan (“for-loop”) sisällä määritelty lähtemistodennäköisyys neljällä vaihtoehtoista muodolla( p[i,1], p[i,2], p[i,3] ja p[i,4]). Eri p-funktioiden indeksöinti (1–4 hakasulkeiden sisällä) tulee mallissa myöhemmin käyttöön, kun vaihtoehtoiset, ainoastaan muuttointensiteetin muodon määrittelyssä toisistaan eroavat neljä osamallia sisällytetään päämalliin. Kaikissa muissa p:n funktioissa, paitsi tasajakauma-muodossa, on käytetty yhteistä odotusarvon parametriä (myy_p), joka on määritelty skriptin lopussa yhdessä eri funktioden hajontaparametrien kanssa.

for(i in 1:60) {

# p tasajakaumalla

p[i,1]<-1/60

# p normaalijakaumalla

pl[i] <- exp(- pow((i-myy_p)/sd_p,2)*0.5)

p[i,2] <- pl[i] / sum(pl[1:60])

# p log-normaalijakaumalla

pn[i] <- (1/i)*exp(- pow(log(i)-location, 2) / scale)

p[i,3] <- pn[i] / sum(pn[1:60])

# p beta-jakaumalla

pbi[i]<-pow(i/60,myy_p*b_eta/60)*pow(1-i/60,(1-myy_p/60)*b_eta)

Koska yllä olevasta beta-jakauman funktiosta puuttuu normalisointivakio B (beta-funktio) nimittäjästä, varmistetaan että beta summautuu yhteen yli päivien, näin:

p[i,4]<-pbi[i]/sum(pbi[])

Lähtevien smolttien määrät

Seuraavaksi mallissa määritellään lähtevien smolttien päivittäiset (i) määrät (n) kokonaismäärästä (N) binomijakaumalla eri lähtemistodennäköisuusmalleille 1,2,3,4, ja varmistetaan, ettei päivittäiset lukumäärät ole nolla, koska se aiheuttaisi ongelmia kokonaispyydystettävyyn määrittelyssä käytettävässä binomijakaumassa.

n2[i]~dbin(p[i,model],N)

n[i]<-n2[i]+1

}

Saalis ja pyydystettävyys

Dataa on tähän mennessä kerätty kahdella eri pyyntivälineellä (smolttiruuvi ja rysä), ja niiden vaikutukset on huomioitu mallissa jo aikaisemmin. Ruuvi on ollut käytössä koko tarkasteluajan tähän asti kun taas rysä otettiin käyttöön vasta myöhemmin, jonka jälkeen se otettiin pois mutta on nyt taas pyynnissä. Periaate kahden eri välineen kokonaispyyntihehon laskemisessa on edelleen sama kuin aiemmin:

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

Muutokset päivittäisen kokonaissaaliin ja pyydystettävyyden arvionnissa ovat malliteknisiä, tehden koodista elegantimman ja mallista mahdollisimman kevyen. Uuden rakenteen myötä myös datan taulukkomuoto on muuttunut, sisältäen seuraavat sarakkeet: ScrewX[] (ruuvin saalis), FykeX[] (rysän saalis), tagged[] (merkityt kalat), Screw_op[] (ruuvin käyttöstatus), Fyke_op[] (rysän käyttöstatus), ScrewR[] (ruuvin takaisinpyyntisaalis), FykeR[] (rysän takaisinpyyntisaalis) ja TotalX[] (kokonaissaalis). Mallissa on tarkasteluajan jokaisen tähänastisen ja seuraavan päivän (1,2,3,…days+1) huomioivan silmukan sisällä määritelty päivittäinen kokonaissaalis (TotalX[d]), päivittyvä kokonaispyyntiteho (Total_q[d]), ruuvin päivittäinen saalis (ScrewX[d]), ruuvin päivittyvä pyyntiteho (Screw_q[d]) ja rysän päivittyvä pyyntiteho (Fyke_q[d]) seuravanlaisesti:

for(d in 1:days+1) {

# Kaikkien havaittujen smolttien yhteenlaskettu määrä päivänä d noudattaa binomijakaumaa parametreilla kokonaispyydystettävyys (Total_q[d]), päivittäiset muuttavat smoltit (n[d])

TotalX[d]~dbin(Total_q[d],n[d])

Total_q[d]<-Screw_q[d]+Fyke_q[d]-Screw_q[d]*Fyke_q[d]

# Ruuvin havainnot noudattavat poisson -jakaumaa parametrilla odotusarvoiset lähtevät smoltit (binomijakaumalla tulee ongelmia, sillä päivittäiset lähtevät määrät saattavat saada arvoja nolla)

ScrewX[d]~dpois(screw_muX[d])

# Odotusarvoiset lähtevät smoltit päivänä d (screw_muX[d]) määrittyvät ruuvin suhteellisella teholla (Screwprop[d])*kuinka moni jää kaiken kaikkiaan kiinni (TotalX[d])

screw_muX[d] <- Screwprop[d]*TotalX[d]

#Ruuvin suhteellinen teho on ruuvin teho jaettuna kokonaisteholla (ruuvi JA rysä)

Screwprop[d]<-Screw_q[d]/Total_q[d]

# Ruuvin suhteellista tehoa päivitetään vain kun ruuvin on käytössä (indikaattorifunktio screw_op saa arvoja 0, kun ei käytössä, arvoja 1 kun käytössä)

Screw_q[d]<-q_S*Screw_op[d]

# Rysän pyydystettävyys päivittyy samalla indikaattorifunktion logiikalla

Fyke_q[d]<-q_F*Fyke_op[d]

# Priori määrälle merkittyjä kaloja (tämä on joku tekninen eli laskennallinen yksityiskohta)

tagged[d]~dbin(1,1000)

}

Seuravaksi mallissa määritellään, kuinka pyydystettävyys päivittyy uudelleenpyydettyjen havainnoilla. Päivät (d) kulkevat toisesta päivästä, koska vasta siitä lähtien on mahdollista pyydystää päivänä 1 merkitty kala.

for(d in 2:days+1) {

# Rysän uudelleenpyytämien kalojen havaintojakauma noudattaa binomijakaumaa parametreilla rysän pyyntiteho(Fyke_q[d]) , edellisenä päivänä merkityt kalat (tagged[d-1])

FykeR[d]~dbin(Fyke_q[d],tagged[d-1])

# Ruuviin voi periaatteessa jäädä kaikki ne kalat, jotka on edellisenä päivänä merkitty, mutta jotka eivät jääneet ruuviin

Screw_nR[d]<-tagged[d-1]-FykeR[d]

# Ruuvin pyydystettävyyden havaintojakauma noudattaa binomijakaumaa parametreilla ruuvin pyydystettävyys (Screw_q[d]), ja edellisessa kaavassa määritellyt mahdolliset määrät (Screw_nR[d])

ScrewR[d]~dbin(Screw_q[d],Screw_nR[d])

}

Rysän pyydystettävyyden prioria (q_F) muutettiin säilyttäen pseudohavaintojen suhteet. Arvoilla 1 ja 4 voidaan tehdä prioriarvio epävarmemmaksi siksi, ettei mandariinikoetta ole toistettu rysälle. Samat suhteet mallintavat kuitenkin ennakkokäsitystä, jonka mukaan rysä ja ruuvi pyydystävät yhtä hyvin. Tämän oletuksen uskottavuudesta voidaan kuitenkin olla eri mieltä, vaikka pyydysten ominaisuuksien eroja ei tiedetäkään, koska rysä pyydystää eri kohdassa jokea ruuviin verrattuna ja koska on epätodennäköistä että smoltit ovat tasaisesti jakautuneita koko virran leveydeltä jo mandariinikokeen tulostenkin perusteella.

q_F~dbeta(1,4)

q_S~dbeta(2,8)

Eri osamallejen sisällyttäminen ja mallikeskiarvoistaminen

Viimeisessä tätä edellisessä versiossa yhdistettiin jo eri lähtemistodennäköisyyksien määrittelemät osamallit yhdeksi malliksi Bayesilaisella mallikeskiarvoistamisella (BMA). Ainoa muutos tämän suhteen tähän malliin, on neljännen mallin, eli lähtemistodennäköisyyden beta-jakaumamallin lisääminen, joka tehtiin seuraavanlaisesti:

# Priori mallin BMA:lle

model~dcat(z[1:4])

# Tässä määritellään prioritodennäköisyys eri mallien totuudelle

z[1]<-1/4 # Tasajakaumamalli

z[2]<-1/4 # Normaalijakaumamalli

z[3]<-1/4 # Log-normaalijakaumamalli

z[4] <-1/4 # Beta-jakaumamalli

# Tästä saadan eri mallien posterioritodennäköisyydet

Z[1]<-equals(model, 1)

Z[2]<-equals(model, 2)

Z[3]<-equals(model, 3)

Z[4] <- equals(model,4)

Lähtemisintensiteettien priorit

myy_p~dnorm(25,tau_myy)I(0.01,)

tau_myy<-pow(7,-2)

sd_p~dnorm(10, tau_sd)I(0.01,)

tau_sd<-pow(5, -2)

scale <- log(1 + pow(sd_p/myy_p,2) )

location <- log(myy_p) -0.5*scale

b_eta~dunif(5,100)

}

Data

list(days=33)

ScrewX[] FykeX[] tagged[] Screw_op[] Fyke_op[] ScrewR[] FykeR[] TotalX[]

0 0 0 1 0 0 0 0

0 0 0 1 0 0 0 0

0 0 0 1 0 0 0 0

1 0 0 1 0 0 0 1

2 0 3 1 0 0 0 2

2 0 2 1 0 0 0 2

0 0 0 1 0 0 0 0

0 0 0 1 0 0 0 0

0 0 0 1 0 0 0 0

1 0 0 1 0 0 0 1

1 0 0 1 0 0 0 1

1 0 3 1 0 0 0 1

0 0 0 1 0 0 0 0

2 0 2 1 0 0 0 2

1 0 1 1 0 0 0 1

0 1 1 1 1 0 0 1

0 1 0 1 1 0 0 1

0 1 0 1 1 0 0 1

1 1 3 1 1 0 0 2

3 2 5 1 1 0 0 5

6 2 8 1 1 0 0 8

4 0 4 1 0 0 0 4

2 0 2 1 0 0 0 2

2 0 2 1 0 0 0 2

1 0 1 1 0 0 0 1

3 0 3 1 0 0 0 3

2 0 2 1 0 0 0 2

1 0 2 1 0 1 0 1

3 0 3 1 0 0 0 3

1 0 1 1 0 0 0 1

1 0 1 1 0 0 0 1

1 0 1 1 0 0 0 1

2 0 2 1 0 0 0 2

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

NA NA NA 1 1 NA NA NA

END

#BUGS:iin ladattavat eri mallien alkuarvot, hyvä käyttää eri alkuarvoja

list(model=1)

list(model=2)

list(model=3)

list(model=4)

Huomioita lopuksi

Vaikka mallin aktiivinen kehittäminen jäänee lähiaikoina hieman vähäisemmäksi, voisi tulevaisuudessa edelleen yrittää sisällyttää tietoa ympäristömuuttujista (lämpötila, vedenkorkeus, ym) ja niiden vaikutuksista, sekä vaellusmatkan pituudesta ja eri smolttituotantoalueiden (koskien) kapasiteetistä ja sijoittumisesta. Lähiaikoina jatkamme mallin sisällyttämisellä päätösanalyysiin, joka koskee kysymystä Vanhankaupunginkosken padon mahdollista purkamista.

 

Asiantuntijatiedon keräämisestä

Jatkan nyt edellisestä postauksestani, jossa mietin Bayes-mallintamista päätöksenteon työkaluna. Bayes-mallinnus on omiaan epävarmuuden ollessa tutkimuskohteen osalta suuri, ja jolloin asiantuntijatiedon rooli on tärkeä. Asiantuntijatietoa voidaan käyttää priorien määrittämiseen, mutta asiantuntijatiedon avulla voidaan myös paremmin hahmoittaa tutkimuksen aihetta sekä suurempia kokonaisuuksia ja kokonaiskuvaa (Uusitalo, 2007).

Kuvailenkin nyt lyhyesti ja pääpiirteittäin, miten asiantuntijatietoa voidaan kerätä. Alla oleva malli kuvaa hyvin yksinkertaisella tasolla asiantuntijatiedon keräämistä. Ensimmäinen askel on tietenkin taustatyö ja ongelmanasettelu, sekä itse asiantuntijoiden valitseminen. Onkin mielenkiintoista pohtia, kenet määritellään asiantuntijaksi? Omasta mielestäni asiantuntija voi olla kuka vaan, joka on jollain tasolla kyseessä olevaan asiaan linkittynyt: Vanhankaupunginkosken meritaimenen kohdalla asiantuntija voi olla niin kalastamisesta vastaava viranomainen, kuin vapaa-ajan kalastajakin. Asiantuntijan voisikin uudelleenmääritellä, ja mielestäni “asiantuntijan”  (expert) sijaan voisi olla hyödyllistä käyttää neutraalimpaa sanavalintaa, ”sidosryhmän jäsentä”(stakeholder), tai pelkästään ”tekijää” (actor).

elicit

Malli 1. Asiantuntijatiedon kerääminen. (O’Hagan et al., 2005).

Seuraava vaihe on itse haastattelu. Haastattelu painottaa asiantuntijan sen hetkistä tietoa aiheesta. Tutkijan roolina haastattelussa on lisätä sujuvuutta ja fasilitoida prosessia: tutkijan tulee valmistella haastateltavaa Bayes-mallintamisesta ja varmistaa, että haastateltava ymmärtää mallintamisen pääpiirteittäin. Jos näin ei tehdä, asiantuntija saattaa suhtautua epäluuloisesti mallintamiseen, eikä haastateltavan ja tutkijan välille synny jo tutkimuksenkin kannalta tärkeää luottamuksen tunnetta (Uusitalo, 2007). Eräänlaista koehaastattelua on hyvä harkita, jota voidaan sitten haastateltavien palautteen kautta muokata sujuvampaan ja ymmärrettävämpään muotoon. Haastattelu on myös hyvä videoida tai dokumentoida, jotta siihen voi palata tutkimuksen myöhemmissä vaiheissa. Haastateltavien mahdolliset henkilökohtaiset intressit tulee myös tiedostaa ja ottaa huomioon. (O’Hagan et al., 2005).

Kolmannessa vaihessa asiantuntijoiden antamat priorit ja arvot voidaan mukauttaa muiden asiantuntijoiden lausuntoihin sekä kerättyyn dataan. Tässä voidaan käyttää myös muissa blogipostauksissa viitattua Bayesilaista keskiarvoistamista, missä eri malleille voidaan antaa eri painoituksia. Viimeisessä vaiheessa mietitään tiedonkeruun kattavuutta ja riittävyyttä, jonka jälkeen voidaan päätyä keräämään vielä täydennystä ja lisää asiantuntijatietoa. (O’Hagan et al., 2005).

Asiantuntijatiedon keräämisessä ja käyttämisessä esintyy erinäisiä ongelmia. Ihmisillä on usein rajallinen kyky ymmärtää todennäköisyyksiä ja asettaa niille tietty numerinen arvo, jolloin erilaisten vääristymien syntyminen on tavanomaista. Haastattelun, eri tehtävien, ja kysymysten asettelu vaikuttavat kaikki haastatteluprosessin tuloksiin. Onnistumista on vaikea määrittää, koska annetut arvot kuvaavat haastateltavan subjektiivistä ”totuutta”, joka vaihtelee haastateltavien kesken. Koehaastattelun järjestäminen ja kysymysten uudellenasettelu voivat auttaa todennäköisyyksien hahmoittamisessa.

Lisäksin, niin kuin jo aikaisemmin mainittu, asiantuntijoilla voi mahdollisesti olla erilaisia henkilökohtaisia intressejä, jotka mallintamisessa pitää tiedostaa ja huomioida. Asiantuntijat voivat niin sanotusti pelata omaa peliään, eli antamalla väärää tietoa määrittelemällä prioreja tai tekijöitä omaksi edukseen. Mallien keskiarvoistamisella vääristymää voidaan korjata: keskiarvoistaessa mallit, jotka eivät tue kerättyä dataa, saavat vähemmän painoarvoa. Sen sijaan, asiantuntijoiden antamia tavoitteita ja arvoja päätösmallintamisessa ei voida, eikä tulekaan, tässä vaiheessa keskiarvoistaa (Mäntyniemi et al., 2013).

Vaikka mallintamisen onnistuminen on tärkeää, ehkä tärkeintä on kuitenkin haastateltavan ja tutkijan välisen dialogin ja yhteistyön mahdollistaminen, sekä mahdollisimman avoin ja läpinäkyvä ongelmanasettelu. Näin voidaan varmistaa, että tutkija ei tutki vaan tutkiakseen, vaan hänen työllään on myös yhteiskunnallinen merkitys ja ulottuvuus.

Mäntyniemi, S., et al., 2013. Incorporating stakeholders’ knowledge to stock assessment: Central Baltic Herring. Can. J. Fish. Aquat. Sci, 70.

O’Hagan, A., et al., 2005. Statistical Methods for Eliciting Probability Distributions. Journal of the American Statistical Association, 100: 470.

Uusitalo, L., 2007. Advantages and challenges of Bayesian networks in environmental modelling. Ecological Modelling, 203.

 

Mietteitä Bayes-mallintamisesta ja päätöksenteosta

Tässä blogipostauksessa on luvassa alustavia mietteitäni Bayes-mallintamisen roolista ja hyödyllisyydestä ympäristöä koskevassa päätöksenteossa, tässä tapauksessa erityisesti Vantaanjoen meritaimenen kohdalla. Toisin sanoen, mihin Bayes-mallintamista tarvitaan, ja mikä mallintamisen hyöty on/voi olla?

Oma taustani on ympäristötieteissä ja kehitysmaatutkimuksessa. Kehitysmaatutkimuksessa painotetaan mahdollisimman yksinkertaisen ja helppokäyttöisen teknologian arvoa, jonka nähdään olevan omiaan lisäämään itse teknologian käyttöä, ja sitä kautta myös hyvinvointia globaalin etelän maiden eri yhteisöissä. Esimerkkinä voisi pitää vaikka aurinkopaneeleiden käyttöä: jos yhteisö ei pysty omaksumaan paneeleiden käyttöä, eikä osaa niitä itse käyttää tai huoltaa, aurinkopaneeleiden edut jäävät pieniksi. Samalla tavalla myös osallistavaa päätöksentekoa lisäävien eri työkalujen ja välineiden nähdään edesauttavan yhteisöiden hyvinvointia ja päätöksenteon oikeudenmukaisuutta. Oletuksena siis on, että nykyisen muotoisen ylhäältä alaspäin toimivan säännöstelyn lisäksi, horisontaalisen päätöksenteon sekä yhteisöstä käsin lähtevä päätöksenteon lisääminen edesauttaisi säännöstelyn legimiteettiä ja toimivuutta.

Myös täällä Itämeren rannikkolla on tarvetta alhaalta ylöspäin suuntautuvalle päätöksenteolle, ja olenkin kiinnostunut osallistavasta päätöksenteosta, jossa tarkoituksena on eri sidosryhmien ja tekijöiden osallistaminen mahdollisimman yksinkertaisin työkaluin. Tälle kurssille tulin miettimään, miten Bayes-mallintaminen voisi toimia yhtenä osallistavan päätöksenteon välineenä. Itselläni ei ole matemaattista taustaa, ja suhtauduinkin mallintamiseen aluksi epäilevästi. Bayes-laskennan hienous piilee kuitenkin siinä, että koko universumia mallintamista ei nähdä välttämättömänä ennen itse tekoihin ryhtymistä, vaan päätöksenteko voi perustua todennäköisyyksille ja niiden mallintamiselle.

Vanhankaupunginkosken meritaimenia ja muita vaelluskaloja koskeva päätöksenteko eri sidosryhmien välillä on hankalaa: kaikki eivät luonnollisestikaan ole samaa mieltä vaelluskaloja koskevista toimenpiteistä. Aiemmissa blogipostauksissa (Juhani ja Mikko) on jo kiinnitetty huomiota Vanhankaupunginkosken padon haittoihin meritaimenien ja muiden vaelluskalojen kohdalla. Vaikka asiaa vaelluskalojen kohdalla on tutkittu vähän, haitoista merkkinä ovat suvannosta kesäaikana löytyneet turbiinien katkomat ankeriaat1. Vanhankaupunginkosken padon purkamisesta on puhuttu viime syksystä, ja tällä hetkellä purkupäätös on saanut jatkoaikaa2. Vastakkain ovat kaksi osapuolta, joista toinen haluaa säilyttää padon sen kulttuurihistoriallisen arvon takia, ja toinen suojella vaelluskaloja.

Vaelluskaloja suojellaan ja säädellään Vanhankaupuginkoskella lähinnä paikallisella tasolla. Vanhankaupunginkosken kalavesien hoidosta vastaa kalastuslain määrämänä Helsingin kalastusalue, jonka tehtävänä on myös kalastuksen säätely ja valvonta, sekä toimiminen yhteistyössä kalavesien omistajien, ammattikalastajien ja vapaa-ajankalastajien kanssa. Helsingin kaupungin liikuntavirasto sen sijaan vastaa kosken kalastusjärjestelyistä. Taimen ja lohi ovat rauhoitettuja ajalla 11.9- 15.11, jolloin niitä ei saa nostaa edes kuvattavaksi. 3

Myös koko Vantaanjoen valuma-alue ja siihen liittyvä päätöksenteko vaikuttavat meritaimeniin. Joka vuotiset kuntien jätevesipuhdistamoilta ja pumppaamoilta tapahtuvat puhdistamattomien jätevesien päästöt, ja maataloudesta huuhtoutuvat ravinteet ja kiintoainekset aiheuttavat pahimillaan kala-ja eliöstökuolemia ja kuormittavat Vantaanjokea, sekä viime kädessä myös Itämerta. Vantaanjoen vedenlaatu olikin 1980-luvulla vielä huono, mutta on siitä lähtien kehittynyt parempaan kuormituksen vähenemisen myötä, ja Vantaanjoen vedenlaadun ekologinen tila on nykyään EU-kriiterein kuvattu ”tyydyttäväksi”. 4

Vantaanjoen varrella ja Vanhankaupunginkoskella on siis monenlaista toimintaa ja monenlaisia toimijoita. Bayes-mallinnus voi toimia yhtenä työkaluna meritaimeniin liittyvässä eri sidosryhmien välisessä keskustelussa ja päätöksenteossa. Eri sidosryhmät ja asiantuntijat voidaan ottaa mukaan jo Bayes-mallin kehittelyvaiheessa, ja heidän avullaan voidaan pohtia eri prioreja (tässä projektissa meritaimenten vuosittaisen määrän priorit on asiantuntijoiden antamia), eri tekijöiden vaikuttavuutta, sekä eri päätösmalleja ja niiden hyödyllisyyttä. Olisikin mielenkiintoista pohtia vielä lisää, miten tarkalleen ottaen asiantuntijat voivat olla mukana kehittelyvaiheessa ja mikä heidän panoksensa mallintamisessa voikaan olla.

Osallistavassa päätöksenteossa työkalujen pitää olla helppokäyttöisiä ja kustannustehokkaita, ja vaikuttaa, että parhaimmillaan Bayes-mallinnus on molempia. Tutkijoiden ja haastateltavien välinen yhteisymmärrys on mallintamisessa tärkeää, ja haasteena on kehittää mallintamisen selkokielisyyttä. Tästä kuitenkin lisää vasta seuraavissa blogipostauksissa….

  1. Virtavesien hoitoyhdistys Ry, 2015. Vantaanjoen vesistö. Available online: http://www.virtavesi.com/index.php?upperCatId=4&catid=4
  2. Vihreä Lanka, 05.05. 2015. Kohuttu patopäätös sai jatkoaikaa. Available online: http://www.vihrealanka.fi/uutiset-kotimaa/kohuttu-patop%C3%A4%C3%A4t%C3%B6s-sai-jatkoaikaa
  3. Vanhankaupunginkoski, 2013. Vanhankaupunginkoski. Available online: http://www.vanhankaupunginkoski.ota.fi/
  4. Vantaanjoen ja Helsingin seudun vesiensuojeluyhdistys ry, 2015. Vedenlaatu. Available online: http://www.vhvsy.fi/content/fi/1007/1048/Jokiluonto.html

 

 

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.

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?

 

Tulokset 30.04-01.05

Image

Ennustukset päiville

22 ja 23

(eli ajalle 30.04-01.05)

 

Eilen (päivä 21) oli saatu yhteensä kuusi taimenta! Vappuluvut ennustetaan tässä nyt uusimmalla mallilla, jossa siis on otettu huomioon edellisessä postauksessa selitetyt tekijät: toinen pyydys (taimenrysä), muutosintensiteetin jakauman uudellenarviointi (käytössä nyt normaali-ja tasamalli), sekä uudet asiantuntijoiden antamat priorit kokonaismäärälle (kokonaislukumäärän odotusarvo 2600, keskihajonta 2000).

Tulokset:

Kuvio 1. Kokonaismäärän posteriori-jakauma normaalimallissa 1

 

 


1

 

Taulukko 1: mallin 1 antamat estimaatit taimenten kokonaismäärän jakauman tunnusluvuille (N Sample size = 147600):

mean
sd MC_error 95% interval
1199.0 834.3   26.68 283-3261

Kokonaismäärän odotusarvon piste-estimaatiksi normaalimallilla on siis laskettu 1199 kappaletta (95% luottamusväli 283-3261), joka, niin kuin saatettiin olettaa, on hieman korkeampi kuin edellisen päivän 1086.

 

Tasajakaumalla (malli 2) kokonaismäärän tulokset näyttävät hieman korkeammilta:

 

Kuvio 2. Kokonaismäärän posteriori-jakauma tasajaukamalla, eli mallissa 2
2

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

mean
sd MC_error 95% interval
1578 1021 36.82 459-4274

 

Ennustukset päiville 22 ja 23 sen sijaan ovat seuraavanlaiset:

Normaalimallilla (malli 1)

Taulukko 3. Posteriorijakauma vaeltavien taimensmolttien määrä torstaina 30.4.2015 (n(22)) ja perjantaina 01.05 (n(23)).

n(i),sample size 147600                       mean sd MC_error 95% interval
n(22) 28.83 18.24 0.5629 7.0-76.0
n(23) 29.83 19.03 0.5872 7.0-79.0

JA

Tasamallilla (malli 2)

Taulukko 3. Posteriorijakauma vaeltavien taimensmolttien määrä torstaina 30.4.2015 (n(22)) ja perjantaina 01.05 (n(23)).

n(i),sample size 156200

  

mean           

 

sd MC_error 95% interval
n(22) 26.30 17.75 0.6137 6.0-73.0
n(23) 26.31 17.76 0.6137 6.0-73.0

 

Viimeiseksi vielä luvut pyydyksiin jäävien taimensmolttien määrästä torstaina 30.4.2015 (x(22)) ja perjantaina 01.05 (x(23)).

 

Normaalimallilla (malli 1):

Kuvio 3 ja 4. Posteriorijakauma ruuviin jäävien taimensmolttien määrästä torstaina 30.4.2015 (x(22)) ja perjantaina 01.05 (x(23)).

4, 5

Mallin mukaan keskiarvo (mean) oli siis:

x(22)= 3,463 ja x(23)= 3,595

Tasamalli (malli 2):

Kuvio 5 ja 6. Posteriorijakauma ruuviin jäävien taimensmolttien määrästä torstaina 30.4.2015 (x(22)) ja perjantaina 01.05 (x(23)).

5,6

Eli tasamallin mukaan keskiarvo (mean) oli:

x(22)= 2,664 ja x(23)= 2,670

Tässä siis esimakua uusimman mallin tuloksista! Mallissa on vielä paljon kehitettävää, ja esim. vieläkin mietimme, miten ruuvin ja rysän vaikutusta voisi parhaiten mallintaa. Nyt kun viime päivinä on satanut paljon ja vedenkorkeus on noussut, rysä otetaan mahdollisesti pois käytöstä, mutta ruuvin pitäisi toimia sitäkin paremmin…

Nyt toivotamme kaikille kuitenkin hauskaa (ja toivottavasti aurinkoista) vappua!!