首页 文章

R中的详尽模型选择是否具有高交互作用项并且可以使用regsubsets（）或其他函数包含主效应？

提问于
浏览 次
1

我想在R中使用7个预测变量（5个连续预测变量和2个分类变量）对数据集执行自动，详尽的模型选择 . 我希望所有连续预测变量都具有交互的可能性（至少最多3个方向的交互）并且还具有不相互作用的平方项 .

我一直在使用 leaps 包中的 regsubsets() 并获得了良好的结果，但是许多模型都包含交互项而不包括主效应（例如， g*h 是包含的模型预测器，但 g 不是） . 由于包含主效应也会影响模型得分（Cp，BIC等），因此将它们包括在与其他模型的比较中是很重要的，即使它们不是强预测因子 .

我可以手动清除结果并交叉包含没有主效应的交互的模型，但是我很可能确定 regsubsets() 或 leaps() 是不可能的，也可能不是 glmulti . 有没有人知道另一个详尽的模型选择功能，允许这样的规范或对脚本建议，将对模型输出进行排序，只找到符合我的规格的模型？

以下是使用 regsubsets() 进行模型搜索的简化输出 . 您可以看到模型3和4确实包含交互术语，而不包括所有相关的主效应 . 如果没有其他功能可以用我的规范运行搜索，那么有关轻松设置此输出以排除没有包含必要主效应的模型的建议将会有所帮助 .

Model adjR2      BIC            CP          n_pred  X.Intercept.    x1      x2      x3      x1.x2   x1.x3   x2.x3   x1.x2.x3
1   0.470344346 -41.26794246    94.82406866 1       TRUE            FALSE   TRUE    FALSE   FALSE   FALSE   FALSE   FALSE
2   0.437034361 -36.5715963     105.3785057 1       TRUE            FALSE   FALSE   TRUE    FALSE   FALSE   FALSE   FALSE
3   0.366989617 -27.54194252    127.5725366 1       TRUE            FALSE   FALSE   FALSE   TRUE    FALSE   FALSE   FALSE
4   0.625478214 -64.64414719    46.08686422 2       TRUE            TRUE    FALSE   FALSE   FALSE   FALSE   FALSE   TRUE
r linear-regression interaction

2 回答

  • 2

    您可以使用MuMIn包中的 dredge() 函数 .

    另见Subsetting in dredge (MuMIn) - must include interaction if main effects are present .

    回复于
  • 2

    使用 dredge 后，我发现我的模型有太多的预测变量和相互作用，可以在合理的时间内运行挖泥机（我计算了40个潜在的预测变量，可能需要300k小时才能在我的计算机上完成搜索） . 但它确实排除了交互与主效应不匹配的模型，所以我想这对许多人来说可能仍然是一个很好的解决方案 .

    根据我的需要，我已经回到 regsubsets 并编写了一些代码来解析搜索输出，以便排除包含未包含在主要效果中的交互中的术语的模型 . 这段代码似乎运行良好，所以我已经有十万个模型来测试你可能想让它更时尚 . （我一直致力于搜索大约50,000个型号和多达40个因素，这些因素需要我的2.4ghz i5核心才能处理几个小时）

    reg.output.search.with.test<- function (search_object) {  ## input an object from a regsubsets search
## First build a df listing model components and metrics of interest
  search_comp<-data.frame(R2=summary(search_object)$rsq,  
                          adjR2=summary(search_object)$adjr2,
                          BIC=summary(search_object)$bic,
                          CP=summary(search_object)$cp,
                          n_predictors=row.names(summary(search_object)$which),
                          summary(search_object)$which)
  ## Categorize different types of predictors based on whether '.' is present
  predictors<-colnames(search_comp)[(match("X.Intercept.",names(search_comp))+1):dim(search_comp)[2]]
  main_pred<-predictors[grep(pattern = ".", x = predictors, invert=T, fixed=T)]
  higher_pred<-predictors[grep(pattern = ".", x = predictors, fixed=T)]
  ##  Define a variable that indicates whether model should be reject, set to FALSE for all models initially.
  search_comp$reject_model<-FALSE  

  for(main_eff_n in 1:length(main_pred)){  ## iterate through main effects
    ## Find column numbers of higher level ters containing the main effect
    search_cols<-grep(pattern=main_pred[main_eff_n],x=higher_pred) 
    ## Subset models that are not yet flagged for rejection, only test these
    valid_model_subs<-search_comp[search_comp$reject_model==FALSE,]  
    ## Subset dfs with only main or higher level predictor columns
    main_pred_df<-valid_model_subs[,colnames(valid_model_subs)%in%main_pred]
    higher_pred_df<-valid_model_subs[,colnames(valid_model_subs)%in%higher_pred]

    if(length(search_cols)>0){  ## If there are higher level pred, test each one
      for(high_eff_n in search_cols){  ## iterate through higher level pred. 
        ##  Test if the intxn effect is present without main effect (working with whole column of models)
        test_responses<-((main_pred_df[,main_eff_n]==FALSE)&(higher_pred_df[,high_eff_n]==TRUE)) 
        valid_model_subs[test_responses,"reject_model"]<-TRUE  ## Set reject to TRUE where appropriate
        } ## End high_eff for
      ## Transfer changes in reject to primary df:
      search_comp[row.names(valid_model_subs),"reject_model"]<-valid_model_subs[,"reject_model"
      } ## End if
    }  ## End main_eff for

  ## Output resulting table of all models named for original search object and current time/date in folder "model_search_reg"
  current_time_date<-format(Sys.time(), "%m_%d_%y at %H_%M_%S")
  write.table(search_comp,file=paste("./model_search_reg/",paste(current_time_date,deparse(substitute(search_object)),
             "regSS_model_search.csv",sep="_"),sep=""),row.names=FALSE, col.names=TRUE, sep=",")
}  ## End reg.output.search.with.test fn
    回复于

相关问题