Number of Horses-Mules Data

We have obtained the data from Banks, 1994 (Chapter 5) and plotted the size profile which suggests that the logistic model is the base model. We first compute the ISRP of K considering the logistic equation as the base model. We observe that the ISRP of K is a decreasing function of time. So, we consider three different variations of K, namely, linear (K(t) = K0 [1 +ct]), hyperbolic (K(t ) = K0 / [1 + ct]), and exponential decay (K(t) = K0 exp[-ct]). Based on the AIC, the linear function has been found to be the best model. Since the analytical solution is not available for this growth equation, we estimate the parameters r, K0, and c by using nonlinear regression, where the empirical RGR has been used as the response variable. The estimates are obtained as r = 0.06046851, K0 = 40.03490149, and c = −0.01062441 and all the parameters are statistically significant. To generate the final fitted growth equation, we utilize the deSolve packages from R. It is important to note that we have also checked for the variation in r using the ISRP of r with the logistic equation as the basic model. A linear regression of ISRP of r on time points gave a statistically insignificant slope parameter. Therefore, we treat r as constant for this particular data.

x = c(7, 8.270, 10.881, 12.170, 14.802, 17.518, 20.557, 20.004,

      22.077, 24.211, 26.490, 25.742,22.569, 19.124, 16.683, 

      14.478, 11.950, 7.781, 4.309, 3.089)

t = seq(0, 95, by = 5)

plot(t, x, lwd=2, col = "magenta", ylab= "N", xlab= "Time (Years)",

     cex=2, type = "p", ylim = c(0,30))  # plot of size profile

n = length(x)

isrp_r = numeric(length = length(t)-2)

h = t[2]-t[1]

num = 1/x[1:(n-2)] - 1/x[2:(n-1)]

den = 1/x[2:(n-1)] - 1/x[3:n]

isrp_r = (1/h)* log(num/den)

ind = is.na(isrp_r)

isrp_r[ind] = 0

plot( isrp_r, type = "p", pch = 19)

fit = lm(isrp_r ~ t[1:(n-2)])

summary(fit)

abline( fit, col = "red", lwd = 2)

denom = exp(-isrp_r*t[1:(n-2)])*(1- exp(-isrp_r*h))

isrp_K = (1/x[1] - num/denom)^(-1)

plot(t[1:(n-2)], isrp_K[1:(n-2)], type = "p", pch = 19, col = "magenta",

     ylab = expression(widehat(K[j](Delta* t))), xlab = "time")

data = data.frame(K = isrp_K, t = t[1:(length(t)-2)])

fit1 = nls(K ~ K0*(1+c*t), data = data, 

           start = list(K0 = 40, c=-0.01) )

K0 = coef(fit1)[1]; c = coef(fit1)[2]

curve(K0*(1+c*x), lwd = 2, add = TRUE, col = "red", lty = 2)

fit2 = nls(K ~ K0/(1+c*t), data = data, 

           start = list(K0 = 40, c=-0.01) )

K0 = coef(fit2)[1]; c = coef(fit2)[2]

curve(K0/(1+c*x), lwd = 2, add = TRUE, lty = 2,

      col = "blue")

fit3 = nls(K ~ K0*exp(c*t), data = data, 

           start = list(K0 = 40, c=-.1) )

K0 = coef(fit2)[1]; c = coef(fit2)[2]

curve(K0*exp(-c*x), lwd = 2, add = TRUE, lty = 2,

      col = "black")

AIC(fit1)

AIC(fit2)

AIC(fit3)

legend("topright", c(expression(K(t)==K[0]*(1+c*t)), expression(K(t)==K[0]/(1+c*t)),

expression(K(t)==K[0]*e^(-c*t))), bty = "n",

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

RGR = (1/h)*log(x[2:length(x)]/x[1:(length(x)-1)])

d = data.frame(RGR = RGR, t = t[1:(length(t)-1)], x = x[1:(length(x)-1)])

fit = nls(RGR ~ a0*(1-x/(K0*(1+c*t))), data = d,

          start = list(a0 = 1, K0 = 40, c = -0.01))

coef(fit)

a0 = coef(fit)[1]

K0 = coef(fit)[2]

c = coef(fit)[3]

library(deSolve)

f = function(t, y, parms){

  list(a0*y*(1-y/(K0*(1+c*t))))

}

yini = x[1]

times = 0:95

parms = c(a0 = a0, K0 = K0, c = c)

out = ode(times = times, y = yini, parms = parms, func = f,

          method = "rk4")

coef(fit)

head(out)

plot(times, out[,2], type = "l", col = "red", xlab = expression(Time), 

     lwd = 2, ylab= "N", ylim = c(0,30))

points(t,x, col = "magenta", pch = 19, cex = 1)