Пошаговая регрессия с использованием P-значения для переменных переменных с несомненными значениями P-значениями
https://stackoverflow.com/questions/3701170
-
02-10-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
Я хочу выполнить Пошаговая линейная регрессия с использованием P-значения В качестве критерия выбора, например: на каждом шагу отбрасывает переменные, которые имеют самые высокие, то есть наиболее незначимые значения p, останавливаясь, когда все значения значимы, определены некоторого порога альфа.
Я полностью осознаю, что должен использовать AIC (например, команда шаг или стеганыйВместо этого или другого критерия, но мой босс не понимает статистики и настаивает на том, чтобы использовать P-значения.
При необходимости я мог бы программировать мою собственную рутину, но мне интересно, есть ли уже реализованная версия этого.
Решение
Покажите свой босс следующее:
set.seed(100)
x1 <- runif(100,0,1)
x2 <- as.factor(sample(letters[1:3],100,replace=T))
y <- x1+x1*(x2=="a")+2*(x2=="b")+rnorm(100)
summary(lm(y~x1*x2))
Который дает :
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.1525 0.3066 -0.498 0.61995
x1 1.8693 0.6045 3.092 0.00261 **
x2b 2.5149 0.4334 5.802 8.77e-08 ***
x2c 0.3089 0.4475 0.690 0.49180
x1:x2b -1.1239 0.8022 -1.401 0.16451
x1:x2c -1.0497 0.7873 -1.333 0.18566
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Теперь, основываясь на P-значениях, вы бы исключали, какой? X2 является наиболее значимым и не значительным одновременно.
Редактировать: Чтобы уточнить: этот экземпль не лучший, как указано в комментариях. Процедура в Stata и SPSS является AFAIK также не основана на P-значениях T-тестирования на коэффициентах, а на F-тесте после удаления одного из переменных.
У меня есть функция, которая имеет именно это. Это выбор «P-значение», но не о T-тесте на коэффициентах или на результатах ANOVA. Ну, не стесняйтесь использовать его, если это выглядит полезно для вас.
#####################################
# 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)
}
Другие советы
Почему бы не попробовать использовать step() Функция Указание вашего метода тестирования?
Например, для устранения обратной удаления вы вводите только команду:
step(FullModel, direction = "backward", test = "F")
А для поэтапного выбора просто:
step(FullModel, direction = "both", test = "F")
Это может отображать как значения AIC, а также значения F и P.
Вот пример. Начните с самой сложной модели: это включает в себя взаимодействие между всеми тремя пояснительными переменными.
model1 <-lm (ozone~temp*wind*rad)
summary(model1)
Coefficients:
Estimate Std.Error t value Pr(>t)
(Intercept) 5.683e+02 2.073e+02 2.741 0.00725 **
temp -1.076e+01 4.303e+00 -2.501 0.01401 *
wind -3.237e+01 1.173e+01 -2.760 0.00687 **
rad -3.117e-01 5.585e-01 -0.558 0.57799
temp:wind 2.377e-01 1.367e-01 1.739 0.08519
temp:rad 8.402e-03 7.512e-03 1.119 0.26602
wind:rad 2.054e-02 4.892e-02 0.420 0.47552
temp:wind:rad -4.324e-04 6.595e-04 -0.656 0.51358
Трехстороннее взаимодействие явно не значимо. Вот как вы его удаляете, чтобы начать процесс упрощения модели:
model2 <- update(model1,~. - temp:wind:rad)
summary(model2)
В зависимости от результатов вы можете продолжать упростить вашу модель:
model3 <- update(model2,~. - temp:rad)
summary(model3)
...
В качестве альтернативы вы можете использовать функцию упрощения автоматической модели step, чтобы увидеть, насколько хорошо он делает:
model_step <- step(model1)
Если вы просто пытаетесь получить лучшую прогнозную модель, то, возможно, это не имеет значения слишком большого количества, но для чего-либо еще не беспокойтесь с таким типом выбора модели. Это неверно.
Используйте методы усадки, такие как регрессия Ridge (в lm.ridge() в упаковке масса например) или лассо или упругость (комбинация хребтов и ограничений лассо). Из них только лассо и эластичная сеть сделает некоторую форму выбора модели, то есть заставить коэффициенты некоторых ковариатов до нуля.
Увидеть раздел регуляризации и усадки Машинное обучение Просмотр заданий на Cran.
Упаковка RMS: Стратегии моделирования регрессии имеет fastbw() Это делает именно то, что вам нужно. Существует даже параметр для ликвидации из AIC до P-значения.
Как упомянуто Гавин Симпсон Функция fastbw от rms Пакет может быть использован для выбора переменных, использующих значение P. Bellow является примером, использующим пример, данный Георгом Донтасом. Используйте опцию rule='p' выбрать критерии P-значения.
require(rms)
model1 <- ols(Ozone ~ Temp * Wind * Solar.R, data=airquality)
model2 <- fastbw(fit=model1, rule="p", sls=0.05)
model2
