używanie predict z listą obiektów lm()

Mam dane, na których regularnie przeprowadzam regresje. Każdy "kawałek" danych dopasowuje się do innej regresji. Każdy stan, na przykład, może mieć inną funkcję, która wyjaśnia wartość zależną. Wydaje się to typowym problemem typu" split-apply-combine", więc używam pakietu plyr. Mogę łatwo utworzyć listę lm() obiektów, która działa dobrze. Nie mogę jednak do końca zrozumieć, w jaki sposób używam tych obiektów do przewidywania wartości w osobnych danych.rama.

Oto totalnie wymyślony przykład ilustrujący to, co próbuję zrobić: 

# 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?

Więc jak wykorzystać lm() obiekty zawarte w modelList do przewidywania wartości za pomocą roku i podać niezależne wartości z newData?

Author: Shreyas Karnik, 2011-12-14
Source

6 answers

Oto moja próba: 

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 edit

Zainspirowałem się na tyle, aby opracować adply() opcję:

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
 9
Author: Joshua Ulrich,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2011-12-13 23:54:12

Rozwiązanie z tylko base R. format wyjścia jest inny, ale wszystkie wartości są tam.

models <- lapply(split(myData, myData$state), 'lm', formula = value ~ year)
pred4  <- mapply('predict', models, split(newData, newData$state))
 6
Author: Ramnath,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2011-12-14 04:43:54

Musisz użyć mdply, aby dostarczyć zarówno model, jak i dane do każdego wywołania funkcji: 

dataList <- dlply(newData, "state")

preds <- mdply(cbind(mod = modelList, df = dataList), function(mod, df) {
  mutate(df, pred = predict(mod, newdata = df))
})
 6
Author: hadley,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2011-12-16 12:56:45

Co jest nie tak z 

lapply(modelList, predict, newData)

?

EDIT:

Dzięki za wyjaśnienie, co w tym złego. A może: 
newData <- data.frame(year)
ldply(modelList, function(model) {
  data.frame(newData, predict=predict(model, newData))
})

Iteruj modele i zastosuj nowe dane(które są takie same dla każdego stanu, ponieważ zrobiłeś expand.grid, aby go utworzyć).

EDIT 2:

Jeśli newData nie ma takich samych wartości dla year dla każdego state, jak w przykładzie, można zastosować bardziej ogólne podejście. Zauważ, że używa się oryginalnej definicji newData, a nie jeden w pierwszej edycji.

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

Pierwsze 15 linii tego wyjścia: 

   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
 4
Author: Brian Diggs,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2011-12-13 23:15:05

Rozumiem, że najtrudniejsze jest dopasowanie każdego stanu w newData do odpowiedniego modelu.

Może coś takiego? 
predList <- dlply(newData, "state", function(x) {
  predict(modelList[[as.character(min(x$state))]], x) 
})

Tutaj użyłem "hakerskiego" sposobu wyodrębnienia odpowiedniego modelu stanu: as.character(min(x$state))

...Jest chyba lepszy sposób?

Wyjście: 

> 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

Lub, jeśli chcesz data.frame jako wyjście: 

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

Wyjście: 

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
 2
Author: Tommy,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2011-12-13 23:11:41

Może czegoś mi brakuje, ale wierzę, że lmList jest idealnym narzędziem tutaj, 

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
 1
Author: baptiste,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2014-07-23 12:59:36