Ho un piccolo pannello N grande T che sto stimando tramite plm (modello di regressione lineare del pannello), con effetti fissi.Esiste una funzione di previsione per PLM in R?
C'è un modo per ottenere valori previsti per un nuovo set di dati? (Voglio i parametri di stima su un sottoinsieme del mio esempio, e quindi usarli per calcolare i valori impliciti del modello per l'intero campione).
Grazie!
Ci sono (almeno) due metodi nel pacchetto di produrre stime dagli oggetti PLM:
- fixef.plm: estrarre gli effetti fissi
- pmodel.response: Una funzione per estrarre the model.response
Mi sembra che l'autore (i) non sia interessato a fornire stime per gli "effetti casuali". Potrebbe essere una questione di "se non sai come farlo da solo, allora non vogliamo darti un coltello affilato per tagliarti troppo profondamente".
ho scritto una funzione chiamata predict.out.plm che può creare le previsioni per i dati originali e per un insieme di dati manipolati(con nomi uguali colonna).
Il predict.out.plm calcola a) l'esito previsto (adattato) dei dati trasformati eb) costruisce il secondo livello di esito. La funzione funziona per le stime di First Difference (FD) e per gli effetti fissi (FE) utilizzando plm. Per FD crea il risultato differenziato nel tempo e per FE genera il risultato decaduto nel tempo.
La funzione è in gran parte non testata e probabilmente funziona solo con frame di dati fortemente bilanciati.
Eventuali suggerimenti e correzioni sono i benvenuti. Aiutare a sviluppare un piccolo pacchetto R sarebbe molto apprezzato.
La funzione predict.out.plm
predict.out.plm<-function(
estimate,
formula,
data,
model="fd",
pname="y",
pindex=NULL,
levelconstr=T
){
# estimate=e.fe
# formula=f
# data=d
# model="within"
# pname="y"
# pindex=NULL
# levelconstr=T
#get index of panel data
if (is.null(pindex) && class(data)[1]=="pdata.frame") {
pindex<-names(attributes(data)$index)
} else {
pindex<-names(data)[1:2]
}
if (class(data)[1]!="pdata.frame") {
data<-pdata.frame(data)
}
#model frame
mf<-model.frame(formula,data=data)
#model matrix - transformed data
mn<-model.matrix(formula,mf,model)
#define variable names
y.t.hat<-paste0(pname,".t.hat")
y.l.hat<-paste0(pname,".l.hat")
y.l<-names(mf)[1]
#transformed data of explanatory variables
#exclude variables that were droped in estimation
n<-names(estimate$aliased[estimate$aliased==F])
i<-match(n,colnames(mn))
X<-mn[,i]
#predict transformed outcome with X * beta
# p<- X %*% coef(estimate)
p<-crossprod(t(X),coef(estimate))
colnames(p)<-y.t.hat
if (levelconstr==T){
#old dataset with original outcome
od<-data.frame(
attributes(mf)$index,
data.frame(mf)[,1]
)
rownames(od)<-rownames(mf) #preserve row names from model.frame
names(od)[3]<-y.l
#merge old dataset with prediciton
nd<-merge(
od,
p,
by="row.names",
all.x=T,
sort=F
)
nd$Row.names<-as.integer(nd$Row.names)
nd<-nd[order(nd$Row.names),]
#construct predicted level outcome for FD estiamtions
if (model=="fd"){
#first observation from real data
i<-which(is.na(nd[,y.t.hat]))
nd[i,y.l.hat]<-NA
nd[i,y.l.hat]<-nd[i,y.l]
#fill values over all years
ylist<-unique(nd[,pindex[2]])[-1]
ylist<-as.integer(as.character(ylist))
for (y in ylist){
nd[nd[,pindex[2]]==y,y.l.hat]<-
nd[nd[,pindex[2]]==(y-1),y.l.hat] +
nd[nd[,pindex[2]]==y,y.t.hat]
}
}
if (model=="within"){
#group means of outcome
gm<-aggregate(nd[, pname], list(nd[,pindex[1]]), mean)
gl<-aggregate(nd[, pname], list(nd[,pindex[1]]), length)
nd<-cbind(nd,groupmeans=rep(gm$x,gl$x))
#predicted values + group means
nd[,y.l.hat]<-nd[,y.t.hat] + nd[,"groupmeans"]
}
if (model!="fd" && model!="within") {
stop('funciton works only for FD and FE estimations')
}
}
#results
results<-p
if (levelconstr==T){
results<-list(results,nd)
names(results)<-c("p","df")
}
return(results)
}
Test della funzione:
##packages
library(plm)
##test dataframe
#data structure
N<-4
G<-2
M<-5
d<-data.frame(
id=rep(1:N,each=M),
year=rep(1:M,N)+2000,
gid=rep(1:G,each=M*2)
)
#explanatory variable
d[,"x"]=runif(N*M,0,1)
#outcome
d[,"y"] = 2 * d[,"x"] + runif(N*M,0,1)
#panel data frame
d<-pdata.frame(d,index=c("id","year"))
##new data frame for out of sample prediction
dn<-d
dn$x<-rnorm(nrow(dn),0,2)
##estimate
#formula
f<- pFormula(y ~ x + factor(year))
#fixed effects or first difffernce estimation
e<-plm(f,data=d,model="within",index=c("id","year"))
e<-plm(f,data=d,model="fd",index=c("id","year"))
summary(e)
##fitted values of estimation
#transformed outcome prediction
predict(e)
c(pmodel.response(e)-residuals(e))
predict.out.plm(e,f,d,"fd")$p
# "level" outcome prediciton
predict.out.plm(e,f,d,"fd")$df$y.l.hat
#both
predict.out.plm(e,f,d,"fd")
##out of sampel prediciton
predict(e,newdata=d)
predict(e,newdata=dn)
# Error in crossprod(beta, t(X)) : non-conformable arguments
# if plm omits variables specified in the formula (e.g. one year in factor(year))
# it tries to multiply two matrices with different length of columns than regressors
# the new funciton avoids this and therefore is able to do out of sample predicitons
predict.out.plm(e,f,dn,"fd")
Sembra che ci sia un nuovo pacchetto di fare previsioni in-campione per una varietà di modelli tra cui plm
https://cran.r-project.org/web/packages/prediction/prediction.pdf
Sembra che stia usando 'lm' sotto il cofano, quindi hai provato a chiamare' predicti.lm'? – James
Sospetto che gli autori sappiano che il rilascio di una funzione 'predict.plm' incoraggerebbe le persone che non capiscono i problemi statistici di applicarla ciecamente quando le ipotesi non sono soddisfatte. IIRC, il pacchetto lme4 non fornisce una funzione di previsione e gli autori del plm notano che stanno valutando sia componenti casuali che fissi. –
predicict.lm non funziona. Suppongo che ci sia un modo per estrarre i coefficienti e le intercettazioni ma immagino che altri abbiano già riscontrato questo problema –