lmtest::bptest 测试 OLS 残差,也将 FE 作为输入

lmtest::bptest tests OLS residuals, also for FE as input

我想用 lmtest::bptest 测试异方差的固定效应回归。我从 bptest 得到的测试结果对于固定效应回归和 OLS 是相同的。

似乎 bptest 自动从输入对象中提取 formula 并且 运行 是具有这些变量的 OLS,而不管输入的回归模型如何。来自 ?bptest:

The Breusch-Pagan test fits a linear regression model to the residuals of a linear regression model (by default the same explanatory variables are taken as in the main regression model) and rejects if too much of the variance is explained by the additional explanatory variables.

没有错误或警告告诉您函数的输出不是基于您用作输入的模型。

问题

首先,有没有办法用 pbtest 来测试我的 plm 对象的异方差性?由于我假设没有,是否有办法 运行 固定效应回归的 Breusch Pagan 测试?

这是一个可重现的例子:

# load the data:
df <- structure(list(country = c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5), year = c(2010, 
         2015, 2010, 2015, 2010, 2015, 2010, 2015, 2010, 2015), dv1 = c(28.61, 
         31.13, 38.87, 39.46, 68.42, 70.39, 79.36, 80.55, 70.14, 71.48
         ), iv1 = c(-20.68, 0, NA, NA, -19.41, -18.73, 24.98, 25.23, 21.23, 
         -21.06), iv2 = c(-4.23, NA, NA, NA, -4.92, -4.22, 9.19, 9.37, 
         4.15, -3.92)), .Names = c("country", "year", "dv1", "iv1", "iv2"
         ), row.names = c(2L, 3L, 5L, 6L, 8L, 9L, 11L, 12L, 14L, 15L),class    
         ="data.frame")
library(plm)
library(lmtest)


# Run the fixed effects regression
FEoutput <- plm(dv1 ~ iv1 + iv2, data = df, 
             model = "within", index = c("country", "year"))
bptest(FEoutput)
#   studentized Breusch-Pagan test
#
#    data:  FEoutput
#    BP = 1.9052, df = 2, p-value = 0.3857


# Run the OLS regression
OLSoutput <- lm(dv1 ~ iv1 + iv2, data = df)
bptest(OLSoutput)
#   studentized Breusch-Pagan test
#
#    data:  OLSoutput
#    BP = 1.9052, df = 2, p-value = 0.3857

将国家虚拟人纳入 OLS 回归以创建固定效应回归也不起作用,因为国家虚拟人增加了检验的自由度:

OLSctry <- lm(dv1 ~ iv1 + iv2 + factor(country), data = df)
bptest(OLSctry)
# studentized Breusch-Pagan test
#
# data:  OLSctry
# BP = 7, df = 5, p-value = 0.2206 

提前致谢!

您需要在 plm 对象的(准)贬低数据上围绕 lmtest::bptest 到 运行 进行包装。这是一个,我称之为 pbptest:

pbptest <-function(x, ...) {
  ## residual heteroskedasticity test based on the residuals of the demeaned
  ## model and the regular bptest() in {lmtest}

  ## structure:
  ## 1: take demeaned data from 'plm' object
  ## 2: est. auxiliary model by OLS on demeaned data
  ## 3: apply bptest() to auxiliary model and return the result

  if (!inherits(x, "plm")) stop("need to supply a panelmodel estimated with plm()")
  model <- plm:::describe(x, "model")
  effect <- plm:::describe(x, "effect")
  theta <- x$ercomp$theta

  ## retrieve demeaned data
  demX <- model.matrix(x, model = model, effect = effect, theta = theta)
  demy <- pmodel.response(model.frame(x), model = model, effect = effect, theta = theta)

  Ti <- pdim(x)$Tint$Ti

  if (is.null(order)) order <- min(Ti)

  ## bgtest on the demeaned model:

  ## check package availability and load if necessary
  lm.ok <- require("lmtest")
  if(!lm.ok) stop("package lmtest is needed but not available")

  ## pbptest is the bptest, exception made for the method attribute
  dots <- match.call(expand.dots=FALSE)[["..."]]
  if (!is.null(dots$type)) type <- dots$type else type <- "Chisq"
  if (!is.null(dots$order.by)) order.by <- dots$order.by else order.by <- NULL

  auxformula <- demy~demX-1
  lm.mod <- lm(auxformula)
  return(lmtest::bptest(lm.mod, ...)) # call and return lmtest::bptest
} # END pbptest()