在模型列表上使用stepAIC

Question

我想在线性模型列表上使用 AIC 进行逐步回归。想法是使用线性模型的 ea 列表，然后对每个列表元素应用 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 中读到太多内容文档。任何更好的理解都会对我有帮助。

旧： 我的数据框 df 中有数据，可以分为几个子组。为此，我创建了一个名为 df$id 的 groupID。 lm() 返回第一个子组的预期系数。我想使用 AIC 作为每个子组的标准分别进行逐步回归。我使用 lmList {lme4} 为每个子组（id）生成一个模型。但如果我对列表元素使用 stepAIC{MASS} ，则会引发错误。见下文。

所以问题是：我的程序/语法中有什么错误？我得到了单个模型的结果，但没有得到使用 lmList 创建的模型的结果。 lmList() 存储的模型信息是否与 lm() 不同？
但在帮助中它指出： 类“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  

嗯，这是我在 Linear.model 上使用 stepAIC 得到的结果。 据我所知，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。好吧，如果有人能简短地描述一下输出，那就太好了。我对逐步回归（假设向后消除）的理解是所有回归量都包含在初始模型中。然后最不重要的一项就被淘汰了。决定的标准是AIC。等等......不知怎的，我在正确解释表格时遇到了问题。如果有人能证实我的解释，那就太好了。 “-”（减号）代表消除的回归量。顶部是“起始”模型，下面的表格中计算了可能消除的 RSS 和 AIC。因此，第一个表中的第一行表示模型 Scherkraft.N~Gap.um+Standoff.um+Voidflaeche.px - Standoff.um 将导致 AIC 293.14 。选择没有 Standoff.um 的：Scherkraft.N~Gap.um+Voidflaeche.px

编辑：
我用 dlply() 替换 lmList{lme4} 以创建模型列表。 stepAIC 仍然无法处理该列表。它抛出另一个错误。实际上，我认为这是AIC需要运行的数据步骤的问题。我想知道它如何仅根据模型数据计算每个步骤的 AIC 值。 我会利用原始数据构建模型，每次都会留下一个回归量。我会计算 AIC 并进行比较。那么如果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

Answer 1

我不确定版本控制中可能发生了什么变化，导致调试变得如此困难，但一种解决方案是使用 do.call，它在执行调用之前评估调用中的表达式。这意味着，update 和 stepAIC 需要去查找 d，而不是在调用中仅存储 d为了完成他们的工作，它存储数据帧本身的完整表示。

也就是说，使用 do

do.call("lm", list(y~x1+x2+x3, data=d))

而不是

lm(y~x1+x2+x3, data=d)

您可以通过查看模型的 call 元素来了解它正在尝试执行的操作，可能如下所示：

dat.lin.model.lst <- lapply(split(dat, dat$id), function(d)
                            do.call("lm", list(y~x1+x2+x3, data=d)) )
dat.lin.model.lst[[1]]$call

还可以在全局环境中创建数据框列表，并且然后构建调用，以便 update 和 stepAIC 依次查找每个数据帧，因为它们的环境链总是返回到全局环境；像这样：

dats <- split(dat, dat$id)
dat.lin.model.list <- lapply(seq_along(dats), function(d)
            do.call("lm", list(y~x1+x2+x3, data=call("[[", quote(dats),i))) )

要查看更改的内容，请再次运行 dat.lin.model.lst[[1]]$call。

Answer 2

由于stepAIC似乎脱离循环环境（即在全局环境中）来查找它需要的数据，我使用分配函数来欺骗它：

    results <- do.call(rbind, lapply(response, function (i) { 
    assign("i", response, envir = .GlobalEnv)
            mdl <- gls(as.formula(paste0(i,"~",paste(expvar, collapse = "+")), data= parevt, correlation = corARMA(p=1,q=1,form= ~as.integer(Year)), weights= varIdent(~1/Linf_var), method="ML")
            mdl <- stepAIC(mdl, direction ="backward")
}))