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.