如何为SITAR拟合得到的均值生长曲线添加95%置信区间

Great question! The default plot() function for sitar models doesn’t include built-in support for overlaying 95% confidence intervals on the mean curves, but we can easily calculate and add them manually using the model’s prediction tools. Here’s a complete, reproducible solution:

Step 1: Fit Your SITAR Model (Your Existing Code) #

First, run your original model fitting code to ensure we have a valid model object:

library(sitar)
data <- na.omit(berkeley[berkeley$sex == 2 & berkeley$age >= 8 & berkeley$age <= 18, c('id', 'age', 'height')])
sitar_model <- sitar(x = age, y = height, id = id, data = data, df = 5)

Step 2: Create a Dense Age Sequence for Smooth Curves #

To generate smooth confidence interval lines, we’ll create a dense sequence of age values spanning your data’s range:

# 100 evenly spaced age points between 8 and 18
new_ages <- seq(min(data$age), max(data$age), length.out = 100)

Step 3: Calculate Mean Predictions & 95% Confidence Intervals #

We’ll use the predict() function to get mean predictions and their standard errors, then compute the 95% CI using the normal approximation (mean ± 1.96 × standard error).

For the Mean Growth Curve (Distance Curve, opt='d') #

# Predict mean height and standard errors
pred_dist <- predict(sitar_model, newdata = data.frame(age = new_ages),
                     type = "response", se.fit = TRUE)
# Compute 95% confidence bounds
ci_dist_lower <- pred_dist$fit - 1.96 * pred_dist$se.fit
ci_dist_upper <- pred_dist$fit + 1.96 * pred_dist$se.fit

For the Mean Velocity Curve (opt='v') #

# Predict mean height velocity and standard errors
pred_vel <- predict(sitar_model, newdata = data.frame(age = new_ages),
                    type = "velocity", se.fit = TRUE)
# Compute 95% confidence bounds
ci_vel_lower <- pred_vel$fit - 1.96 * pred_vel$se.fit
ci_vel_upper <- pred_vel$fit + 1.96 * pred_vel$se.fit

Step 4: Plot Curves with Confidence Intervals #

Now we’ll plot the original mean curves, then overlay the confidence interval lines.

Distance Curve with 95% CI #

par(mar = c(4,4,1,1) + 0.1, cex = 0.8)
# Plot the base mean growth curve
plot(sitar_model, opt = 'd', las = 1, apv = TRUE)
# Add lower and upper CI lines (gray dashed for visibility)
lines(new_ages, ci_dist_lower, col = "gray50", lty = 3)
lines(new_ages, ci_dist_upper, col = "gray50", lty = 3)
# Optional: Add a legend for clarity
legend("topleft", legend = c("Mean Growth Curve", "95% Confidence Interval"),
       col = c("black", "gray50"), lty = c(1, 3), cex = 0.8)

Velocity Curve with 95% CI #

# Plot the base mean velocity curve
plot(sitar_model, opt = 'v', las = 1, apv = TRUE, lty = 2)
# Add lower and upper CI lines
lines(new_ages, ci_vel_lower, col = "gray50", lty = 3)
lines(new_ages, ci_vel_upper, col = "gray50", lty = 3)
# Optional: Add a legend
legend("topleft", legend = c("Mean Velocity Curve", "95% Confidence Interval"),
       col = c("black", "gray50"), lty = c(2, 3), cex = 0.8)

Quick Notes #

• The 95% CI here uses the normal approximation (1.96 × SE), which is standard for large datasets. If you need more robust intervals (e.g., bootstrap-based), you can use the boot() function with your sitar model to generate resampled predictions, but this method works well for most visualization purposes.
• Feel free to tweak the col (color) and lty (line type) parameters to match your preferred style.

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