2

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

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

OLD： 我在数据帧 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 仍然无法处理该列表。它引发另一个错误。实际上，我认为这是 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
4

2 回答 2

4

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

也就是说，做

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再次运行。

于 2012-04-23T14:32:33.217 回答
3

由于 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")
}))
于 2015-02-03T15:53:23.490 回答