Я хочу сделать ступенчатую регрессию с использованием AIC в списке линейных моделей. Идея состоит в том, чтобы использовать список линейных моделей и затем применять stepAIC к каждому элементу списка. Не удается.
Привет, ребята. Я пытался отследить проблему.
Я думаю, что нашел проблему. Однако я не понимаю причину.
Попробуйте код, чтобы увидеть разницу между тремя случаями.
require(MASS)
n<-30
x1<-rnorm(n, mean=0, sd=1) #create rv x1
x2<-rnorm(n, mean=1, sd=1)
x3<-rnorm(n, mean=2, sd=1)
epsilon<-rnorm(n,mean=0,sd=1) # random error variable
dat<-as.data.frame(cbind(x1,x2,x3,epsilon)) # combine to a data frame
dat$id<-c(rep(1,10),rep(2,10),rep(3,10))
# y is combination from all three x and a random uniform variable
dat$y<-x1+x2+x3+epsilon
# apply lm() only resulting in a list of models
dat.lin.model.lst<-lapply(split(dat,dat$id),function(d) lm(y~x1+x2+x3,data=d))
stepAIC(dat.lin.model.lst[[1]]) # FAIL!!!
# apply function stepAIC(lm())- works
dat.lin.model.stepAIC.lst<-lapply(split(dat,dat$id),function(d) stepAIC(lm(y~x1+x2+x3,data=d)))
# create model for particular group with id==1
k<-which(dat$id==1) # manually select records with id==1
lin.model.id1<-lm(dat$y[k]~dat$x1[k]+dat$x2[k]+dat$x3[k])
stepAIC(lin.model.id1) # check stepAIC - works!
Я почти уверен, что stepAIC () нужны исходные данные из data.frame "dat". Это то, о чем я думал раньше. (Надеюсь, я прав в этом)
Но в stepAIC () нет параметра, по которому я могу передать исходный фрейм данных. Очевидно, что для простых моделей, не включенных в список, достаточно пройти модель. (последние три строки в коде) Так что мне интересно:
Q1: Как stepAIC знает, где найти исходные данные "dat" (не только данные модели, которые передаются в качестве параметра)?
В2: Как я могу узнать, что в stepAIC () есть другой параметр, который явно не указан на страницах справки? (может быть, мой английский слишком плох для поиска)
Q3: Как я могу передать этот параметр в stepAIC ()?
Должно быть где-то в среде функции apply и передачи данных. Либо lm () или stepAIC () и указатель / ссылка на необработанные данные должны где-то потеряться. Я не очень хорошо понимаю, что окружение в R делает. Для меня это было своего рода изоляция локальных от глобальных переменных. Но, может быть, это сложнее. Кто-нибудь, кто может объяснить это мне в связи с проблемой выше? Честно говоря, я не слишком много читаю из R документации . Любое лучшее понимание поможет мне. Спасибо.
OLD:
У меня есть данные в фрейме данных df, которые можно разбить на несколько подгрупп. Для этого я создал groupID с именем df $ id. lm () возвращает коэффициент, как и ожидалось для первой подгруппы. Я хочу сделать ступенчатую регрессию, используя AIC в качестве критерия для каждой подгруппы в отдельности. Я использую lmList {lme4}, который приводит к модели для каждой подгруппы (id). Но если я использую stepAIC {MASS} для элементов списка, он выдаст ошибку. увидеть ниже.
Итак, вопрос: какая ошибка в моей процедуре / синтаксисе? Я получаю результаты для отдельных моделей, но не те, которые созданы с помощью lmList. Хранит ли lmList () информацию о модели, отличную от lm ()?
Но в справке говорится:
class "lmList": список объектов класса lm с общей моделью.
>lme4.list.lm<-lmList(formula=Scherkraft.N~Gap.um+Standoff.um+Voidflaeche.px |df$id,data = df)
>lme4.list.lm[[1]]
Call: lm(formula = formula, data = data)
Coefficients:
(Intercept) Gap.um Standoff.um Voidflaeche.px
62.306133 -0.009878 0.026317 -0.015048
>stepAIC(lme4.list.lm[[1]], direction="backward")
#stepAIC on first element on the list of linear models
Start: AIC=295.12
Scherkraft.N ~ Gap.um + Standoff.um + Voidflaeche.px
Df Sum of Sq RSS AIC
- Standoff.um 1 2.81 7187.3 293.14
- Gap.um 1 29.55 7214.0 293.37
<none> 7184.4 295.12
- Voidflaeche.px 1 604.38 7788.8 297.97
Error in terms.formula(formula, data = data) :
'data' argument is of the wrong type
Очевидно, что-то не работает со списком. Но я понятия не имею, что это может быть.
Так как я пытался сделать то же самое с базовым пакетом, который создает ту же модель (по крайней мере, те же коэффициенты). Результаты ниже:
>lin.model<-lm(Scherkraft.N ~ Gap.um + Standoff.um + Voidflaeche.px,df[which(df$id==1),])
# id is in order, so should be the same subgroup as for the first list element in lmList
Coefficients:
(Intercept) Gap.um Standoff.um Voidflaeche.px
62.306133 -0.009878 0.026317 -0.015048
Хорошо, это то, что я получаю, используя stepAIC на моей linear.model.
Насколько мне известно, информационный критерий akaike можно использовать для оценки того, какая модель лучше балансирует между подгонкой и обобщением с учетом некоторых данных.
>stepAIC(lin.model,direction="backward")
Start: AIC=295.12
Scherkraft.N ~ Gap.um + Standoff.um + Voidflaeche.px
Df Sum of Sq RSS AIC
- Standoff.um 1 2.81 7187.3 293.14
- Gap.um 1 29.55 7214.0 293.37
<none> 7184.4 295.12
- Voidflaeche.px 1 604.38 7788.8 297.97
Step: AIC=293.14
Scherkraft.N ~ Gap.um + Voidflaeche.px
Df Sum of Sq RSS AIC
- Gap.um 1 28.51 7215.8 291.38
<none> 7187.3 293.14
- Voidflaeche.px 1 717.63 7904.9 296.85
Step: AIC=291.38
Scherkraft.N ~ Voidflaeche.px
Df Sum of Sq RSS AIC
<none> 7215.8 291.38
- Voidflaeche.px 1 795.46 8011.2 295.65
Call: lm(formula = Scherkraft.N ~ Voidflaeche.px, data = df[which(df$id == 1), ])
Coefficients:
(Intercept) Voidflaeche.px
71.7183 -0.0151
Я прочитал из вывода, что я должен использовать модель: Scherkraft.N ~ Voidflaeche.px , потому что это минимальный AIC. Что ж, было бы неплохо, если бы кто-то смог кратко описать результаты. Мое понимание ступенчатой регрессии (при условии обратного исключения) заключается в том, что все регрессоры включены в исходную модель. Тогда наименее важный исключается. Критерием принятия решения является АПК. и так далее ... Почему-то у меня возникают проблемы с правильной интерпретацией таблиц. Было бы хорошо, если бы кто-то мог подтвердить мою интерпретацию. «-» (минус) обозначает исключенный регрессор. Сверху находится «стартовая» модель, а в таблице под таблицами RSS и AIC рассчитаны на возможные исключения. Итак, первая строка в первой таблице говорит, что модель Scherkraft.N ~ Gap.um + Standoff.um + Voidflaeche.px - Standoff.um приведет к AIC 293.14. Выберите тот без Standoff.um: Scherkraft.N ~ Gap.um + Voidflaeche.px
EDIT:
Я заменил lmList {lme4} на dlply (), чтобы создать список моделей.Все же stepAIC не справляется со списком.Выдает еще одну ошибку.На самом деле, я считаю, что это проблема с данными, которые необходимо пройти через AIC.Мне было интересно, как он вычисляет значение AIC для каждого шага только из данных модели. I будет использовать исходные данные для построения моделей, каждый раз оставляя по одному регрессору.Вследствие этого я бы рассчитал АПК и сравнил.Итак, как работает stepAIC, если он не имеет доступа к исходным данным.(Я не вижу параметр, в котором я передаю исходные данные в stepAIC).Тем не менее, я понятия не имею, почему он работает с простой моделью, но не с моделью, включенной в список.
>model.list.all <- dlply(df, .id, function(x)
{return(lm(Scherkraft.N~Gap.um+Standoff.um+Voidflaeche.px,data=x)) })
>stepAIC(model.list.all[[1]])
Start: AIC=295.12
Scherkraft.N ~ Gap.um + Standoff.um + Voidflaeche.px
Df Sum of Sq RSS AIC
- Standoff.um 1 2.81 7187.3 293.14
- Gap.um 1 29.55 7214.0 293.37
<none> 7184.4 295.12
- Voidflaeche.px 1 604.38 7788.8 297.97
Error in is.data.frame(data) : object 'x' not found