在模型列表上使用 stepAIC

Question

我想在线性模型列表上使用 AIC 进行逐步回归。想法是使用 e 线性模型列表，然后在每个列表元素上应用 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”（不仅仅是作为参数传递的模型数据）？
• Q2：我怎么可能知道 stepAIC() 中有另一个参数在帮助页面中没有明确说明？ （也许我的英语太糟糕了，找不到）
• Q3：如何将该参数传递给 stepAIC()？

它必须在 apply 函数的环境中的某个地方并传递数据。 lm() 或 stepAIC() 以及指向原始数据的指针/链接必须在某处丢失。我不太了解 R 中的环境是做什么的。对我来说，这是一种将局部变量与全局变量隔离开来。但也许它更复杂。任何人都可以就上述问题向我解释一下吗？老实说，我没有从R documentation 中读到太多内容。任何更好的理解都会对我有所帮助。

旧： 我在数据帧 df 中有数据，可以分成几个子组。为此，我创建了一个名为 df$id 的 groupID。 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  

嗯，这就是我在 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 仍然无法处理该列表。它抛出另一个错误。实际上，我认为这是 stepAIC 需要运行的数据的问题。我想知道它如何仅根据模型数据计算每个步骤的 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，它会在执行调用之前评估调用中的表达式。这意味着，与其在调用中仅存储d，以便update 和stepAIC 需要找到d 来完成它们的工作，它还存储了数据帧本身的完整表示。

也就是说，做

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 似乎跳出循环环境（即在全局环境中）来寻找它需要的数据，我使用 assign 函数来欺骗它：

    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")
}))