1.写在前面

以下代码记录了立地指数的计算过程，包括了优势树筛选、误差清理、非线性拟合以及结果成图。
优势树木确定以及数据清理过程：

相关导向函数：

2.实现代码

##*******************************************************************************----
##*******************************************************************************
## @ author:JAckson Zhao
#  @ time: 2024年8月23日17:34:07
# @ description：立地指数数据拟合
library(tidyverse)
library(mgcv)
library(dplyr)

setwd("C:\\Users\\YP\\Desktop\\Site index")

data <- read.csv("Final_data.csv", sep = ",", fileEncoding = "GBK")

# 获取前431行数据
data <- head(data, 430)
nrow(data)
# nihe <- data %>%
#   select(Hight, Age) %>%
#   rename(height = Hight, age = Age)
# nrow(nihe)

## 样地优势树高获取------------------------------------------------------------
# 处理数据
result <- data %>%
  group_by(Long, Lat, Site, PlotsID) %>%  # 根据 site 和 PlotsID 进行分组
  arrange(desc(height)) %>%    # 根据 height 降序排列
  slice(1:5) %>%              # 选取每组中最大的三个 height 值
  ungroup() %>%                # 取消分组
  group_by(Long, Lat, Site, PlotsID) %>%  # 再次分组
  summarise(
    avg_height = mean(height, na.rm = TRUE),  # 计算高度的均值
    avg_age = mean(Age, na.rm = TRUE),        # 计算对应的年龄的均值
    .groups = "drop"                        # 汇总后取消分组
  )

# 查看结果
head(result)
summary(result)

# 如果存在负值或零值，可能需要进行数据过滤
nihe <- result %>% filter(avg_height > 0 & avg_height < 23, avg_age > 0 & avg_age < 200) %>%
  rename(height = avg_height, age = avg_age)
head(nihe)
nrow(nihe)
summary(nihe)

# 对每个age组计算高度的均值和3倍标准差，并过滤掉超出这个范围的数据
nihe_clean <- nihe %>%
  group_by(age) %>%
  mutate(mean_height = mean(height, na.rm = TRUE),
         sd_height = sd(height, na.rm = TRUE)) %>%
  filter(height > (mean_height - 3 * sd_height), height < (mean_height + 3 * sd_height)) %>%
  ungroup()  # 移除分组，以便进行后续操作

# 查看清理后的数据
nrow(nihe_clean)
head(nihe_clean)
summary(nihe_clean)


# 绘制散点图，并添加拟合线
ggplot(nihe, aes(x = age, y = height)) + 
  geom_point() +  # 添加散点图层
  geom_smooth(method = "gam", formula = y ~ s(x),
              method.args = list(family = "gaussian"),
              color = "blue") + # 添加GAM拟合线
  labs(x = "林龄（年）", y = "群落高度（m）", title = "林龄与群落高度的关系") + 
  theme_minimal()  # 使用简洁主题

# 定义不同的非线性模型方程和初始参数---------------------------------------------
# 定义不同的非线性模型方程和初始参数
models <- list(
  list(formula = log(height) ~ a + b / log(age), start = list(a = 0.01, b = 1)),
  list(formula = log(height) ~ a + b / age, start = list(a = 0.01, b = 1)),
  list(formula = height ~ a * (1 - b * exp(-c * age)) ^ (1 / (1-d)), start = list(a = 15.618, b = 13.312, k = 1.255, d = 1)),
  list(formula = height ~ a * (1 - exp(-b * age)), start = list(a = 13.934, b = 0.114)),
  list(formula = height ~ a * (1 - exp(-b * age)^c), start = list(a = 14.531, b = 0.056, c = 1.304)),
  list(formula = height ~ a + b * age + I(age^2), start = list(a = 0.01, b = 1)),
  list(formula = height ~ a * exp(-b * exp( -c * age)), start = list(a = 13.668, b = 1.785, c = 0.182)),
  list(formula = height ~ a + b / log(age), start = list(a = 0.01, b = 1)),
  list(formula = height ~ a / (1 + b * exp(-c * age)), start = list(a = 16.848, b = 8.068, c = 0.182))
)

# 定义计算拟合优度的函数
calculate_fit_metrics <- function(fit, actual_values) {
  fitted_values <- fitted(fit)  # 计算预测值
  
  # 1、计算 MAE
  MAE <- mean(abs(actual_values - fitted_values))
  
  # 2、计算 RMSE
  RMSE <- sqrt(mean((actual_values - fitted_values)^2))
  
  # 3、计算普通的 R²
  SST <- sum((actual_values - mean(actual_values))^2)
  SSE <- sum((actual_values - fitted_values)^2)
  R_squared <- 1 - (SSE / SST)
  
  # 4、计算 Adjusted R²
  n <- length(actual_values)
  p <- length(coef(fit))
  Adjusted_R_squared <- 1 - ((1 - R_squared) * (n - 1) / (n - p - 1))
  
  return(list(MAE = MAE, RMSE = RMSE, R_squared = R_squared, Adjusted_R_squared = Adjusted_R_squared))
}

# 拟合每个模型并计算拟合优度
results <- lapply(models, function(model) {
  tryCatch({
    fit <- nls(
      model$formula,
      data = nihe,
      start = model$start,
      control = nls.control(maxiter = 100, minFactor = 1e-3)
    )
    actual_values <- if (grepl("log", deparse(model$formula))) log(nihe$height) else nihe$height
    metrics <- calculate_fit_metrics(fit, actual_values)
    list(fit = fit, metrics = metrics)
  }, error = function(e) {
    message("Error in fitting model: ", deparse(model$formula))
    NULL
  })
})

# 绘制模型拟合曲线的函数，并在图上显示 R2、MAE 和 RMSE
plot_model_fit <- function(model, data, actual_values, fitted_values, metrics, model_name) {
  p <- ggplot(data, aes(x = age)) +
    geom_point(aes(y = actual_values), color = "blue", size = 1.5) +
    geom_line(aes(y = fitted_values), color = "red", size = 1) +
    labs(
      title = paste("Model:", model_name),
      x = "Age",
      y = "Height"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(size = 14, family = "Times New Roman", face = "bold"),  # 修改标题的字体大小和字体样式
      axis.title.x = element_text(size = 12, family = "Times New Roman"),  # 修改 x 轴标签的字体大小和字体样式
      axis.title.y = element_text(size = 12, family = "Times New Roman")  # 修改 y 轴标签的字体大小和字体样式
    )
  
  # 添加 R², MAE, RMSE 到图上
  p <- p + annotate("text", x = Inf, y = Inf, label = sprintf("R²: %.2f\nMAE: %.2f\nRMSE: %.2f", metrics$R_squared, metrics$MAE, metrics$RMSE),
                    hjust = 1, vjust = 1, size = 3.5, color = "black", fontface = "bold")
  
  return(p)
}


# 遍历results列表，绘制每个成功拟合的模型，并显示指标
for (i in 1:length(results)) {
  if (!is.null(results[[i]])) {
    fit <- results[[i]]$fit
    metrics <- results[[i]]$metrics
    
    # 提取拟合值
    fitted_values <- fitted(fit)
    actual_values <- if (grepl("log", deparse(models[[i]]$formula))) log(nihe$height) else nihe$height
    
    # 绘制图形并显示指标
    p <- plot_model_fit(fit, nihe, actual_values, fitted_values, metrics, deparse(models[[i]]$formula))
    print(p)
  }
}