第一季度

你在这里做错了两件事。第一个通常是坏事；一般不要深入研究模型对象并删除组件。在这种情况下，学习使用提取器功能resid()。在这种情况下，您将获得一些有用的东西，但如果您有不同类型的模型对象，例如来自 的 GLM glm()，那么mod$residuals将包含上一次 IRLS 迭代的工作残差，这是您通常不想要的东西！

您做错的第二件事也让我感到困惑。您提取的残差（如果您使用过，也会提取resid()）是原始残差或响应残差。本质上，这是响应的拟合值和观测值之间的差异，仅考虑固定效应项。这些值将包含与 相同的残差自相关，m1因为两个模型中的固定效应（或者如果您愿意，线性预测变量）是相同的 ( ~ time + x)。

要获得包含您指定的相关项的残差，您需要标准化残差。您可以通过以下方式获得这些：

resid(m1, type = "normalized")

这（和其他类型的可用残差）在以下内容中进行了描述?residuals.gls：

type: an optional character string specifying the type of residuals
      to be used. If ‘"response"’, the "raw" residuals (observed -
      fitted) are used; else, if ‘"pearson"’, the standardized
      residuals (raw residuals divided by the corresponding
      standard errors) are used; else, if ‘"normalized"’, the
      normalized residuals (standardized residuals pre-multiplied
      by the inverse square-root factor of the estimated error
      correlation matrix) are used. Partial matching of arguments
      is used, so only the first character needs to be provided.
      Defaults to ‘"response"’.

通过比较，这里是原始（响应）和归一化残差的 ACF

layout(matrix(1:2))
acf(resid(m2))
acf(resid(m2, type = "normalized"))
layout(1)

要了解为什么会发生这种情况，以及原始残差不包括相关项的地方，请考虑您拟合的模型

$$y = \beta_0 + \beta_1 \mathrm{time} + \beta_2 \mathrm{x} + \varepsilon$$

在哪里

$$\varepsilon \sim N(0, \sigma^2 \Lambda)$$

和是由带有参数的 AR(1) 定义的相关矩阵，其中矩阵的非对角元素用值填充，其中是正整数残差对的时间单位分离。$\Lambda$$\hat{\rho}$$\rho^{|d|}$$d$

原始残差，返回的默认resid(m2)值仅来自线性预测器部分，因此来自该位

$$\beta_0 + \beta_1 \mathrm{time} + \beta_2 \mathrm{x}$$

中包含的相关项的任何信息。$\Lambda$

第二季度

似乎您正试图用线性函数拟合非线性趋势，time并解释与 AR(1)（或其他结构）的“趋势”不匹配。如果您的数据与您在此处提供的示例数据类似，我会拟合 GAM 以实现协变量的平滑函数。这个模型将是

$$y = \beta_0 + f_1(\mathrm{time}) + f_2(\mathrm{x}) + \varepsilon$$

最初我们将假设与 GLS 相同的分布，除了最初我们假设 （单位矩阵，因此是独立的残差）。该模型可以使用$\Lambda = \mathbf{I}$

library("mgcv")
m3 <- gam(y ~ s(time) + s(x), select = TRUE, method = "REML")

whereselect = TRUE应用一些额外的收缩以允许模型从模型中删除任何一项。

该模型给出

> summary(m3)

Family: gaussian 
Link function: identity 

Formula:
y ~ s(time) + s(x)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  23.1532     0.7104   32.59   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
          edf Ref.df      F  p-value    
s(time) 8.041      9 26.364  < 2e-16 ***
s(x)    1.922      9  9.749 1.09e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

并具有如下所示的平滑术语：

该模型的残差也表现得更好（原始残差）

acf(resid(m3))

现在提醒一句；平滑时间序列存在一个问题，因为决定函数平滑或摆动程度的方法假设数据是独立的。这实际上意味着时间 ( s(time)) 的平滑函数可以拟合真正随机自相关误差的信息，而不仅仅是潜在趋势。因此，在将平滑器拟合到时间序列数据时应该非常小心。

有很多方法可以解决这个问题，但一种方法是切换到拟合模型，通过gamm()它在内部调用lme()并允许您correlation使用用于gls()模型的参数。这是一个例子

mm1 <- gamm(y ~ s(time, k = 6, fx = TRUE) + s(x), select = TRUE,
            method = "REML")
mm2 <- gamm(y ~ s(time, k = 6, fx = TRUE) + s(x), select = TRUE,
            method = "REML", correlation = corAR1(form = ~ time))

请注意，我必须修复自由度，s(time)因为这些数据存在可识别性问题。该模型可能是摇摆不定的s(time)，没有 AR(1) ($\rho = 0$) 或线性s(time)(1 个自由度) 和强 AR(1) ($\rho >> .5$）。因此，我对潜在趋势的复杂性做出了有根据的猜测。（注意我没有在这个虚拟数据上花费太多时间，但是您应该查看数据并使用您对时间可变性的现有知识来确定样条曲线的适当自由度数。）

与没有 AR(1) 的模型相比，具有 AR(1) 的模型没有显着改进：

> anova(mm1$lme, mm2$lme)
        Model df      AIC      BIC    logLik   Test   L.Ratio p-value
mm1$lme     1  9 301.5986 317.4494 -141.7993                         
mm2$lme     2 10 303.4168 321.0288 -141.7084 1 vs 2 0.1817652  0.6699

如果我们查看 $\hat{\rho}} 的估计值，我们会看到

> intervals(mm2$lme)
....

 Correlation structure:
         lower      est.     upper
Phi -0.2696671 0.0756494 0.4037265
attr(,"label")
[1] "Correlation structure:"

Phi我叫什么$\rho$. 因此，0 是可能的值$\rho$. 该估计值略大于零，因此对模型拟合的影响可以忽略不计，因此如果有很强的先验理由假设残差自相关，您可能希望将其留在模型中。