-
Notifications
You must be signed in to change notification settings - Fork 0
/
13_Calculating_p-tail.R
91 lines (65 loc) · 2.3 KB
/
13_Calculating_p-tail.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
library(tidyverse)
options(scipen=999)
#### ODBA ####
# Set up matrices and data frames
mintemp <- matrix(NA, ncol=35, nrow=120) # 35 birds x 120 days longest migration
mintemp <- as.data.frame(mintemp)
prcp <- matrix(NA, ncol=35, nrow=120)
prcp <- as.data.frame(prcp)
# Load file list
files <- list.files("E:/ResearchProjects/GeeseBehavior-Weather/results/dlmODBA", pattern=".Rdata", all.files=TRUE, full.names=TRUE)
# Loop over each file, calculate number of posterior samples >0
for (i in 1:length(files)) {
# Load model output
load(files[i])
# Extract bird ID
x <- str_split(files[i], "[:punct:]")
x <- unlist(x[[1]][11])
# save sims list
beta1 <- out$sims.list$beta1
beta2 <- out$sims.list$beta2
for (j in 1:dim(beta1)[2]) {
prcp[j,i] <- ifelse(unique(is.na(beta1[,j])),NA,sum(beta1[,j]>0)/out$mcmc.info$n.samples)
names(prcp)[i] <- x
mintemp[j,i] <- ifelse(unique(is.na(beta2[,j])),NA,sum(beta2[,j]>0)/out$mcmc.info$n.samples)
names(mintemp)[i] <- x
}
print(x)
print(dim(beta1)[2])
print(dim(beta2)[2])
}
# Save to file for plotting
write_csv(prcp, "results/ODBAprcp_ptail.csv")
write_csv(mintemp, "results/ODBAmintemp_ptail.csv")
#### PTF ####
# Set up matrices and data frames
mintemp <- matrix(NA, ncol=35, nrow=120)
mintemp <- as.data.frame(mintemp)
prcp <- matrix(NA, ncol=35, nrow=120)
prcp <- as.data.frame(prcp)
# Load file list
files <- list.files("E:/ResearchProjects/GeeseBehavior-Weather/results/dlmPTF",
pattern=".Rdata", all.files=TRUE, full.names=TRUE)
# Loop over each file, calculate number of posterior samples >0
for (i in 1:length(files)) {
# Load model output
load(files[i])
# Extract bird ID
x <- str_split(files[i], "[:punct:]")
x <- unlist(x[[1]][11])
# save sims list
beta1 <- out$sims.list$beta1
beta2 <- out$sims.list$beta2
for (j in 1:dim(beta1)[2]) {
prcp[j,i] <- ifelse(unique(is.na(beta1[,j])),NA,sum(beta1[,j]>0)/out$mcmc.info$n.samples)
names(prcp)[i] <- x
mintemp[j,i] <- ifelse(unique(is.na(beta2[,j])),NA,sum(beta2[,j]>0)/out$mcmc.info$n.samples)
names(mintemp)[i] <- x
}
print(x)
print(dim(beta1)[2])
print(dim(beta2)[2])
}
# Save to file for plotting
write_csv(prcp, "results/PTFprcp_ptail.csv")
write_csv(mintemp, "results/PTFmintemp_ptail.csv")