Есть ли функция прогнозирования для PLM в R?

У меня есть маленькая N большая T панель, которую я оцениваю с помощью plm (модель линейной регрессии панели), с фиксированными эффектами.

Есть ли способ получить прогнозируемые значения для нового набора данных? (Я хочу оценить параметры в подмножестве моей выборки, а затем использовать их для расчета значений, подразумеваемых моделью, для всей выборки).

Спасибо!

Источник

4 ответа

Решение

В пакете есть (как минимум) два метода для получения оценок из объектов plm:

- fixef.plm: извлечь фиксированные эффекты

- pmodel.response: функция для извлечения model.response

Мне кажется, что автор (ы) не заинтересованы в предоставлении оценок для "случайных эффектов". Это может быть вопросом "если вы не знаете, как это сделать самостоятельно, то мы не хотим давать вам острый нож, чтобы порезаться слишком глубоко".

Источник

Я написал функцию под названием predict.out.plm который может создавать прогнозы для исходных данных и для манипулируемого набора данных (с одинаковыми именами столбцов).

predict.out.plm вычисляет а) прогнозируемый (подобранный) результат преобразованных данных и б) строит в соответствии с уровнем результата. Функция работает для оценок First Difference (FD) и Fixed Effects (FE), используя plm, Для 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")
Источник

plm теперь имеет predict.plm() функция, хотя она не документирована / не экспортирована.

Обратите внимание, что predict работает с преобразованной моделью (т. е. после выполнения внутри / между /fd преобразования), а не с исходной. Я предполагаю, что причина этого в том, что прогнозировать в каркасных данных панели сложнее. В самом деле, вам нужно подумать, предсказываете ли вы:

  • новые периоды времени, для существующего человека, и вы использовали индивидуальный FE? Затем вы можете добавить прогноз к существующему индивидуальному среднему значению.
  • новые периоды времени, для нового человека? Тогда вам нужно выяснить, какого человека вы собираетесь использовать?
  • то же самое еще сложнее, если вы используете модель со случайным эффектом, так как эффекты не легко получить

В приведенном ниже коде я иллюстрирую, как использовать подогнанные значения в существующем примере:

library(plm)
#> Loading required package: Formula
library(tidyverse)

data("Produc", package = "plm")
zz <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp,
          data = Produc, index = c("state","year"))


## produce a dataset of prediction, added to the group means
Produc_means <- Produc %>% 
  mutate(y = log(gsp)) %>% 
  group_by(state) %>% 
  transmute(y_mean = mean(y),
            y = y, 
            year = year) %>% 
  ungroup() %>% 
  mutate(y_pred = predict(zz) + y_mean) %>% 
  select(-y_mean)

## plot it
Produc_means %>% 
  gather(type, value, y, y_pred) %>% 
  filter(state %in% toupper(state.name[1:5])) %>% 
  ggplot(aes(x = year, y = value, linetype = type))+
  geom_line() +
  facet_wrap(~state) +
  ggtitle("Visualising in-sample prediction, for 4 states")
#> Warning: attributes are not identical across measure variables;
#> they will be dropped

Создано в 2018-11-20 пакетом представлением (v0.2.1)

Источник

Похоже, что есть новый пакет для прогнозирования в различных образцах, включая plm.

https://cran.r-project.org/web/packages/prediction/prediction.pdf

Источник

Вы можете рассчитать остатки через residuals(reg_name). Отсюда вы можете вычесть их из переменной ответа и получить прогнозируемые значения.

Источник
Другие вопросы по тегам