2

我正在尝试将 MuMIn 的 model.avg 函数与粘贴的模型公式一起使用，并使用索引而不是直接输入，例如：

m1<-gls(as.formula(paste(response,"~",paste(combns[,j], collapse="+"))), data=dat)

'combns' 是由 combn() 创建的包含预测变量组合的二维数组。如果 gls 函数直接包含公式，则会产生模型平均系数和 AICc 值，例如：

m1<-gls(median_Ta ~ day_of_season + hour_of_day + pct_grey_cover + 
    foliage_height_diversity + tree_shannon_diversity + median_patch_size, data=dat)

但是，相对变量重要性不是计算，我认为这与使用 for 循环或使用变量访问存储模型的列表的索引有关，导致组件模型项不正确“阅读”（参见模型的术语代码）：

Component models: 
         df  logLik   AICc delta weight
1234567b  7 -233.08 481.43  0.00   0.59
1234567f  3 -237.97 482.21  0.78   0.40
1234567e  4 -241.32 491.08  9.65   0.00
1234567a  9 -241.15 502.39 20.96   0.00
1234567c  6 -248.37 509.68 28.25   0.00
1234567d  5 -250.22 511.11 29.68   0.00

Term codes: 
           day_of_season foliage_height_diversity              hour_of_day 
                       1                        2                        3 
       median_patch_size           pct_grey_cover   tree_shannon_diversity 
                       4                        5                        6 
 urban_boundary_distance 
                       7

这导致相对变量重要性被给出为：

Relative variable importance: 
                     day_of_season foliage_height_diversity hour_of_day
Importance:          1             1                        1          
N containing models: 6             6                        6          
                     median_patch_size pct_grey_cover     tree_shannon_diversity
Importance:          1                 1              1                     
N containing models: 6                 6              6                     
                     urban_boundary_distance
Importance:          1                      
N containing models: 6

而如果我在单独键入公式的相同模型上使用 model.avg，我会得到以下正确的输出：

Component models: 
        df  logLik   AICc delta weight
23456    7 -233.08 481.43  0.00   0.59
1        3 -237.97 482.21  0.78   0.40
57       4 -241.32 491.08  9.65   0.00
1234567  9 -241.15 502.39 20.96   0.00
1467     6 -248.37 509.68 28.25   0.00
147      5 -250.22 511.11 29.68   0.00

Relative variable importance: 
                     pct_grey_cover median_patch_size tree_shannon_diversity
Importance:            0.6           0.59              0.59                 
N containing models:     3              4                 3                 
                      foliage_height_diversity hour_of_day day_of_season
Importance:           0.59                     0.59         0.4        
N containing models:     2                        2           4        
                     urban_boundary_distance
Importance:          <0.01                  
N containing models:     4

如何让 model.avg 正确读取公式中的预测变量？我在这里只包含了六个模型作为示例，但我想比较完整的 128 个模型（我还有其他响应变量和更多的预测变量），因此单独输入它们是不可行的。

提前致谢。

编辑：可重现的例子

我花了一段时间来缩小问题的范围。第一个示例 m.ave 显示了使用 for 循环的问题。第二个示例 m.ave2 显示它使用键入的索引而不是使用变量。显然，这只是预测变量的一小部分。

require(nlme)
require(MuMIn)

dat<-data.frame(median_Ta=rnorm(100), day_of_season=runif(100), hour_of_day=runif(100), pct_grey_cover=rnorm(100),
        foliage_height_diversity=rnorm(100), urban_boundary_distance=runif(100), tree_shannon_diversity=rnorm(100), 
        median_patch_size=rnorm(100))

f1<-"median_Ta ~ day_of_season + hour_of_day + pct_grey_cover + foliage_height_diversity + 
        urban_boundary_distance + tree_shannon_diversity + median_patch_size"

f1<-gsub("\\s", "", f1) # remove whitespace
f1split <- strsplit(f1, split="~") # split predictors and response
response <- f1split[[1]][1]
predictors <- strsplit(f1split[[1]][2], split="+", fixed=TRUE)[[1]]

modelslist<-list()

combns <- combn(predictors, 6)
for (j in 1:7) {
    modelslist[[j]]<-gls(as.formula(paste(response,"~",paste(combns[,j], collapse="+"))), data=dat)
}

m.ave<-model.avg(modelslist[[2]], modelslist[[3]], modelslist[[4]],
        modelslist[[5]], modelslist[[6]], modelslist[[7]], modelslist[[8]])

summary(m.ave)

#compare....

modelslist2<-list()
modelslist2[[1]]<-gls(as.formula(paste(response,"~",paste(combns[,1], collapse="+"))), data=dat)
modelslist2[[2]]<-gls(as.formula(paste(response,"~",paste(combns[,2], collapse="+"))), data=dat)
modelslist2[[3]]<-gls(as.formula(paste(response,"~",paste(combns[,3], collapse="+"))), data=dat)
modelslist2[[4]]<-gls(as.formula(paste(response,"~",paste(combns[,4], collapse="+"))), data=dat)
modelslist2[[5]]<-gls(as.formula(paste(response,"~",paste(combns[,5], collapse="+"))), data=dat)
modelslist2[[6]]<-gls(as.formula(paste(response,"~",paste(combns[,6], collapse="+"))), data=dat)
modelslist2[[7]]<-gls(as.formula(paste(response,"~",paste(combns[,7], collapse="+"))), data=dat)

m.ave2<-model.avg(modelslist2[[1]], modelslist2[[2]], modelslist2[[3]], modelslist2[[4]],
        modelslist2[[5]], modelslist2[[6]], modelslist2[[7]])

summary(m.ave2)
4

1 回答 1

1

这是（在包中）formula方法中的一个错误。由于实际公式未存储在对象中的任何位置，因此它计算函数调用中的参数。如果它们的元素都是相同的，例如：glsnlme"model"modellist

modelslist[[1]]$call$model
modelslist[[7]]$call$model

两者都返回

> formula(paste(response, "~", paste(combns[, j], collapse = "+")))

其中，当evaluated 使用 的当前（最后）值时j，所有formula(modellist[[N]])返回最后一个模型公式。

 all.equal(formula(modelslist[[1]]), formula(modelslist[[7]]))

返回

> TRUE

不用说，所有这些都混淆model.avg了使用公式来构建模型选择表（这是一个后备，因为也gls缺乏terms）。

编辑：可能的解决方法

获得您想要的东西的更简单方法：

model.avg(dredge(..., m.lim = c(6,6)))

或者，如果您想进行预测：

modellist <- lapply(dredge(..., m.lim = c(6,6), evaluate = FALSE), eval)

但是，如果您想使用任意一组模型，请将$call$model每个模型对象中的元素替换gls为适当的公式，例如

combns <- combn(1:7, 6)
modellist <- vector("list", 7)
for (j in 1:7) {
    f <- reformulate(predictors[combns[, j]], response = response)
    fm <- gls(f, data = dat)
    fm$call$model <- f # assign the actual formula
    modellist[[j]] <- fm
}
于 2017-08-05T06:19:37.077 回答