【发布时间】：2021-04-12 20:58:53
【问题描述】：
我想在国家年份数据的几对不同变量之间进行格兰杰因果关系。我似乎能够让它在循环/函数之外工作，但是在将它与我的其余代码集成时遇到了一些麻烦。我在下面提供了一个最小工作示例和所需输出的代码。任何帮助将不胜感激。提前谢谢你！
编辑：我应该在原帖中更清楚地说明。此生成数据中的每一列都包含多个国家/地区的时间序列数据。我想在各个国家/地区进行平均，然后使用下面详述的相同方法对这些变量执行 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()
在数据框中提取 pearson 相关性和格式的代码 感谢 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
再次感谢您的阅读。
【问题讨论】：
标签： r dataframe dplyr time-series correlation