将 predict 与 lm() 对象列表一起使用

Posted 2023-02-16

【中文标题】将 predict 与 lm() 对象列表一起使用

【英文标题】：using predict with a list of lm() objects

【发布时间】：2012-01-19 19:05:24

【问题描述】：

我有一些数据，我会定期对其进行回归分析。每个“块”数据都适合不同的回归。例如，每个状态可能有不同的函数来解释依赖值。这似乎是一个典型的“拆分应用组合”类型的问题，所以我使用的是 plyr 包。我可以轻松创建一个运行良好的lm() 对象列表。但是，我不能完全理解以后如何使用这些对象来预测单独 data.frame 中的值。

这是一个完全人为的例子，说明了我正在尝试做的事情：

# 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 的年份和状态独立值的值？

【参考方案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

京东长编辑

我受到了足够的启发，想出了一个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

【讨论】：

这完全可以！非常感谢。您能解释一下 data.frame piece 的来源吗？它是由 ddply 自动生成的吗？

@JDLong: .fun 最终在名为 piece 的数据帧上调用。但是，正如@BrianDiggs 在聊天中指出的那样，这不应该被依赖。最好封装在匿名函数中（请参阅我的更新）。

嗨，如果你能看看我的问题，那就太好了***.com/questions/43427392/…。谢谢！

@JDLong 我能用这种方法得到标准错误吗？

@juliamm2011 我认为您所要做的就是根据这个问题转se.fit=TRUE：***.com/a/33660779/37751 请注意，在这个问题得到回答后我们现在已经 8 年了，我将不再使用 @ 987654330@ 不再是，但可能会使用broom

【参考方案2】：

只有base R 的解决方案。输出的格式不同，但所有值都在那里。

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

【讨论】：

感谢@ramnath。我真的很喜欢将基本 R 解决方案与使用包完成的解决方案进行比较。它既可以帮助我提高对 R 基础的理解，也可以理解我在使用 plyr 之类的抽象时所做的妥协。

这就是我通常解决问题的方式 - 但使用 dlply 和 mdply

@hadley 您能否为这个案例展示一个工作示例？我尝试用mdply 构建一个，但不知道怎么做，因为.data 必须是矩阵或data.frame，predict 的两个参数是lm 对象和data.frame .我无法将lm 对象列表填充为data.frame 中的一列。我尝试的另一种方法是将.data 设为列表列表（.data=list(object=modelList, newData=newDataList) 其中newDataList <- dlply(newData, .(state), identity)）不起作用，因为.data 不是矩阵或data.frame（根据文档）。跨度>

简而言之，将两个列表绑定在一起

【参考方案3】：

您需要使用mdply 为每个函数调用提供模型和数据：

dataList <- dlply(newData, "state")

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

【参考方案4】：

有什么问题

lapply(modelList, predict, newData)

?

编辑：

感谢您解释这有什么问题。怎么样：

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

迭代模型，并应用新数据（这对于每个状态都是相同的，因为您刚刚使用expand.grid 来创建它）。

编辑 2：

如果newData 对于每个state 的year 值与示例中的不同，则可以使用更通用的方法。请注意，这使用了 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

【讨论】：

这正是我一直在做的事情，但这并不是我真正想要的。这将每个模型应用于每个州。我只希望将 state==50 的模型应用于 state==50

的数据

【参考方案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

或者，如果您想要 data.frame 作为输出：

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

【参考方案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

【讨论】：

嗯，是的，这似乎是最好的！ lmList 有自己的 predict() 方法真是太好了。