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、计算 MAEMAE <- mean(abs(actual_values - fitted_values))# 2、计算 RMSERMSE <- 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$heightmetrics <- 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]]$fitmetrics <- 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)}
}