Derivation of ISRP

# we simulate data from a multivariate normal setup with logistic mean and then will the estimates of r and K.

# Initial parameter set up

K = 100 # carrying capacity     

x0 = 10 # initial population size

r0 = 0.3 # initial intrinsic growth rate

t = seq(from = 1, to = 20, by = 1) # time

q= length(t)

h = numeric(length = length(t)-1) # difference between two consecutive time points

for(i in 1:(length(t)-1)){

  h[i] = t[i+1] - t[i]

}

mu = K/(1 + (K/x0 - 1)*exp(-r0*t))   # logistic mean function

par(mfrow = c(1,1))

par(mar = c(5.1, 5.1, 4.1, 2.1))

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.1

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

library(mvtnorm) # to generate data from multivariate setup

data = mvtnorm::rmvnorm(n, mean = mu, sigma = cov_mat)

matplot(t(data), type = "l", ylab = "Size",  xlab = "Time")

legend("topleft", legend = paste("n = ", n), bty = "n")

# section fro derivation of ISRP

d1 = matrix(data = NA, nrow = n, ncol = length(t)-2)

d2 = matrix(data = NA, nrow = n, ncol = length(t)-2)

for(j in 1:(length(t)-2)){

  d1[,j] = 1/(data[,j]) - 1/(data[,j+1])

  d2[,j] = 1/(data[,j+1]) - 1/(data[,j+2])

}

# Computation for r_hat

r_hat_vals = matrix(data = NA, nrow = n, ncol = length(t)-2)

for (j in 1:(length(t)-2)) {

  r_hat_vals[,j] = (1/h[j])*log(d1[,j]/d2[,j])

}

View(r_hat_vals)

# as sigma value increases in the derivation of ISRP estimates of r, higher NAN value will be produced. So now we will only consider those values where NAN is not present

selected_index = rep(FALSE, n)

for (i in 1:n) {

  if(sum(is.na(r_hat_vals[i,])) ==0)

    selected_index[i] = TRUE

}

selected_index    # List of selected rows for further calculations

sum(selected_index)

# Selected final series

final_data = data[selected_index, ]

matplot(t(final_data), type = "l", ylab = "Size", xlab = "Time")

legend("topleft", paste("n = ", sum(selected_index)), bty = "n")

par(mfrow=c(1,1))

par(mar = c(5.1, 5.1, 4.1, 2.1))

# Plot of r_hat vals

r_hat_vals = r_hat_vals[selected_index, ]

plot(t[1:(length(t)-2)], r_hat_vals[1,], col = "red", type = "p", 

     ylab = expression(widehat(r(Delta(t)))), xlab = "Time", 

     ylim = c(min(r_hat_vals),max(r_hat_vals)), lwd= 2)

for(i in 2:nrow(r_hat_vals)){

  lines(t[1:(length(t)-2)], r_hat_vals[i,], col = "red", type = "p", lwd = 2)

}

boxplot(r_hat_vals, col = "red")

# Plot of mean, median of ISRP of r

r_mean_vals = colMeans(r_hat_vals)

r_median_vals = apply(r_hat_vals, 2, median) 

lines(t[1:(length(t)-2)], r_mean_vals, type = "p", col = "black", lwd = 2)

lines(t[1:(length(t)-2)], r_median_vals, type = "p", col = "blue", lwd = 2)

legend("topleft", legend = c("ISRP of r", "mean of r", "median of r"), 

       col = c("red","black", "blue"), lwd = c(2,2,2))

# View of the density functions together

par(mfrow=c(1,1))

plot(density(r_hat_vals[,1]), col = 1, lwd = 1, lty = 1, 

     xlim=c(-5, 5), ylim= c(0,4.1),main = "Density plots of ISRP of r")

for(j in 2:(length(t)-2)){

  lines(density(r_hat_vals[,j]), col = j, lwd= 1, lty = j)

}

# Computation of K_hat_vals

K_hat_vals = matrix(data = NA, nrow = nrow(final_data), ncol = ncol(r_hat_vals))

for(j in 1:ncol(K_hat_vals)){

  K_hat_vals[,j] = (1/final_data[,1] - d1[selected_index, j]/(exp(-r_hat_vals[,j]*t[j])*

                          (1- exp(-r_hat_vals[,j]*h[j]))))^(-1)

}

# Plot of K_hat values

plot(t[1:(length(t)-2)], K_hat_vals[1,], col = "red", type = "p", 

     ylab = expression(hat(K)), xla = "time", 

     ylim = c(-20000, 20000))

for(i in 2:nrow(K_hat_vals)){

  lines(t[1:(length(t)-2)], K_hat_vals[i,], col = "red",type= "p")

}

boxplot(K_hat_vals, col = "red", ylim = c(-20000, 20000))

# Removing the outlier estimates

par(mfrow=c(1,2))

plot(t[1:(length(t)-2)], colMeans(r_hat_vals), col = "red", lwd=2, pch= 19, type = "b",

     ylab = expression(hat(r)), xlab = "time")

plot(t[1:(length(t)-2)], colMeans(K_hat_vals), col = "red", lwd=2, pch= 19, type = "b",

     ylab = expression(hat(K)), xlab = "time")

apply(r_hat_vals, 2, median)

apply(K_hat_vals, 2, median)

colMeans(K_hat_vals)

colMeans(r_hat_vals)

par(mfrow=c(3,3))

for (j in 1:(length(t)-2)) {

  plot(density(r_hat_vals[,j]), col = "red", main = paste("t = ", t[j]))

}

for (j in 1:(length(t)-2)) {

  plot(density(K_hat_vals[,j]), col = "red", main = paste("t = ", t[j]))

}