如何在R中为Kozak(2004)可变指数树削度模型添加随机效应

Kozak(2004)的可变指数削度函数公式如下：

diameter = b₀ × (dbh^b₁) × (height^b₂) × [(1 - hih^(1/3)) / (1 - p^(1/3))]x

其中指数x的计算式为：

x = b3 * hih^4 + 
    b4 * (1 / exp(dbh / height)) + 
    b5 * [(1 - hih^(1/3)) / (1 - p^(1/3))]^0.1 + 
    b6 * (1 / dbh) + 
    b7 * (height^(1 - hih^(1/3))) + 
    b8 * [(1 - hih^(1/3)) / (1 - p^(1/3))]

注：hih为相对高度（测量高度/树高），p为胸高相对比例（通常取1.3/height，对应胸径测量高度1.3m）

1. 加载所需包 #

使用nlme包拟合非线性混合效应模型，tidyverse用于数据处理（可选）：

# 首次运行时安装包
# install.packages(c("nlme", "tidyverse"))

# 加载包
library(nlme)
library(tidyverse)

2. 准备数据集 #

这里模拟一组符合模型结构的示例数据，可替换为你的实际数据：

set.seed(123)

# 生成100棵树，每棵树5个不同高度的直径测量值
n_trees <- 100
n_measurements <- 5

tree_data <- tibble(
  tree_id = rep(1:n_trees, each = n_measurements),
  dbh = rnorm(n_trees, mean = 25, sd = 5) %>% rep(each = n_measurements),
  height = rnorm(n_trees, mean = 20, sd = 3) %>% rep(each = n_measurements)
) %>%
  mutate(
    hih = runif(n(), min = 0.1, max = 1),
    p = 1.3 / height,
    # 用真实参数模拟x值
    x = 0.2 * hih^4 + 
        0.5 * (1 / exp(dbh / height)) + 
        0.3 * ((1 - hih^(1/3)) / (1 - p^(1/3)))^0.1 + 
        0.1 * (1 / dbh) + 
        0.4 * (height^(1 - hih^(1/3))) + 
        0.6 * ((1 - hih^(1/3)) / (1 - p^(1/3))),
    # 模拟带误差的直径数据
    diameter = 1.2 * (dbh^0.8) * (height^0.3) * ((1 - hih^(1/3)) / (1 - p^(1/3)))^x + 
               rnorm(n(), mean = 0, sd = 0.5)
  )

# 查看数据结构
glimpse(tree_data)

3. 定义Kozak削度模型函数 #

将非线性模型封装为nlme可识别的函数形式：

kozak_model <- function(hih, dbh, height, p, b0, b1, b2, b3, b4, b5, b6, b7, b8) {
  # 计算指数x
  x <- b3 * hih^4 + 
       b4 * (1 / exp(dbh / height)) + 
       b5 * ((1 - hih^(1/3)) / (1 - p^(1/3)))^0.1 + 
       b6 * (1 / dbh) + 
       b7 * (height^(1 - hih^(1/3))) + 
       b8 * ((1 - hih^(1/3)) / (1 - p^(1/3)))
  
  # 计算预测直径
  diameter_pred <- b0 * (dbh^b1) * (height^b2) * ((1 - hih^(1/3)) / (1 - p^(1/3)))^x
  return(diameter_pred)
}

4. 拟合非线性混合效应模型 #

示例中将b0设为随tree_id变化的随机效应，可根据需求调整随机效应参数：

mixed_model <- nlme(
  # 固定效应公式
  fixed = diameter ~ kozak_model(hih, dbh, height, p, b0, b1, b2, b3, b4, b5, b6, b7, b8),
  # 随机效应：b0随单棵树变化
  random = b0 ~ 1 | tree_id,
  # 初始参数值（可通过nls拟合固定效应模型获取更优初始值）
  start = c(b0 = 1.0, b1 = 0.8, b2 = 0.3, b3 = 0.2, b4 = 0.5, b5 = 0.3, b6 = 0.1, b7 = 0.4, b8 = 0.6),
  data = tree_data,
  # 调整迭代次数确保收敛
  control = nlmeControl(maxIter = 1000, pnlsMaxIter = 500)
)

# 查看模型结果
summary(mixed_model)

5. 模型诊断与预测 #

# 绘制残差拟合图，评估模型拟合效果
plot(mixed_model, resid(., type = "p") ~ fitted(.), abline = 0)

# 对新数据进行预测
new_data <- tibble(
  tree_id = 1,
  dbh = 25,
  height = 20,
  hih = seq(0.1, 1, by = 0.1),
  p = 1.3 / 20
) %>%
  mutate(pred_diameter = predict(mixed_model, newdata = .))

# 查看预测结果
print(new_data)

关键说明 #

• 随机效应调整：若需多个参数随分组变化，可修改random参数，例如random = b0 + b1 ~ 1 | tree_id表示b0和b1均随单棵树变化。
• 初始参数优化：非线性模型对初始值敏感，建议先用nls函数拟合固定效应模型，将得到的参数作为混合效应模型的初始值。
• 收敛问题：若模型不收敛，可尝试简化随机效应结构、调整初始参数或增加迭代次数。

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