2

I would like to extract the variance-covariance matrix from a simple plm fixed effects model. For example:

library(plm)
data("Grunfeld")

M1 <- plm(inv ~ lag(inv) + value + capital, index = 'firm',
          data = Grunfeld)

The usual vcov function gives me:

vcov(M1)

              lag(inv)         value       capital
lag(inv)  3.561238e-03 -7.461897e-05 -1.064497e-03
value    -7.461897e-05  9.005814e-05 -1.806683e-05
capital  -1.064497e-03 -1.806683e-05  4.957097e-04

plm's fixef function only gives:

fixef(M1)
          1           2           3           4           5           6           7 
-286.876375  -97.190009 -209.999074  -53.808241  -59.348086  -34.136422  -34.397967 
      8           9          10 
-65.116699  -54.384488   -6.836448

Any help extracting the variance-covariance matrix that includes the fixed effects would be much appreciated.

Christopher Gandrud
  • 88
  • 1
  • 8
  • Anyway, maybe there is method for the `plm` class, or something strange here because I get different results with you code. For me `vcov(M1)` gives the same result as `M1$vcov` but you result is very different. – SabDeM Jul 17 '15 at 18:54

2 Answers2

1

Using names sometimes is very useful:

 names(M1)
[1] "coefficients" "vcov"         "residuals"    "df.residual" 
[5] "formula"      "model"        "args"         "call"        
 M1$vcov
              lag(inv)         value       capital
lag(inv)  1.265321e-03  3.484274e-05 -3.395901e-04
value     3.484274e-05  1.336768e-04 -7.463365e-05
capital  -3.395901e-04 -7.463365e-05  3.662395e-04
SabDeM
  • 7,050
  • 2
  • 25
  • 38
  • Right, the difference I think just has to do with the different methods available in `plm` for finding the covariance matrix. There are the the Beck and Katz, Heteroskedasticity-consistent, and Driscoll and Kraay approaches at least. My issue is that I would like to find them for the fixed effect parameters as well as the three main covariates. – Christopher Gandrud Jul 20 '15 at 13:11
1

Picking up your example, do the following to get the standard errors (if that is what you are interested in; it is not the whole variance-covariance matrix):

library(plm)
data("Grunfeld")

M1 <- plm(inv ~ lag(inv) + value + capital, index = 'firm',
          data = Grunfeld)

fix <- fixef(M1)
fix_se <- attr(fix, "se")
fix_se
        1         2         3         4         5         6         7         8         9        10 
43.453642 25.948160 20.294977 11.245009 12.472005  9.934159 10.554240 11.083221 10.642589  9.164694

You can also use the summary function for more info:

 summary(fix)
    Estimate Std. Error  t-value  Pr(>|t|)    
1  -286.8764    43.4536  -6.6019 4.059e-11 ***
2   -97.1900    25.9482  -3.7455 0.0001800 ***
3  -209.9991    20.2950 -10.3473 < 2.2e-16 ***
4   -53.8082    11.2450  -4.7851 1.709e-06 ***
5   -59.3481    12.4720  -4.7585 1.950e-06 ***
6   -34.1364     9.9342  -3.4363 0.0005898 ***
7   -34.3980    10.5542  -3.2592 0.0011174 ** 
8   -65.1167    11.0832  -5.8753 4.222e-09 ***
9   -54.3845    10.6426  -5.1101 3.220e-07 ***
10   -6.8364     9.1647  -0.7460 0.4556947

Btw, the documentation expains the "se" attribute:

Value An object of class "fixef". It is a numeric vector containing the fixed effects with attribute se which contains the standard errors. [...]"

Note: You might need the latest development version for that because much has improved there about fixef: https://r-forge.r-project.org/R/?group_id=406

Helix123
  • 3,502
  • 2
  • 16
  • 36