Bayesianske fiskevann

Oslomarkas fiskeadministrasjon administrerer hundrevis av fiskevann i Oslo og omegn, men uten erfaring er det vanskelig å vite hvor man kan forvente å få fisk. Heldigvis har de også et datasett med løpende oversikt over fangst. Kan vi finne ut hvilke vann som har store fisk, selv når det er registrert svært få fangster i de fleste vannene?

Kartet viser noen av fiskevannene i Oslomarka.

Datasettet ser slik ut:

art vekt lengde vann
Gjedde nan 45 Nøklevann
Røye 400 37 Langvann
Ørret 430 nan Store sandungen
Ørret 700 40 Nordvannet
Ørret 352 36 Elvann

Filtrering og antagelser

Vi begynner med å filtrere dataene:

  • Vi ser kun på abbor, gjedde og ørret. Disse artene representerer \(95 \, \%\) av observasjonene i datasettet.
  • Vi ser kun på vann med minst \(4\) observasjoner. Da gjenstår det \(21\) vann.

Disse to filtreringene reduserer datasettet fra \(413\) til \(167\) observasjoner. I \(106\) av disse er fiskens vekt registrert. Vi bør også være klar over to potensielt store svakheter med datasettet:

  • Det er mulig at fangst ikke blir registrert, og det er nok ikke “missing at random”. Det er mulig at store eller små fangster har større sannsynlighet for å mangle data.
  • Vi har bare observasjoner av fanget fisk, ikke observasjoner av hvor mange timer man har prøvd å få fisk i de ulike vannene.

Figuren nedenfor viser fordelingen til fiskeartenes vekt i datasettet (til venstre), samt forholdet mellom lengde og vekt (til høyre). Legg merke til at forhold mellom lengde og vekt ikke er lineært–vi kommer tilbake til dette senere.

Modell 1 - ett nivå

Det er alltid greit å begynne med en enkel modell. La \(S\) være student-t fordelingen med \(\nu\) frihetsgrader. Vi trener følgende bayesianske modell på \(n=106\) datapunkter der vi har vekta på fisken.

\begin{align*} \text{vekt}_i &\sim \text{S}(\mu_i, \sigma, \nu) \\ \mu_i &= \alpha + \beta_{\text{art}[i]} + \gamma_{\text{vann}[i]} \\ \alpha &\sim N(500, 150) \\ \beta_{\text{art}}, \gamma_{\text{vann}} &\sim N(0, 250) \\ \sigma &\sim \text{cauchy}(0, 1) \\ \nu &\sim \text{gamma}(2, 0.1) \end{align*}

Vi er i utgangspunktet interessert i å estimere parametrene \(\gamma_{\text{vann}}\), men legger også til \(\beta_{\text{art}}\) slik at modellen kontrollerer for art og vann. For å forklare dette, anta at vi har observert følgende data fra to ulike vann A og B:

Art / Vann A B
Abbor 10 fisk med snitt lik 300g 1 fisk med snitt lik 150g
Gjedde 1 fisk med snitt lik 2000g 10 fisk med snitt lik 1000g

Hvilket vann har størst fisk? Tar vi snitt over hvert vann vinner B. Men ser vi på hver art individuelt, så har vann A både større abbor og gjedde. At modellen kontrollerer for artene, betyr enkelt forklart at den ville svart at vann A har størst fisk.

La oss plotte sannsynlighetsfordelingene til \(\alpha + \gamma_{\text{vann}}\) for hvert vann. En måte å tolke figuren på er at den viser fordelingen til gjennomsnittsfisken i hvert vann, ikke fiskepopulasjonen. Eksempelvis har Lutvann har et gjennomsnitt på \(745\), og \(50 \, \%\) av fordelingen ligger mellom \(657\) og \(833\). Det betyr ikke at halvparten av fiskepopulasjonen ligger i dette intervallet, men at modellen er \(50 \, \%\) sikker på at gjennomsnittsfisken ligger i dette intervallet.

\(x\)-aksen er det oppgitt gram, men dette må tolkes med forsiktighet. For å tolke tallene helt bokstavelig må vi anta at det er \(1/3\) av abbor, gjedde og ørret i hvert vann. Det er mest hensiktmessig å bruke figuren til å sammenligne vannene med hverandre, ikke så tolke de absolutte verdiene til enkelt vann.

De blå prikkene viser gjennomsnittet til fordelingen. Vi ser at Store skillingen ser ut som et bra fiskevann, mens Paddetjern ikke er fullt så bra.

Modell 2 - hierarkisk modell

Vi bygger videre på forrige moell ved å gjøre to endringer:

Hyper-priors. I forrige modell satte vi \(\gamma_{\text{vann}} \sim N(0, 250)\). Hvert vann fikk en separat prior-fordeling, uten mulighet til for å dele informasjon. For å dele informasjon og regularisere estimatene, setter vi opp en hierarkisk modell. Dette vil i praksis vekte ned få ekstreme observasjoner. Den teoretiske forkjellen er at vi bytter ut \(N(0, 250)\) med \(N(0, \sigma_{\gamma})\), der \(\sigma_{\gamma}\) er en parameter som modellen selv lærer. Vi velger å sette hyper-prioren til \(\sigma_{\gamma} \sim N(250, 50)\) fordi da prior-distribusjonene \(N(250, 50)\) og \(N(0, \sigma_{\gamma})\) veldig like om man integrerer ut \(\sigma_{\gamma}\).

Ulike feil per fiskeart. I forrige modell hadde vi en overordnet \(\sigma\) for alle kombinasjoner av vann og art. Å anta at \(\sigma\) er lik for hver art virker urealistisk. Det er mye større variasjon mellom gjedder enn mellom ørret. Vi endrer modellen slik at hver art får hver sin egen \(\sigma\).

\begin{align*} \text{vekt}_i &\sim \text{S}(\mu_i, \sigma_{\text{art}[i]}, \nu) \\ \mu_i &= \alpha + \beta_{\text{art}[i]} + \gamma_{\text{vann}[i]} \\ \alpha &\sim N(500, 150) \\ \beta_{\text{art}} &\sim N(0, \sigma_{\beta}) \\ \sigma_{\beta}&\sim N(250, 50) \\ \gamma_{\text{vann}} &\sim N(0, \sigma_{\gamma}) \\ \sigma_{\gamma}&\sim N(250, 50) \\ \sigma_{\text{art}} &\sim \text{cauchy}(0, 1) \\ \nu &\sim \text{gamma}(2, 0.1) \end{align*}

Resultatet er rimelig likt, men det er noen forskjeller:

  • Fordelingene er noe regularisert (dratt inn mot midten).
  • Estimatene blir noe lavere fordi modellen klarer å plassere variasjonen i gjeddene på \(\sigma_{\text{gjedde}}\).
  • Enkelte vann endrer rangering.
  • Fordelingen til Nøklevann har blitt smalere, og ser mer fornuftig ut.

Modell 3 - imputering

I datasettet er det enten registrert lengde eller vekt på fisken, eller begge deler. Til nå har vi trent på \(106\) observasjoner der vi har vekten. Men vi har også \(61\) observasjoner der vi mangler vekt, men har lengde. Å kaste bort en så stor andel av dataene virker lite hensiktsmessig, så la oss modellere vekt som en funksjon av lengde og lære vekten der vi ikke har data. Å imputere data vil si at vi estimerer de ukjente verdiene istedet for å fjerne observasjonene.

Vi kunne modellert forholdet mellom vekt og lengde som lineært, men det ser ikke ut til å stemme godt med dataene (se første figur i artikkelen). La oss modellere en fisk som en sylinder, slik at \(\text{vekt} = \rho \pi r^2 \ell\). Her er \(\rho\) massetettheten, \(r\) er fiskens radius og \(\ell\) er lengden. Anta at forholdet mellom radius og lengde er konstant, slik at \(r = k \ell\). Da får vi

\begin{align*} \text{vekt} = \rho \pi (k \ell)^2 \ell = (\rho \pi k^2) \ell^3 = \eta \ell^3, \end{align*}

der vi har byttet ut den ukjente konstanten \(\rho \pi k^2\) med \(\eta\). En kubisk likning passer godt til dataene, så vi tar utgangspunkt i denne modellen og sier at i de tilfellene der vi mangler vekt, men har lengde, så er

\begin{align*} \text{vekt}_i &\sim \text{N}(a + \text{lengde}_i^p, \epsilon) \\ a &\sim N(0, 1) \\ p &\sim N(3, 0.2) \\ \epsilon &\sim \text{cauchy}(0, 1) \end{align*}

Denne modellen for imputering brukes i kombinasjon med den hierarkiske modellen for å få litt mer informasjon ut av dataene:

Her en en figur som viser sannsynlighetsfordelingene til \(\alpha + \beta_{\text{art}}\) for hvert art.

Oppsummering

  • Vi modellerte forholdet mellom vekt og lengde ved å se for oss en fisk som en sylinder. Dette er et eksempel på å bruke en vitenskapelig modell til å lage statististiske modeller, istedet for å anta en lineær modell. Kapittel 16 i boka “Statistical Rethinking” av McElreath gir noen flere eksempler. Artikkelen Model building and expansion for golf putting er også god.

Notater og referanser

  • Vi modellerte forholdet mellom vekt og lengde ved å se for oss en fisk som en sylinder. Dette er et eksempel på å bruke en vitenskapelig modell til å lage statististiske modeller, istedet for å anta en lineær modell. Kapittel 16 i boka “Statistical Rethinking” av McElreath gir noen flere eksempler. Artikkelen Model building and expansion for golf putting er også god.