Testing of Multivariate Delta Method for the variance of r in the extended logistic model. It is crucial to note that we do not need to derive further interval estimators from the extended growth model, however, the simulation was carried out using the extended growth model. Based on the simulated data we shall evaluate the interval estimators for the logistic model only. The profile would no longer be the same as we observed in the logistic model. This gives us the indication that we need to consider an extended growth model.
# we simulate data from a multivariate normal setup with logistic mean and then will the estimates of r and K.
# Initial parameter setup
K = 100 # carrying capacity
x0 = 10 # initial population size
r0 = 0.3 # initial intrinsic growth rate
t = seq(from = 0, to = 19, by = 1) # time
h = numeric(length = length(t)-1) # difference between two consequtive time pionts
for(i in 1:(length(t)-1)){
h[i] = t[i+1] - t[i]
}
# extended logistic Model
c=1.5
mu = K/(1 + (K/x0 - 1)*exp(-(r0/c)*t^c)) # extended logistic mean function
par(mfrow = c(1,1))
par(mar = c(5.1, 5.1, 4.1, 2.1))
# Visualization of the size profile of extended logistic model.
plot(t, mu, col = "red", lwd=2, type = "b", pch = 19, ylab = "size, X(t)", xlab = "time (t)")
# Now we will create the covariance matrix which follows Koopmans structure
sigma2 = 0.001
rho = 0.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 = 1000 # number of trajectories
rep = 1000 # number of replication
# As the variation is present in r only, we are only interested to see the boxplot of r over time.
est_interval_r = matrix(data = NA, nrow = rep, ncol = length(t)-2)
for (i in 1:rep) {
library(mvtnorm)
data = mvtnorm::rmvnorm(n, mean = mu, sigma = cov_mat)
data_mean = colMeans(data)
for(j in 1:(length(t)-2)){
num = 1/data_mean[j] - 1/data_mean[j+1]
den = 1/data_mean[j+1] - 1/data_mean[j+2]
est_interval_r[i,j] = (1/h[j])*log(num/den)
}
}
# Visualization of the boxplot of r over time
boxplot(est_interval_r[,(1:15)], col = "grey", xlab=expression(bold(Time(t))), lwd=2,
ylab = expression(bold(widehat(r[j](Delta~t)))))
axis(1,font = 2); axis(2,font = 2)
# Logistic model with linearly growing r
c=0.5
mu = K/(1 + (K/x0 - 1)*exp(-(r0*t)*(1+c*t/2))) # linearly growing r mean function
par(mfrow = c(1,1))
par(mar = c(5.1, 5.1, 4.1, 2.1))
# Visualization of size profile of Logistic model with linearly growing r
plot(t, mu, col = "red", lwd=2, type = "b", pch = 19,
ylab = "size, X(t)", xlab = "time (t)")
# Now we will create the covariance matrix which follows Koopmans structure
sigma2 = 0.001
rho = 0.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 = 1000 # number of trajectories
rep = 1000 # number of replications
# As the variation is present in r only, we are only interested to see the boxplot of r over time.
est_interval_r = matrix(data = NA, nrow = rep, ncol = length(t)-2)
for (i in 1:rep) {
library(mvtnorm)
data = mvtnorm::rmvnorm(n, mean = mu, sigma = cov_mat)
data_mean = colMeans(data)
for(j in 1:(length(t)-2)){
num = 1/data_mean[j] - 1/data_mean[j+1]
den = 1/data_mean[j+1] - 1/data_mean[j+2]
est_interval_r[i,j] = (1/h[j])*log(num/den)
}
}
# Visualization of the boxplot of r over time
boxplot(est_interval_r[,(1:10)], col = "grey", xlab=expression(bold(Time(t))), lwd=2,
ylab = expression(bold(widehat(r[j](Delta~t)))))
axis(1,font = 2); axis(2,font = 2)