即使在具有適當相關結構的模型中，自相關殘差模式是否仍然存在，以及如何選擇最佳模型？
語境
這個問題使用 R，但是關於一般統計問題。
我正在分析死亡率因素（疾病和寄生引起的死亡率百分比）對飛蛾種群增長率的影響，其中每年一次從 12 個地點對幼蟲種群進行抽樣，為期 8 年。人口增長率數據隨著時間的推移呈現出明顯但不規則的周期性趨勢。
來自簡單廣義線性模型的殘差（增長率 ~ 疾病百分比 + 寄生百分比 + 年）隨著時間的推移顯示出同樣清晰但不規則的周期性趨勢。因此，同樣形式的廣義最小二乘模型也適用於具有適當相關結構的數據，以處理時間自相關，例如復合對稱、自回歸過程階數1和自回歸移動平均相關結構。
模型都包含相同的固定效應，使用 AIC 進行比較，並通過 REML 擬合（允許通過 AIC 比較不同的相關結構）。我正在使用 R 包 nlme 和 gls 函數。
問題 1
GLS 模型的殘差在隨時間繪製時仍顯示幾乎相同的周期性模式。即使在準確解釋自相關結構的模型中，這種模式是否會一直存在？
我在第二個問題下方的 R 中模擬了一些簡化但相似的數據，這基於我目前對評估模型殘差中時間自相關模式所需的方法的理解顯示了這個問題，我現在知道這是錯誤的（見答案）。
問題2
我已經為我的數據擬合了具有所有可能合理相關結構的 GLS 模型，但實際上沒有一個模型比沒有任何相關結構的 GLM 更適合：只有一個 GLS 模型稍微好一點（AIC 得分 = 低 1.8），而所有其他模型都有更高的 AIC 值。但是，只有當所有模型都由 REML 擬合時才會出現這種情況，而不是 ML，其中 GLS 模型顯然要好得多，但我從統計書籍中了解到，您只能使用 REML 來比較具有不同相關結構和相同固定效應的模型，原因是我不會在這裡詳細說明。
鑑於數據明顯的時間自相關性質，如果沒有模型比簡單的 GLM 稍微好一點，假設我使用的是適當的方法（我最終想使用AIC 比較不同的變量組合）？
Q1“模擬”探索具有和不具有適當相關結構的模型中的殘差模式
生成具有“時間”循環效應和“x”正線性效應的模擬響應變量：
time <- 1:50 x <- sample(rep(1:25,each=2),50) y <- rnorm(50,5,5) + (5 + 15*sin(2*pi*time/25)) + (x/1)y 應顯示隨“時間”隨機變化的周期性趨勢：
plot(time,y)與隨機變化的“x”呈正線性關係：
plot(x,y)創建一個“y ~ time + x”的簡單線性加法模型：
require(nlme) m1 <- gls(y ~ time + x, method="REML")當針對“時間”繪製時，該模型在殘差中顯示出清晰的周期性模式，正如預期的那樣：
plot(time, m1$residuals)當針對“x”進行繪製時，殘差中應該沒有任何模式或趨勢，這應該是一個很好的、明顯的缺失：
plot(x, m1$residuals)當使用 AIC 評估時，由於自相關結構，包含 1 階自回歸相關結構的簡單模型“y ~ time + x”應該比之前的模型更好地擬合數據：
m2 <- gls(y ~ time + x, correlation = corAR1(form=~time), method="REML") AIC(m1,m2)但是，模型仍應顯示幾乎相同的“時間”自相關殘差：
plot(time, m2$residuals)非常感謝您的任何建議。
第一季度
你在這裡做錯了兩件事。第一個通常是壞事；一般不要深入研究模型對象並刪除組件。在這種情況下，學習使用提取器功能
resid()。在這種情況下，您會得到一些有用的東西，但如果您有不同類型的模型對象，例如來自 的 GLMglm()，那麼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)
要了解為什麼會發生這種情況，以及原始殘差不包括相關項的情況，請考慮您擬合的模型
在哪裡
和是由帶有參數的 AR(1) 定義的相關矩陣其中矩陣的非對角元素填充有值， 在哪裡是殘差對的時間單位的正整數分離。
原始殘差，返回的默認
resid(m2)值僅來自線性預測器部分，因此來自該位因此它們不包含有關包含在.
第二季度
似乎您正試圖用線性函數擬合非線性趨勢，
time並解釋與 AR(1)（或其他結構）的“趨勢”不匹配。如果您的數據與您在此處提供的示例數據類似，我會擬合 GAM 以實現協變量的平滑函數。這個模型將是並且最初我們將假設與 GLS 相同的分佈，除了最初我們將假設 （單位矩陣，所以獨立殘差）。該模型可以使用
library("mgcv") m3 <- gam(y ~ s(time) + s(x), select = TRUE, method = "REML")where
select = 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) () 或線性s(time)(1 個自由度) 和強 AR(1) (）。因此，我對潛在趨勢的複雜性做出了有根據的猜測。（注意我沒有在這個虛擬數據上花費太多時間，但是您應該查看數據並使用您對時間可變性的現有知識來確定樣條曲線的適當數量的自由度。）與沒有 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我叫什麼. 因此，0 是可能的值. 該估計值略大於零，因此對模型擬合的影響可以忽略不計，因此如果有很強的先驗理由假設殘差自相關，您可能希望將其留在模型中。
