Does the Granger Causality test in the "vars" package make sense? We all understand that the Granger Causality test entails constructing two models. The first one is simply an autoregressive model with $y_{t-1}$ being the single independent variable. The second one is adding an $x_{t-1}$ to the autoregressive model. Next, you check whether the residuals of those two models are different using an $F$-test. If they are, you can say that $x$ Granger-causes $y$. 
However, within the "vars" package in R, when you use the causality function it works differently.  That's because you can only specify what is $x$. The package's author calls it a "cause" variable. But you can't specify what is $y$, or the response variable. So when the R output comes out, you get that $x$ Granger-causes several different variables simultaneously with one single given $p$-value. 
Within the framework of Granger Causality as depicted above, this does not seem to make much sense. I can understand variable $A$ Granger-causes variable $B$. But I don't understand how variable $A$ Granger-causes variables $B$, $C$, $D$ with a single $p$-value estimate just as if you had done $A$ Granger-causes $B$ alone.    
 A: When I first used the causality function of the vars package I had the same doubt. Here is what I thought.
Imagine a trivariate VAR(1) model:
$$
Y_t = a_0 + a_1 Y_{t-1} + a_2 X_{t-1} + a_3 Z_{t-1} + \epsilon_{y,t} \\
X_t = b_0 + b_1 Y_{t-1} + b_2 X_{t-1} + b_3 Z_{t-1} + \epsilon_{x,t} \\
Z_t = c_0 + c_1 Y_{t-1} + c_2 X_{t-1} + c_3 Z_{t-1} + \epsilon_{z,t}
$$
In order to $X_t$ not Granger cause $Y_t$ you need to make sure that $H_0: a_2 = a_3 = 0$ or $H_0: a_2 = c_2 = 0$. You can build an $F$-test for that, but what the causality function seems to be doing is checking if $X_t$ Granger causes all other variables in the model. Why would they do that? Maybe because it is simpler. That would be the same as to just check if $H_0: a_2 = c_2 = 0$. This is a much more simple $F$-test, and it extends easily to more than three variables. So, if I am correct, this is not a bug, they simple took the easiest path. So if the null is not rejected, this means that $X_t$ does not Granger cause $Y_t$ and $Z_t$. But if the null is rejected, the test doesn't say much, $X_t$ may be Granger causing either $Y_t$, $Z_t$ or both.
Note that in the help of the causality function they only show a bivariate case, but from that example you can infer that the trivariate case would be as I described. To make sure that this is the case, one can build the corresponding $F$-test and check if the values are equal. Also, if you want to specify the response variable you can use the grangertest function of the lmtest package.
A: As Regis suggested, using the function grangertest from the package lmtest is a way to produce the pairwise test result. This should be equivalent to using the causality function from vars on a bivariate VAR model. In the following, the F-statistics are indentical, and the p-values are slightly different:
data(ChickEgg)
res1 = grangertest(egg ~ chicken, data = ChickEgg, order = 3)
res2 = causality(VAR(ChickEgg, p = 3), cause = 'chicken')
res1$F[2]
res1[4][2,]
res2$Granger$statistic
res2$Granger$p.value
For multivariate VAR, we can examine each pair of variables separately.
A: Another simple solution is the following:
    NAMES = colnames(df)
    k = ncol(df)
    for (j in 1:k) {
      for (i in 1:k) {
        if (i != j) {
          print(paste(NAMES[j], "->", NAMES[i]))
          VARest = vars::VAR(df[,c(j,i)], p=1)
          print(causality(VARest, cause=NAMES[j]))
        }
      }
    }

A: An workaround to test cause and response within the vars pkg with the help of the "exogen" settings in the VAR function:
granger_bivariate <- function(varest, causal, dep){
  dtmat <- varest$datamat
  mat_target <- dtmat[, c(causal, dep)]
  other_as_exo <- dtmat[, setdiff(names(dtmat),c(causal, dep,'const',names(dtmat)[grepl(paste0('^',causal),names(dtmat)) | grepl(paste0('^',dep),names(dtmat)) ]))]
  var_target <- VAR(mat_target, p = varest$p, exogen = other_as_exo)
  gr_target <- causality(var_target, cause = causal)
  g1 <- gr_target$Granger
  result <- cbind(g1$statistic[1,1], g1$p.value)
  return(result)
}

This works since the causality test ignore the variables in the exogen setting.
