-
Notifications
You must be signed in to change notification settings - Fork 0
/
partialCoh.R
133 lines (98 loc) · 3.78 KB
/
partialCoh.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
## partial coherence
## only significant coherence within a station (for now) but set up framework just in case
## same examples for now
#D22 D10 do
#but not D4
#lag 6
test1=lm(D22$residGAMdo~D4$residGAMdo)
partial=resid(test1)
station1Data=D10
station2Data=D4
varName="residGAMdo"
varName2="residGAMdo"
empP<-c()
for(i in 2:nrow(D22) ){
test1=partial[c((i):length(partial),1:(i-1))]
test=spectrum( cbind(station1Data[,varName],test1),taper=.2,log="no",spans=c(16,16),demean=T,detrend=F,plot=F)
probCoh=df(test$coh,2,2*test$df-2)
f=qf(0.95,2,2*test$df-2)
C = f/(16-1+f) ## what about a double pass filter?
#test=ccf(station1Data[,varName],test1,lag.max=12,plot=F)
#ccfOfInterest=test$acf[14:25]
empP<-c( empP,probCoh[which(probCoh>C)][which.max(probCoh[which(probCoh>C)])])
# print(i)
}
test=spectrum( cbind(station1Data[,varName],partial),taper=.2,log="no",spans=c(16,16),demean=T,detrend=F,plot=F)
probCoh=df(test$coh,2,2*test$df-2)
f=qf(0.95,2,2*test$df-2)
C = f/(16-1+f)
maxPhase=test$phase[which(probCoh>C)][which.max(probCoh[which(probCoh>C)])]
maxCoh=test$coh[which.max(probCoh[which(probCoh>C)])]
maxFreq=1/test$freq[which(probCoh>C)][which.max(probCoh[which(probCoh>C)])]
pVal=length(which(empP>probCoh[which(probCoh>C)][which.max(probCoh[which(probCoh>C)])] ))/nrow(station1Data)
minPval=1/nrow(station2Data)
maxPhase ## 1.96693
maxCoh ## 0.0006343933
maxFreq ## 24
pVal ## 0.4677419
minPval ## 0.005376344
#D26 pheo
#D4 sal
#but not D22
#lag 3
test1=lm(D26$residGAMpheoTransform~D22$residGAMsal)
test2=lm(D26$residGAMpheoTransform~D22$residGAMpheo)
partial=resid(test1)
partial2=resid(test2)
station1Data=D4
varName="residGAMsal"
empP<-c()
for(i in 2:nrow(D26) ){
test1=partial[c((i):length(partial),1:(i-1))]
test=spectrum( cbind(station1Data[,varName],test1),taper=.2,log="no",spans=c(16,16),demean=T,detrend=F,plot=F)
#test=ccf(station1Data[,varName],test1,lag.max=12,plot=F)
#ccfOfInterest=test$acf[14:25]
probCoh=df(test$coh,2,2*test$df-2)
f=qf(0.95,2,2*test$df-2)
C = f/(16-1+f)
empP<-c( empP,probCoh[which(probCoh>C)][which.max(probCoh[which(probCoh>C)])] )
# print(i)
}
test=spectrum( cbind(station1Data[,varName],partial),taper=.2,log="no",spans=c(16,16),demean=T,detrend=F,plot=F)
probCoh=df(test$coh,2,2*test$df-2)
f=qf(0.95,2,2*test$df-2)
C = f/(16-1+f)
maxPhase=test$phase[which(probCoh>C)][which.max(probCoh[which(probCoh>C)])]
maxCoh=test$coh[which.max(probCoh[which(probCoh>C)])]
maxFreq=1/test$freq[which(probCoh>C)][which.max(probCoh[which(probCoh>C)])]
pVal=length(which(empP>probCoh[which(probCoh>C)][which.max(probCoh[which(probCoh>C)])] ))/nrow(station1Data)
minPval=1/nrow(station2Data)
maxPhase ## 9.625096e-16
maxCoh ## 0.01996144
maxFreq ## 2
pVal ## 0.9731183
minPval ## 0.005376344
empP<-c()
for(i in 2:nrow(D26) ){
test1=partial2[c((i):length(partial2),1:(i-1))]
test=spectrum( cbind(station1Data[,varName],test1),taper=.2,log="no",spans=c(16,16),demean=T,detrend=F,plot=F)
probCoh=df(test$coh,2,2*test$df-2)
f=qf(0.95,2,2*test$df-2)
C = f/(16-1+f) ## what about a double pass filter?
empP<-c( empP,probCoh[which(probCoh>C)][which.max(probCoh[which(probCoh>C)])] )
# print(i)
}
test=spectrum( cbind(station1Data[,varName],partial2),taper=.2,log="no",spans=c(16,16),demean=T,detrend=F,plot=F)
probCoh=df(test$coh,2,2*test$df-2)
f=qf(0.95,2,2*test$df-2)
C = f/(16-1+f)
maxPhase=test$phase[which(probCoh>C)][which.max(probCoh[which(probCoh>C)])]
maxCoh=test$coh[which.max(probCoh[which(probCoh>C)])]
maxFreq=1/test$freq[which(probCoh>C)][which.max(probCoh[which(probCoh>C)])]
pVal=length(which(empP>probCoh[which(probCoh>C)][which.max(probCoh[which(probCoh>C)])] ))/nrow(station1Data)
minPval=1/nrow(station2Data)
maxPhase ##-1.183674
maxCoh ## 0.0008370846
maxFreq ## 8.347826
pVal ## 0.4623656
minPval ## 0.005376344