Click here for other works of the author on RPubs

The von Bertalanffy growth function is used to describe body length of fish as a function of its age. It is defined as \(L(t) = L_{\infty}(1 - e^{- K(t - t_0)})\) where \(L_{\infty}\) is the hypothesized (mean) maximum length, \(K\) is growth rate, \(t\) is current age and \(t_0\) is the hypothesized age when body length is 0.

In this assignment, we try to model body length of fish using the following regression model: \(L(t) = L_{\infty}(1 - e^{-K(t - t_0)}) * e^{\epsilon}\) with \(\epsilon\) ~ \(N(0, \sigma^2)\). Parameters will be estimated using the Bayesian grid search with priors \(L_{\infty}\) ~ \(U(40, 100)\) and \(K\) ~ \(U(0.1, 0.6)\).

Q. Assume \(t_0\)=0 (two parameters). Please estimate the posterior distributions of parameters using the method Bayesian grid search for the average length-at-age for the female Pacific Hake, respectively, considering a multiplicative error model with log-normal distribution.

Load package

library(knitr)
library(ggplot2)

Load data

age = c(1, 2, 3.3, 4.3, 5.3, 6.3, 7.3, 8.3, 9.3, 10.3, 11.3, 12.3, 13.3)
length = c(15.4, 28.03, 41.18, 46.2, 48.23, 50.26, 51.82, 54.27, 56.98, 58.93, 59, 60.91, 61.83)

#show data
kable(cbind(age, length), col.names = c("Age (year)", "Female mean length (cm)"))
Age (year) Female mean length (cm)
1.0 15.40
2.0 28.03
3.3 41.18
4.3 46.20
5.3 48.23
6.3 50.26
7.3 51.82
8.3 54.27
9.3 56.98
10.3 58.93
11.3 59.00
12.3 60.91
13.3 61.83

Explore the relationship between body length and age

#prepare data for plotting
dataset = data.frame(age, length)
colnames(dataset) <- c("Age", "Body length")

#draw graph
ggplot(dataset, aes(x = Age, y = `Body length`)) +
    geom_point()

It seems that growth rate declines as age increases and the von Bertalanffy growth function is appropriate.

The von Bertalanffy growth function

Define function VBGF that calculates estimated length using the von Bertalanffy growth function, which would be used in our regression model. The parameter \(t_0\) would be set to 0.

VBGF <-function(x, Linf, K){
    y = Linf * (1 - exp(- K * (x - 0)))
    return(y)
}

Log-likelihood function for the regression

Define function lognormal_like which calculates the total likelihood of the regression model. The parameter sigma would be set to 0.2.

lognormal_like = function(Linf, K){
    
    like = numeric(length(length))
    LL = numeric(length(length))
    ypred = VBGF(age, Linf, K)
    dev2 =(log(length) - log(ypred)) ^ 2
    sigma = 0.2
    
    for (i in 1:length(length)){
        like[i] = (1 / (age[i] * sqrt(2 * pi) * sigma)) * exp(-dev2[i] / (2 * sigma ^ 2))
        LL[i] = log(like[i])
    }
    
    tot_like = exp(sum(LL))
    return(tot_like)
}

Assign the prior distribution values in a grid

Prior distributions for our parameters are: \(L_{\infty}\) ~ \(U(40, 100)\) and \(K\) ~ \(U(0.1, 0.6)\)

Linf_grid = seq(40, 100, length.out = 1000)
K_grid = seq(0.1, 0.6, length.out = 1000)

# prior porbability
pLinf_grid = rep(1, length(Linf_grid))
pK_grid = rep(1, length(K_grid))
pTheta_matrix = expand.grid(Linf = pLinf_grid, K = pK_grid)

# expand the combination of grid 
Theta_matrix = expand.grid(Linf = Linf_grid, K = K_grid)

Calculate posterior probabilities

# negtive likelihood values of each grid
pDataGivenTheta = mapply(lognormal_like, Theta_matrix[,1], Theta_matrix[,2])

# Compute the posterior density
pData = sum(pDataGivenTheta * pTheta_matrix[, 1] * pTheta_matrix[, 2])
pThetaGivenData = (pDataGivenTheta * pTheta_matrix[, 1] * pTheta_matrix[, 2]) / pData

# arrange the results
result = cbind(Theta_matrix, pThetaGivenData)

Plot posterior distributions of parameters \(L_{\infty}\) and \(K\)

par(mfrow=c(2,1),mar=c(4,4,2,2))

# calculate the marginal probability
marginal_Linf = with(result, tapply(pThetaGivenData, Linf, sum))
plot(Linf_grid, marginal_Linf, type = "h", lwd = 3, xlab = expression(L[infinity]), ylab = bquote(paste("p(", theta, "|D)")))

# calculate the marginal probability
marginal_K = with(result, tapply(pThetaGivenData, K, sum))
plot(K_grid, marginal_K, type = "h", lwd = 3, xlab = "K", ylab = bquote(paste("p(", theta, "|D)")))