基于nlme拟合的线性混合效应模型，按太阳时间绘制预测移动速率

模型背景与绘图需求 #

已使用nlme包拟合带有corCAR1自相关结构的线性混合效应模型：

• 响应变量：对数转换的动物移动速率（单位：m/h）
• 核心解释变量：
  • 二水平因子：study_area（研究区域）、season（季节）
  • 连续变量：dist_road（距最近道路距离，单位：米）
  • 弧度制太阳时间sun_time：采用Richter等人（2020）的三角函数结构（sin(sun_time)+cos(sun_time)+sin(2*sun_time)+cos(2*sun_time)）建模时间效应
• 模型结构：包含个体ID的随机截距，以及study_area与dist_road、太阳时间项、season的双向交互项

绘图需求：

• 绘制太阳时间（关键刻度为0、π/2、π、3π/2、2π）对原始尺度移动速率的预测曲线（需将对数预测值反变换）
• 分别展示dist_road=0米和500米、各季节下的预测结果

示例数据与模型拟合代码 #

library(nlme)

df <- data.frame(
  ID = rep(1:3, each = 30),
  timestep = rep(0:29, 3),
  sun_time = c(
    1.23, 1.64, 1.81, 1.98, 2.15, 2.31, 2.48, 2.65, 2.81, 2.98, 3.15, 3.32, 3.48, 3.65, 4.15, 4.32, 4.49, 4.66, 5.11, 5.72,
    0.03, 0.64, 1.25, 1.65, 1.82, 1.98, 2.15, 2.32, 2.48, 2.65, 0.49, 1.00, 1.52, 1.73, 1.90, 2.08, 2.26, 2.43, 2.60, 2.78, 2.96, 3.13, 3.31, 3.49, 3.66, 3.83, 4.01, 4.36, 4.54,
    4.71, 5.22, 5.74, 6.25, 0.49, 0.99, 1.50, 1.72, 1.90, 2.07, 2.25, 5.27, 5.95, 0.51, 1.21, 1.69, 1.82, 1.97, 2.15, 2.60, 2.77, 2.92, 3.54, 3.72, 3.86, 4.00, 4.18, 4.33, 4.49, 4.62, 5.25,
    5.92, 0.44, 1.24, 1.97, 2.15, 2.28, 2.77, 2.92, 3.08, 3.24
  ),
  mov_rate_log = c(3.54, 3.53, 3.52, 0.00, 5.86, 4.65, 0.00, 4.36, 3.74, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 3.77, 5.81, 7.15, 6.78, 0.00,
                   0.00, 5.22, 6.56, 6.88, 7.31, 0.00, 0.00, 0.00, 0.00, 5.33, 5.91, 7.20, 4.60, 3.00, 4.53, 7.58, 6.95, 4.43, 4.70, 6.51, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2.90, 0.00, 3.59, 7.30,
                   5.59, 5.71, 6.59, 7.23, 6.97, 6.91, 7.01, 0.00, 0.00, 0.00, 5.05, 7.14, 4.60, 7.32, 7.10, 3.58, 3.95, 4.44, 2.78, 0.00, 5.57, 0.00, 0.00, 0.00, 6.71, 6.42, 7.40, 8.07, 6.32, 6.51,
                   0.00, 5.84, 7.57, 0.00, 4.41, 3.00, 0.00, 4.08, 4.19, 0.00),
  season = rbinom(90, 1, 0.5),
  study_area = rbinom(90, 1, 0.5),
  dist_road = c(329, 334, 315, 290, 286, 521, 455, 457, 408, 389, 397, 396, 391, 392, 393, 384, 342, 389, 203, 73, 78, 71,  31, 143, 154, 204,
                204, 204, 201, 202, 310, 159, 175, 242, 233, 231, 268, 101, 124,  78, 298, 291, 298, 295, 301, 299, 304, 287, 291, 293, 107, 117, 199, 98, 140, 189,
                173, 121, 110, 114, 137,  90, 548, 602, 246,  88, 122,  73, 199, 196, 198, 278, 269, 272, 261, 607,   4, 653, 262, 359, 390, 387, 272, 273, 271, 350,
                340, 332, 307, 344)
)

df$season <- as.factor(df$season)
df$study_area <- as.factor(df$study_area)

model <- lme(
  mov_rate_log ~ study_area * (
    sin(sun_time) + cos(sun_time) + sin(2 * sun_time) + cos(2 * sun_time)
  ) + study_area * dist_road + study_area * season,
  random = list(ID = ~ 1),
  control = lmeControl(msMaxIter = 100),
  data = df,
  method = "REML",
  correlation = corCAR1(form = ~ timestep | ID)
)

summary(model)

解决方案：绘制预测曲线 #

步骤1：加载依赖包并构建预测数据集 #

生成包含所有要展示的变量组合的数据集，同时生成足够密集的sun_time序列以绘制平滑曲线：

library(ggplot2)

# 生成太阳时间序列：0到2π，间隔0.05弧度
sun_seq <- seq(0, 2*pi, by = 0.05)
# 定义要展示的分组变量组合
pred_df <- expand.grid(
  sun_time = sun_seq,
  dist_road = c(0, 500),
  season = levels(df$season),
  study_area = levels(df$study_area),
  timestep = 0,  # 固定timestep，不影响边际预测
  ID = 1         # 虚拟ID，用于predict函数识别分组
)

步骤2：生成预测值并反变换 #

使用predict()函数获取群体水平（不含随机效应）的对数预测值，再通过指数变换还原为原始移动速率：

# 生成对数尺度预测值，level=0表示群体水平
pred_df$mov_rate_log_pred <- predict(model, newdata = pred_df, level = 0)
# 反变换为原始速率（指数函数）
pred_df$mov_rate_pred <- exp(pred_df$mov_rate_log_pred)

步骤3：可视化预测曲线 #

使用ggplot2绘制曲线，通过分面区分研究区域，颜色区分季节，x轴设置指定的太阳时间刻度：

# 定义x轴刻度与标签
sun_ticks <- c(0, pi/2, pi, 3*pi/2, 2*pi)
sun_labels <- c("0", "π/2\n(日出)", "π", "3π/2\n(日落)", "2π")

ggplot(pred_df, aes(x = sun_time, y = mov_rate_pred, color = season)) +
  geom_line(linewidth = 1) +
  # 按研究区域和距路距离分面
  facet_grid(study_area ~ dist_road, 
             labeller = labeller(dist_road = function(x) paste0("距路距离: ", x, "米"))) +
  # 设置x轴刻度
  scale_x_continuous(breaks = sun_ticks, labels = sun_labels) +
  # 自定义主题与标签
  labs(
    x = "太阳时间（弧度）",
    y = "预测移动速率（m/h）",
    color = "季节",
    title = "太阳时间对动物移动速率的预测曲线"
  ) +
  theme_bw() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title = element_text(size = 12),
    legend.title = element_text(size = 12),
    strip.text = element_text(size = 11)
  )

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