usando predire con un elenco di oggetti lm()

Question

Ho dati su cui gestisco regolarmente le regressioni. Ogni "frammento" di dati si adatta a una regressione diversa. Ad esempio, ogni stato potrebbe avere una funzione diversa che spiega il valore dipendente. Questo sembra un tipico tipo di problema "split-apply-combine", quindi sto usando il pacchetto plyr. Posso facilmente creare un elenco di oggetti lm() che funziona bene. Tuttavia non riesco a capire come utilizzare questi oggetti in un secondo momento per prevedere i valori in un data.frame separato.

Ecco un esempio del tutto artificiosa che illustra quello che sto cercando di fare:

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

Quindi, come si usa il lm() oggetti contenuti in modelList per prevedere i valori utilizzando l'anno e statali valori indipendenti dalla newData?

Answer 1

Ecco il mio tentativo:

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 lungo modificare

mi sono ispirato abbastanza per elaborare un'opzione 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

Cosa c'è di sbagliato con

lapply(modelList, predict, newData) 

?

EDIT:

Grazie per spiegare che cosa è sbagliato in questo. Che ne dite:

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

iterare i modelli, e applicare i nuovi dati (che è lo stesso per ogni stato dato che hai appena fatto un expand.grid per crearlo).

EDIT 2:

Se newData non ha gli stessi valori per year indipendentemente dai state come nell'esempio, un approccio più generale può essere utilizzato. Si noti che questo utilizza la definizione originale di newData, non quella nella prima modifica.

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

Primi 15 linee di questa uscita:

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

mi prendono la parte difficile sta abbinando ogni stato newData al modello corrispondente.

Forse qualcosa del genere?

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

Qui ho usato un modo "hacky" di estrarre il modello di stato corrispondente: as.character(min(x$state))

... C'è probabilmente un modo migliore?

uscita:

> 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 

Oppure, se si vuole un data.frame come output:

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

uscita:

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

Una soluzione con un solo base R. Il formato dell'output è diverso, ma tutti i valori sono proprio qui.

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

Answer 5

è necessario utilizzare mdply per fornire sia il modello che i dati a ogni chiamata di funzione:

dataList <- dlply(newData, "state") 

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

Answer 6

forse mi manca qualcosa, ma credo lmList è lo strumento ideale qui,

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