(MuMIn)当全局混合效应模型缺乏秩时

Question

我试图使用glmer()和dredge()对泊松混合效应模型进行变量选择.由于多个变量是共线的，所以我使用dredge的子集函数来避免相关变量。然而，要有效地使用dredge()，就需要有一个包含所有术语的完整模型--这可能导致完全模型缺乏等级。

为了给出一个可重复的例子，让我们生成一个随机数据集：

dfdat<-data.frame(replicate(6, round(rnorm(6),2)))
dfdat$group<-factor(sample(1:2,nrow(dfdat),replace=T))
dfdat$Y<-rpois(nrow(dfdat),10)+rpois(nrow(dfdat),as.numeric(dfdat$group))
dfdat
     X1    X2    X3    X4    X5    X6 group  Y
1 -0.88  0.05  1.33 -1.51  0.61 -0.09     2  8
2 -0.12 -0.57  0.05 -1.12  0.60 -0.41     1  7
3  0.14 -0.97 -1.04  0.40  0.87  0.27     1  9
4 -1.04 -0.26 -1.33  0.77 -1.84  1.67     1 11
5 -1.06  1.10 -0.09  0.50 -2.62  2.15     1 10
6 -1.74 -0.61  0.72 -0.29 -0.30 -0.93     1  8

试图运行一个包含所有6个术语的模型是行不通的，因为该模型缺乏等级：

 #library(MuMIn) # not run
 #library(lme4) # not run
 vars<-names(dfdat)[1:6]
 form<-formula(paste0('Y~',paste0(vars,collapse='+'),'+(1|group)'))
 fmod<-glmer(form,data=dfdat,family='poisson')
 fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
 Error: pwrssUpdate did not converge in (maxit) iterations

在dredge上使用fmod将导致glmer删除的一个变量始终被排除在外。

解决方案似乎是: 1.运行一个收敛的模型，2.技巧通过更改收敛模型中的公式来考虑变量的完整列表。

## full model is rank deficient, so use smaller subset
vars.red<-vars[1:3]
form.red<-formula(paste0('Y~',paste0(vars.red,collapse='+'),'+(1|group)'))
fmod.red<-glmer(form.red,data=dfdat,family='poisson')

这种新模型fmod.red收敛，但只包含变量X1、X2和X3。

现在是“欺骗疏浚”部分。上面链接的页面上提出的解决方案并不适用于glmer，因为mermod的结构不同于gamms。所以我试着用：

fom.red@call$formula<-form

其中form有我所有的协变量(将被子设置)。这是行不通的，但是使用框架元素中的公式(如Kamil Bartoń下面所建议的那样)确实有效：

# replace formula in the frame element of fmod.red
attr(fmod.red@frame,"formula")<-form
# check
formula(fmod.red)
# now apply dredge function with covariates
# exclude variable combinations (randomly chosen for the sake of example)
sexpr<-expression(!((X1 && X3) || (X1&&X6) || (X4 && X6) || (X4 && X5)))
# run dredge()
options(na.action = na.fail)
ms<-dredge(fmod.red,subset=sexpr)

更新

而ms似乎包含所有变量，如下所示：

names(ms)
[1] "(Intercept)" "X1"          "X2"          "X3"          "X4"          "X5"          "X6"         
[8] "df"          "logLik"      "AICc"        "delta"       "weight"  

新变量(X4、X5、X6)实际上从未包括在内(到处都是NAs)：

summary(ms)
  (Intercept)          X1                X2                X3              X4        X5        X6    
 Min.   :2.407   Min.   :0.09698   Min.   :-0.4026   Min.   :-0.42078   +   : 0   +   : 0   +   : 0  
 1st Qu.:2.443   1st Qu.:0.22688   1st Qu.:-0.3204   1st Qu.:-0.35303   NA's:26   NA's:26   NA's:26  
 Median :2.474   Median :0.27361   Median :-0.2980   Median :-0.22444                                
 Mean   :2.535   Mean   :0.27539   Mean   :-0.3059   Mean   :-0.23517                                
 3rd Qu.:2.515   3rd Qu.:0.32357   3rd Qu.:-0.2718   3rd Qu.:-0.17472                                
 Max.   :3.009   Max.   :0.45664   Max.   :-0.2177   Max.   : 0.08802                                
             NA's   :20        NA's   :13        NA's   :16                                      

这是怎么回事？

Answer 1

在"merMod"对象中，首先在attr(<object>@frame, "formula")中查找公式(参见getS3method("formula", "merMod")的函数代码)。因此，在call元素中替换它是无效的，可以用formula()或getAllTerms()进行测试。替换"formula"属性@frame。

编辑：事实证明，欺骗dredge并不那么容易，因为在构建表时，它也会查看coef (在本例中是fixef )。要解决这个问题，首先生成调用evaluate，然后用model.sel构建表

model.sel(lapply(dredge(..., evaluate = FALSE), eval), ...)