Distribution of ISRP
# First, we simulate data from a multivariate normal setup with logistic mean and then # will the estimates of r and K and then their distribution. We also check the accuracy # of the delta method in approximating the sampling distribution of the estimates of the # parameters.
# Initial parameters 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 consecutive time pionts
for(i in 1:(length(t)-1)){
h[i] = t[i+1] - t[i]
}
mu = K/(1 + (K/x0 - 1)*exp(-r0*t)) # Mean of the logistic model
par(mfrow = c(1,1))
par(mar = c(5.1, 5.1, 4.1, 2.1)) # for resizing the plot
# Plot of the size profile for 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.01 # variance of X(t)
rho = 0.1 # Correlation between X(t) and X(t+1)
cov_mat = matrix(data = NA, nrow = length(t), ncol = length(t)) # covariance matrix
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
# Define the Matrix to store estmated values of r
est_interval_r = matrix(data = NA, nrow = rep, ncol = length(t)-2)
# simulation of artificial data sets for n trajectories and 1000 replication with multivariate normal setup
for (i in 1:rep) {
library(mvtnorm)
data = mvtnorm::rmvnorm(n, mean = mu, sigma = cov_mat)
data_mean = colMeans(data)
# Inteval estimates of r
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 sampling distribution of the interval estimates of r by histogram
par(mfrow=c(3,3))
for (j in 1:(length(t)-2)) {
hist(est_interval_r[,j], probability = TRUE, col = "grey",
xlab = expression(r), main = paste("t = ", j), breaks = 30)
points(r0, 0, col = "red", pch = 19, cex = 2)
}
# Calculation of the mean and variance of the interval estimators of r.
# The analytical expression is given in Karim et al. (2013), which are calculated by uing the delta method.
# Calculation of Mean
delta_mean_r = numeric(length = length(t)-2)
for(i in 1:(length(t)-2)){
num = 1/mu[i] - 1/mu[i+1]
den = 1/mu[i+1] - 1/mu[i+2]
delta_mean_r[i] = (1/h[i])*log(num/den)
}
print(delta_mean_r)
# Calculation of Variance
delta_var_r = numeric(length = length(t)-2)
grad_phi = matrix(data = NA, nrow = 3, ncol = length(t)-2)
for (j in 1:(length(t)-2)) {
partial_1 = (-mu[j+1])/(mu[j]*(mu[j+1]-mu[j]))
partial_2 = (mu[j+2]-mu[j])/((mu[j+1]-mu[j])*(mu[j+2]-mu[j+1]))
partial_3 = (-mu[j+1])/(mu[j+2]*(mu[j+2]-mu[j+1]))
grad_phi[,j] = (1/h[j])*c(partial_1, partial_2, partial_3)
}
sub_cov_mat = cov_mat[1:3, 1:3]
for (j in 1:(length(t)-2)) {
delta_var_r[j] = t(grad_phi[,j])%*%sub_cov_mat%*%grad_phi[,j]
}
print(delta_var_r)
# Section for the histogram of estimates of r and calculation of estimated variance
centered_r_vals = matrix(data = NA, nrow = rep, ncol = length(t)-2)
for(j in 1:(length(t)-2)){
centered_r_vals[,j] = sqrt(n)*(est_interval_r[,j] - delta_mean_r[j])
}
for (j in 1:(length(t)-2)) {
hist(centered_r_vals[,j], probability = TRUE, col = "grey",
main = paste("t = ",j), xlab = expression(r), breaks = 30)
curve(dnorm(x, mean = 0, sd = sqrt(delta_var_r[j])), lwd=2,
add= TRUE, col = "red" )
}
# Now we will compare the estimated variance with the variance derived by using the delta method
par(mfrow=c(1,1))
par(mar = c(5.1, 5.1, 4.1, 2.1))
plot(apply(centered_r_vals[,(1:18)], 2, var),
type = "b", col = "red", lwd= 2, xlab=expression(bold(Time(t))),
ylab = expression(bold(widehat(r[j](Delta~t)))))
lines(delta_var_r, lwd=2, col = "grey")
legend("topleft", c("estimated variance", "delta method"),
col = c("red", "grey"), lwd= c(2,2), bty = "n")
axis(1,font = 2); axis(2,font = 2)
# Boxplot of estimates of r for logistic model
boxplot(est_interval_r[,(1:18)], 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)
# Define the Matrix to store estmated values of K
est_interval_K = matrix(data = NA, nrow = rep, ncol = length(t)-2)
# simulation of artificial data sets for n trajectories and 1000 replication with multivariate normal setup for K.
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 = ((data_mean[j+1] - data_mean[j])^((t[j]/h[j])+2))*((data_mean[j+2])^((t[j]/h[j])+1))
den = ((data_mean[j+2] - data_mean[j+1])^(t[j]/h[j]))*((data_mean[j])^((t[j]/h[j])+1))*
(data_mean[j+1])
den2= ((data_mean[j+2]*(data_mean[j+1] - data_mean[j])))-
(data_mean[j]*(data_mean[j+2]-data_mean[j+1]))
est_interval_K[i,j] = ((1/data_mean[1])-(num/(den*den2)))^(-1)
}
}
# Visualization of the sampling distribution of the interval estimates of K by histogram
par(mfrow=c(3,3))
for (j in 1:(length(t)-2)) {
hist(est_interval_K[,j], probability = TRUE, col = "grey",
xlab = expression(K), main = paste("t = ", j), breaks = 30)
points(K, 0, col = "red", pch = 19, cex = 2)
}
# Calculation of the mean and variance of the interval estimators of K
# The analytical expression is given in Karim et al. (2013), which are calculated by uing the delta method.
# Section for mean
delta_mean_K = numeric(length = length(t)-2)
for(i in 1:(length(t)-2)){
num = ((mu[i+1] - mu[i])^((t[i]/h[i])+2))*((mu[i+2])^((t[i]/h[i])+1))
den = ((mu[i+2] - mu[i+1])^(t[i]/h[i]))*((mu[i])^((t[i]/h[i])+1))*(mu[i+1])
den2= ((mu[i+2]*(mu[i+1] - mu[i])))-(mu[i]*(mu[i+2]-mu[i+1]))
delta_mean_K[i] = ((1/mu[1])-(num/(den*den2)))^(-1)
}
print(delta_mean_K)
# Section for variance
delta_var_K = numeric(length = length(t)-2)
grad_phi = matrix(data = NA, nrow = 3, ncol = length(t)-2)
for (j in 1:(length(t)-2)) {
num = ((mu[i+1] - mu[i])^((t[i]/h[i])+2))*((mu[i+2])^((t[i]/h[i])+1))
den = ((mu[i+2] - mu[i+1])^(t[i]/h[i]))*((mu[i])^((t[i]/h[i])+1))*(mu[i+1])
den2= ((mu[i+2]*(mu[i+1] - mu[i])))-(mu[i]*(mu[i+2]-mu[i+1]))
partial_1 = -(((t[j]/h[j])+2)/(mu[j+1]-mu[j]))-(((t[j]/h[j])+1)/mu[j])+
((2*mu[j+2]-mu[j+1])/den2)
partial_2 = (((t[j]/h[j])+2)/(mu[j+1]-mu[j]))+((t[j]/h[j])/(mu[j+2]-mu[j+1]))-(1/mu[j+1])-
((mu[j+2]+mu[j])/den2)
partial_3 = (((t[j]/h[j])+1)/(mu[j+2]))-((t[j]/h[j])/(mu[j+2]-mu[j+1]))+
((2*mu[j]-mu[j+1])/den2)
zi=num/(den*den2)
phi=((1/mu[1])-(num/(den*den2)))^(2)
grad_phi[,j] = (zi/phi)*c(partial_1, partial_2, partial_3)
}
sub_cov_mat = cov_mat[1:3, 1:3]
for (j in 1:(length(t)-2)) {
delta_var_K[j] = t(grad_phi[,j])%*%sub_cov_mat%*%grad_phi[,j]
}
print(delta_var_K)
# Section for the histogram of estimates of K and calculation of estimated variance
centered_K_vals = matrix(data = NA, nrow = rep, ncol = length(t)-2)
for(j in 1:(length(t)-2)){
centered_K_vals[,j] = sqrt(n)*(est_interval_K[,j] - delta_mean_K[j])
}
for (j in 1:(length(t)-2)) {
hist(centered_K_vals[,j], probability = TRUE, col = "grey",
main = paste("t = ",j), xlab = expression(r), breaks = 30)
curve(dnorm(x, mean = 0, sd = sqrt(delta_var_K[j])), lwd=2,
add= TRUE, col = "red" )
}
# Now we will compare the estimated variance with the variance derived by using the delta method
par(mfrow=c(1,1))
par(mar = c(5.1, 5.1, 4.1, 2.1))
plot(apply(centered_K_vals, 2, var),
type = "b", col = "red", lwd= 2, xlab=expression(bold(Time(t))),
ylab = expression(bold(widehat(K[j](Delta~t)))))
lines(delta_var_K, lwd=2, col = "grey")
legend("topleft", c("estimated variance", "delta method"),
col = c("red", "grey"), lwd= c(2,2), bty = "n")
axis(1,font = 2); axis(2,font = 2)
# Boxplot of estimates of K for logistic model
boxplot(est_interval_K, col = "grey", xlab=expression(bold(Time(t))), lwd=2,
ylab = expression(bold(widehat(K[j](Delta~t)))))
axis(1,font = 2); axis(2,font = 2)