Alometry law for the Nile Tilapia

(image from Fishsource)

We use the alometry law $y = A x^B$ to correlate the length $x$ of the fish tilapia with its mass $y$. The following table for the length-weight relationship of growing-finishing cage-farmed Nile tilapia (Oreocromis niloticus)is provided in the article T. S. de Castro Silva, L. D. dos Santos, L. C. R. da Silva, M. Michelato, V. R. B. Furuya, W. M. Furuya, Length-weight relationship and prediction equations of body composition for growing-finishing cage-farmed Nile tilapia, R. Bras. Zootec. vol.44 no.4 Viçosa Apr. 2015):

We use Turing.jl with the compound model

\[ \begin{align*} A & \sim \mathrm{InverseGamma}(\alpha_A, \beta_A) \\ B & \sim \mathrm{InverseGamma}(\alpha_B, \beta_B) \\ \sigma^2 & \sim \mathrm{InverseGamma}(2, 1) \\ y & \sim \mathrm{Normal}(A x^B, \sigma^2) \end{align*}\]

with suitable hyperparameters $(\alpha_A, \beta_A)$ and $(\alpha_B, \beta_B)$ to be chosen shortly.

We start by loading the necessary packages:

using Distributions, Turing, StatsPlots

Then we collect the adimensionalized data in vector form and plot it for the fun of it:

xx = [10.9, 15.3, 19.1, 22.8, 26.3, 31.3]
yy = [28.6, 88.6, 177.6, 313.8, 423.7, 774.4]

scatter(xx, yy, xlabel="Body length (cm)", ylabel="Body weight (g)", xlims=(0.0, 35.0), ylims=(0.0, 1000.0), title="Length-weight relationship of growing-finishing cage-farmed Nile tilapia", titlefont=10, legend=nothing)

We define the Turing.jl model with parameters $A$ and $B$ as Inverse Gamma functions with hyperparameters Ah = (Ah.α, Ah.β) and Bh=(Bh.α, Bh.β).

@model function alometry(x, q; Ah, Bh)
    A ~ InverseGamma(Ah.α, Ah.β)
    B ~ InverseGamma(Bh.α, Bh.β)
    σ² ~ InverseGamma(2, 1)
    σ = sqrt(σ²)

    for i in eachindex(x)
        y = A * x[i] ^ B
        q[i] ~ Normal(y, σ)
    end
end

alometry (generic function with 2 methods)

As in any nonlinear optimization problem, the starting point is crucial. Here, the starting point is our prior. The following prior does not work properly, as we can see.

model = alometry(xx, yy; Ah=(α=1,β=1), Bh=(α=1,β=1))
chain = sample(model, NUTS(0.65), 1_000)

Chains MCMC chain (1000×15×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 3.69 seconds
Compute duration  = 3.69 seconds
parameters        = A, B, σ²
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Use `describe(chains)` for summary statistics and quantiles.

plt = scatter(xx, yy, xlabel="Body length (cm)", ylabel="Body weight (g)", xlims=(0.0, 35.0), ylims=(0.0, 1000.0), title="Length-weight relationship of growing-finishing cage-farmed Nile tilapia", titlefont=10, legend=nothing)
xxx = range(0.9*first(xx), 1.1*last(xx), length=200)
yyy = mean(chain, :A) * xxx .^ mean(chain, :B)
plot!(plt, xxx, yyy)

If we start with a more informative prior, we get a suitable result. If we pick two data points $(x_1, y_1)$ and $(x_2, y_2)$, assuming $y \approx Ax^B$, we have $y_2/y_1 = (x_2/x_1)^B$, so that

\[ B = \frac{\ln(y_2/y_1)}{\ln(x_2/x_1)}.\]

If we choose the second and third points, we get

    B = log(yy[3]/yy[2])/log(xx[3]/xx[2])

3.1347640575458167

From that, we can also estimate $A$ from $A = y/x^B$, so that, choosing the second point

    A = yy[2]/xx[2]^B

0.017127958628527618

With that in mind, we choose the hyperparameters for the prior as $(\alpha_A, \beta_A) = (57, 1)$, with mean $1/(57+1) \approx 0.017$, and $(\alpha_B, \beta_B) = (1, 6)$, with mean $6/(1+1) = 3.0$.

model = alometry(xx, yy; Ah=(α=57,β=1), Bh=(α=1,β=6))

DynamicPPL.Model{typeof(Main.alometry), (:x, :q), (:Ah, :Bh), (), Tuple{Vector{Float64}, Vector{Float64}}, Tuple{@NamedTuple{α::Int64, β::Int64}, @NamedTuple{α::Int64, β::Int64}}, DynamicPPL.DefaultContext}(Main.alometry, (x = [10.9, 15.3, 19.1, 22.8, 26.3, 31.3], q = [28.6, 88.6, 177.6, 313.8, 423.7, 774.4]), (Ah = (α = 57, β = 1), Bh = (α = 1, β = 6)), DynamicPPL.DefaultContext())

With this prior, we attempt again to fit the model.

# chain = sample(model, HMC(0.05, 10), 4_000) # HMC seems quite unstable here
chain = sample(model, NUTS(0.65), MCMCSerial(), 1000, 3; progress=false)

Chains MCMC chain (1000×15×3 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 3
Samples per chain = 1000
Wall duration     = 1.16 seconds
Compute duration  = 1.0 seconds
parameters        = A, B, σ²
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Use `describe(chains)` for summary statistics and quantiles.

Here is the result of the chain.

plot(chain)

Taking the mean of the parameters $A$ and $B$ we plot the fitted curve.

plt = scatter(xx, yy, xlabel="Body length (cm)", ylabel="Body weight (g)", xlims=(0.0, 35.0), ylims=(0.0, 1000.0), title="Length-weight relationship of growing-finishing cage-farmed Nile tilapia", titlefont=10, legend=nothing)
xxx = range(0.9*first(xx), 1.1*last(xx), length=200)
yyy = mean(chain, :A) * xxx .^ mean(chain, :B)
plot!(plt, xxx, yyy)

This seems successful. Now we compute the 95% credible interval

quantiles = reduce(
    hcat,
    quantile(
        [
            A * x^B for (A, B) in eachrow(view(chain.value.data, :, 1:2, 1))
        ],
        [0.025, 0.975]
        )
    for x in xxx
)

2×200 Matrix{Float64}:
 19.7533  20.5525  21.3735  22.2091  23.0665  …   978.839   989.694  1000.63
 23.0062  23.8972  24.8111  25.7481  26.7085     1046.48   1058.3    1070.21

and plot it along the data:

plt = plot(xlabel="Body length (cm)", ylabel="Body weight (g)", xlims=(0.0, 35.0), ylims=(0.0, 1000.0), title="Length-weight relationship of growing-finishing cage-farmed Nile tilapia", titlefont=10, legend=nothing)
plot!(plt, xxx, yyy, ribbon=(yyy .- view(quantiles, 1, :), view(quantiles, 2, :) .- yyy), label="Bayesian fitted line", color=2)
scatter!(plt, xx, yy, color=1)

We end this section plotting an ensemble of lines generated with the chain.

plt = plot(xlabel="Body length (cm)", ylabel="Body weight (g)", xlims=(0.0, 35.0), ylims=(0.0, 1000.0), title="Length-weight relationship of growing-finishing cage-farmed Nile tilapia", titlefont=10, legend=nothing)
plot!(plt, xxx, yyy, label="Bayesian fitted line", color=2)
for (a, b) in eachrow(view(chain.value.data, :, 1:2, 1))
    plot!(plt, xxx, a .* xxx .^b, alpha=0.01, color=2, label=false)
end
scatter!(plt, xx, yy, color=1)