r - 使用 R 对数据框中成对变量进行格兰杰因果关系分析

Question

我想在国家年数据的几个不同变量对之间进行格兰杰因果关系。我似乎能够让它在循环/函数之外工作，但是在将它与我的其余代码集成时遇到了一些麻烦。我在下面提供了一个最小工作示例和所需输出的代码。任何帮助将不胜感激。先感谢您!
编辑 : 我应该在原帖中更清楚地说明。此生成数据中的每一列都包含多个国家/地区的时间序列数据。我想对各个国家/地区进行平均，然后使用下面详述的相同方法对这些变量执行 granger 操作
模拟时间序列数据的代码

library(tidyindexR)
library(tidyverse)
library(dplyr)
library(reshape2)
library(lmtest)

simulateCountryData = function(N=90, NEACH = 30, SEED=100){
        
        # Set seed for reproducibility
        set.seed(SEED)
        
        #""" 
        # This function simulates 
        # N is the total number of observations you want across all subjecs/geocodes
        # NEACH is the number of observations per geocode/sample
        # see https://blog.revolutionanalytics.com/2009/02/how-to-choose-a-random-number-in-r.html for more reources
        #"""
        
        variableOne<-rnorm(N,sample(1:100, NEACH),0.5)
        variableOne[variableOne<0]<-0

        variableTwo<-rnorm(N,sample(100:1, NEACH),0.5)
        variableTwo[variableTwo<0]<-0
        
        variableThree<-rnorm(N,sample(1:100, NEACH),0.5)
        variableThree[variableTwo<0]<-0
        
        variableFour<-rnorm(N,sample(1:100, NEACH),0.5)
        variableFour[variableFour<0]<-0
        
        variableFive<-rnorm(N,sample(1:100, NEACH),0.5)
        variableFive[variableFive<0]<-0
        
        # Defining the subjects/geocodes
        # Makes a sequence of 1:N/NEACH such that the geocode appears NEACH times (i.e. 1,60/3,each=3)
        geocodeNum<-factor(rep(seq(1,N/NEACH),each=NEACH))
        
        # Defining the years/time periods to be repeated 
        # Lists the tear 3x each (3 years per subject/geocode)
        year<-rep(seq(2000,2000+NEACH-1,1),N/NEACH)
        
        # Putting it all together
        AllData<-data.frame(geocodeNum,
                            year,
                            variableOne,
                            variableTwo,
                            variableThree,
                            variableFour,
                            variableFive)
        
        return(AllData)
}

mySimData = simulateCountryData()

在数据帧中提取皮尔逊相关性和格式的代码
信用转至 the-mad-statter谁好心提供了答案 here

# matrix of unique pairs coded as numeric
mx_combos <- combn(1:length(myVariables), 2)
# list of unique pairs coded as numeric
ls_combos <- split(mx_combos, rep(1:ncol(mx_combos), each = nrow(mx_combos)))
# for each pair in the list, create a 1 x 4 dataframe
ls_rows <- lapply(ls_combos, function(p) {
        # lookup names of variables
        v1 <- myVariables[p[1]]
        v2 <- myVariables[p[2]]
        # perform the cor.test()
        htest <- cor.test(mySimData[[v1]], mySimData[[v2]])
        # record pertinent info in a dataframe
        data.frame(Var1 = v1,
                   Var2 = v2,
                   Pval = htest$p.value,
                   Rval = unname(htest$estimate))
})
# row bind the list of dataframes

my_output = dplyr::bind_rows(ls_rows)

此函数生成以下数据帧:

        Var1          Var2        Pval        Rval
variableFive  variableFour 0.749283286 -0.03415477
variableFive   variableOne 0.865119584 -0.01815724
variableFive variableThree 0.004120881 -0.29960240
variableFive   variableTwo 0.024713897  0.23666723
variableFour   variableOne 0.249864859  0.12254729
variableFour variableThree 0.587474758  0.05794634
variableFour   variableTwo 0.624329543  0.05231733
 variableOne variableThree 0.056216554 -0.20200708
 variableOne   variableTwo 0.368598424 -0.09589547

变量三变量二 0.056192121 -0.20202672
我想在 my_output 中的每个成对变量之间添加格兰杰因果关系分析的 p 值。数据框作为附加列。目前，我可以为特定的成对比较提取 pvalues，但我无法弄清楚如何以简单的方式对所有成对变量执行此操作(因为在实际示例中还有更多变量)。如果有人能指出我正确的方向或提供解决方案，那就太棒了。
格兰杰分析代码
此代码似乎对各国每年的数据进行平均，然后在变量之间执行格兰杰因果关系。我只需要一些帮助来弄清楚如何改变它以使用 my_output数据框。

a= mySimData %>% 
        select(geocodeNum, year,variableFive) %>%
        group_by(year) %>%
        summarize(mean(variableFive))

b = mySimData %>% 
        select(geocodeNum, year,variableTwo) %>%
        group_by(year) %>%
        summarize(mean(variableTwo))

c = left_join(a,b)

d = grangertest(c$`mean(variableFive)` ~ c$`mean(variableTwo)`)
d[2,4]

所需输出
理想情况下，输出应如下所示(尽管格兰杰的 P 值正确):

        Var1          Var2        Pval        Rval Granger
variableFive  variableFour 0.749283286 -0.03415477   0.050
variableFive   variableOne 0.865119584 -0.01815724   0.021
variableFive variableThree 0.004120881 -0.29960240   0.090
variableFive   variableTwo 0.024713897  0.23666723   0.042
variableFour   variableOne 0.249864859  0.12254729   0.050
variableFour variableThree 0.587474758  0.05794634   0.021
variableFour   variableTwo 0.624329543  0.05231733   0.092
variableOne variableThree 0.056216554 -0.20200708   0.046
variableOne   variableTwo 0.368598424 -0.09589547   0.072
variableThree   variableTwo 0.056192121 -0.20202672   0.033

再次感谢阅读。

Answer 1

您可以将该函数应用于 combn 中的所有组合。函数本身，所以不需要单独的 lapply称呼。

cols <- grep('variable', names(mySimData), value = TRUE)

result <- do.call(rbind, combn(cols, 2, function(p) {
  v1 <- p[1]
  v2 <- p[2]
  # perform the cor.test()
  htest <- cor.test(mySimData[[v1]], mySimData[[v2]])
  val1 <- aggregate(mySimData[[v1]], list(mySimData$year), mean)[[2]]
  val2 <- aggregate(mySimData[[v2]], list(mySimData$year), mean)[[2]]
  #perform grangertest
  mod <- lmtest::grangertest(val1, val2)
  # record pertinent info in a dataframe
  data.frame(Var1 = v1,
             Var2 = v2,
             Pval = htest$p.value,
             Rval = unname(htest$estimate), 
             Granger = mod[2, 4])
}, simplify = FALSE))

result

#            Var1          Var2     Pval     Rval Granger
#1    variableOne   variableTwo 0.368598 -0.09590  0.8309
#2    variableOne variableThree 0.056217 -0.20201  0.2831
#3    variableOne  variableFour 0.249865  0.12255  0.8860
#4    variableOne  variableFive 0.865120 -0.01816  0.8842
#5    variableTwo variableThree 0.056192 -0.20203  0.9049
#6    variableTwo  variableFour 0.624330  0.05232  0.7993
#7    variableTwo  variableFive 0.024714  0.23667  0.7060
#8  variableThree  variableFour 0.587475  0.05795  0.2510
#9  variableThree  variableFive 0.004121 -0.29960  0.2813
#10  variableFour  variableFive 0.749283 -0.03415  0.8396