2

MuMIN の model.avg 関数を、直接入力するのではなく貼り付けてインデックスを使用するモデル式で使用しようとしています。たとえば、次のようになります。

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

'combns' は、combn() によって作成された 2D 配列で、予測子変数の組み合わせが含まれています。これにより、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 に数式内の予測変数を適切に読み取らせるにはどうすればよいですか? ここでは例として 6 つのモデルのみを含めましたが、128 のモデルの完全なセットを比較したいので (さらに多数の予測変数を持つ他の応答変数があります)、それらを個別に入力することは現実的ではありません。

前もって感謝します。

編集:再現可能な例

問題を絞り込むのに時間がかかりました。最初の例 m.ave は、for ループを使用した問題の動作を示しています。2 番目の例 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