预测测试数据,使用 R 中的 plm 包,并计算测试数据的 RMSE

Predict on test data, using plm package in R, and calculate RMSE for test data

我使用 plm 包构建了一个模型。示例数据集是 here.

我正在尝试预测测试数据并计算指标。

# Import package
library(plm)
library(tidyverse)
library(prediction)
library(nlme)

# Import data 
df <- read_csv('Panel data sample.csv')

# Convert author to character
df$Author <- as.character(df$Author) 

# Split data into train and test
df_train <- df %>% filter(Year != 2020) # 2017, 2018, 2019
df_test <- df %>% filter(Year == 2020) # 2020

# Convert data
panel_df_train <- pdata.frame(df_train, index = c("Author", "Year"), drop.index = TRUE, row.names = TRUE)
panel_df_test <- pdata.frame(df_train, index = c("Author", "Year"), drop.index = TRUE, row.names = TRUE)

# Create the first model
plmFit1 <- plm(Score ~ Articles, data = panel_df_train)

# Print
summary(plmFit1)

# Get the RMSE for train data
sqrt(mean(plmFit1$residuals^2))

# Get the MSE for train data
mean(plmFit1$residuals^2)

现在我正在尝试计算测试数据的指标

首先,我尝试使用 prediction package 中的 prediction(),它有一个 plm.

的选项
predictions <- prediction(plmFit1, panel_df_test)

出现错误:

Error in crossprod(beta, t(X)) : non-conformable arguments

我阅读了以下问题:

我也读了this question,但是

fitted <- as.numeric(plmFit1$model[[1]] - plmFit1$residuals) 给出了与我的训练或测试数字不同数量的值。

关于out-of-sample固定效应模型的预测,不清楚如何处理原始模型中不存在的与固定效应相关的数据,例如原始数据集中不包含的个体数据该模型估计在。 (与其说这是一个编程问题,不如说这是一个方法论问题)。

plm (https://github.com/ycroissant/plm) 的开发版本现在允许 predict 具有原始数据和 out-of-sample 数据的固定效应模型(参见 ?predict.plm).

下面找一个例子,有10家公司做模型估计,预测的数据中有一个公司没有包含在原始数据集中(除了那个公司,还有一些年份没有包含在原始模型对象中,但是这些在这里无关紧要,因为它是一个 one-way 个人模型)。目前尚不清楚 out-of-sample 公司的固定效应是什么。因此,默认情况下,没有给出预测值(NA 值)。如果参数 na.fill 设置为 TRUE,则使用原始模型对象中包含的固定效应的(加权)平均值作为最佳猜测。

library(plm)
data("Grunfeld", package = "plm")

# fit a fixed effect model
fit.fe <- plm(inv ~ value + capital, data = Grunfeld, model = "within")

# generate 55 new observations of three firms used for prediction:
#  * firm 1 with years 1935:1964 (has out-of-sample years 1955:1964), 
#  * firm 2 with years 1935:1949 (all in sample),
#  * firm 11 with years 1935:1944 (firm 11 is out-of-sample)
set.seed(42L)

new.value2   <- runif(55, min = min(Grunfeld$value),   max = max(Grunfeld$value))
new.capital2 <- runif(55, min = min(Grunfeld$capital), max = max(Grunfeld$capital))

newdata <- data.frame(firm = c(rep(1, 30), rep(2, 15), rep(11, 10)),
                      year = c(1935:(1935+29), 1935:(1935+14), 1935:(1935+9)),
                      value = new.value2, capital = new.capital2)
# make pdata.frame
newdata.p <- pdata.frame(newdata, index = c("firm", "year"))

## predict from fixed effect model with new data as pdata.frame
predict(fit.fe, newdata = newdata.p) # has NA values for the 11'th firm

## set na.fill = TRUE to have the weighted mean used to for fixed effects -> no NA values
predict(fit.fe, newdata = newdata.p, na.fill = TRUE)

注意:当您输入 data.frame 为 newdata 时,不清楚数据与个人和时间段的关系,这就是为什么固定效应的加权平均值来自原始模型对象用于 newdata 中的所有观察,并打印一条警告。对于固定效应模型预测,假设用户可以提供信息(通过 pdata.frame)用户想要用于预测的数据如何与面板数据的个体和时间维度相关。