2011-08-19 28 views
10

我有一個小N大T面板,我估計通過plm(面板線性迴歸模型),具有固定的影響。R中是否有PLM的預測函數?

有什麼辦法可以獲得新數據集的預測值嗎? (我想 估計我樣本子集上的參數,然後使用這些來計算整個樣本的模型隱含值)。

謝謝!

+0

它似乎在引擎蓋下使用'lm',所以你嘗試調用'predict.lm'嗎? – James

+2

我懷疑作者知道釋放一個'predict.plm'函數會鼓勵那些不瞭解統計問題的人盲目地應用它,當這些假設沒有得到滿足時。 IIRC,lme4軟件包也沒有提供預測功能,而且plm的作者指出他們正在評估隨機和固定組件。 –

+0

predict.lm不起作用。我想有一種方法可以提取係數和攔截,但我想其他人已經遇到過這個問題 –

回答

7

有(至少)兩種方法中的包,以產生從PLM對象估計:

- fixef.plm:提取固定效應

- pmodel.response:一個函數來提取模型。響應

在我看來,作者(s)沒有興趣提供「隨機效應」的估計。這可能是一個問題,「如果你不知道如何自己做,那麼我們不希望給你一把鋒利的刀把自己切得太深。」

2

我寫了一個叫predict.out.plm功能,可以創建爲數據和一個操縱數據集的預測(以相等的列名)。

predict.out.plm計算a)轉換數據的預測(擬合)結果和b)構建根據水平結果的結果。該函數適用於使用plm的一階差分(FD)估算和固定效應(FE)估算。對於FD來說,它隨着時間的推移創造了不同的結果,對於FE來說它創造了時間貶低的結果。

該功能很大程度上未經測試,可能只適用於強平衡的數據幀。

任何建議和更正都非常受歡迎。幫助開發一個小R包將非常感激。

功能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) 
} 

測試的功能:

##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")