library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(forecast)
library(lubridate)
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
##
## date
Everyone
load("track_all.rdata")
track_all = track_all %>%
select(name = MovDataID, x=X, y=Y, timestamp=Fixtime) %>%
mutate_if(is.factor, as.character) %>%
mutate(timestamp = ymd_hms(timestamp)) %>%
group_by(name) %>%
arrange(timestamp) %>%
distinct()
track_all %>%
count(name)
## # A tibble: 57 x 2
## # Groups: name [57]
## name n
## <chr> <int>
## 1 Abu 17864
## 2 Amelia 11217
## 3 Angelia 6735
## 4 Annabelle_ww 10778
## 5 Beti 10352
## 6 Boniface 16582
## 7 BraBrou 17811
## 8 David 10301
## 9 Dedi 10199
## 10 Doudou 10483
## # … with 47 more rows
Interesting
interesting_elephants = c("Amelia", "Mboumba", "Annabelle_ww", "Beti", "Nzamba")
track_interesting = track_all %>%
filter(name %in% interesting_elephants)
all_lags = purrr::map_dfr(
1:500,
function(i) {
track_interesting %>%
transmute(
time_diff = difftime(timestamp, lag(timestamp, i), units="hours"),
x_diff = x - lag(x, i),
y_diff = y - lag(y, i),
lag = i
) %>%
na.omit()
}
)
all_lags %>% count(name)
## # A tibble: 5 x 2
## # Groups: name [5]
## name n
## <chr> <int>
## 1 Amelia 5483250
## 2 Annabelle_ww 5263750
## 3 Beti 5050750
## 4 Mboumba 5168250
## 5 Nzamba 4929250
all_lags = all_lags %>%
mutate(time_diff = as.numeric(time_diff)) %>%
filter(time_diff <= 120) %>%
arrange(name, time_diff) %>%
mutate(bin = ceiling(time_diff))
all_lags %>% count(name)
## # A tibble: 5 x 2
## # Groups: name [5]
## name n
## <chr> <int>
## 1 Amelia 1615517
## 2 Annabelle_ww 1206441
## 3 Beti 1106102
## 4 Mboumba 1179699
## 5 Nzamba 1102792
all_binned_params = all_lags %>%
group_by(name, bin) %>%
summarize(
var_x = var(x_diff),
var_y = var(y_diff),
cor = cor(x_diff, y_diff)
)
all_binned_params %>%
tidyr::gather(var, value, -bin, -name) %>%
ggplot(aes(x=bin, y=value, color=name)) +
#geom_smooth(method = "lm", formula=y~1+x+I(x^2), alpha=0.2) +
geom_line() +
facet_wrap(~var, scale="free_y")
Naieve Model
models = all_binned_params %>%
rename(t = bin) %>%
group_by(name) %>%
tidyr::nest() %>%
transmute(
name = name,
model = purrr::map(
data,
function(d) {
list(
var_x = lm(var_x ~ 1 + t + I(t^2), data=d),
var_y = lm(var_y ~ 1 + t + I(t^2), data=d),
cor = lm(cor ~ 1 + t + I(t^2), data=d)
)
}
)
) %>%
tidyr::unnest() %>%
mutate(
param = purrr::map_chr(model, ~ .x$terms[[2]] %>% as.character()),
coefs = purrr::map(model, ~ broom::tidy(.x, quick=TRUE) %>% tidyr::spread(term, estimate))
) %>%
tidyr::unnest(coefs) %>%
rename(
intercept = `(Intercept)`,
b_t = t,
b_t2 = `I(t^2)`
)
models
## # A tibble: 15 x 6
## name model param intercept b_t2 b_t
## <chr> <list> <chr> <dbl> <dbl> <dbl>
## 1 Amelia <S3: lm> var_x -0.0000551 -9.33e-9 0.0000109
## 2 Amelia <S3: lm> var_y -0.000105 -5.64e-8 0.0000175
## 3 Amelia <S3: lm> cor -0.00269 5.02e-5 -0.0103
## 4 Annabelle_ww <S3: lm> var_x 0.0000282 -5.45e-9 0.000000885
## 5 Annabelle_ww <S3: lm> var_y 0.0000677 -3.24e-8 0.00000544
## 6 Annabelle_ww <S3: lm> cor -0.0103 2.79e-5 -0.00331
## 7 Beti <S3: lm> var_x 0.0000327 -2.34e-8 0.00000520
## 8 Beti <S3: lm> var_y -0.000209 8.78e-8 0.0000389
## 9 Beti <S3: lm> cor 0.115 5.03e-7 -0.00456
## 10 Mboumba <S3: lm> var_x -0.0000267 -3.17e-8 0.0000311
## 11 Mboumba <S3: lm> var_y -0.0000678 -5.75e-8 0.0000347
## 12 Mboumba <S3: lm> cor 0.0656 6.39e-6 0.00124
## 13 Nzamba <S3: lm> var_x -0.0000531 -2.00e-8 0.0000165
## 14 Nzamba <S3: lm> var_y -0.0000351 -6.08e-8 0.0000133
## 15 Nzamba <S3: lm> cor 0.115 -2.39e-5 0.00327
Seed Shadow
library(brms, )
d = readr::read_csv("Gut passage data.csv") %>%
filter(beadsFound == "Y")
## Warning: Missing column names filled in: 'X1' [1]
n = 10000
nthin = 5
nburn = 10000
fit_gamma = brm(
bf(GPT~1), data = d, family = Gamma(link="log"),
chains = 1, cores = 1, iter = nburn+nthin*n, thin = nthin
)
##
## SAMPLING FOR MODEL 'dafe6b2fe1e53b13a456a79b2df05d7f' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 3.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 60000 [ 0%] (Warmup)
## Chain 1: Iteration: 6000 / 60000 [ 10%] (Warmup)
## Chain 1: Iteration: 12000 / 60000 [ 20%] (Warmup)
## Chain 1: Iteration: 18000 / 60000 [ 30%] (Warmup)
## Chain 1: Iteration: 24000 / 60000 [ 40%] (Warmup)
## Chain 1: Iteration: 30000 / 60000 [ 50%] (Warmup)
## Chain 1: Iteration: 30001 / 60000 [ 50%] (Sampling)
## Chain 1: Iteration: 36000 / 60000 [ 60%] (Sampling)
## Chain 1: Iteration: 42000 / 60000 [ 70%] (Sampling)
## Chain 1: Iteration: 48000 / 60000 [ 80%] (Sampling)
## Chain 1: Iteration: 54000 / 60000 [ 90%] (Sampling)
## Chain 1: Iteration: 60000 / 60000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.45281 seconds (Warm-up)
## Chain 1: 1.77175 seconds (Sampling)
## Chain 1: 3.22456 seconds (Total)
## Chain 1:
chain = fit_gamma$fit@sim$samples[[1]]
gamma_sim = data_frame(
mu = exp(chain$b_Intercept),
shape = chain$shape
) %>%
slice((nburn/nthin+1):(n+nburn/nthin)) %>%
mutate(
scale = mu / shape,
t = rgamma(n(), shape=shape, scale=scale)
)
sims = models %>%
group_by(name) %>%
summarize(
models = list(model %>% setNames(param))
) %>%
mutate(
shadow = purrr::map(
models,
function(m) {
data_frame(
t = gamma_sim$t,
var_x = predict(m$var_x, newdata = gamma_sim),
var_y = predict(m$var_y, newdata = gamma_sim),
cor = predict(m$cor, newdata = gamma_sim)
)
}
)
) %>%
tidyr::unnest(shadow) %>%
mutate(
x = rnorm(n(), 0, sqrt(var_x)),
y = rnorm(n(), sqrt(var_y/var_x)*cor*x, sqrt( (1-cor^2)*var_y ))
)
## Warning in sqrt(var_x): NaNs produced
## Warning in rnorm(n(), 0, sqrt(var_x)): NAs produced
## Warning in sqrt(var_y/var_x): NaNs produced
## Warning in sqrt((1 - cor^2) * var_y): NaNs produced
## Warning in rnorm(n(), sqrt(var_y/var_x) * cor * x, sqrt((1 - cor^2) *
## var_y)): NAs produced
sims %>% filter(is.nan(x) | is.nan(y))
## # A tibble: 7 x 7
## name t var_x var_y cor x y
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Amelia 5.94 0.00000922 -0.00000250 -0.0619 -0.000336 NaN
## 2 Amelia 1.91 -0.0000343 -0.0000712 -0.0222 NaN NaN
## 3 Amelia 4.51 -0.00000615 -0.0000266 -0.0480 NaN NaN
## 4 Beti 1.91 0.0000426 -0.000135 0.106 -0.00200 NaN
## 5 Beti 4.51 0.0000557 -0.0000323 0.0941 -0.000942 NaN
## 6 Mboumba 1.91 0.0000326 -0.00000173 0.0680 -0.00492 NaN
## 7 Nzamba 1.91 -0.0000216 -0.00000994 0.121 NaN NaN
sims %>%
ggplot(aes(x=x,y=y)) +
stat_density_2d(geom = "polygon", aes(fill = stat(nlevel))) +
scale_fill_viridis_c() +
facet_wrap(~name)
## Warning: Removed 7 rows containing non-finite values (stat_density2d).
