The parameters of a growth model are an integral and important part of that model and ecologically the study of the dynamical properties of parameters is meaningful as it dictates the dynamical behavior of the studied model. We also conduct the sensitivity analysis using the Von Bertalanffy model as the test bed model and plot the parameter behavior in various time points t = 1, 4, 7, 10, 13, 16, and 18. From these figures, we observe that initially, the likelihood surface is more intrinsic growth rate (r) dominant but as time increases it becomes more K dominant. Therefore, estimates of ISRP may provide information about the sensitivity of the parameters for real data collected from the natural growth process. This approach may be useful when estimating models with a high degree of non-linearity and sufficient expert knowledge about the growth process is not available.
# Initial parameters setups
K = 100 # Carring capacity
x0 = 10 # Initial population size
r = 0.8 # Intrinsic growth rate
t = seq(from = 1, to = 20, by = 1) # Time points
q = length(t) # length of time points
# Required Packages
library(mvtnorm)
h = numeric(length = length(t)-1)
for(i in 1:(length(t)-1)){
h[i] = t[i+1] - t[i]
}
# Von Bertalanffy mean function
mu = K*((1 + ((x0/K)^(1/3) - 1)*exp(-r*(t/3)))^3)
#Visualization of Size profile over time
par(mfrow = c(1,1))
plot(t, mu, col = "red", lwd=2, type = "b", pch = 19,
ylab = "size, X(t)", xlab = "time (t)")
# Formation of the variance Co-variance matrix
sigma2 = 0.5 # variance of X_t
rho = 0.5 # corr(X_t, X_{t+1})
cov_mat = matrix(data = NA, nrow = length(t), ncol = length(t))
for (i in 1:length(t)) {
for (j in 1:length(t)) {
cov_mat[i,j] = sigma2*rho^(abs(i-j))
}
}
n = 100 # number of independent trajectories
# Generation of Multivariate normal data for n trajectories
data = mvtnorm::rmvnorm(n, mean = mu, sigma = cov_mat)
matplot(t, t(data), type = "l")
# Computation of likelihood function for overall estimates that will give us a single
# estimate of the parameters using the whole data at all time points 1 to q.
fun_likelihood = function(param){
r = param[1]
K = param[2]
sigma2 = param[3]
rho = param[4]
mu = K*((1 + ((x0/K)^(1/3) - 1)*exp(-r*(t/3)))^3)
cov_mat = matrix(data = NA, nrow = length(t), ncol = length(t))
for (i in 1:q) {
for (j in 1:q) {
cov_mat[i,j] = sigma2*rho^(abs(i-j))
}
}
exponent = 0
for(i in 1:n){
exponent = exponent + t(data[i,]-mu)%*%(solve(cov_mat))%*%(data[i,]-mu)
}
log_lik = - (n*q/2)*log(2*pi) - (n/2)*log(det(cov_mat)) - (1/2)*exponent
return(-log_lik)
}
param_0 = c(0.8, 100, 0.5, 0.5)
out = optim(param_0, fun_likelihood, method = "L-BFGS-B",
hessian = TRUE, lower = c(0.001, 50, 0.0001, -0.001),
upper = c(3, 150, 10, 0.999))
out$par
solve(out$hessian)
# Local likelihood surface
fun_local_likelihood = function(r,K){
sigma2 = 0.5
rho = 0.5
mu = K*((1 + ((x0/K)^(1/3) - 1)*exp(-r*(t/3)))^3) # VB Model
cov_mat = matrix(data = NA, nrow = length(t), ncol = length(t))
for (i in 1:q) {
for (j in 1:q) {
cov_mat[i,j] = sigma2*rho^(abs(i-j))
}
}
# Local localized estimates
mu_local = mu[ind:(ind + 2)]
data_local = data[,ind:(ind + 2)]
cov_mat_local = cov_mat[ind:(ind+2), ind:(ind+2)]
exponent = 0
for(i in 1:n){
exponent = exponent + t(data_local[i,] - mu_local)%*%(solve(cov_mat_local))%*%(data_local[i,] - mu_local)
}
log_lik = - (n*length(ind:(ind+2))/2)*log(2*pi) - (n/2)*log(det(cov_mat_local)) - (1/2)*exponent
return(-log_lik)
}
# Range of parameters for 3D plot
r = seq(0.01, 3, length = 10)
K = seq(50, 150, length = 10)
# 3D ploting
par(mfrow = c(1,1))
for (k in 1:(q-2)) {
ind = k
val = matrix(data = NA, nrow = length(r), ncol = length(K))
for(i in 1:length(r)){
for (j in 1:length(K)) {
val[i,j] = fun_local_likelihood(r[i], K[j])
}
}
persp(r,K,val, col = "yellow", main = paste("t=",ind), theta = 50, phi = 30)
}
# Likelihood surface plot at different time points which will indicate the sensitivity
# of the surface on parameters
# Surface plot for t=1
ind = 1
val = matrix(data = NA, nrow = length(r), ncol = length(K))
for(i in 1:length(r)){
for (j in 1:length(K)) {
val[i,j] = fun_local_likelihood(r[i], K[j])
}
}
persp(r,K,val, col = "grey", theta = 50, phi = 30, main = paste("t=",ind))
# Surface plot for t=4
ind = 4
val = matrix(data = NA, nrow = length(r), ncol = length(K))
for(i in 1:length(r)){
for (j in 1:length(K)) {
val[i,j] = fun_local_likelihood(r[i], K[j])
}
}
persp(r,K,val, col = "grey", main = paste("t=",ind),theta = 50, phi = 30)
# Surface plot for t=7
ind = 7
val = matrix(data = NA, nrow = length(r), ncol = length(K))
for(i in 1:length(r)){
for (j in 1:length(K)) {
val[i,j] = fun_local_likelihood(r[i], K[j])
}
}
persp(r,K,val, col = "grey", main = paste("t=",ind), theta = 50, phi = 30)
# Surface plot for t=10
ind = 10
val = matrix(data = NA, nrow = length(r), ncol = length(K))
for(i in 1:length(r)){
for (j in 1:length(K)) {
val[i,j] = fun_local_likelihood(r[i], K[j])
}
}
persp(r,K,val, col = "grey", main = paste("t=",ind),theta = 50, phi = 30)
# Surface plot for t=13
ind = 13
val = matrix(data = NA, nrow = length(r), ncol = length(K))
for(i in 1:length(r)){
for (j in 1:length(K)) {
val[i,j] = fun_local_likelihood(r[i], K[j])
}
}
persp(r,K,val, col = "grey", main = paste("t=",ind), theta = 50, phi = 30)
# Surface plot for t=16
ind = 16
val = matrix(data = NA, nrow = length(r), ncol = length(K))
for(i in 1:length(r)){
for (j in 1:length(K)) {
val[i,j] = fun_local_likelihood(r[i], K[j])
}
}
persp(r,K,val, col = "grey", main = paste("t=",ind), theta = 50, phi = 30)
# Surface plot for t=18
ind = 18
val = matrix(data = NA, nrow = length(r), ncol = length(K))
for(i in 1:length(r)){
for (j in 1:length(K)) {
val[i,j] = fun_local_likelihood(r[i], K[j])
}
}
persp(r,K,val, col = "grey", main = paste("t=",ind), theta = 50, phi = 30)