используя прогноз со списком объектов lm()

Question

У меня есть данные, которые я регулярно запускаю регрессии. Каждый «кусок» данных подходит для другой регрессии. Каждое состояние, например, может иметь другую функцию, которая объясняет зависимое значение. Это похоже на типичную проблему типа «split-apply-comb», поэтому я использую пакет 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? 

Так как я использую lm() объекты, содержащиеся в modelList для прогнозирования значений, используя год и государственные независимые значения от newData?

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 

Answer 2

Что с

lapply(modelList, predict, newData) 

?

EDIT:

Спасибо за объяснения, что случилось с этим. Как насчет:

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

перебирать модели, и применять новые данные (что то же самое для каждого государства, так как вы просто сделали expand.grid создать его).

EDIT 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 

Answer 3

Я принимаю это твердая часть соответствие каждого состояния, в 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 

Answer 4

Раствор только с base R. Формат выходных данных отличается, но все значения прямо там.

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

Answer 5

Вы должны использовать mdply поставлять как модель и данные для каждого вызова функции:

dataList <- dlply(newData, "state") 

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

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