方差成分的测试 – 混合模型

huangapple
go评论59阅读模式
英文:

Tests for Variance Components - mixed model

问题

以下是您要翻译的内容：

看下面的情况：

方差成分的测试 – 混合模型

好的，基于这个情况，我在R中拟合了上面的以下模型（不过，我不确定这些模型是否正确）：

library(nlme)
model1 <- lm(Y ~ Treatm * VarT, data = datarats)
model2 <- lme(Y ~ Treatm * VarT, data = datarats, random = ~ 1|RAT, method = "ML")
model3 <- lme(Y ~ Treatm * VarT, data = datarats, random = list(RAT = pdDiag(~ VarT)), method = "ML")
model4 <- lme(Y ~ Treatm * VarT, data = datarats, random = ~ 1 + VarT | RAT, method = "ML")

好的，现在我的兴趣是进行涉及以下比较的方差分量测试：

模型2 vs. 模型1
模型3 vs. 模型2
模型4 vs. 模型2

请看下面的代码：

library(varTestnlme)
library(EnvStats)
library(lmeInfo)
Comp1 <- varCompTest(model2, model1, pval.comp = "both")
summary(Comp1)
print.htestEnvStats(Comp1)
------
Comp2 <- varCompTest(model3, model2, pval.comp = "both")
summary(Comp2)
print.htestEnvStats(Comp2)
------
Comp3 <- varCompTest(model4, model2, pval.comp = "both")
summary(Comp3)
print.htestEnvStats(Comp3)

然而，出现了一些问题（我怀疑主要是在我的模型4的规范中）。这里的问题是：基于我希望得出的结果，模型2和模型4之间的比较的p值应该是p-value = 0.85（((1 - pchisq(0.1,1))/2 + (1 - pchisq(0.1,2))/2，其中0.1是两个模型的loglik之间的差异）。

请看：

方差成分的测试 – 混合模型

因此，当使用varTestnlme包进行测试时，我得到的是p-value = 1，即模型2和模型4是相等的（两个模型的logLik()之间的差异等于零）。请看结果：

summary(Comp3)
混合效应模型中的方差分量测试
测试：
 与VarT相关联的随机效应的方差是否等于0
对立：
 与VarT相关联的随机效应的方差是否> 0

似然比检验统计量：
  LRT =  -1.759843e-07

限制分布：
 2个自由度的两个卡方-平方分布的混合物
相关（精确）权重：0.5 0.5

测试的p值：
从精确权重：1

我在下面放置了随机效应协方差结构的结果：

> VarCorr(model2)
RAT = pdLogChol(1) 
            Variance StdDev  
(Intercept) 3.289539 1.813709
Residual    1.423777 1.193221
> VarCorr(model4)
RAT = pdLogChol(1 + VarT) 
            Variance     StdDev      Corr  
(Intercept) 3.289539e+00 1.813708744 (Intr)
VarT        1.492819e-08 0.000122181 0     
Residual    1.423777e+00 1.193221261

在这个主题中：https://stackoverflow.com/questions/32407631/mixed-effects-model-with-negative-variances中，提到了“nlme”包中关于模型中负方差的限制，可以看到D矩阵中的斜率变化有一个负的条目。然而，在同一个帖子：https://stackoverflow.com/questions/32407631/mixed-effects-model-with-negative-variances中，也提到解决问题的一种方法是考虑“复合对称方差结构”，它允许组内出现负相关。因此，我重新调整了模型如下：

model2 <- lme(Y ~ Treatm * VarT, data = datarats, random = ~ 1|RAT, 
              correlation = corSymm(form = ~1), 
              method = "ML")

model3 <- lme(Y ~ Treatm * VarT, data = datarats, random = list(RAT = pdDiag(~ VarT)), 
              correlation = corSymm(form = ~1), 
              control = lmeControl(opt = "optim", optCtrl = list(method = "BFGS")),
              method = "ML")

model4 <- lme(Y ~ Treatm * VarT, data = datarats, random = ~ 1 + VarT | RAT,
              correlation = corSymm(form = ~1), 
              control = lmeControl(opt = "optim", optCtrl = list(method = "BFGS")), 
              method = "ML")

请看随机效应协方差结构值：

> VarCorr(model2)
RAT = pdLogChol(1) 
            Variance StdDev  
(Intercept) 3.338813 1.827242
Residual    1.417336 1.190519
> VarCorr(model3)
RAT = pdDiag(VarT) 
            Variance    StdDev    
(Intercept) 3.370298164 1.83583718
VarT        0.001353776 0.03679369
Residual    1.384231329 1.17653361
&gt

<details>
<summary>英文:</summary>

See the following situation:

[![enter image description here](https://i.stack.imgur.com/AUsq1.png)](https://i.stack.imgur.com/AUsq1.png)

Ok, based on this, I have fitted the following models above in ````R```` (however, I am not sure if these models are right):

library(nlme)
model1 <- lm(Y ~ Treatm * VarT, data = datarats)
model2 <- lme(Y ~ Treatm * VarT, data = datarats, random = ~ 1|RAT, method = "ML")
model3 <- lme(Y ~ Treatm * VarT, data = datarats, random = list(RAT = pdDiag(~ VarT)), method = "ML")
model4 <- lme(Y ~ Treatm * VarT, data = datarats, random = ~ 1 + VarT | RAT, method = "ML")

Ok, now my interest is to do the test for variance components involving the following comparisons:

Model 2 vs. Model 1
Model 3 vs. Model 2
Model 4 vs. Model 2

See below the code:

library(varTestnlme)
library(EnvStats)
library(lmeInfo)
Comp1 <- varCompTest(model2, model1, pval.comp = "both")
summary(Comp1)
print.htestEnvStats(Comp1)

Comp2 <- varCompTest(model3, model2, pval.comp = "both")
summary(Comp2)
print.htestEnvStats(Comp2)

Comp3 <- varCompTest(model4, model2, pval.comp = "both")
summary(Comp3)
print.htestEnvStats(Comp3)


However, there is something wrong going on (I suspect it is in the specification of my model 4 - mainly).
Here&#39;s the thing: based on the results I would like to arrive at, the p-value involving the comparison between models 2 and 4 should be ````p-value = 0.85````
````((1 - pchisq(0.1,1))/2 + (1 - pchisq(0.1,2))/2````, where 0.1 is the difference  between the loglik of both models).

See:

[![enter image description here](https://i.stack.imgur.com/B4k02.png)](https://i.stack.imgur.com/B4k02.png)

Therefore, when performing the test with the ````varTestnlme package````, I am getting ````p-value = 1````, i.e. models 2 and 4 are equal (difference between the ````logLik()```` of both models is equal to zero).
See the results:

summary(Comp3)
Variance components testing in mixed effects models
Testing that:
variance of the random effect associated to VarT is equal to 0
against the alternative that:
variance of the random effect associated to VarT > 0

Likelihood ratio test statistic:
LRT = -1.759843e-07

Limiting distribution:
mixture of 2 chi-bar-square distributions with degrees of freedom 1 2
associated (exact) weights: 0.5 0.5

p-value of the test:
from exact weights: 1


I put below the results of the random effects covariance structure:

> VarCorr(model2)
RAT = pdLogChol(1)
Variance StdDev
(Intercept) 3.289539 1.813709
Residual 1.423777 1.193221
> VarCorr(model4)
RAT = pdLogChol(1 + VarT)
Variance StdDev Corr
(Intercept) 3.289539e+00 1.813708744 (Intr)
VarT 1.492819e-08 0.000122181 0
Residual 1.423777e+00 1.193221261


In this topic: https://stackoverflow.com/questions/32407631/mixed-effects-model-with-negative-variances it is said that there is a limitation with the ````nlme```` package referring to negative variances in the model and see in the print above that the D matrix has a negative entry for the slope variation. However, in this same post: https://stackoverflow.com/questions/32407631/mixed-effects-model-with-negative-variances it is said that one of the ways to solve the problem is to consider the *compound symmetry variance structure*, which will allow negative correlations within groups. Therefore, I readjusted the models as follows:

model2 <- lme(Y ~ Treatm * VarT, data = datarats, random = ~ 1|RAT,
correlation = corSymm(form = ~1),
method = "ML")

model3 <- lme(Y ~ Treatm * VarT, data = datarats, random = list(RAT = pdDiag(~ VarT)),
correlation = corSymm(form = ~1),
control = lmeControl(opt = "optim", optCtrl = list(method = "BFGS")),
method = "ML")

model4 <- lme(Y ~ Treatm * VarT, data = datarats, random = ~ 1 + VarT | RAT,
correlation = corSymm(form = ~1),
control = lmeControl(opt = "optim", optCtrl = list(method = "BFGS")),
method = "ML")


See random effects covariance structure values:

> VarCorr(model2)
RAT = pdLogChol(1)
Variance StdDev
(Intercept) 3.338813 1.827242
Residual 1.417336 1.190519
> VarCorr(model3)
RAT = pdDiag(VarT)
Variance StdDev
(Intercept) 3.370298164 1.83583718
VarT 0.001353776 0.03679369
Residual 1.384231329 1.17653361
> VarCorr(model4)
RAT = pdLogChol(1 + VarT)
Variance StdDev Corr
(Intercept) 3.372970029 1.83656474 (Intr)
VarT 0.002481283 0.04981247 0.019
Residual 1.375129804 1.17265929


However, the problem still persists. So I have a feeling that my models are not implemented correctly in ````R````, based on the initial information in this post, since the resulting test values for variance components are not consistent with what they should be. Is there any way around this?

To exemplify, I am putting the same situation, however considering ````Orthodont data```` from ````R````:

data(Orthodont)
library(nlme)
library(varTestnlme)
library(EnvStats)
library(lmeInfo)
model11 <- lm(distance ~ Sex * age, data = Orthodont)
model22 <- lme(distance ~ Sex * age, data = Orthodont, random = ~ 1|Subject,
method = "ML")
model33 <- lme(distance ~ Sex * age, data = Orthodont, random = list(Subject = pdDiag(~ age)),
method = "ML")
model44 <- lme(distance ~ Sex * age, data = Orthodont, random = ~ 1 + age | Subject,
method = "ML")

Comp11 <- varCompTest(model22, model11, pval.comp = "both")
summary(Comp11)
print.htestEnvStats(Comp11)

Comp22 <- varCompTest(model33, model22, pval.comp = "both")
summary(Comp22)
print.htestEnvStats(Comp22)

Comp33 <- varCompTest(model44, model22, pval.comp = "both")
summary(Comp33)
print.htestEnvStats(Comp33)


</details>


# 答案1
**得分**: 3

模型已经正确指定。然而，您的参考结果在模型3和4中包含了随机斜率方差分量的负估计。

根据[2015年的这个回答](https://stackoverflow.com/a/32407816/4550695)，`nlme`或`lme4`无法拟合这样的模型。[最近（2023年6月）的讨论](https://stat.ethz.ch/pipermail/r-sig-mixed-models/2023q2/030425.html)在[R-sig-mixed-models](https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models)邮件列表上也没有找到任何能够实现这一点的免费可用的R软件包。

目前在R中似乎很难复制这些结果。

<details>
<summary>英文:</summary>

The models are specified correctly. However, your reference result contains negative estimates of the random slope variance component in models 3 and 4.

According to [this answer from 2015](https://stackoverflow.com/a/32407816/4550695) such models cannot be fitted with `nlme` or `lme4`. A [recent (June 2023) discussion](https://stat.ethz.ch/pipermail/r-sig-mixed-models/2023q2/030425.html) on the [R-sig-mixed-models](https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models) mailing list did not identify any freely available R packages that are able to do so, either.

It seems you cannot easily replicate these results in R at present.

</details>

huangapple
  • 本文由 发表于 2023年7月20日 18:19:27
  • 转载请务必保留本文链接：https://go.coder-hub.com/76728884.html
  • mixed-models
  • nlme
  • r
  • regression
  • statistics
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定