library(CTMCdive)
# Simulate data -----------------------------------------------------------
# total observation time
T <- 24 * 60 * 7
# exposure time
exp_T <- T / 2
# time step
dt <- 0.01
# time-varying intensities
tgr <- seq(0, T, by = dt)
dive_I <- function(t) {
eff <- 0.06 * (sin(2 * pi * t / T - pi / 2)) + 6 * 0.05
return(eff)
}
surf_I <- function(t, exp = TRUE) {
eff <- 0.02 * (sin(2 * pi * t / T) + 1.2) + 0.06
a <- exp_T - 24 * 60
b <- exp_T + 12 * 60
if(exp) eff <- ifelse(t > exp_T, eff - 0.03 * (t - a)/(exp_T - a) * (t - b)/(exp_T - b) * (t > a) * (t < b), eff)
return(eff)
}
# kappa
kappa <- list(dive = 3, surf = 3)
# plot truth
plot(tgr, dive_I(tgr), type = "l", lwd = 1.5, xlab = "Time", ylab = "Dive Intensity")
lines(tgr, dive_I(tgr), col = "steelblue")
plot(tgr, surf_I(tgr), type = "l", lwd = 1.5, xlab = "Time", ylab = "Surface Intensity")
abline(v = c(exp_T, exp_T + 24 * 60), col = "firebrick", lty = "dashed")
# simulate data
seed <- sample(1:65555, size = 1)
set.seed(seed)
dat <- simulateCTMC2(dive_I, surf_I, T, dt, kappa = kappa)
# add exposure data
window_size <- 24 * 60
dat$expt <- ifelse(dat$time >= exp_T & dat$time < exp_T + window_size, dat$time - exp_T, 0)
dat$exp <- ifelse(dat$time >= exp_T & dat$time < exp_T + window_size, 1, 0)
dat$exp <- factor(dat$exp, ordered = TRUE)
# plot data
plot(dat$time, dat$dive, pch = 19, xlab = "Time of Dive Start", ylab = "Dive Duration")
abline(v = c(exp_T, exp_T + window_size), col = "firebrick", lty = "dashed")
plot(dat$time, dat$surf, pch = 19, xlab = "Time of Dive Start", ylab = "Surface Duration")
abline(v = c(exp_T, exp_T + window_size), col = "firebrick", lty = "dashed")
# Fit Model ---------------------------------------------------------------
# setup model
forms <- list(surface ~ s(time, bs = "ts"),
dive ~ s(time, bs = "ts"))
# fit model
m0 <- FitCTMCdive(forms, dat, dt = 1, print = TRUE)
# fit exposure effect
mexp <- update(m0, ~.+ exp + s(expt, by = exp, bs = "cs", k = 10, m = 1), print = TRUE)
# select model
mod <- mexp$both
# see results
summary(mod)
exp(mod$res$surface[,1])
exp(mod$res$dive[,1])
# plot fitted model
plot(mod)
# get predicted values
pred <- predict(mod)
# plot residuals
rdive <- pred$rdive
plot(dat$time, rdive, pch = 19)
hist(rdive)
qqnorm(rdive); qqline(rdive)
rsurf <- pred$rsurf
plot(dat$time, rsurf, pch = 19)
hist(rsurf)
qqnorm(rsurf); qqline(rsurf)
# estimated exposure effect
expeff <- GetExposureEff(mod, exp_var = "exp")
plotExposureEffect(expeff, pick = "dive")
abline(v = c(exp_T, exp_T + window_size), lty = "dotted", col = "firebrick")
# plot exposure against baseline prediction
limits <- c(exp_T - 24*60, exp_T + 3*24*60)
plot(dat$time, dat$dive, pch = 19, xlab = "Time of Dive Start", ylab = "Surface Duration", xlim = limits)
abline(v = c(exp_T, exp_T + 24 * 60), col = "firebrick", lty = "dashed")
lines(dat$time[dat$time < exp_T], pred$dive[dat$time < exp_T], col = "steelblue", lwd = 1.5)
lines(dat$time[dat$time > exp_T + 24 * 60], pred$dive[dat$time > exp_T + 24 * 60], col = "steelblue", lwd = 1.5)
matlines(expeff$time, expeff$dive$pred[,1:20], col = "firebrick", lwd = 1.5, alpha = 0.7)
matlines(expeff$time, expeff$dive$base[,1:20], col = "steelblue")
# plot estimated intensities against truth
plot(tgr, dive_I(tgr), type = "l", lwd = 1.5, xlab = "Time", ylab = "Dive Intensity")
lines(mod$sm$ints, pred$diveI, lwd = 1.5, col = "firebrick")
plot(tgr, surf_I(tgr), type = "l", lwd = 1.5, xlab = "Time", ylab = "Surface Intensity")
lines(mod$sm$ints, pred$surfI, lwd = 1.5, col = "firebrick")
lines(tgr, surf_I(tgr, exp = FALSE), col = "steelblue")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.