用 R 中的重复测量/时间序列数据解释多级线性模型

问题描述 投票：0回答：2

我正在阅读 Andy Field 教科书 使用 R 发现统计，在第 899 页，我的 R 代码没有返回与书中所示相同的结果。

具体来说，该示例涉及将多级线性模型拟合到重复测量数据，并通过一次一个参数创建更简单的模型来慢慢构建最终模型来实现这一点。

我的问题源于将时间变量从一个因子强制转换为数字，这很简单，但我的解决方案涉及将时间整数输入为月份（例如，数据文件使用“Satisfaction_6_Months”和“Satisfaction_12_Months”），因此我对这四个进行了编码时间级别为 c(0,6,12,18)。

文本（实际上是在线勘误表）使用了不同的方法，将四个时间级别编码为 c(0,1,2,3)。我本以为时间代码之间的差异对数据影响不大，但我错了。

下面的代码演示了 anova 中的最后一行（包含 ARModel）在这两种情况下有何不同。

我试图理解 1) 为什么会观察到这种差异，2) 为什么仅在模型构建练习中的最后一点观察到这种差异，以及 3) 当我将这些内容付诸实践时，它通常意味着什么可能有充分的理由采用按月编码结构。

首先，运行完整脚本来复制结果。然后，取消注释指示的行并再次运行。

library(reshape2)

data_url <- "https://studysites.sagepub.com/dsur/study/DSUR%20Data%20Files/Chapter%2019/Honeymoon%20Period.dat"

satisfactionData <- read.delim(data_url, header = TRUE)

restructuredData <- melt(satisfactionData, 
                         id = c("Person", "Gender"), 
                         measured = c("Satisfaction_Base", "Satisfaction_6_Months", "Satisfaction_12_Months", "Satisfaction_18_Months"))

names(restructuredData) <- c("Person", "Gender", "Time", "Life_Satisfaction")

restructuredData$Time<-as.numeric(restructuredData$Time)-1

# ***************************************************************
# On second pass, uncomment the line below and rerun whole script
# ***************************************************************
# restructuredData$Time = restructuredData$Time*6


# create models one parameter at a time
intercept <- gls(Life_Satisfaction~1, 
                 data = restructuredData, 
                 method = "ML", 
                 na.action = na.exclude)

randomIntercept <- lme(Life_Satisfaction ~1, 
                       data = restructuredData, 
                       random = ~1|Person, 
                       method = "ML",  
                       na.action = na.exclude, 
                       control = list(opt="optim"))

timeRI <- update(randomIntercept, .~. + Time)

timeRS <- update(timeRI, random = ~Time|Person)

ARModel <- update(timeRS, correlation = corAR1(0, form = ~Time|Person))

anova(intercept, randomIntercept, timeRI, timeRS, ARModel)

# Pass #1 - with Time as 0,1,2,3

#                 Model df      AIC  BIC         logLik   Test   L.Ratio p-value
# intercept           1  2 2064.053 2072.217 -1030.0263                         
# randomIntercept     2  3 1991.396 2003.642  -992.6978 1 vs 2  74.65704  <.0001
# timeRI              3  4 1871.728 1888.057  -931.8642 2 vs 3 121.66714  <.0001
# timeRS              4  6 1874.626 1899.120  -931.3131 3 vs 4   1.10224  0.5763
# ARModel             5  7 1872.891 1901.466  -929.4453 4 vs 5   3.73564  0.0533


# Pass #2 - with Time as 0,6,12,18

#                 Model df      AIC  BIC         logLik   Test   L.Ratio p-value
# intercept           1  2 2064.053 2072.217 -1030.0263                         
# randomIntercept     2  3 1991.396 2003.642  -992.6978 1 vs 2  74.65704  <.0001
# timeRI              3  4 1871.728 1888.057  -931.8642 2 vs 3 121.66714  <.0001
# timeRS              4  6 1874.627 1899.120  -931.3135 3 vs 4   1.10151  0.5765
# ARModel             5  7 1876.627 1905.203  -931.3135 4 vs 5   0.00001  0.9978
r time-series multilevel-analysis
2个回答
1
投票

这并不是一个真正的编程问题，但我们开始吧。

一阶自回归协方差背后的假设是，每个时间点通过单个参数（

Phi
输出中的
nlme::lme
）与其他时间点相关，并且这种相关性与时间距离成比例回归。具体来说，如果 
t
和 
t+1
之间的相关性是 
Phi
，则 
t
和 
t+2
之间的相关性是 
Phi^2
（较小，因为 Phi 几乎总是 <1). This generalizes to 
cor(t, t+x) = Phi^x
。

如果您查看这两个模型的（初始）相关矩阵，您实际上可以看到这一点。我在这里没有使用默认的

Phi=0
，因为那样你最终会得到一个充满零的矩阵：

restructuredData$Time6 <- restructuredData$Time * 6

## correlation matrix with Time = 0, 1, 2, 3
cor1 <- nlme::corAR1(.8, form = ~Time|Person) |> nlme::Initialize(data=restructuredData)
nlme::corMatrix(cor1)[[1]]

#>      [,1] [,2] [,3] [,4]
#> [1,] 1.00 0.80 0.64 0.51
#> [2,] 0.80 1.00 0.80 0.64
#> [3,] 0.64 0.80 1.00 0.80
#> [4,] 0.51 0.64 0.80 1.00

## correlation matrix with Time = 0, 6, 12, 18
cor6 <- nlme::corAR1(.8, form = ~Time6|Person) |> nlme::Initialize(data=restructuredData)
nlme::corMatrix(cor6)[[1]]

#>      [,1] [,2] [,3] [,4]
#> [1,] 1.00 0.26 0.07 0.02
#> [2,] 0.26 1.00 0.26 0.07
#> [3,] 0.07 0.26 1.00 0.26
#> [4,] 0.02 0.07 0.26 1.00

请注意，对于相同的 

Phi
，将时间尺度增加 6 会如何将其他相同观测值之间的（假设）相关性降低 6 次幂或 
0.80^6 ~ 0.26
。这与您在前面步骤中添加的效果不同，因为您的拟合参数及其方差将简单地增大六倍（本质上是单位的变化，t/F统计数据将保持相同）。

在 0:3 比例中拟合的 

Phi
为 0.215，因此在六倍时间尺度中对应于 
0.215^6 ~ 1e-4
，非常接近于零。与实际为零的默认初始化值相结合，六倍时间的模型fits相关性为零。这就是为什么只有最后一步不同，以及为什么这两个模型在六倍尺度的方差分析中没有不同。

我还想在这里强调，这对

corAR1
中的初始化值非常敏感：使其非零将允许第二个模型拟合非零相关性，尽管它似乎对您的值非常敏感请指定。所有这些可能表明相关性实际上相当低和/或 AR1 不是最合适的协方差结构。

最后，请注意，这也是相关矩阵指定方式的结果；另一种方法是将 

Phi
指定为两个相邻时间点之间的相关性，而不管它们的距离如何（即仅假设它们都是等距的，就像原始模型基本上所做的那样）。您可以通过将 
form = Time|Person
更改为 
form = 1|Person
来完成此操作。拟合例程现在假设您的数据是按组内的时间顺序提供的，并且所有组也可能包含相同数量的观察值或至少还包含中间缺失值。然后，组内的每个记录都会占用相关矩阵中的下一行/列，而不考虑它们可能相距多远（除了通过排序设置的时间顺序）。这样，两个模型将再次相同。

0
投票

按照 PBulls 的建议，我稍微修改了代码，现在得到的结果与书上的结果几乎相同：

library(reshape2)

data_url <- "https://studysites.sagepub.com/dsur/study/DSUR%20Data%20Files/Chapter%2019/Honeymoon%20Period.dat"

satisfactionData <- read.delim(data_url, header = TRUE)

restructuredData <- melt(satisfactionData, 
                         id = c("Person", "Gender"), 
                         measured = c("Satisfaction_Base", "Satisfaction_6_Months", "Satisfaction_12_Months", "Satisfaction_18_Months"))

names(restructuredData) <- c("Person", "Gender", "Time", "Life_Satisfaction")

restructuredData$Time<-as.numeric(restructuredData$Time)-1
# ***************************************************************
# Leave my original c(0,6,12,18) codes in
# ***************************************************************
restructuredData$Time = restructuredData$Time*6


# create models one parameter at a time
intercept <- gls(Life_Satisfaction~1, 
                 data = restructuredData, 
                 method = "ML", 
                 na.action = na.exclude)

randomIntercept <- lme(Life_Satisfaction ~1, 
                       data = restructuredData, 
                       random = ~1|Person, 
                       method = "ML",  
                       na.action = na.exclude, 
                       control = list(opt="optim"))

timeRI <- update(randomIntercept, .~. + Time)

timeRS <- update(timeRI, random = ~Time|Person)

ARModel_corAR1 <- update(timeRS, correlation = corAR1(form = ~Time|Person))
anova(intercept, randomIntercept, timeRI, timeRS, ARModel_corAR1)

#                 Model df      AIC  BIC         logLik   Test   L.Ratio p-value
# intercept           1  2 2064.053 2072.217 -1030.0263                         
# randomIntercept     2  3 1991.396 2003.642  -992.6978 1 vs 2  74.65704  <.0001
# timeRI              3  4 1871.728 1888.057  -931.8642 2 vs 3 121.66714  <.0001
# timeRS              4  6 1874.627 1899.120  -931.3135 3 vs 4   1.10151  0.5765
# ARModel_corAR1      5  7 1876.627 1905.203  -931.3135 4 vs 5   0.00001  0.9978

ARModel_corCAR1 <- update(timeRS, correlation = corCAR1(form = ~Time|Person))
anova(intercept, randomIntercept, timeRI, timeRS, ARModel_corCAR1)

#                 Model df      AIC  BIC         logLik   Test   L.Ratio p-value
# intercept           1  2 2064.053 2072.217 -1030.0263                         
# randomIntercept     2  3 1991.396 2003.642  -992.6978 1 vs 2  74.65704  <.0001
# timeRI              3  4 1871.728 1888.057  -931.8642 2 vs 3 121.66714  <.0001
# timeRS              4  6 1874.627 1899.120  -931.3135 3 vs 4   1.10151  0.5765
# ARModel_corCAR1     5  7 1872.884 1901.459  -929.4419 4 vs 5   3.74319  0.0530

最后一个模型与这本书非常接近， 其中 L.Ratio = 3.736 且 p = .0533。

简单地说，我的时间代码更好地建模为连续 (CAR1)，而不是离散 (AR1)。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.