regsubsets() 或其他函数是否可以在 R 中选择具有高交互项和包含主要效果的详尽模型？
Is exhaustive model selection in R with high interaction terms and inclusion of main effects possible with regsubsets() or other functions?
我想在 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
您可以使用 MuMIn 包中的 dredge() 函数。
另见 Subsetting in dredge (MuMIn) - must include interaction if main effects are present。
在使用 dredge 之后，我发现我的模型有太多的预测变量和交互作用，无法在合理的时间内 运行 挖掘（我计算出 40 多个潜在预测变量可能需要 30 万小时在我的电脑上完成搜索）。但它确实排除了交互与主要效果不匹配的模型，所以我想这对许多人来说可能仍然是一个很好的解决方案。
为了满足我的需要，我已经回到 regsubsets 并编写了一些代码来解析搜索输出，以便排除包含未作为主效应包含的交互项的模型。这段代码似乎运行良好，所以我将在这里分享。警告：它是为人类的权宜之计而不是计算而编写的，因此它可能会被重新编码以更快。如果您有 100,000 个模型要测试，您可能希望使其更时尚。 （我一直在使用约 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
我想在 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
您可以使用 MuMIn 包中的 dredge() 函数。
另见 Subsetting in dredge (MuMIn) - must include interaction if main effects are present。
在使用 dredge 之后，我发现我的模型有太多的预测变量和交互作用，无法在合理的时间内 运行 挖掘（我计算出 40 多个潜在预测变量可能需要 30 万小时在我的电脑上完成搜索）。但它确实排除了交互与主要效果不匹配的模型，所以我想这对许多人来说可能仍然是一个很好的解决方案。
为了满足我的需要，我已经回到 regsubsets 并编写了一些代码来解析搜索输出，以便排除包含未作为主效应包含的交互项的模型。这段代码似乎运行良好，所以我将在这里分享。警告：它是为人类的权宜之计而不是计算而编写的，因此它可能会被重新编码以更快。如果您有 100,000 个模型要测试，您可能希望使其更时尚。 （我一直在使用约 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