使用lm()对象列表进行预测

Question

我有一些数据，我经常运行回归分析。每个“块”数据都会配合不同的回归方程。例如每个州可能有一个不同的函数来解释因变量。这似乎是一个典型的“拆分-应用-组合”问题，因此我正在使用 plyr 包。我可以轻松创建 lm() 对象的列表，这很有效。但我还不太确定如何稍后使用这些对象来预测另一个数据框中的值。

下面是一个完全人为制造的例子，说明我想要做什么：

# setting up some fake data
set.seed(1)
funct <- function(myState, myYear){
   rnorm(1, 100, 500) +  myState + (100 * myYear) 
}
state <- 50:60
year <- 10:40
myData <- expand.grid( year, state)
names(myData) <- c("year","state")
myData$value <- apply(myData, 1, function(x) funct(x[2], x[1]))
## ok, done with the fake data generation. 

require(plyr)

modelList <- dlply(myData, "state", function(x) lm(value ~ year, data=x))
## if you want to see the summaries of the lm() do this:  
    # lapply(modelList, summary)

state <- 50:60
year <- 50:60
newData <- expand.grid( year, state)
names(newData) <- c("year","state") 
## now how do I predict the values for newData$value 
   # using the regressions in modelList? 

我怎样使用在modelList中包含的lm()对象，利用来自newData的年份和州独立值进行值预测呢？

- JD Long

Answer 1

这里是我的尝试：

predNaughty <- ddply(newData, "state", transform,
  value=predict(modelList[[paste(piece$state[1])]], newdata=piece))
head(predNaughty)
#   year state    value
# 1   50    50 5176.326
# 2   51    50 5274.907
# 3   52    50 5373.487
# 4   53    50 5472.068
# 5   54    50 5570.649
# 6   55    50 5669.229
predDiggsApproved <- ddply(newData, "state", function(x)
  transform(x, value=predict(modelList[[paste(x$state[1])]], newdata=x)))
head(predDiggsApproved)
#   year state    value
# 1   50    50 5176.326
# 2   51    50 5274.907
# 3   52    50 5373.487
# 4   53    50 5472.068
# 5   54    50 5570.649
# 6   55    50 5669.229

JD Long编辑

我受到启发，为adply()编写了一个选项：

pred3 <- adply(newData, 1,  function(x)
    predict(modelList[[paste(x$state)]], newdata=x))
head(pred3)
#   year state        1
# 1   50    50 5176.326
# 2   51    50 5274.907
# 3   52    50 5373.487
# 4   53    50 5472.068
# 5   54    50 5570.649
# 6   55    50 5669.229

- Joshua Ulrich

---

完美！非常感谢。您能解释一下数据框“piece”是从哪里来的吗？它是由ddply自动生成的吗？ - JD Long

@JDLong：.fun 最终在名为 piece 的数据框上调用。但是，正如 @BrianDiggs 在聊天中指出的那样，不应该依赖于此。最好将其包装在匿名函数中（请参见我的更新）。 - Joshua Ulrich

嗨，如果您能看一下我的问题，那就太好了。http://stackoverflow.com/questions/43427392/apply-predict-between-data-frames-within-two-lists。谢谢！ - aaaaa

@JDLong 我能用这种方法获取标准误吗？ - juliamm2011

@juliamm2011 我认为你所要做的就是像这个问题中所说的那样打开 se.fit=TRUE。请注意，现在距离这个问题被回答已经过去了8年，我不再使用 adply，而是可能会使用 broom。 - JD Long

Answer 2

您需要使用mdply在每次函数调用中提供模型和数据：

dataList <- dlply(newData, "state")

preds <- mdply(cbind(mod = modelList, df = dataList), function(mod, df) {
  mutate(df, pred = predict(mod, newdata = df))
})

- hadley

Answer 3

一个只使用base R的解决方案。输出格式不同，但所有值都在其中。

models <- lapply(split(myData, myData$state), 'lm', formula = value ~ year)
pred4  <- mapply('predict', models, split(newData, newData$state))

- Ramnath

---

感谢@ramnath。我非常喜欢将基本的R解决方案与使用软件包完成的解决方案进行比较。这有助于我提高对基本R的理解，同时也能够了解在使用像plyr这样的抽象时所做出的妥协。 - JD Long

这通常是我解决问题的方式 - 使用dlply和mdply函数 - hadley

@hadley 能否给这个案例展示一个实际的例子吗？我尝试使用 mdply 构建了一个，但是因为 .data 必须是矩阵或数据框架，而 predict 的两个参数是一个 lm 对象和一个 data.frame，所以我无法弄清楚如何做。我无法将 lm 对象的列表作为列放入 data.frame 中。我尝试的另一种方法是将 .data 设为列表的列表（.data=list(object=modelList, newData=newDataList)，其中 newDataList <- dlply(newData, .(state), identity)），但是由于文档规定 .data 不是矩阵或数据框架，所以它不起作用。 - Brian Diggs

简而言之，将这两个列表 cbind 在一起。 - hadley

Answer 4

这句话的意思是“出了什么问题”。

lapply(modelList, predict, newData)

这段文字的意思是：

？

编辑：

谢谢您解释了这个问题。那么下面这个怎么翻译呢：

newData <- data.frame(year)
ldply(modelList, function(model) {
  data.frame(newData, predict=predict(model, newData))
})

遍历模型，并应用新数据（由于您刚刚使用 expand.grid 创建它，因此对于每个状态，新数据相同）。

编辑2：

如果newData的year值在每个state中不同于示例中的值，则可以使用更通用的方法。请注意，这使用newData的原始定义，而不是第一个编辑中的定义。

ldply(state, function(s) {
  nd <- newData[newData$state==s,]
  data.frame(nd, predict=predict(modelList[[as.character(s)]], nd))
})

这个输出的前15行是什么？

   year state  predict
1    50    50 5176.326
2    51    50 5274.907
3    52    50 5373.487
4    53    50 5472.068
5    54    50 5570.649
6    55    50 5669.229
7    56    50 5767.810
8    57    50 5866.390
9    58    50 5964.971
10   59    50 6063.551
11   60    50 6162.132
12   50    51 5514.825
13   51    51 5626.160
14   52    51 5737.496
15   53    51 5848.832

- Brian Diggs

---

这正是我一直在想的事情，但这并不是我真正想要的。它将每个模型应用于每个状态。我只想将状态==50的模型应用于状态==50的数据。 - JD Long

Answer 5

我认为难点在于将newData中的每个状态与相应的模型匹配。也许可以尝试以下方法？

predList <- dlply(newData, "state", function(x) {
  predict(modelList[[as.character(min(x$state))]], x) 
})

在这里，我使用了一种“hacky”的方法来提取相应的状态模型：as.character(min(x$state))

......可能有更好的方法吗？

输出：

> predList[1:2]
$`50`
       1        2        3        4        5        6        7        8        9       10       11 
5176.326 5274.907 5373.487 5472.068 5570.649 5669.229 5767.810 5866.390 5964.971 6063.551 6162.132 

$`51`
      12       13       14       15       16       17       18       19       20       21       22 
5514.825 5626.160 5737.496 5848.832 5960.167 6071.503 6182.838 6294.174 6405.510 6516.845 6628.181

或者，如果你需要一个数据框作为输出：

predData <- ddply(newData, "state", function(x) {
  y <-predict(modelList[[as.character(min(x$state))]], x)
  data.frame(id=names(y), value=c(y))
})

输出：

head(predData)
  state id    value
1    50  1 5176.326
2    50  2 5274.907
3    50  3 5373.487
4    50  4 5472.068
5    50  5 5570.649
6    50  6 5669.229

- Tommy

Answer 6

也许我有所遗漏，但我认为lmList是这里的理想工具。

library(nlme)
ll = lmList(value ~ year | state, data=myData)
predict(ll, newData)


## Or, to show that it produces the same results as the other proposed methods...
newData[["value"]] <- predict(ll, newData)
head(newData)
#   year state    value
# 1   50    50 5176.326
# 2   51    50 5274.907
# 3   52    50 5373.487
# 4   53    50 5472.068
# 5   54    50 5570.649
# 6   55    50 5669.229

- baptiste

---

嗯，是的，这似乎是最好的选择！lmList 有自己的 predict() 方法真的很不错。 - Josh O'Brien