The ANOVA in linear regression models is equivalent to the Wald test (and the likelihood ratio test) of the corresponding nested models. So when you want to conduct the corresponding test using heteroskedasticity-consistent (HC) standard errors, this cannot be obtained from a decomposition of the sums of squares but you can carry out the Wald test using a HC covariance estimate. This idea is used in both Anova() and linearHypothesis() from the car package and coeftest() and waldtest() from the lmtest package. The latter three can also be used with plm objects.
A simple (albeit not very interesting/meaningful) example is the following. We use the standard model from the ?plm manual page and want to carry out a Wald test for the significance of both log(pcap) and unemp. We need these packages:
library("plm")
library("sandwich")
library("car")
library("lmtest")
The model (under the alternative) is:
data("Produc", package = "plm")
mod <- plm(log(gsp) ~ log(pc) + log(emp) + log(pcap) + unemp,
data = Produc, index = c("state", "year"))
First, let's look at the marginal Wald tests with HC standard errors for all individual coefficients:
coeftest(mod, vcov = vcovHC)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
log(pc) 0.2920069 0.0617425 4.7294 2.681e-06 ***
log(emp) 0.7681595 0.0816652 9.4062 < 2.2e-16 ***
log(pcap) -0.0261497 0.0603262 -0.4335 0.66480
unemp -0.0052977 0.0024958 -2.1226 0.03411 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
And then we carry out a Wald test for both log(pcap) and unemp:
linearHypothesis(mod, c("log(pcap)", "unemp"), vcov = vcovHC)
Linear hypothesis test
Hypothesis:
log(pcap) = 0
unemp = 0
Model 1: restricted model
Model 2: log(gsp) ~ log(pc) + log(emp) + log(pcap) + unemp
Note: Coefficient covariance matrix supplied.
Res.Df Df Chisq Pr(>Chisq)
1 766
2 764 2 7.2934 0.02608 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Alternatively, we can also fit the model under the null hypothesis (mod0 say) without the two coefficients and then call waldtest():
mod0 <- plm(log(gsp) ~ log(pc) + log(emp),
data = Produc, index = c("state", "year"))
waldtest(mod0, mod, vcov = vcovHC)
Wald test
Model 1: log(gsp) ~ log(pc) + log(emp)
Model 2: log(gsp) ~ log(pc) + log(emp) + log(pcap) + unemp
Res.Df Df Chisq Pr(>Chisq)
1 766
2 764 2 7.2934 0.02608 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The test statistic and p-value computed by linearHypothesis() and waldtest() is exactly the same. Just the interface and output formatting is somewhat different. In some cases one or the other is simpler or more intuitive.
Note: If you supply a covariance matrix estimate (i.e., a matrix like vocvHC(mod)) instead of a covariance matrix estimator (i.e., a function like vocvHC), make sure that you supply the HC covariance matrix estimate of the model under the alternative, i.e., the non-restricted model.
It's hard to comment on your overall goal without data. But, the reason your random effects go down with adding more predictors is a consequence of explaining away doctor variance. For example, suppose that one of your doctors treats a large number of end-stage renal disease patients, who are very sick. Without end-stage renal disease as a fixed effect, that doctors random effect will likely be large to explain the outcomes of treating those sick patients. Whereas, if you add the fixed effect, you'll correctly adjust for these sick patients, and therefore reduce that doctors random effect.
Most likely you just don't have enough data to form more accurate predictions of your random effects.
Best Answer
Look at % variance farther up in the tutorial where the R-squared values are reported. You look at the model for the quantity of variance explained and the test (
anova) for whether the increase across model complexity is significant.