0

我一直在使用glmulti为我感兴趣的变量获取模型平均估计值和相对重要性值。在运行中，glmulti我指定了一个候选模型，其中所有变量和交互都基于先验知识包含在内（参见下面的代码）。

运行glmutli模型后，我使用函数summary()和weightable(). 结果似乎发生了许多奇怪的事情，我不明白。

首先，当我使用lme4 glmer()函数运行我的候选模型时，我获得了 2086 的 AIC 值。在 glmulti 输出中，这个候选模型（具有完全相同的公式）具有较低的 AIC 值（2107），因此它出现在所有潜在模型列表中的 26 个位置中的第 8 位（通过 weigtable() 函数获得）。

似乎导致此问题的原因是 logArea:Habitat 交互从候选模型中删除，尽管level=2已指定。与通过 提供的公式相比，该函数summary(output_new@objects[[8]]) 提供了不同的公式（没有 logArea:Habitat 交互变量）weightable()。这解释了为什么候选模型 AIC 值与通过 获得的不一样lme4，但我不明白为什么公式中缺少交互变量 logArea:Habitat。其他可能的模型也是如此。似乎对于具有 2 个或更多交互的所有模型，一个交互被丢弃。

有人对发生的事情有解释吗？任何帮助将非常感激！

最好的，罗伯特

注意：我创建了我的数据子集 ( https://drive.google.com/open?id=1rc0Gkp7TPdnhW6Bw87FskL5SSNp21qxl ) 并通过删除变量来简化候选模型以减少模型运行时间。（问题依旧）

     newdat <- Data_ommited2[, c("Presabs","logBodymass", "logIsolation", "Matrix", "logArea", "Protection","Migration", "Habitat", "Guild", "Study","Species", "SpeciesStudy")] 

     glmer.glmulti <- function (formula, data, random, ...) {
     glmer(paste(deparse(formula), random), data = data, family=binomial(link="logit"),contrasts=list(Matrix=contr.sum, Habitat=contr.treatment, Protection=contr.treatment, Guild=contr.sum),glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
     }

   output_new <- glmulti(y = Presabs ~  Matrix  + logArea*Protection + logArea*Habitat,
    data = sampledata,
    random = '+(1|Study)+(1|Species)+(1|SpeciesStudy)',
    family = binomial,
    method = 'h',
    level=2,
    marginality=TRUE,
    crit = 'aic',
    fitfunc = glmer.glmulti,
    confsetsize = 26)

    print(output_new)

    summary(output_new)

    weightable(output_new)
4

1 回答 1

0

我发现了一个遇到同样问题的人的帖子（https://stats.stackexchange.com/questions/341356/glmulti-package-in-r-reporting-incorrect-aicc-values），看来问题是由通过这行代码：

glmer.glmulti <- function (formula, data, random, ...) {
  glmer(paste(deparse(formula), random), data = data, family=binomial(link="logit"))
}

通过将这部分代码更改为以下内容，问题得到解决：

glmer.glmulti<-function(formula,data,random,...) {
  newf <- formula
  newf[[3]] <- substitute(f+r,
                          list(f=newf[[3]],
                               r=reformulate(random)[[2]]))
  glmer(newf,data=data,
        family=binomial(link="logit"))
}
于 2019-05-19T08:41:16.923 回答