线性混合模型置信区间问题

Question

希望你能澄清我脑子里的困惑。

用lmerTest构造线性混合模型

MODEL <- lmer(Ca content ~ SYSTEM +(1 | YEAR/replicate) + 
               (1 | YEAR:SYSTEM), data = IOSDV1)

当我试图获得主要效果的特定水平的置信区间时，乐趣就开始发生了。

命令emmeans和lsmeans产生相同的间隔(例如；SYSTEM A3: 23.9-128.9, mean 76.4, SE:8.96)。

但是，命令as.data.frame(effect("SYSTEM", MODEL))产生不同的、更窄的置信区间(例如；SYSTEM A3: 58.0-94.9, mean 76.4, SE:8.96)。

我少了什么，我应该报告什么号码？

总之，对于钙的含量，我有6次总测量(每年3次，每个来自不同的复制)。我将把代码中的名称保留在我使用的语言中。其目的是检验某些生产实践是否会影响谷物中特定矿物的含量。模型中保留了不带残差方差的随机效应。

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: CA ~ SISTEM + (1 | LETO/ponovitev) + (1 | LETO:SISTEM)
   Data: IOSDV1

REML criterion at convergence: 202.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.60767 -0.74339  0.04665  0.73152  1.50519 

Random effects:
 Groups         Name        Variance Std.Dev.
 LETO:SISTEM    (Intercept)   0.0     0.0    
 ponovitev:LETO (Intercept)   0.0     0.0    
 LETO           (Intercept) 120.9    11.0    
 Residual                   118.7    10.9    
Number of obs: 30, groups:  LETO:SISTEM, 10; ponovitev:LETO, 8; LETO, 2

Fixed effects:
               Estimate Std. Error      df t value Pr(>|t|)  
(Intercept)      76.417      8.959   1.548   8.530   0.0276 *
SISTEM[T.C0]     -5.183      6.291  24.000  -0.824   0.4181  
SISTEM[T.C110]  -13.433      6.291  24.000  -2.135   0.0431 *
SISTEM[T.C165]   -7.617      6.291  24.000  -1.211   0.2378  
SISTEM[T.C55]   -10.883      6.291  24.000  -1.730   0.0965 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
             (Intr) SISTEM[T.C0 SISTEM[T.C11 SISTEM[T.C16
SISTEM[T.C0  -0.351                                      
SISTEM[T.C11 -0.351  0.500                               
SISTEM[T.C16 -0.351  0.500       0.500                   
SISTEM[T.C5  -0.351  0.500       0.500        0.500      
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see ?isSingular

> ls_means(MODEL, ddf="Kenward-Roger")
Least Squares Means table:

           Estimate Std. Error  df t value    lower    upper Pr(>|t|)  
SISTEMA3    76.4167     8.9586 1.5  8.5299  23.9091 128.9243  0.02853 *
SISTEMC0    71.2333     8.9586 1.5  7.9514  18.7257 123.7409  0.03171 *
SISTEMC110  62.9833     8.9586 1.5  7.0305  10.4757 115.4909  0.03813 *
SISTEMC165  68.8000     8.9586 1.5  7.6797  16.2924 121.3076  0.03341 *
SISTEMC55   65.5333     8.9586 1.5  7.3151  13.0257 118.0409  0.03594 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Confidence level: 95%
  Degrees of freedom method: Kenward-Roger

> emmeans(MODEL, spec = c("SISTEM"))
 SISTEM emmean   SE   df lower.CL upper.CL
 A3       76.4 8.96 1.53     23.9      129
 C0       71.2 8.96 1.53     18.7      124
 C110     63.0 8.96 1.53     10.5      115
 C165     68.8 8.96 1.53     16.3      121
 C55      65.5 8.96 1.53     13.0      118

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95

> as.data.frame(effect("SISTEM", MODEL))
  SISTEM      fit       se    lower    upper
1     A3 76.41667 8.958643 57.96600 94.86734
2     C0 71.23333 8.958643 52.78266 89.68400
3   C110 62.98333 8.958643 44.53266 81.43400
4   C165 68.80000 8.958643 50.34933 87.25067
5    C55 65.53333 8.958643 47.08266 83.98400

非常感谢。

Answer 1

我非常肯定，这与可怕的“分母自由度”问题有关，也就是采用了什么样的有限样本修正(如果有的话)。tl； emmeans博士正在使用肯沃德-罗杰修正，这或多或少是最准确的可用选项--唯一不使用for的原因是，如果你有一个大的数据集，它会变得非常慢。

装载包，模拟数据，拟合模型

library(lmerTest)
library(emmeans)
library(effects)
dd <- expand.grid(f=factor(letters[1:3]),g=factor(1:20),rep=1:10)
set.seed(101)
dd$y <- simulate(~f+(1|g), newdata=dd, newparams=list(beta=rep(1,3),theta=1,sigma=1))[[1]]
m <- lmer(y~f+(1|g), data=dd)

比较默认的emmeans方法和效果

emmeans(m, ~f)
##  f emmean    SE   df lower.CL upper.CL
##  a  0.848 0.212 21.9    0.409     1.29
##  b  1.853 0.212 21.9    1.414     2.29
##  c  1.863 0.212 21.9    1.424     2.30

## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 

as.data.frame(effect("f",m))
##   f       fit        se     lower    upper
## 1 a 0.8480161 0.2117093 0.4322306 1.263802
## 2 b 1.8531805 0.2117093 1.4373950 2.268966
## 3 c 1.8632228 0.2117093 1.4474373 2.279008

effects并没有明确告诉我们它使用的是什么/是否使用有限样本更正:我们可以在文档或代码中进行挖掘以试图找到答案。或者，我们可以告诉emmeans不要使用有限样本校正：

emmeans(m, ~f, lmer.df="asymptotic")
##  f emmean    SE  df asymp.LCL asymp.UCL
##  a  0.848 0.212 Inf     0.433      1.26
##  b  1.853 0.212 Inf     1.438      2.27
##  c  1.863 0.212 Inf     1.448      2.28

## Degrees-of-freedom method: asymptotic 
## Confidence level used: 0.95 

测试表明，它们相当于大约0.001的公差(可能足够接近)。原则上，我们应该能够指定KR=TRUE来让effects使用Kenward-Roger校正，但是我还没能让它起作用。

不过，我也要说，你的例子有点古怪。如果我们以标准误差单位计算平均与较低CI之间的距离，对于emmeans，我们得到了(76.4-23.9)/8.96 = 5.86，这意味着一个很小的影响自由度(例如约1.55)。这对我来说很可疑除非你的数据集非常小..。

从你更新的帖子来看，肯沃德-罗杰确实只估计了1.5分母df。

一般来说，如果分组变量有少量的水平(尽管反参数参见这里 )，那么尝试拟合随机效应是不可取的。我试着把LETO (只有两个层次)当作一个固定的效果，即

CA ~ SISTEM + LETO + (1 | LETO:ponovitev) + (1 | LETO:SISTEM)

看看这是否有帮助。(我希望你得到的是7df，这将使你的CIs±2.4se而不是±6se.)