(MuMIn) 当全局混合效应模型秩亏时疏通

Question

我正在尝试运行使用 glmer() 和 dredge() 在泊松混合效应模型上进行变量选择。由于几个变量是共线的，我使用 dredge 的子集函数来避免相关变量。但是，要有效地使用 dredge()，需要一个包含所有项的完整模型 - 这可能导致完整模型排名不足。

[2016 年 2 月 15 日编辑]为了给出一个可重现的例子，让我们生成一个随机数据集：

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

在 fmod 上使用 dredge 会导致 glmer 删除的一个变量总是被排除在外。

解法， 1. 运行一个收敛的模型，2.trick 通过改变收敛模型中的公式疏通到考虑完整的变量列表。

## 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。

现在进入 "tricking dredge" 部分。上面链接的页面上提出的解决方案不适用于 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）（到处都是 NA）：

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")的函数代码）。因此，在调用元素中替换它是无效的，可以用 formula() 或 getAllTerms() 进行测试。替换 @frame.

的 "formula" 属性

编辑：事实证明要欺骗dredge并不那么容易，因为它还会查看coef（或fixef 在这种情况下）构建 table 时。要解决这个问题，首先生成调用，evaluate，然后使用 model.sel:

构建 table

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

(MuMIn) 当全局混合效应模型秩亏时疏通

(MuMIn) Dredge when global mixed-effects model is rank deficient

r

lme4

更新