如何在R的betareg包中处理beta回归的时间自相关问题？

我有一组不同时间间隔的植物病害严重度时间序列数据，响应变量取值在(0,1)区间（不含端点），因此选用beta回归模型最合适。预测变量为气象变量，研究目标是明确气象变量对病害严重度的影响。

可复现数据代码：

df <- structure(
  list(
    year = structure(
      c(1L, 2L, 3L, 5L, 4L, 6L, 7L,
        8L, 9L, 10L),
      .Label = c(
        "2007",
        "2012",
        "2013",
        "2014",
        "2014.1",
        "2015.1",
        "2015.2",
        "2016",
        "2017",
        "2020"
      ),
      class = "factor"
    ),
    mean_rh = c(
      83.9025107032967,
      86.3309364921875,
      78.7225209283154,
      82.3598611111111,
      77.8125490392157,
      77.5694460507813,
      77.340364207483,
      78.601888359589,
      77.9234626042403,
      79.1268228283582
    ),
    mean_temp = c(
      8.06137667087912,
      3.75335204210526,
      8.39571386344086,
      7.57235900444444,
      10.5797098501961,
      8.52468121914062,
      9.63416591190476,
      12.2749429150685,
      9.65211886219081,
      10.1620900981343
    ),
    mean_ws = c(
      1.84656288406593,
      2.0747284924812,
      1.92455694623656,
      2.12702791944444,
      1.91802733215686,
      1.77915314179687,
      1.71631475340136,
      1.78024162876712,
      1.45554503710247,
      1.60409589440299
    ),
    total_rain = c(
      469.5,
      367.1,
      509.3,
      358.7,
      562.3,
      547.6,
      756.4,
      789.5,
      640.5,
      665.1
    ),
    severity = c(
      0.81667,
      0.01325,
      0.06125,
      0.81667,
      0.0198,
      0.0623,
      0.0035,
      0.00475,
      0.885,
      0.0348
    )
  ),
  class = c("grouped_df", "tbl_df", "tbl", "data.frame"),
  row.names = c(NA,-10L),
  groups = structure(
    list(
      year = structure(
        1:10,
        .Label = c(
          "2007",
          "2012",
          "2013",
          "2014",
          "2014.1",
          "2015.1",
          "2015.2",
          "2016",
          "2017",
          "2020"
        ),
        class = "factor"
      ),
      .rows = structure(
        list(1L,
             2L, 3L, 5L, 4L, 6L, 7L, 8L, 9L, 10L),
        ptype = integer(0),
        class = c("vctrs_list_of",
                  "vctrs_vctr", "list")
      )
    ),
    row.names = c(NA,-10L),
    class = c("tbl_df",
              "tbl", "data.frame"),
    .drop = TRUE
  )
)

拟合的beta回归模型代码：

mod05 <- 
  betareg(severity_ps ~ mean_temp +  mean_ws + total_rain + mean_rh,
          data = dat_preplanting)
summary(mod05)

用acf()检查模型残差自相关时，绘图显示存在显著的时间自相关；但用pacf(mod05$residuals)检查偏自相关时，未检测到自相关。

我有两个问题：

1. 若ACF图显示自相关但PACF图无，是否存在自相关问题？
2. 如何在R的betareg包中处理时间自相关？我查阅了文档但未找到相关方法。

另外，我曾尝试用glmmTMB和gam包拟合模型，将year作为随机效应，但模型无法收敛，说明自由度非常有限，因此急需找到在betareg包中处理时间自相关的方法。

问题1：ACF有自相关但PACF无，是否存在自相关问题？ #

是，这确实存在自相关问题。ACF反映的是残差与所有滞后项的总自相关，而PACF是控制了中间滞后项影响后的偏自相关。当ACF显示显著自相关但PACF无显著值时，通常意味着残差存在移动平均(MA)过程（比如MA(1)：当前残差仅与前一期残差的误差项相关），或者是数据中存在某种平滑的趋势性自相关（比如随机游走的残差）。无论哪种情况，自相关的存在都会导致模型的标准误被低估，进而影响系数的显著性检验结果，必须加以处理。

问题2：betareg包中处理时间自相关的方法 #

betareg包本身没有内置的时间序列自相关处理功能，但可以通过以下几种方式解决：

1. 手动加入滞后项作为预测变量 #

如果数据的时间间隔是固定的，可以尝试将响应变量的滞后项（比如severity_lag1）加入模型，直接控制自相关：

# 先创建滞后项（假设数据已按时间排序）
dat_preplanting <- dat_preplanting %>%
  arrange(year) %>%
  mutate(severity_lag1 = lag(severity_ps, 1))

# 拟合包含滞后项的beta回归
mod_lag <- betareg(severity_ps ~ mean_temp + mean_ws + total_rain + mean_rh + severity_lag1,
                   data = dat_preplanting)
summary(mod_lag)

这种方法的优势是简单直接，不需要额外的包，适合小样本数据。注意要处理滞后项带来的缺失值（可以删除第一行或用合适方法填充）。

2. 使用glmmTMB拟合带自相关结构的beta模型 #

虽然你之前尝试加入year随机效应时模型不收敛，但可以尝试直接在glmmTMB中指定自相关协方差结构（比如AR1），而不是用随机效应：

library(glmmTMB)

# 确保数据按时间排序
dat_preplanting <- dat_preplanting %>% arrange(year)

# 拟合带AR1自相关结构的beta模型
mod_ar1 <- glmmTMB(severity_ps ~ mean_temp + mean_ws + total_rain + mean_rh,
                   data = dat_preplanting,
                   family = beta_family(),
                   autocor = corAR1(form = ~1 | year)) # 这里的form指定时间分组

summary(mod_ar1)

如果模型仍不收敛，可以尝试调整初始值，或者简化模型（比如先去掉不显著的气象变量）。glmmTMB的beta族支持指定多种自相关结构，包括AR、MA等。

3. 使用nlme包结合beta变换 #

另一种思路是对响应变量进行logit变换（因为beta回归的链接函数通常是logit），然后用nlme包拟合带自相关结构的线性模型：

library(nlme)

# logit变换响应变量
dat_preplanting <- dat_preplanting %>%
  arrange(year) %>%
  mutate(severity_logit = log(severity_ps / (1 - severity_ps)))

# 拟合带AR1自相关的线性模型
mod_nlme <- gls(severity_logit ~ mean_temp + mean_ws + total_rain + mean_rh,
                data = dat_preplanting,
                correlation = corAR1(form = ~1))

summary(mod_nlme)

这种方法的缺点是变换后的模型假设误差服从正态分布，不如beta回归贴合原始数据的分布，但在小样本下可能更容易收敛，且能有效处理自相关。

4. 尝试betareg的稳健标准误 #

如果上述方法都无法实现，可以退而求其次，使用稳健标准误来修正自相关带来的标准误偏误，不需要改变模型结构：

library(sandwich)
library(lmtest)

# 计算稳健标准误（适用于时间序列的Newey-West标准误）
coeftest(mod05, vcov = NeweyWest(mod05, lag = 1))

Newey-West标准误可以同时处理自相关和异方差，虽然没有从模型层面消除自相关，但能得到更可靠的系数显著性检验结果。

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