基于R语言的高斯函数模型：藤壶密度与深度、珊瑚盖度关联分析

Hey there! Let's break down how to implement and validate a Gaussian function model for your barnacle density-depth relationship, since you’ve already tested linear and negative binomial models (m2 <- glm.nb(dens.cm ~ depth + coral.cover+location+depth:location, data =data)) and suspect a Gaussian pattern is a better fit.

1. First: Clarify the Gaussian Model Form #

A standard Gaussian curve tailored to your data will take this form:
dens.m2 = a * exp(-((depth - b)^2)/(2*c^2)) + d
Where:

• a: The peak height of barnacle density (difference between peak and baseline levels)
• b: The depth where barnacle density is highest (the curve's peak position)
• c: Controls the width of the curve (smaller values mean a narrower, sharper peak)
• d: Baseline barnacle density at depths far from the peak

If you want to incorporate your other variables (coral cover, location), you can extend this model to let parameters like b (peak depth) or a (peak height) vary with these covariates—we’ll cover that later.

2. Fit the Gaussian Model in R #

First, make sure your density units are consistent (you mentioned "每平方米藤壶密度" but your existing model uses dens.cm—convert if needed):

# Convert density from per cm² to per m² (1 m² = 10000 cm²)
data$dens.m2 <- data$dens.cm * 10000

Start with a simple univariate Gaussian model (focused on depth vs density) before adding covariates. You’ll need initial parameter guesses—this is critical for non-linear model convergence:

# Plot your data to estimate starting values visually
plot(dens.m2 ~ depth, data = data, pch = 16, main = "Barnacle Density vs Depth")

# Set initial values based on the plot's patterns
start_vals <- list(
  a = max(data$dens.m2, na.rm = TRUE),  # Estimate of peak density
  b = median(data$depth, na.rm = TRUE), # Midpoint depth as initial peak position
  c = sd(data$depth, na.rm = TRUE),     # Spread of depth values
  d = min(data$dens.m2, na.rm = TRUE)   # Baseline density at extreme depths
)

# Fit the model with base R's nls()
gauss_model <- nls(
  dens.m2 ~ a * exp(-((depth - b)^2)/(2*c^2)) + d,
  data = data,
  start = start_vals
)

# Check the model output
summary(gauss_model)

If you hit convergence issues, use nlsLM() from the minpack.lm package—it’s more robust for tricky non-linear fits:

library(minpack.lm)

gauss_model_robust <- nlsLM(
  dens.m2 ~ a * exp(-((depth - b)^2)/(2*c^2)) + d,
  data = data,
  start = start_vals
)
summary(gauss_model_robust)

Extend to Include Covariates

To incorporate location and coral cover, modify the model to let parameters vary with these variables. For example, let peak depth differ by location and baseline density depend on coral cover:

# Updated starting values for the extended model
start_vals_extended <- list(
  a = max(data$dens.m2, na.rm = TRUE),
  b_site1 = median(data$depth[data$location == "Site1"], na.rm = TRUE),
  b_site2_diff = median(data$depth[data$location == "Site2"], na.rm = TRUE) - median(data$depth[data$location == "Site1"], na.rm = TRUE),
  c = sd(data$depth, na.rm = TRUE),
  d_baseline = min(data$dens.m2, na.rm = TRUE),
  d_coral_effect = 0.1  # Initial guess for coral cover's impact on baseline density
)

# Fit the extended model
gauss_model_extended <- nlsLM(
  dens.m2 ~ a * exp(-((depth - (b_site1 + b_site2_diff*(location == "Site2")))^2)/(2*c^2)) + (d_baseline + d_coral_effect*coral.cover),
  data = data,
  start = start_vals_extended
)

summary(gauss_model_extended)

3. Compare Gaussian Model to Your Existing Negative Binomial Model #

To confirm the Gaussian model is a better fit, use AIC (lower values mean better fit) and visualize both models:

# Refit your negative binomial model with m² units for fair comparison
m2_m2 <- glm.nb(dens.m2 ~ depth + coral.cover + location + depth:location, data = data)

# Compare AIC values
AIC(m2_m2, gauss_model)

# Visualize both model fits
depth_seq <- seq(min(data$depth), max(data$depth), length.out = 100)

# Predict from negative binomial model (using Site1 and mean coral cover)
pred_nb <- predict(m2_m2, newdata = data.frame(
  depth = depth_seq,
  coral.cover = mean(data$coral.cover),
  location = "Site1"
), type = "response")

# Predict from Gaussian model
pred_gauss <- predict(gauss_model, newdata = data.frame(depth = depth_seq))

# Plot the raw data and both fits
plot(dens.m2 ~ depth, data = data, col = ifelse(data$location == "Site1", "blue", "red"), pch = 16)
lines(depth_seq, pred_nb, col = "blue", lwd = 2)
lines(depth_seq, pred_gauss, col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("Negative Binomial (Site1)", "Gaussian Model"), col = c("blue", "red"), lty = c(1,2), lwd = 2)

This visual comparison will make it clear how well the Gaussian curve captures the density-depth pattern compared to the linear-like fit of the negative binomial model.

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