I'm running a random effects model using the plm package and now I need to test for the presence of heteroscedasticity, but I'm not sure how to process it in the mentioned package.
My model:
random <- plm(Y ~ X, data=panel_data, model= "random", effect = "twoways")
One can test for heteroskedasticity and cross-sectional dependence using the plm::pcdtest() function, as documented on page 50 of the plm package vignette. A comprehensive walkthrough illustrating how to interpret the results from plm random and fixed effect models is Getting Started with Fixed and Random Effects Models in R and is available on the Princeton University's Data and Statistical Services website.
Using an example from the plm vignette:
library(plm)
data("Grunfeld", package = "plm")
g <- plm(inv ~ value + capital, data = Grunfeld, index = c("firm", "year"))
pcdtest(g)
...and the results:
> pcdtest(g)
Pesaran CD test for cross-sectional dependence in panels
data: inv ~ value + capital
z = 4.6612, p-value = 3.144e-06
alternative hypothesis: cross-sectional dependence