Varying r in Logistic Model

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