#####################################
# Automated model selection
# Author : Joris Meys
# version : 0.2
# date : 12/01/09
#####################################
#CHANGE LOG
# 0.2 : check for empty scopevar vector
#####################################
# Function has.interaction checks whether x is part of a term in terms
# terms is a vector with names of terms from a model
has.interaction <- function(x,terms){
out <- sapply(terms,function(i){
sum(1-(strsplit(x,":")[[1]] %in% strsplit(i,":")[[1]]))==0
})
return(sum(out)>0)
}
# Function Model.select
# model is the lm object of the full model
# keep is a list of model terms to keep in the model at all times
# sig gives the significance for removal of a variable. Can be 0.1 too (see SPSS)
# verbose=T gives the F-tests, dropped var and resulting model after
model.select <- function(model,keep,sig=0.05,verbose=F){
counter=1
# check input
if(!is(model,"lm")) stop(paste(deparse(substitute(model)),"is not an lm object\n"))
# calculate scope for drop1 function
terms <- attr(model$terms,"term.labels")
if(missing(keep)){ # set scopevars to all terms
scopevars <- terms
} else{ # select the scopevars if keep is used
index <- match(keep,terms)
# check if all is specified correctly
if(sum(is.na(index))>0){
novar <- keep[is.na(index)]
warning(paste(
c(novar,"cannot be found in the model",
"\nThese terms are ignored in the model selection."),
collapse=" "))
index <- as.vector(na.omit(index))
}
scopevars <- terms[-index]
}
# Backward model selection :
while(T){
# extract the test statistics from drop.
test <- drop1(model, scope=scopevars,test="F")
if(verbose){
cat("-------------STEP ",counter,"-------------\n",
"The drop statistics : \n")
print(test)
}
pval <- test[,dim(test)[2]]
names(pval) <- rownames(test)
pval <- sort(pval,decreasing=T)
if(sum(is.na(pval))>0) stop(paste("Model",
deparse(substitute(model)),"is invalid. Check if all coefficients are estimated."))
# check if all significant
if(pval[1]<sig) break # stops the loop if all remaining vars are sign.
# select var to drop
i=1
while(T){
dropvar <- names(pval)[i]
check.terms <- terms[-match(dropvar,terms)]
x <- has.interaction(dropvar,check.terms)
if(x){i=i+1;next} else {break}
} # end while(T) drop var
if(pval[i]<sig) break # stops the loop if var to remove is significant
if(verbose){
cat("\n--------\nTerm dropped in step",counter,":",dropvar,"\n--------\n\n")
}
#update terms, scopevars and model
scopevars <- scopevars[-match(dropvar,scopevars)]
terms <- terms[-match(dropvar,terms)]
formul <- as.formula(paste(".~.-",dropvar))
model <- update(model,formul)
if(length(scopevars)==0) {
warning("All variables are thrown out of the model.\n",
"No model could be specified.")
return()
}
counter=counter+1
} # end while(T) main loop
return(model)
}
7条答案
按热度按时间voase2hg1#
向老板展示以下内容:
其中:
现在,基于p值,您将排除哪一个?x2同时最显著和最不显著。
编辑:澄清:这个例子不是最好的,正如评论中所指出的。Stata和SPSS中的程序是AFAIK也不是基于系数的T检验的P值,而是基于删除其中一个变量后的F检验。
我有一个函数可以做到这一点。这是一个关于"p值"的选择,但不是关于系数或方差分析结果的T检验。好吧,如果它对你有用,你可以随意使用它。
js81xvg62#
为什么不尝试使用
step()
函数指定测试方法呢?例如,对于向后消除,只需键入一个命令:
对于逐步选择,只需:
这可以显示AIC值以及F和P值。
ui7jx7zq3#
下面是一个例子,从最复杂的模型开始:这包括所有三个解释变量之间的相互作用。
三向交互作用显然并不重要,这就是如何将其移除,以开始模型简化过程:
根据结果,您可以继续简化模型:
或者,您可以使用自动模型简化函数
step
,看看它的效果如何:vfwfrxfs4#
软件包rms: Regression Modeling Strategies中的
fastbw()
完全满足您的需求,甚至还有一个参数可以从AIC转换为基于p值的消除。7qhs6swi5#
如果你只是想得到最好的预测模型,那么也许这并不太重要,但对于其他任何事情,不要为这种模型选择而烦恼,这是错误的。
使用收缩方法,如岭回归(例如MASS软件包中的
lm.ridge()
)、套索或弹性网络(岭约束和套索约束的组合)。在这些方法中,只有套索和弹性网络会执行某种形式的模型选择,即强制某些协变量的系数为零。请参见CRAN上Machine Learning任务视图的正则化和收缩部分。
mec1mxoz6#
如Gavin Simpson所述,
rms
程序包中的函数fastbw
可用于选择使用p值的变量。Bellow是使用乔治Dontas给出的示例的示例。使用选项rule='p'
选择p值标准。xu3bshqb7#
OLSRR封装可能是有用。
您可以定义pent(进入模型的p值)和prem(移除的p值)
输出提供了您需要的所有指标,甚至更多。