Ken S. - 3 months ago 10

R Question

I want to test a fixed effects regression for heteroscedasticity with

`lmtest::bptest`

`bptest`

It seems as if

`bptest`

`formula`

`?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.

No error or warning is produced to tell you that the output of the function is not based on the model you used as input.

First of all, is there a way to test my plm-object for heteroscedasticity with

`pbtest`

Here is a reproducible example:

`# 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

Including country dummies into the OLS regression to create a fixed effects regression also didn't work, because the country dummies increase the degrees of freedom of the test:

`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

Thanks in advance!

Answer

You will need a wrapper around `lmtest::bptest`

to run it on plm object's (quasi-)demeaned data. Here is one, I called it `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()
```