libraries = c("nleqslv","dplyr","magrittr","tidyr","rEDM","ggplot2","readxl",
"gridExtra","zoo","scales","humidity")
for(x in libraries) { library(x,character.only=TRUE,warn.conflicts=FALSE,quietly=TRUE) }
### to show the plots as svg-graphics in Jupyter
### although this has some problems when you export the notebook
### to html or other formats
options(jupyter.plot_mimetypes = "image/svg+xml")
packageVersion('rEDM')
set.seed(599213)
nsurr = 500
'%&%' <- function(x,y)paste0(x,y)[1] ‘0.7.2’cntr = "Denmark"
df.long = read.csv('Sugihara2012_dataset.txt',header = TRUE)
df.long %<>%
filter(country == cntr) %>%
filter(year >= 1996)
df.long %>% head| country | variable | year | month | day | value |
|---|---|---|---|---|---|
| Denmark | flu | 1996 | 1 | 1 | NaN |
| Denmark | flu | 1996 | 1 | 8 | NaN |
| Denmark | flu | 1996 | 1 | 15 | NaN |
| Denmark | flu | 1996 | 1 | 22 | NaN |
| Denmark | flu | 1996 | 1 | 29 | NaN |
| Denmark | flu | 1996 | 2 | 5 | NaN |
df.in = df.long %>%
select(-country) %>%
spread(variable,value) %>%
mutate(date = ISOdate(year,month,day)) %>%
select(-year,-month,-day) %>%
select(date,flu,everything())
df.in %>% head| date | flu | AH | PRCP | RH | T |
|---|---|---|---|---|---|
| 1996-01-01 12:00:00 | NaN | 3.314742 | 6.60 | 0.827500 | 27.1 |
| 1996-01-08 12:00:00 | NaN | 4.991437 | 9.54 | 0.954702 | 34.0 |
| 1996-01-15 12:00:00 | NaN | 4.310679 | 8.50 | 0.928307 | 30.9 |
| 1996-01-22 12:00:00 | NaN | 3.044399 | 7.71 | 0.790498 | 26.1 |
| 1996-01-29 12:00:00 | NaN | 3.531248 | 7.20 | 0.878098 | 27.2 |
| 1996-02-05 12:00:00 | NaN | 2.367563 | 5.40 | 0.729602 | 21.8 |
make_block = function(data,cols,delays,lib=c(1,NROW(data))){
lib = matrix(lib,ncol = 2)
data = as.matrix(data)
ncol = length(cols)
nrow = dim(data)[1]
block = array(NA,dim = c(nrow,ncol))
colnames(block) = 1:ncol
for (i in 1:ncol){
I = 1:nrow
I_delay <- intersect(I,I+delays[i])
block[I_delay-delays[i],i] = data[I_delay,cols[i]]
if (delays[i] < 0){
# remove data points that fall at start of lib segments
block[lib[,1] - (0:(delays[i]+1)),i] = NA
colnames(block)[i] = paste(colnames(data)[cols[i]],'_t-',abs(delays[i]),sep="")
} else if (delays[i] > 0) {
# remove data points that fall at end of lib segments
block[lib[,2] - (0:(delays[i]+1)),i] = NA
colnames(block)[i] <- paste(colnames(data)[cols[i]],'_t+',abs(delays[i]),sep="")
} else {
colnames(block)[i] <- paste(colnames(data)[cols[i]],'_t',sep="")
}
}
return(block)
}We define a function yearday_anom() that computes the seasonal cycle and day-of-year anomaly for a time-series. As described in the main text, this is done with smoothing splines.
Anomaly is the difference between the observed value and the smoothed spline.
yearday_anom = function(t,x){
# t: date formatted with POSIXt
# x: time-series values to compute seasonal mean and anomaly
doy = as.numeric(strftime(t, format = "%j"))
I_use = which(!is.na(x))
# create time indices to use for smoothing, replicating data to "wrap around"
doy_sm = rep(doy[I_use],3) + rep(c(-366,0,366),each=length(I_use))
x_sm = rep(x[I_use],3)
xsp = smooth.spline(doy_sm, y = x_sm, w = NULL, spar = 0.8, cv = NA,
all.knots = TRUE,keep.data = TRUE, df.offset = 0)
xbar = data.frame(t=t,doy=doy) %>%
left_join(data.frame(doy=xsp$x,xbar=xsp$y),by='doy') %>%
select(xbar)
out = data.frame(t=t,mean=xbar,anomaly=(x - xbar))
names(out) = c('date','mean','anomaly')
return(out)
}out = yearday_anom(df.in$date,df.in$AH)
AH.bar = out$mean
AH.tilde = out$anomalyoptions(repr.plot.width=7,repr.plot.height=4.5)
out %>%
ggplot(aes(x=date)) +
geom_line(aes(y=mean)) +
geom_point(aes(y=mean+anomaly),color="blue",shape=1) +
labs(x="date",y="Absolute humidity",title=cntr) +
theme_bw() -> p1_Denmark
p1_DenmarkWarning message:
“Removed 4 rows containing missing values (geom_point).”
png
AH.surrs = do.call(cbind,
lapply(1:500, function(i) {
I_na = is.na(AH.tilde)
out = AH.bar
out[I_na] = NA
out[!I_na] = out[!I_na] + sample(AH.tilde[!I_na],sum(!I_na),replace = FALSE)
return(out)
})
)options(repr.plot.width=7,repr.plot.height=5)
plot(df.in$date,df.in$AH,type='l',col='blue',xlab='date',ylab='Absolute Humidity')
lines(df.in$date,AH.surrs[,1],col='red')
legend( x="topright", legend=c("Observed","Seasonal Surrogate"), col=c("blue","red"), lwd=1, lty=c(1,1),bg='white')
png
3. CCM
For set up, we need to determine the set of points to be predicted (pred input) to deal with the fact that many of the influenza time-series have strings of 0 counts. These strings will contain little or no causal information, and more practically can cause problems with how nearest neighbors are calculated for simplex projection when there are many 0 count vectors. We write a function that eliminates all vectors from the prediction that contain a string of at least \(E\) 0’s.
make_pred_nozero = function(time_series,E) {
I_zero_strings = which(time_series==0)
I_zero_strings = Reduce(intersect, lapply((0:E),function(offset) I_zero_strings-offset))
I_zero_strings = c(0,I_zero_strings,NROW(df.in)) # for convenience in next step
N_zero_strings = length(I_zero_strings)
lib_nozeroes = cbind(I_zero_strings[1:(N_zero_strings-1)]+1,I_zero_strings[2:(N_zero_strings)]-1)
lib_out = lib_nozeroes[which(lib_nozeroes[,2] > lib_nozeroes[,1]),]
return(lib_out)
}df.in %>% mutate(index=1:n()) %>% filter(index>=600&index<615)| date | flu | AH | PRCP | RH | T | index |
|---|---|---|---|---|---|---|
| 2007-06-25 12:00:00 | NaN | 10.598504 | 55.70 | 0.867344 | 57.60 | 600 |
| 2007-07-02 12:00:00 | NaN | 11.874325 | 69.20 | 0.898884 | 59.90 | 601 |
| 2007-07-09 12:00:00 | NaN | 10.996224 | 59.64 | 0.846580 | 59.40 | 602 |
| 2007-07-16 12:00:00 | NaN | 12.403229 | 19.48 | 0.835184 | 63.40 | 603 |
| 2007-07-23 12:00:00 | NaN | 11.702477 | 40.23 | 0.847995 | 61.20 | 604 |
| 2007-07-30 12:00:00 | NaN | 10.683401 | 42.80 | 0.779359 | 61.00 | 605 |
| 2007-08-06 12:00:00 | NaN | 14.497789 | 12.50 | 0.862870 | 67.15 | 606 |
| 2007-08-13 12:00:00 | NaN | 11.939355 | 52.16 | 0.806622 | 63.30 | 607 |
| 2007-08-20 12:00:00 | NaN | 14.267609 | 11.34 | 0.915726 | 64.85 | 608 |
| 2007-08-27 12:00:00 | NaN | 9.605718 | 26.29 | 0.818989 | 56.40 | 609 |
| 2007-09-03 12:00:00 | NaN | 9.804182 | 21.30 | 0.813357 | 57.20 | 610 |
| 2007-09-10 12:00:00 | NaN | 9.695090 | 17.28 | 0.835153 | 56.10 | 611 |
| 2007-09-17 12:00:00 | NaN | 10.795756 | 7.14 | 0.914172 | 56.60 | 612 |
| 2007-09-24 12:00:00 | NaN | 10.025608 | 18.80 | 0.922036 | 54.20 | 613 |
| 2007-10-01 12:00:00 | NaN | 8.960859 | 40.93 | 0.899003 | 51.70 | 614 |
I_zero_strings = which(df.in$flu==0)
I_zero_strings- 59
- 61
- 106
- 107
- 109
- 118
- 119
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 147
- 148
- 149
- 150
- 154
- 167
- 168
- 170
- 191
- 192
- 193
- 194
- 195
- 196
- 197
- 198
- 199
- 200
- 201
- 203
- 204
- 205
- 216
- 218
- 219
- 220
- 221
- 222
- 248
- 249
- 250
- 251
- 252
- 253
- 254
- 255
- 256
- 257
- 258
- 260
- 262
- 275
- 276
- 277
- 278
- 279
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 310
- 311
- 312
- 313
- 314
- 319
- 326
- 327
- 328
- 329
- 330
- 331
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 363
- 364
- 365
- 366
- 367
- 379
- 380
- 381
- 382
- 383
- 384
- 404
- 405
- 406
- 407
- 408
- 409
- 410
- 411
- 422
- 423
- 426
- 428
- 430
- 431
- 432
- 433
- 434
- 435
- 436
- 457
- 458
- 459
- 460
- 461
- 462
- 463
- 464
- 465
- 466
- 467
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 519
- 520
- 521
- 522
- 523
- 525
- 538
- 539
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 571
- 572
- 573
- 574
- 575
- 578
- 591
- 592
- 593
- 594
- 615
- 616
- 617
- 618
- 619
- 620
- 622
- 625
- 643
- 644
- 665
- 666
- 667
- 668
- 669
- 671
- 740
- 751
- 752
- 753
- 754
- 755
- 756
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 771
- 772
- 773
- 774
- 775
- 776
- 802
- 804
- 805
- 806
- 808
- 809
- 810
- 811
- 812
- 814
- 815
- 816
- 817
- 820
- 823
- 824
- 825
- 826
- 827
- 828
- 830
- 832
- 833
- 834
- 835
- 836
- 837
- 838
- 839
- 855
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 864
- 865
- 866
- 867
- 868
- 869
- 870
- 872
- 873
- 876
- 908
- 909
- 911
- 913
- 916
- 918
- 919
- 920
- 921
- 922
- 923
- 924
- 925
- 926
I_zero_strings2 = Reduce(intersect, lapply((0:5),function(offset) I_zero_strings-offset))
I_zero_strings2- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 191
- 192
- 193
- 194
- 195
- 196
- 248
- 249
- 250
- 251
- 252
- 253
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 326
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 379
- 404
- 405
- 406
- 430
- 431
- 457
- 458
- 459
- 460
- 461
- 462
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 615
- 751
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 771
- 823
- 832
- 833
- 834
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 864
- 865
- 918
- 919
- 920
- 921
lapply((0:8),function(offset) I_zero_strings-offset)-
- 59
- 61
- 106
- 107
- 109
- 118
- 119
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 147
- 148
- 149
- 150
- 154
- 167
- 168
- 170
- 191
- 192
- 193
- 194
- 195
- 196
- 197
- 198
- 199
- 200
- 201
- 203
- 204
- 205
- 216
- 218
- 219
- 220
- 221
- 222
- 248
- 249
- 250
- 251
- 252
- 253
- 254
- 255
- 256
- 257
- 258
- 260
- 262
- 275
- 276
- 277
- 278
- 279
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 310
- 311
- 312
- 313
- 314
- 319
- 326
- 327
- 328
- 329
- 330
- 331
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 363
- 364
- 365
- 366
- 367
- 379
- 380
- 381
- 382
- 383
- 384
- 404
- 405
- 406
- 407
- 408
- 409
- 410
- 411
- 422
- 423
- 426
- 428
- 430
- 431
- 432
- 433
- 434
- 435
- 436
- 457
- 458
- 459
- 460
- 461
- 462
- 463
- 464
- 465
- 466
- 467
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 519
- 520
- 521
- 522
- 523
- 525
- 538
- 539
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 571
- 572
- 573
- 574
- 575
- 578
- 591
- 592
- 593
- 594
- 615
- 616
- 617
- 618
- 619
- 620
- 622
- 625
- 643
- 644
- 665
- 666
- 667
- 668
- 669
- 671
- 740
- 751
- 752
- 753
- 754
- 755
- 756
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 771
- 772
- 773
- 774
- 775
- 776
- 802
- 804
- 805
- 806
- 808
- 809
- 810
- 811
- 812
- 814
- 815
- 816
- 817
- 820
- 823
- 824
- 825
- 826
- 827
- 828
- 830
- 832
- 833
- 834
- 835
- 836
- 837
- 838
- 839
- 855
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 864
- 865
- 866
- 867
- 868
- 869
- 870
- 872
- 873
- 876
- 908
- 909
- 911
- 913
- 916
- 918
- 919
- 920
- 921
- 922
- 923
- 924
- 925
- 926
-
- 58
- 60
- 105
- 106
- 108
- 117
- 118
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 147
- 148
- 149
- 153
- 166
- 167
- 169
- 190
- 191
- 192
- 193
- 194
- 195
- 196
- 197
- 198
- 199
- 200
- 202
- 203
- 204
- 215
- 217
- 218
- 219
- 220
- 221
- 247
- 248
- 249
- 250
- 251
- 252
- 253
- 254
- 255
- 256
- 257
- 259
- 261
- 274
- 275
- 276
- 277
- 278
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 310
- 311
- 312
- 313
- 318
- 325
- 326
- 327
- 328
- 329
- 330
- 351
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 363
- 364
- 365
- 366
- 378
- 379
- 380
- 381
- 382
- 383
- 403
- 404
- 405
- 406
- 407
- 408
- 409
- 410
- 421
- 422
- 425
- 427
- 429
- 430
- 431
- 432
- 433
- 434
- 435
- 456
- 457
- 458
- 459
- 460
- 461
- 462
- 463
- 464
- 465
- 466
- 509
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 519
- 520
- 521
- 522
- 524
- 537
- 538
- 561
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 571
- 572
- 573
- 574
- 577
- 590
- 591
- 592
- 593
- 614
- 615
- 616
- 617
- 618
- 619
- 621
- 624
- 642
- 643
- 664
- 665
- 666
- 667
- 668
- 670
- 739
- 750
- 751
- 752
- 753
- 754
- 755
- 757
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 771
- 772
- 773
- 774
- 775
- 801
- 803
- 804
- 805
- 807
- 808
- 809
- 810
- 811
- 813
- 814
- 815
- 816
- 819
- 822
- 823
- 824
- 825
- 826
- 827
- 829
- 831
- 832
- 833
- 834
- 835
- 836
- 837
- 838
- 854
- 856
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 864
- 865
- 866
- 867
- 868
- 869
- 871
- 872
- 875
- 907
- 908
- 910
- 912
- 915
- 917
- 918
- 919
- 920
- 921
- 922
- 923
- 924
- 925
-
- 57
- 59
- 104
- 105
- 107
- 116
- 117
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 147
- 148
- 152
- 165
- 166
- 168
- 189
- 190
- 191
- 192
- 193
- 194
- 195
- 196
- 197
- 198
- 199
- 201
- 202
- 203
- 214
- 216
- 217
- 218
- 219
- 220
- 246
- 247
- 248
- 249
- 250
- 251
- 252
- 253
- 254
- 255
- 256
- 258
- 260
- 273
- 274
- 275
- 276
- 277
- 298
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 310
- 311
- 312
- 317
- 324
- 325
- 326
- 327
- 328
- 329
- 350
- 351
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 363
- 364
- 365
- 377
- 378
- 379
- 380
- 381
- 382
- 402
- 403
- 404
- 405
- 406
- 407
- 408
- 409
- 420
- 421
- 424
- 426
- 428
- 429
- 430
- 431
- 432
- 433
- 434
- 455
- 456
- 457
- 458
- 459
- 460
- 461
- 462
- 463
- 464
- 465
- 508
- 509
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 519
- 520
- 521
- 523
- 536
- 537
- 560
- 561
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 571
- 572
- 573
- 576
- 589
- 590
- 591
- 592
- 613
- 614
- 615
- 616
- 617
- 618
- 620
- 623
- 641
- 642
- 663
- 664
- 665
- 666
- 667
- 669
- 738
- 749
- 750
- 751
- 752
- 753
- 754
- 756
- 757
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 771
- 772
- 773
- 774
- 800
- 802
- 803
- 804
- 806
- 807
- 808
- 809
- 810
- 812
- 813
- 814
- 815
- 818
- 821
- 822
- 823
- 824
- 825
- 826
- 828
- 830
- 831
- 832
- 833
- 834
- 835
- 836
- 837
- 853
- 855
- 856
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 864
- 865
- 866
- 867
- 868
- 870
- 871
- 874
- 906
- 907
- 909
- 911
- 914
- 916
- 917
- 918
- 919
- 920
- 921
- 922
- 923
- 924
-
- 56
- 58
- 103
- 104
- 106
- 115
- 116
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 147
- 151
- 164
- 165
- 167
- 188
- 189
- 190
- 191
- 192
- 193
- 194
- 195
- 196
- 197
- 198
- 200
- 201
- 202
- 213
- 215
- 216
- 217
- 218
- 219
- 245
- 246
- 247
- 248
- 249
- 250
- 251
- 252
- 253
- 254
- 255
- 257
- 259
- 272
- 273
- 274
- 275
- 276
- 297
- 298
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 310
- 311
- 316
- 323
- 324
- 325
- 326
- 327
- 328
- 349
- 350
- 351
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 363
- 364
- 376
- 377
- 378
- 379
- 380
- 381
- 401
- 402
- 403
- 404
- 405
- 406
- 407
- 408
- 419
- 420
- 423
- 425
- 427
- 428
- 429
- 430
- 431
- 432
- 433
- 454
- 455
- 456
- 457
- 458
- 459
- 460
- 461
- 462
- 463
- 464
- 507
- 508
- 509
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 519
- 520
- 522
- 535
- 536
- 559
- 560
- 561
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 571
- 572
- 575
- 588
- 589
- 590
- 591
- 612
- 613
- 614
- 615
- 616
- 617
- 619
- 622
- 640
- 641
- 662
- 663
- 664
- 665
- 666
- 668
- 737
- 748
- 749
- 750
- 751
- 752
- 753
- 755
- 756
- 757
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 771
- 772
- 773
- 799
- 801
- 802
- 803
- 805
- 806
- 807
- 808
- 809
- 811
- 812
- 813
- 814
- 817
- 820
- 821
- 822
- 823
- 824
- 825
- 827
- 829
- 830
- 831
- 832
- 833
- 834
- 835
- 836
- 852
- 854
- 855
- 856
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 864
- 865
- 866
- 867
- 869
- 870
- 873
- 905
- 906
- 908
- 910
- 913
- 915
- 916
- 917
- 918
- 919
- 920
- 921
- 922
- 923
-
- 55
- 57
- 102
- 103
- 105
- 114
- 115
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 150
- 163
- 164
- 166
- 187
- 188
- 189
- 190
- 191
- 192
- 193
- 194
- 195
- 196
- 197
- 199
- 200
- 201
- 212
- 214
- 215
- 216
- 217
- 218
- 244
- 245
- 246
- 247
- 248
- 249
- 250
- 251
- 252
- 253
- 254
- 256
- 258
- 271
- 272
- 273
- 274
- 275
- 296
- 297
- 298
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 310
- 315
- 322
- 323
- 324
- 325
- 326
- 327
- 348
- 349
- 350
- 351
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 363
- 375
- 376
- 377
- 378
- 379
- 380
- 400
- 401
- 402
- 403
- 404
- 405
- 406
- 407
- 418
- 419
- 422
- 424
- 426
- 427
- 428
- 429
- 430
- 431
- 432
- 453
- 454
- 455
- 456
- 457
- 458
- 459
- 460
- 461
- 462
- 463
- 506
- 507
- 508
- 509
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 519
- 521
- 534
- 535
- 558
- 559
- 560
- 561
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 571
- 574
- 587
- 588
- 589
- 590
- 611
- 612
- 613
- 614
- 615
- 616
- 618
- 621
- 639
- 640
- 661
- 662
- 663
- 664
- 665
- 667
- 736
- 747
- 748
- 749
- 750
- 751
- 752
- 754
- 755
- 756
- 757
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 771
- 772
- 798
- 800
- 801
- 802
- 804
- 805
- 806
- 807
- 808
- 810
- 811
- 812
- 813
- 816
- 819
- 820
- 821
- 822
- 823
- 824
- 826
- 828
- 829
- 830
- 831
- 832
- 833
- 834
- 835
- 851
- 853
- 854
- 855
- 856
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 864
- 865
- 866
- 868
- 869
- 872
- 904
- 905
- 907
- 909
- 912
- 914
- 915
- 916
- 917
- 918
- 919
- 920
- 921
- 922
-
- 54
- 56
- 101
- 102
- 104
- 113
- 114
- 134
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 149
- 162
- 163
- 165
- 186
- 187
- 188
- 189
- 190
- 191
- 192
- 193
- 194
- 195
- 196
- 198
- 199
- 200
- 211
- 213
- 214
- 215
- 216
- 217
- 243
- 244
- 245
- 246
- 247
- 248
- 249
- 250
- 251
- 252
- 253
- 255
- 257
- 270
- 271
- 272
- 273
- 274
- 295
- 296
- 297
- 298
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 314
- 321
- 322
- 323
- 324
- 325
- 326
- 347
- 348
- 349
- 350
- 351
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 374
- 375
- 376
- 377
- 378
- 379
- 399
- 400
- 401
- 402
- 403
- 404
- 405
- 406
- 417
- 418
- 421
- 423
- 425
- 426
- 427
- 428
- 429
- 430
- 431
- 452
- 453
- 454
- 455
- 456
- 457
- 458
- 459
- 460
- 461
- 462
- 505
- 506
- 507
- 508
- 509
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 520
- 533
- 534
- 557
- 558
- 559
- 560
- 561
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 573
- 586
- 587
- 588
- 589
- 610
- 611
- 612
- 613
- 614
- 615
- 617
- 620
- 638
- 639
- 660
- 661
- 662
- 663
- 664
- 666
- 735
- 746
- 747
- 748
- 749
- 750
- 751
- 753
- 754
- 755
- 756
- 757
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 771
- 797
- 799
- 800
- 801
- 803
- 804
- 805
- 806
- 807
- 809
- 810
- 811
- 812
- 815
- 818
- 819
- 820
- 821
- 822
- 823
- 825
- 827
- 828
- 829
- 830
- 831
- 832
- 833
- 834
- 850
- 852
- 853
- 854
- 855
- 856
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 864
- 865
- 867
- 868
- 871
- 903
- 904
- 906
- 908
- 911
- 913
- 914
- 915
- 916
- 917
- 918
- 919
- 920
- 921
-
- 53
- 55
- 100
- 101
- 103
- 112
- 113
- 133
- 134
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 148
- 161
- 162
- 164
- 185
- 186
- 187
- 188
- 189
- 190
- 191
- 192
- 193
- 194
- 195
- 197
- 198
- 199
- 210
- 212
- 213
- 214
- 215
- 216
- 242
- 243
- 244
- 245
- 246
- 247
- 248
- 249
- 250
- 251
- 252
- 254
- 256
- 269
- 270
- 271
- 272
- 273
- 294
- 295
- 296
- 297
- 298
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 313
- 320
- 321
- 322
- 323
- 324
- 325
- 346
- 347
- 348
- 349
- 350
- 351
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 373
- 374
- 375
- 376
- 377
- 378
- 398
- 399
- 400
- 401
- 402
- 403
- 404
- 405
- 416
- 417
- 420
- 422
- 424
- 425
- 426
- 427
- 428
- 429
- 430
- 451
- 452
- 453
- 454
- 455
- 456
- 457
- 458
- 459
- 460
- 461
- 504
- 505
- 506
- 507
- 508
- 509
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 519
- 532
- 533
- 556
- 557
- 558
- 559
- 560
- 561
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 572
- 585
- 586
- 587
- 588
- 609
- 610
- 611
- 612
- 613
- 614
- 616
- 619
- 637
- 638
- 659
- 660
- 661
- 662
- 663
- 665
- 734
- 745
- 746
- 747
- 748
- 749
- 750
- 752
- 753
- 754
- 755
- 756
- 757
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 796
- 798
- 799
- 800
- 802
- 803
- 804
- 805
- 806
- 808
- 809
- 810
- 811
- 814
- 817
- 818
- 819
- 820
- 821
- 822
- 824
- 826
- 827
- 828
- 829
- 830
- 831
- 832
- 833
- 849
- 851
- 852
- 853
- 854
- 855
- 856
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 864
- 866
- 867
- 870
- 902
- 903
- 905
- 907
- 910
- 912
- 913
- 914
- 915
- 916
- 917
- 918
- 919
- 920
-
- 52
- 54
- 99
- 100
- 102
- 111
- 112
- 132
- 133
- 134
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 147
- 160
- 161
- 163
- 184
- 185
- 186
- 187
- 188
- 189
- 190
- 191
- 192
- 193
- 194
- 196
- 197
- 198
- 209
- 211
- 212
- 213
- 214
- 215
- 241
- 242
- 243
- 244
- 245
- 246
- 247
- 248
- 249
- 250
- 251
- 253
- 255
- 268
- 269
- 270
- 271
- 272
- 293
- 294
- 295
- 296
- 297
- 298
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 312
- 319
- 320
- 321
- 322
- 323
- 324
- 345
- 346
- 347
- 348
- 349
- 350
- 351
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 372
- 373
- 374
- 375
- 376
- 377
- 397
- 398
- 399
- 400
- 401
- 402
- 403
- 404
- 415
- 416
- 419
- 421
- 423
- 424
- 425
- 426
- 427
- 428
- 429
- 450
- 451
- 452
- 453
- 454
- 455
- 456
- 457
- 458
- 459
- 460
- 503
- 504
- 505
- 506
- 507
- 508
- 509
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 518
- 531
- 532
- 555
- 556
- 557
- 558
- 559
- 560
- 561
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 571
- 584
- 585
- 586
- 587
- 608
- 609
- 610
- 611
- 612
- 613
- 615
- 618
- 636
- 637
- 658
- 659
- 660
- 661
- 662
- 664
- 733
- 744
- 745
- 746
- 747
- 748
- 749
- 751
- 752
- 753
- 754
- 755
- 756
- 757
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 795
- 797
- 798
- 799
- 801
- 802
- 803
- 804
- 805
- 807
- 808
- 809
- 810
- 813
- 816
- 817
- 818
- 819
- 820
- 821
- 823
- 825
- 826
- 827
- 828
- 829
- 830
- 831
- 832
- 848
- 850
- 851
- 852
- 853
- 854
- 855
- 856
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 865
- 866
- 869
- 901
- 902
- 904
- 906
- 909
- 911
- 912
- 913
- 914
- 915
- 916
- 917
- 918
- 919
-
- 51
- 53
- 98
- 99
- 101
- 110
- 111
- 131
- 132
- 133
- 134
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 146
- 159
- 160
- 162
- 183
- 184
- 185
- 186
- 187
- 188
- 189
- 190
- 191
- 192
- 193
- 195
- 196
- 197
- 208
- 210
- 211
- 212
- 213
- 214
- 240
- 241
- 242
- 243
- 244
- 245
- 246
- 247
- 248
- 249
- 250
- 252
- 254
- 267
- 268
- 269
- 270
- 271
- 292
- 293
- 294
- 295
- 296
- 297
- 298
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 311
- 318
- 319
- 320
- 321
- 322
- 323
- 344
- 345
- 346
- 347
- 348
- 349
- 350
- 351
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 371
- 372
- 373
- 374
- 375
- 376
- 396
- 397
- 398
- 399
- 400
- 401
- 402
- 403
- 414
- 415
- 418
- 420
- 422
- 423
- 424
- 425
- 426
- 427
- 428
- 449
- 450
- 451
- 452
- 453
- 454
- 455
- 456
- 457
- 458
- 459
- 502
- 503
- 504
- 505
- 506
- 507
- 508
- 509
- 510
- 511
- 512
- 513
- 514
- 515
- 517
- 530
- 531
- 554
- 555
- 556
- 557
- 558
- 559
- 560
- 561
- 562
- 563
- 564
- 565
- 566
- 567
- 570
- 583
- 584
- 585
- 586
- 607
- 608
- 609
- 610
- 611
- 612
- 614
- 617
- 635
- 636
- 657
- 658
- 659
- 660
- 661
- 663
- 732
- 743
- 744
- 745
- 746
- 747
- 748
- 750
- 751
- 752
- 753
- 754
- 755
- 756
- 757
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 794
- 796
- 797
- 798
- 800
- 801
- 802
- 803
- 804
- 806
- 807
- 808
- 809
- 812
- 815
- 816
- 817
- 818
- 819
- 820
- 822
- 824
- 825
- 826
- 827
- 828
- 829
- 830
- 831
- 847
- 849
- 850
- 851
- 852
- 853
- 854
- 855
- 856
- 857
- 858
- 859
- 860
- 861
- 862
- 864
- 865
- 868
- 900
- 901
- 903
- 905
- 908
- 910
- 911
- 912
- 913
- 914
- 915
- 916
- 917
- 918
Reduce is from the functional programming module of R ## Iterative function application: Funcall = function(f, …) f(…) ## Compute log(exp(acos(cos(0)) Reduce(Funcall, list(log, exp, acos, cos), 0, right = TRUE)
I_zero_strings = I_zero_strings2
I_zero_strings = c(0,I_zero_strings,NROW(df.in)) # for convenience in next step
I_zero_strings- 0
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 191
- 192
- 193
- 194
- 195
- 196
- 248
- 249
- 250
- 251
- 252
- 253
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 326
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 379
- 404
- 405
- 406
- 430
- 431
- 457
- 458
- 459
- 460
- 461
- 462
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 615
- 751
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 771
- 823
- 832
- 833
- 834
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 864
- 865
- 918
- 919
- 920
- 921
- 970
N_zero_strings = length(I_zero_strings)
N_zero_strings106
lib_nozeroes = cbind(I_zero_strings[1:(N_zero_strings-1)]+1,I_zero_strings[2:(N_zero_strings)]-1)
lib_nozeroes %>% head(16)| 1 | 138 |
| 140 | 139 |
| 141 | 140 |
| 142 | 141 |
| 143 | 142 |
| 144 | 143 |
| 145 | 144 |
| 146 | 190 |
| 192 | 191 |
| 193 | 192 |
| 194 | 193 |
| 195 | 194 |
| 196 | 195 |
| 197 | 247 |
| 249 | 248 |
| 250 | 249 |
lib_out = lib_nozeroes[which(lib_nozeroes[,2] > lib_nozeroes[,1]),]
lib_out| 1 | 138 |
| 146 | 190 |
| 197 | 247 |
| 254 | 299 |
| 310 | 325 |
| 327 | 351 |
| 363 | 378 |
| 380 | 403 |
| 407 | 429 |
| 432 | 456 |
| 463 | 509 |
| 519 | 561 |
| 571 | 614 |
| 616 | 750 |
| 752 | 757 |
| 772 | 822 |
| 824 | 831 |
| 835 | 856 |
| 866 | 917 |
| 922 | 969 |
3.a Determine optimal \(E\) for CCM
We now determine the proper embedding dimension for cross-mapping following the procedure of Deyle et al. (2016). We determine \(E^*\) that gives optimal prediction of cross-mapping lagged 1 week. Note that we specificy the lib_size to be the largest possible number of points in the library (even though the actual number will be smaller because of NAs and the zero strings). This is for convenience so that the code will automatically use the largest possible lib_size.
block_temp = df.in %>% select(flu,AH)
lib_ccm = c(1,NROW(df.in))
Emax=8
out.temp = do.call(
rbind,
lapply(1:Emax, function(E_i) {
pred_ccm = make_pred_nozero(block_temp$flu,E_i)
ccm(block=block_temp,
E=E_i,
lib=lib_ccm,
pred=pred_ccm,
lib_sizes = NROW(block_temp),
exclusion_radius=0,
random_libs = FALSE,
num_sample=1,
tp = -1,
lib_column = 1,
target_column = 2,
silent = TRUE) }
)
)
out.temp| E | tau | tp | num_neighbors | lib_column | target_column | lib_size | num_pred | rho | mae | rmse |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 1 | -1 | 2 | 1 | 2 | 614 | 400 | 0.17910471 | 1.843047 | 2.501677 |
| 2 | 1 | -1 | 3 | 1 | 2 | 594 | 400 | 0.05135658 | 1.951606 | 2.597812 |
| 3 | 1 | -1 | 4 | 1 | 2 | 577 | 399 | 0.22749008 | 1.751233 | 2.342287 |
| 4 | 1 | -1 | 5 | 1 | 2 | 559 | 399 | 0.04313842 | 1.883103 | 2.497428 |
| 5 | 1 | -1 | 6 | 1 | 2 | 542 | 406 | 0.06881939 | 2.056252 | 2.737585 |
| 6 | 1 | -1 | 7 | 1 | 2 | 526 | 416 | 0.29636337 | 1.700103 | 2.285878 |
| 7 | 1 | -1 | 8 | 1 | 2 | 510 | 415 | 0.34598041 | 1.645424 | 2.229613 |
| 8 | 1 | -1 | 9 | 1 | 2 | 494 | 415 | 0.32404071 | 1.673133 | 2.290292 |
E_star = out.temp$E[which.max(out.temp$rho)]
E_star7
3.b Calculate CCM for real data
We then measure the un-lagged cross-map skill with this value of \(E^*\).
pred_ccm = make_pred_nozero(block_temp$flu,E_star)
df.out.ccm = ccm(block=block_temp,
E=E_star,
lib=lib_ccm,
pred = pred_ccm,
lib_sizes = NROW(block_temp),
exclusion_radius=0,
random_libs = FALSE,
num_sample=1,
tp = 0,
silent = TRUE)
df.out.ccm| E | tau | tp | num_neighbors | lib_column | target_column | lib_size | num_pred | rho | mae | rmse |
|---|---|---|---|---|---|---|---|---|---|---|
| 7 | 1 | 0 | 8 | 1 | 2 | 510 | 415 | 0.3955595 | 1.560958 | 2.186334 |
df.out.ccm.surr = data.frame()
for (i_surr in 1:dim(AH.surrs)[2]) {
block_temp = data.frame(flu=df.in$flu,AH=AH.surrs[,i_surr])
out.temp = do.call(
rbind,
lapply(1:8, function(E_i){
pred_ccm = make_pred_nozero(block_temp$flu,E_i)
ccm(block=block_temp,
E=E_i,
lib=lib_ccm,
pred=pred_ccm,
lib_sizes = NROW(block_temp),
exclusion_radius=0,
random_libs = FALSE,
num_sample=1,
tp = -1,
lib_column = 1,
target_column = 2,
silent=TRUE)
}
)
)
E_star = out.temp$E[which.max(out.temp$rho)]
pred_ccm = make_pred_nozero(block_temp$flu,E_star)
out.temp = ccm(block = block_temp,
E = E_star,
lib = lib_ccm,
pred = pred_ccm,
lib_sizes = NROW(block_temp),
exclusion_radius = 0,
random_libs = FALSE,
num_sample=1,
tp = 0,
silent=TRUE)
df.out.ccm.surr <- df.out.ccm.surr %>% bind_rows(out.temp)
} # for i_surroptions(repr.plot.width=3.5,repr.plot.height=5)
boxplot(df.out.ccm.surr$rho, ylab='rho (CCM)', xlab=cntr, main='Flu cross-map AH')
points(1,df.out.ccm$rho, pch=8, col='red', cex=2, lwd=2)
png
df.out.ccm| E | tau | tp | num_neighbors | lib_column | target_column | lib_size | num_pred | rho | mae | rmse |
|---|---|---|---|---|---|---|---|---|---|---|
| 7 | 1 | 0 | 8 | 1 | 2 | 510 | 415 | 0.3955595 | 1.560958 | 2.186334 |