重复测量逻辑增长曲线建模及2018目标拟合技术咨询

Great question! You absolutely should combine your 2016 and 2017 data into a single model instead of fitting separate ones. Here's why, plus a practical implementation to hit all your goals (incorporating repeated measures, target asymptotes, and month-to-month variability confidence intervals):

一、Why merging is better #

• More robust estimates: Combining two years of data gives you 24 observations instead of 12, which reduces parameter uncertainty and makes your model more reliable.
• Capture year-to-year variability: A merged model lets you explicitly account for differences between 2016 and 2017 (via random effects or fixed group terms), which is exactly what you need to generate meaningful confidence intervals for month-to-month trends.
• Aligns with your repeated measures framing: Treating each year as a "group" of repeated monthly observations fits perfectly with your goal of modeling cumulative growth across time.

二、Model choice: Mixed-effects logistic growth model #

Your cumulative population data follows a classic S-shaped (logistic) growth curve, which matches the model you started with. For the merged dataset, we'll use a mixed-effects nonlinear model (via the nlme package) to balance average trends and year-to-year variability.

2.1 Fixed asymptote (tied to 2018 target) #

If you want to force the model's asymptote to your 2018 year-end target (e.g., let's say target_pop = 16000), use this code:

library(nlme)

# Define your 2018 year-end target population
target_pop <- 16000

# Fit mixed-effects logistic model
mixed_fit <- nlme(
  Total ~ target_pop / (1 + exp(-(beta0 + beta1*Time + (b0 + b1*Time)|Year))),
  data = df_tot,
  fixed = beta0 + beta1 ~ 1,  # Average intercept and slope across years
  random = b0 + b1 ~ 1,       # Random variation in intercept/slope by year
  start = list(fixed = c(beta0 = -1.5, beta1 = 0.4)),  # Initial guesses (use your logit model results here!)
  control = nlmeControl(maxIter = 100)
)

• beta0/beta1: The average starting point and growth rate across 2016-2017
• (b0 + b1*Time)|Year: Captures how 2016 and 2017 differ from the average trend (this drives your variability-based confidence intervals)
• The asymptote is locked to target_pop, so the model directly reflects your 2018 goal.

2.2 Unconstrained asymptote (estimate first, adjust later) #

If you want to first estimate the natural asymptote from 2016-2017 data before setting it to your 2018 target, use this unconstrained version:

unconstrained_fit <- nlme(
  Total ~ phi1 / (1 + exp(-(beta0 + beta1*Time + (b0 + b1*Time)|Year))),
  data = df_tot,
  fixed = phi1 + beta0 + beta1 ~ 1,  # Estimate asymptote, intercept, slope
  random = b0 + b1 ~ 1,
  start = list(fixed = c(phi1 = 15000, beta0 = -1.5, beta1 = 0.4)),
  control = nlmeControl(maxIter = 100)
)

# Extract the estimated asymptote to inform your 2018 target
estimated_asymptote <- coef(unconstrained_fit)[["phi1"]]

三、Generating confidence intervals for month-to-month variability #

The mixed-effects model automatically accounts for year-to-year variation, so you can generate prediction intervals that reflect this uncertainty:

# Create a dataset for predictions (include 2018 to test your target)
pred_data <- expand.grid(Year = c(2016, 2017, 2018), Time = 1:12)

# Get predicted values + 95% confidence intervals
preds <- predict(mixed_fit, newdata = pred_data, interval = "confidence", level = 0.95)
pred_data <- cbind(pred_data, preds)

# Check 2018's predicted range to compare with your actual 2018 monthly data
subset(pred_data, Year == 2018)

These intervals will show you the expected range of cumulative population for each month, accounting for how much 2016 and 2017 varied from each other.

四、Comparing to your single-year model #

Your original 2016-only nls fit works, but it has key limitations:

• It only uses 12 observations, so parameter estimates are much less precise.
• It can't capture year-to-year variability, so you can't generate meaningful intervals to judge 2018's performance.
• It doesn't let you directly compare 2016 and 2017 trends.

The merged model fixes all these issues while staying true to your repeated-measures framing.

五、Assessing if 2018 can hit a higher target #

Once you have partial 2018 data, here's how to use the model to evaluate progress:

1. Add your 2018 monthly data to df_tot (create a new dataframe like df_combined <- rbind(df_tot, df_2018)).
2. Refit the model with your higher target (e.g., higher_target = 16500).
3. Compare your actual 2018 data points to the model's predicted confidence intervals:
   • If your actual values fall within or above the interval, you're on track to hit the higher target.
   • You can also predict the 2018 year-end value directly:

     higher_target <- 16500
     
     # Refit model with higher target and 2018 data
     mixed_fit_higher <- nlme(
       Total ~ higher_target / (1 + exp(-(beta0 + beta1*Time + (b0 + b1*Time)|Year))),
       data = df_combined,
       fixed = beta0 + beta1 ~ 1,
       random = b0 + b1 ~ 1,
       start = list(fixed = c(beta0 = -1.5, beta1 = 0.4)),
       control = nlmeControl(maxIter = 100)
     )
     
     # Predict 2018 year-end population
     predict(mixed_fit_higher, newdata = data.frame(Year = 2018, Time = 12))
     

If the predicted year-end value is close to or exceeds higher_target, you're in good shape.

Quick tip: Use your logit model for initial values #

Your lm(logit(Norm) ~ Time) approach is a great way to get reliable starting values for the nonlinear model. Just extract the coefficients:

logit_fit <- glm(Norm ~ Time, data = df_tot, family = binomial(link = "logit"))
start_beta0 <- coef(logit_fit)[1]
start_beta1 <- coef(logit_fit)[2]

Plug these into the start parameter of your nlme model to improve convergence.

内容的提问来源于stack exchange，提问作者zcs