Итак, я хочу исследовать линейные модели с взаимодействиями, используя функцию "regsubsets" из пакета "прыжки". Я легко могу сделать модель без взаимодействий, но когда я пытаюсь включить взаимодействия, вот где она становится сложной. Похоже, что не вполне возможно запустить «regsubsets» в модели взаимодействия, так как у меня есть 8 переменных-предикторов и 388 наблюдений для каждой, и я видел в других постах, что ищет взаимодействия с таким большим количеством переменных и таким количеством наблюдений займет невероятно много времени, чтобы найти (как не в моей жизни). Я довольно новичок в этой функции, так как обычно я нашел другие достойные способы найти полезную модель, поэтому я не знаю, упускаю ли я какой-то очевидный способ вычисления "regsubsets" с взаимодействиями достаточно быстрым способом. Поэтому, если у кого-то есть какие-либо другие предложения, кроме следующего, который я нашел, пожалуйста, go впереди.
Fortunatley, я нашел в Интернете предложение о том, как обойти мою проблему. Один человек предоставил код, который «проанализировал бы результаты поиска, чтобы исключить модели, которые содержат термины во взаимодействиях, которые не включены в качестве основных эффектов». К сожалению, похоже, что код содержит некоторые ошибки и не работает вообще. Я пытался это исправить, но безрезультатно. Мне было интересно, может ли кто-то еще взглянуть на код и помочь исправить его, или помочь предложить лучший способ достижения того же желаемого результата.
Код с собственными комментариями следующий:
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
Код без комментариев следующий:
reg.output.search.with.test<- function (search_object) {
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)
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)]
search_comp$reject_model<-FALSE
for(main_eff_n in 1:length(main_pred)){
search_cols<-grep(pattern=main_pred[main_eff_n],x=higher_pred)
valid_model_subs<-search_comp[search_comp$reject_model==FALSE,]
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){
for(high_eff_n in search_cols){
test_responses<-((main_pred_df[,main_eff_n]==FALSE)&(higher_pred_df[,high_eff_n]==TRUE))
valid_model_subs[test_responses,"reject_model"]<-TRUE
}
search_comp[row.names(valid_model_subs),"reject_model"]<-valid_model_subs[,"reject_model"
}
}
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=",")
}
Любая помощь с этим будет высоко ценится! Спасибо!