Logistic model prediction

Isle Royale National Park has been the site of a long-term ecological study tracking the population of moose and wolves. The data spans several decades and shows how the moose population fluctuates in response to predation, food availability, and other environmental factors.

Here’s a simplified version of the Moose population data on Isle Royale, which I’ll use for this example. The actual data spans several decades, but for simplicity, let’s use a subset of the data:

# Year (from 1959 to 1995)
year <- c(1959, 1965, 1970, 1975, 1980, 1985, 1990, 1995)

# Moose Population (in hundreds)
moose_population <- c(6.3, 8.0, 15.2, 18.2, 22.7, 14.2, 12.1, 16.5)

data <- data.frame(
  year,
  moose_population
)


logistic_model <- function(t, K, r, P0) {
  K / (1 + ((K - P0) / P0) * exp(-r * t))
}


# Time (in years since 1959)
time <- year - 1959

# Initial guesses for K, r, and P0
K_start <- 25  # an initial guess for the carrying capacity in hundreds
r_start <- 0.1
P0_start <- moose_population[1]

# Fitting the logistic model
fit <- nls(moose_population ~ logistic_model(time, K, r, P0_start),
           start = list(K = K_start, r = r_start))

# Summary of the fitted model
summary(fit)
## 
## Formula: moose_population ~ logistic_model(time, K, r, P0_start)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## K  16.7726     1.9989   8.391 0.000156 ***
## r   0.2125     0.1318   1.612 0.158056    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.75 on 6 degrees of freedom
## 
## Number of iterations to convergence: 12 
## Achieved convergence tolerance: 3.774e-06
# Extracting estimated parameters
K_est <- coef(fit)["K"]
r_est <- coef(fit)["r"]

K_est
##        K 
## 16.77263
r_est
##         r 
## 0.2124951
# Plotting the data and the fitted model
plot(time, moose_population, main="Logistic Growth Model Fit for Moose on Isle Royale", 
     xlab="Time (Years since 1959)", ylab="Moose Population (Hundreds)", pch=19)
curve(logistic_model(x, K_est, r_est, P0_start), add=TRUE, col="blue", lwd=2)

model = nls(moose_population ~ SSlogis(time, a, b, c), data = data)
lines(time, predict(model), col="red", lwd=2)
summary(model)
## 
## Formula: moose_population ~ SSlogis(time, a, b, c)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a   16.817      2.034   8.269 0.000422 ***
## b    3.965      3.136   1.264 0.261804    
## c    3.889      3.039   1.279 0.256894    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.02 on 5 degrees of freedom
## 
## Number of iterations to convergence: 2 
## Achieved convergence tolerance: 1.785e-06
legend("bottomright", legend=c("Own function", "nls model"),
       col=c("blue", "red"), lwd=2, lty=c(1, 2))