Effect size for r lmtest granger causality?

Question

Is it possible to calculate the effect size from the output of grangertest in package lmtest? I could calculate it by hand, but it only gives the F value, and not the sums of squares.

Granger causality test

Model 1: apwbc ~ Lags(apwbc, 1:1) + Lags(other, 1:1)
Model 2: apwbc ~ Lags(apwbc, 1:1)
  Res.Df Df      F  Pr(>F)  
1    163                    
2    164 -1 4.8495 0.02906 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Answer 1

We say a variable X Granger causes a variable Y if predictions of Y based on lagged values of both X and Y are significantly better than predictions of Y based on lagged values of Y alone. This means that a Granger causality test is a test of nested models; the model with lags of both variables is the full model, and the model with lags of Y alone is the nested or restricted model. We can calculate partial η² from a test of nested models as (SS_restricted - SS_full) / SS_restricted (see for example Wright and London (2009), p. 21).

This is pretty easy to code up:

nested_model_partial_eta_sq <- function(full_mod, rest_mod) {
    # full_mod is the full model (an lm object)
    # rest_mod is the model that omits one or more explanatory variables
    SS_full <- sum(full_mod$residuals^2)
    SS_rest <- sum(rest_mod$residuals^2)
    return((SS_rest - SS_full) / SS_rest)
}

However, we might want a function that generates the regressions for us also, as well as runs the Granger causality test:

custom_granger <- function(x, y, p = 1) {
    # Does x Granger cause y? What is the effect size?
    # We want to test Granger causality,
    # but then also get partial eta squared to measure effect size
    # First it will be convenient to store the variable names
    varnames <- c(deparse(substitute(x)), deparse(substitute(y)))
    lagnames <- paste0(varnames, "_lag", rep(1:p, each = 2))
    # Then we created the lagged variables / data for models
    VAR_data <- embed(as.matrix(cbind(x, y)), p + 1)
    colnames(VAR_data) <- c(varnames, lagnames)
    VAR_data <- VAR_data[, -1]
    # Run the full model
    model_formula <- paste(colnames(VAR_data)[-1], collapse = " + ")
    model_formula <- formula(paste0(varnames[2], " ~ ", model_formula))
    full_mod <- lm(model_formula, data = as.data.frame(VAR_data))
    # Take out the lags of x and run the nested model
    VAR_data <- VAR_data[ , seq(from = 1, to = p * 2 + 1, by = 2)]
    model_formula <- paste(colnames(VAR_data)[-1], collapse = " + ")
    model_formula <- formula(paste0(varnames[2], " ~ ", model_formula))
    rest_mod <- lm(model_formula, data = as.data.frame(VAR_data))
    # Then we can do the Granger test
    granger_test <- anova(full_mod, rest_mod)
    # and get partial eta squared
    SS_full <- granger_test$RSS[1]
    SS_rest <- granger_test$RSS[2]
    partial_eta_squared <- (SS_rest - SS_full) / SS_rest
    # And return all of it
    return(list(VAR_result = full_mod,
                granger_test = granger_test,
                partial_eta_squared = partial_eta_squared))
}

With your data this results in

df <- read.csv('cormanaz-data.txt')
with(df, custom_granger(other, apwbc))

$VAR_result

Call:
lm(formula = model_formula, data = as.data.frame(VAR_data))

Coefficients:
(Intercept)   other_lag1   apwbc_lag1  
    1.97732     -0.01671      0.48997  


$granger_test
Analysis of Variance Table

Model 1: apwbc ~ other_lag1 + apwbc_lag1
Model 2: apwbc ~ apwbc_lag1
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1    163 6267.9                              
2    164 6454.4 -1   -186.48 4.8495 0.02906 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$partial_eta_squared
[1] 0.02889191

You can see the F statistic for the Granger causality test is the same as from lmtest::grangertest(), but this will also give you the VAR coefficients and the partial η², which for your particular case (I guess the number you've been after all along) is about 0.029.