-
Notifications
You must be signed in to change notification settings - Fork 0
/
model2.R
109 lines (95 loc) · 4.32 KB
/
model2.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
library(dplyr)
library(ggplot2)
library(tidyr)
library(rstan)
d <- readRDS("data/t2.rds")
d <- d %>% mutate(itime = 12*year+imonth) %>% mutate(itime = as.integer(1+itime-min(itime)))
if (F) {
m1 <- lm(t ~ name*month + month/decade + I(lat-62):decade + I(lon-29):decade, data=d)
summary(m1)
}
itimes <- 1:max(d$itime)
names <- unique(d$name)
d.complete <- expand.grid(name=names, itime=itimes) %>%
left_join(d, by=c("name", "itime"))
temp <- matrix(d.complete$t, nrow=length(itimes), ncol=length(names), byrow=T)
if (F)
image(temp)
time.vars <- data.frame(itime=itimes) %>% left_join(d %>%
select(itime, year, month, imonth, decade) %>% distinct(), by="itime")
station.vars <- data.frame(name=names) %>% left_join(d %>%
select(name, sid, sub, lat, lon, subsid) %>% distinct(), by="name")
darks <- which(is.na(temp), T)
temp[is.na(temp)] <- 10000
stan.data <- list(
S=length(names), T=length(itimes),
temp=temp, lat=station.vars$lat-60, lon=station.vars$lon-25,
month=time.vars$imonth,
darkN = nrow(darks))
m <- stan_model(file="model2.stan")
s <- sampling(m, data=stan.data, warmup=250, iter=1000, thin=5, init=0, chains=4, refresh=1)
saveRDS(s, "s.rds")
s <- readRDS("s.rds")
plot(s)
plot(s, par="trend")
plot(s, par="rho")
plot(s, par="tau_month")
traceplot(s, "trend", inc_warmup=F, ask=T)
trend.samples <- extract(s, "trend")[[1]]
plot(density(apply(trend.samples, 1, mean)))
hist(apply(trend.samples, 1, mean), n=100)
Omega <- matrix(apply(apply(extract(s, "LOmega")[[1]], 1, function (m) m %*% t(m)), 1, mean), rep(length(names), 2))
rownames(Omega) <- names
colnames(Omega) <- names
heatmap(Omega, symm=T)
traceplot(s, "tau_month", inc_warmup=F, ask=F)
traceplot(s, "trend_lat", inc_warmup=F)
#traceplot(s, "LOmega", inc_warmup=F, ask=T)
sort(setNames(svd(Omega)[[2]][,1], names), decr=T)
sort(setNames(svd(Omega)[[2]][,2], names), decr=T)
trend.samples %>% apply(., 2, function (v) quantile(v, c(.025, .25, .5, .75, .975))/10) %>%
t %>% data.frame %>% setNames(c("ll", "l", "m", "u", "uu")) %>%
data.frame(., month=reorder(levels(d$month), 1:12)) %>%
ggplot(., aes(x=month, y=m)) + #geom_point(size=2) +
geom_pointrange(aes(ymax=u, ymin=l), size=1) +
geom_pointrange(aes(ymax=uu, ymin=ll), size=.5) +
ylab("°C / decade") + xlab("") + ggtitle("1980-2014 (bars 50% and 95%)") +
expand_limits(y=0) +
theme_bw(20)
ggsave("figs/monthly.png", scale=.7)
trend.samples %>% apply(., 1, mean) %>% data.frame(trend=./10) %>%
ggplot(., aes(x=trend)) + geom_density(fill="grey80", color="white") + theme_bw(20) +
xlab("°C / decade") + ylab("") + scale_y_continuous(breaks=NULL) +
geom_vline(aes(xintercept=mean(trend.samples)/10), color="red") +
geom_vline(aes(xintercept=quantile(apply(trend.samples, 1,
function (x) mean(x)/10), c(.025, .25, .75, .975))), color="red", linetype=2) +
ggtitle("1980-2014 (lines 50% and 95%)")
ggsave("figs/trend.png", scale=.7)
ggplot(d.complete, aes(x=itime, y=name, fill=t)) + geom_tile() + ggtitle("Raw data")
ggsave("figs/data.png", scale=.7)
ggplot(reshape2::melt(Omega), aes(x=Var1, y=Var2, fill=value)) + geom_raster() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + xlab("") + ylab("") +
ggtitle("Station correlations")
ggsave("figs/statcorr.png", scale=.7)
extract(s, "tau_month")[[1]] %>%
apply(., 2, function (v) quantile(v, c(.025, .25, .5, .75, .975))) %>%
t %>% data.frame %>% setNames(c("ll", "l", "m", "u", "uu")) %>%
data.frame(., month=reorder(levels(d$month), 1:12)) %>%
ggplot(., aes(x=month, y=m)) + #geom_point(size=2) +
geom_pointrange(aes(ymax=u, ymin=l), size=1) +
geom_pointrange(aes(ymax=uu, ymin=ll), size=.5) +
ylab("") + xlab("") + ggtitle("Monthly sd (relative)") +
expand_limits(y=0) +
theme_bw()
ggsave("figs/tau_month.png", scale=.7)
extract(s, "tau")[[1]] %>%
apply(., 2, function (v) quantile(v, c(.025, .25, .5, .75, .975))) %>%
t %>% data.frame %>% setNames(c("ll", "l", "m", "u", "uu")) %>%
data.frame(., month=names) %>%
ggplot(., aes(x=month, y=m)) + #geom_point(size=2) +
geom_pointrange(aes(ymax=u, ymin=l), size=1) +
geom_pointrange(aes(ymax=uu, ymin=ll), size=.5) +
ylab("") + xlab("") + ggtitle("Station sd (relative)") +
expand_limits(y=0) + coord_flip() +
theme_bw()
ggsave("figs/tau.png", scale=.7)