# vcovHC (heteroskedasticity) in pooled and panel probit

I run two types of regressions in R:

1) Using a panel probit of the following form glm_panel <- pglm(formula=Formula,data=Regress_Sample,na.action=na.exclude,family=binomial(probit),index=c("ID","Date")))

2) A pooled probit of the follwing form glm_pooled <- glm(formula=Formula,data=Regress_Sample,na.action=na.exclude,family=binomial(probit)))

I then computed my standard errors in order to account for potential heteroskedasticity in my Regress_Sample

coeftest(glm_pooled, function(x) vcovHC(x, method="arellano",type="HC1"))

It works perfecly with the pooled probit however with the panel probit I get

Error in terms.default(object) : no terms component nor attribute

Any idea why?

Unfortunately it seems like the documentation for pglm is incorrect, and the function does not work as advertised.

The error message is quite simply because the generic function terms does not know how to extract terms in a format it understands from the result of pglm; it results to the default method (terms.default) before failing there.

The return object of pglm is NOT of class pglm (which is what is suggested in ?pglm), and does not contain all of the components listed therein (typically a regression object, e.g. from lm, returns a terms component)

Following the example in ?pglm:

data('Fairness', package = 'pglm')
Parking <- subset(Fairness, good == 'parking')
op <- pglm(as.numeric(answer) ~ education + rule,
Parking[1:105, ],
family = ordinal('probit'), R = 5, print.level = 3,
method = 'bfgs', index = 'id',  model = "random")
class(op)
# [1] "maxLik" "maxim"
names(op)
#  [1] "maximum"     "estimate"    "gradient"    "hessian"     "code"
#  [6] "message"     "last.step"   "fixed"       "iterations"  "type"
# [11] "constraints" "gradientObs" "control"     "call"        "args"
# [16] "model"