source("~/Dropbox/R_package_testing/brmsish/R/html_summary.R")
clim_dat_2020 <- read_excel("./data/datos_metereologicos/Estacion_met/Dat_met_estacion_2022-01-11.xlsx",
    sheet = "2020")

clim_dat_2019 <- read_excel("./data/datos_metereologicos/Estacion_met/Dat_met_estacion_2022-01-11.xlsx",
    sheet = "2019")

clim_dat <- rbind(clim_dat_2019, clim_dat_2020)

clim_dat <- clim_dat[, c("filename", "Año", "Mes", "Día", "Hora", "Temp (°C)",
    "Humedad Relat.", "Precipitación")]

names(clim_dat) <- c("filename", "year", "month", "day", "hour", "temp", "HR", "rain")

clim_dat <- aggregate(cbind(rain, temp, HR) ~ filename + year + month + day + hour,
    clim_dat, mean)

clim_dat$year <- clim_dat$year + 2000

clim_dat$date <- as.Date(paste(clim_dat$year, clim_dat$month, clim_dat$day, sep = "-"))


clim_dat$date_hour <- paste(gsub("-", "", clim_dat$date), clim_dat$hour, sep = "-")

call_dat <- read.csv("./data/processed/call_rate_per_date_time_and_site.csv")

# remove 5 pm data and keep only form Sukia
call_dat <- call_dat[call_dat$hour != 17, ]
call_dat <- call_dat[call_dat$site != "LAGCHIMU", ]


call_dat$date_hour <- paste(sapply(as.character(call_dat$site_date_hour), function(x) strsplit(x,
    "_")[[1]][2]), call_dat$hour, sep = "-")


call_dat$temp <- sapply(1:nrow(call_dat), function(x) {
    y <- clim_dat$temp[clim_dat$date_hour == call_dat$date_hour[x]]

    if (length(y) < 1)
        y <- NA

    return(y)
})


call_dat$HR <- sapply(1:nrow(call_dat), function(x) {
    y <- clim_dat$HR[clim_dat$date_hour == call_dat$date_hour[x]]

    if (length(y) < 1)
        y <- NA

    return(y)
})


call_dat$rain <- sapply(1:nrow(call_dat), function(x) {
    y <- clim_dat$rain[clim_dat$date_hour == call_dat$date_hour[x]]

    if (length(y) < 1)
        y <- NA

    return(y)
})

# proportion of acoustic data with climatic data
sum(call_dat$date_hour %in% clim_dat$date_hour)/nrow(call_dat)
sum(!is.na(call_dat$temp))/nrow(call_dat)


call_dat <- call_dat[!is.na(call_dat$temp), ]

call_dat$day <- as.numeric(substr(sapply(as.character(call_dat$site_date_hour), function(x) strsplit(x,
    "_")[[1]][2]), 7, 8))

call_dat$date <- as.Date(paste(call_dat$year, call_dat$month, call_dat$day, sep = "-"))

call_dat$moon.date <- ifelse(call_dat$hour < 12, as.Date(call_dat$date - 1), as.Date(call_dat$date))

call_dat$moon.date <- as.Date(call_dat$moon.date, origin = "1970-01-02")

## add moon
call_dat$moonlight <- lunar.illumination(call_dat$moon.date, shift = -6)


call_dat$date_hour_min <- strptime(paste(paste(call_dat$year, call_dat$month, call_dat$day,
    sep = "-"), paste(call_dat$hour, "00", sep = ":")), format = "%Y-%m-%d  %H:%M")

call_dat$hour_diff <- as.numeric(call_dat$date_hour_min - min(call_dat$date_hour_min))/3600


call_dat$rain_24 <- sapply(1:nrow(call_dat), function(x) sum(clim_dat$rain[strptime(clim_dat$date,
    format = "%Y-%m-%d") == (strptime(call_dat$date[x], format = "%Y-%m-%d") - 60 *
    60 * 24)]))

call_dat$rain_48 <- sapply(1:nrow(call_dat), function(x) sum(clim_dat$rain[strptime(clim_dat$date,
    format = "%Y-%m-%d") == (strptime(call_dat$date[x], format = "%Y-%m-%d") - 60 *
    60 * 48)]))

clim_dat$date_hour_min <- strptime(paste(paste(clim_dat$year, clim_dat$month, clim_dat$day,
    sep = "-"), paste(clim_dat$hour, "00", sep = ":")), format = "%Y-%m-%d %H:%M")

clim_dat$hour_diff <- as.numeric(clim_dat$date_hour_min - min(call_dat$date_hour_min))/3600

call_dat$prev_temp <- sapply(1:nrow(call_dat), function(x) {

    # if(call_dat$hour_diff[x] < 48) pt <- NA else
    pt <- mean(clim_dat$temp[clim_dat$hour_diff %in% (call_dat$hour_diff[x] - 48):(call_dat$hour_diff[x] -
        24)])

    return(pt)
})

write.csv(call_dat, "./data/processed/acoustic_and_climatic_data_by_hour.csv")

 

Purpose

  • Evaluate effect of enviromental factors on vocal activity of A. lemur

 

Prepare data

Descriptive stats

call_dat_site <- read.csv("./data/processed/call_rate_per_date_time_and_site.csv")
  • Total number of recordings: 9681
  • Total recordings per site:
agg_recs <- aggregate(rec_time ~ site, data = call_dat_site, length)
names(agg_recs)[1:2] <- c("site", "rec_count")

# print table as kable
kb <- kable(agg_recs, row.names = TRUE, digits = 3)

kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

print(kb)
site rec_count
1 LAGCHIMU 5127
2 SUKIA 4554

  • Total recording time: 3224

  • Total recording time per site:

agg_recs <- aggregate(rec_time ~ site, data = call_dat_site, sum)
names(agg_recs)[1:2] <- c("site", "recording_time")

agg_recs$recording_time <- round((agg_recs$recording_time)/60)

# print table as kable
kb <- kable(agg_recs, row.names = TRUE, digits = 3)

kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

print(kb)
site recording_time
1 LAGCHIMU 1707
2 SUKIA 1517

  • Total detections: 540203

  • Total detections per site:

agg_recs <- aggregate(n_call ~ site, data = call_dat_site, sum)
names(agg_recs)[1:2] <- c("calls", "count")

# print table as kable
kb <- kable(agg_recs, row.names = TRUE, digits = 3)

kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

print(kb)
calls count
1 LAGCHIMU 169485
2 SUKIA 370718

  • Call rate: 18.598587

  • Call rate per site:

agg_recs <- aggregate(call_rate ~ site, data = call_dat_site, mean)
agg_recs$sd <- aggregate(call_rate ~ site, data = call_dat_site, sd)[, 2]

names(agg_recs)[1:3] <- c("site", "call_rate", "sd")

# print table as kable
kb <- kable(agg_recs, row.names = TRUE, digits = 3)

kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

print(kb)
site call_rate sd
1 LAGCHIMU 11.018 18.039
2 SUKIA 27.133 87.778 
call_rate_hour <- read.csv("./data/processed/acoustic_and_climatic_data_by_hour.csv")


agg <- aggregate(cbind(temp, prev_temp, HR, rain, rain_24, rain_48, moonlight) ~
    1, call_rate_hour, function(x) round(c(mean(x), sd(x), min(x), max(x)), 3))

agg <- as.data.frame(matrix(unlist(agg), ncol = 4, byrow = TRUE, dimnames = list(c("Temperature",
    "Previous temperature", "Relative humidity", "Night rain", "Rain 24 hours", "Rain 48 hours",
    "Moonlight"), c("mean", "sd", "min", "max"))))
# print table as kable
kb <- kable(agg, row.names = TRUE, digits = 3)

kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

print(kb)
mean sd min max
Temperature 24.229 2.126 19.472 31.763
Previous temperature 24.097 0.988 20.005 27.111
Relative humidity 90.056 8.464 55.158 99.940
Night rain 0.079 0.457 0.000 10.417
Rain 24 hours 1.441 3.161 0.000 25.620
Rain 48 hours 1.453 3.170 0.000 25.620
Moonlight 0.495 0.355 0.000 1.000

Mean and sd temperature: 24.229179

Mean previous temperature: 24.097073

Mean temperature: 24.229179

Mean cumulative rain per hour: 0.079232

Mean cumulative rain per hour previous 24 hours: 1.44064

Mean daily cumulative rain per hour previous 48 hours: 1.45317

Bayesian regression models

Scale variables and set model parameters

call_rate_hour <- read.csv("./data/processed/acoustic_and_climatic_data_by_hour.csv")

# make hour a factor
call_rate_hour$hour <- factor(call_rate_hour$hour)

# scale and mean-center
call_rate_hour$sc_temp <- scale(call_rate_hour$temp)
call_rate_hour$sc_HR <- scale(call_rate_hour$HR)
call_rate_hour$sc_rain <- scale(call_rate_hour$rain)
call_rate_hour$sc_rain_24 <- scale(call_rate_hour$rain_24)
call_rate_hour$sc_rain_48 <- scale(call_rate_hour$rain_48)
call_rate_hour$sc_moonlight <- scale(call_rate_hour$moonlight)
call_rate_hour$sc_prev_temp <- scale(call_rate_hour$prev_temp)

priors <- c(prior(normal(0, 4), class = "b"))
chains <- 4
iter <- 10000

Global models

The increase in relative humidity, decrease in temperature, increase in the previous accumulated rain, decrease in the night rain, decrease in the percentage of the moon illuminated cause the activity of A. lemur to increase?

fit_glob1 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_moonlight +
    sc_rain + sc_rain_24 + sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour),
    data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(),
    prior = priors, file = "./data/processed/regression_models/global_rain_24", file_refit = "always")


fit_glob2 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_moonlight +
    sc_rain + sc_rain_48 + sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour),
    data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(),
    prior = priors, file = "./data/processed/regression_models/global_rain_48", file_refit = "always")
html_summary(read.file = "./data/processed/regression_models/global_rain_24.rds",
    n.posterior = 1000)
global_rain_24
n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_moonlight + sc_rain + sc_rain_24 + sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 0 0 9377.5 12503.8 536473866
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.698 1 9377.5 12503.8 0.636 0.760
sc_temp 0.573 1 10930.3 12951.7 0.493 0.654
sc_HR 0.307 1 14333.7 14034.8 0.233 0.381
sc_moonlight -0.090 1 12403.9 13102.6 -0.149 -0.031
sc_rain 0.038 1 18321.1 15080.0 0.002 0.075
sc_rain_24 0.091 1 18798.2 14951.0 0.053 0.130
sc_prev_temp -0.024 1 14546.5 13679.8 -0.069 0.022 

html_summary(read.file = "./data/processed/regression_models/global_rain_48.rds",
    n.posterior = 1000)
global_rain_48
n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_moonlight + sc_rain + sc_rain_48 + sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 0 0 12375.1 13473.4 58053599
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.703 1 12375.1 13699.9 0.640 0.766
sc_temp 0.566 1 13593.6 13473.4 0.485 0.649
sc_HR 0.303 1 17757.9 15101.0 0.228 0.379
sc_moonlight -0.098 1 18857.2 14680.5 -0.158 -0.039
sc_rain 0.031 1 22865.9 15592.2 -0.006 0.068
sc_rain_48 -0.006 1 24427.4 16403.0 -0.044 0.033
sc_prev_temp -0.037 1 20592.6 15806.9 -0.084 0.010

Temperature

The increase in temperature at night causes that activity of A. lemur to decrease?

fit2.1 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2,
    time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
    cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/temp2.1",
    file_refit = "always")

fit2.2 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2,
    time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
    cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/temp2.2",
    file_refit = "always")

fit2.3 <- brm(n_call | resp_rate(rec_time) ~ sc_prev_temp + ar(p = 2, time = hour_diff,
    gr = hour), data = call_rate_hour, iter = iter, chains = chains, cores = chains,
    family = negbinomial(), prior = priors, file = "./data/processed/regression_models/temp2.3",
    file_refit = "always")
html_summary(read.file = "./data/processed/regression_models/temp2.1.rds", n.posterior = 1000)
temp2.1
n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 1 0 10328.4 12111.8 1846130948
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.698 1 10328.4 12111.8 0.634 0.762
sc_temp 0.312 1 12331.5 14350.8 0.261 0.361
sc_rain 0.044 1 20658.2 15981.9 0.008 0.083
sc_rain_24 0.090 1 20542.9 16470.9 0.052 0.128 

html_summary(read.file = "./data/processed/regression_models/temp2.2.rds", n.posterior = 1000)
temp2.2
n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 0 0 8873.8 11982.2 2046380001
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.701 1 8873.8 11982.2 0.635 0.766
sc_temp 0.305 1 10391.0 13564.9 0.255 0.355
sc_rain 0.036 1 17117.4 15242.3 0.000 0.073
sc_rain_48 0.006 1 17732.8 15884.1 -0.031 0.044 

html_summary(read.file = "./data/processed/regression_models/temp2.3.rds", n.posterior = 1000)
temp2.3
n_call | resp_rate(rec_time) ~ sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 0 0 20200 15620.2 1962664000
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.675 1 20200.0 15620.2 0.605 0.745
sc_prev_temp -0.002 1 36045.7 16766.2 -0.048 0.044

Relative humidity

Does an increase in relative humidity cause the activity of A. lemur to increase?

fit.3.1 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_24 +
    ar(p = 2, time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
    cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/RH3.1",
    file_refit = "always")

fit.3.2 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_48 +
    ar(p = 2, time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
    cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/RH3.2",
    file_refit = "always")  # Este tengo que doblemente confirmarlo
html_summary(read.file = "./data/processed/regression_models/RH3.1.rds", n.posterior = 1000)
RH3.1
n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 0 0 10000 12606.5 1854758574
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.699 1 10000.0 12606.5 0.637 0.760
sc_temp 0.578 1 13756.3 13884.5 0.499 0.656
sc_HR 0.314 1 17498.1 14651.1 0.242 0.385
sc_rain 0.037 1 24287.1 16414.2 0.002 0.075
sc_rain_24 0.096 1 22809.6 15579.0 0.058 0.135 

html_summary(read.file = "./data/processed/regression_models/RH3.2.rds", n.posterior = 1000)
RH3.2
n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 0 0 8492.95 12550.7 1444793977
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.703 1 8492.95 12697.7 0.641 0.765
sc_temp 0.567 1 9626.63 12550.7 0.488 0.648
sc_HR 0.308 1 12079.96 13794.6 0.235 0.382
sc_rain 0.029 1 18138.82 15107.9 -0.007 0.066
sc_rain_48 0.005 1 16926.99 14915.2 -0.032 0.043

Moon

Decreasing the percentage of the moon illuminated causes an increase in A. lemur activity?

fit.4 <- brm(n_call | resp_rate(rec_time) ~ sc_moonlight + ar(p = 2, time = hour_diff,
    gr = hour), data = call_rate_hour, iter = iter, chains = chains, cores = chains,
    family = negbinomial(), prior = priors, file = "./data/processed/regression_models/moon4",
    file_refit = "always")
html_summary(read.file = "./data/processed/regression_models/moon4.rds", n.posterior = 1000)
moon4
n_call | resp_rate(rec_time) ~ sc_moonlight + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 0 0 19357 14632 1734290538
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.674 1 19357.0 14632.0 0.606 0.743
sc_moonlight -0.128 1 31443.5 15209.9 -0.194 -0.063

Night rain

If the night rain decreases can affect the A. lemur activity to increase?

fit.5.1 <- brm(n_call | resp_rate(rec_time) ~ sc_prev_temp + sc_rain + sc_rain_48 +
    ar(p = 2, time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
    cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/night_rain5.1",
    file_refit = "always")

fit.5.2 <- brm(n_call | resp_rate(rec_time) ~ sc_prev_temp + sc_rain + sc_rain_24 +
    ar(p = 2, time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
    cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/night_rain5.2",
    file_refit = "always")
html_summary(read.file = "./data/processed/regression_models/night_rain5.1.rds",
    n.posterior = 1000)
night_rain5.1
n_call | resp_rate(rec_time) ~ sc_prev_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 0 0 15449.5 13446.5 1903747053
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.676 1.001 15449.5 13446.5 0.607 0.745
sc_prev_temp 0.003 1 27583.8 14753.3 -0.046 0.052
sc_rain -0.006 1 29045.4 16584.0 -0.042 0.030
sc_rain_48 0.014 1 27750.2 15632.5 -0.027 0.054 

html_summary(read.file = "./data/processed/regression_models/night_rain5.2.rds",
    n.posterior = 1000)
night_rain5.2
n_call | resp_rate(rec_time) ~ sc_prev_temp + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 0 0 21570 14952.6 1800703315
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.673 1 21570.0 14952.6 0.604 0.742
sc_prev_temp 0.007 1 39064.3 16543.7 -0.039 0.054
sc_rain 0.001 1 36294.6 16198.6 -0.035 0.038
sc_rain_24 0.077 1 34404.5 17412.9 0.038 0.117

Previous Rain

Does an increase in the accumulated previous rain causes that A. lemur activity to increase?

fit.6.1 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2,
    time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
    cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/previous_rain6.1",
    file_refit = "always")

fit.6.2 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2,
    time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
    cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/previous_rain6.2",
    file_refit = "always")

fit.6.3 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain_24 + ar(p = 2,
    time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
    cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/night_rain6.3",
    file_refit = "always")
html_summary(read.file = "./data/processed/regression_models/previous_rain6.1.rds",
    n.posterior = 1000)
previous_rain6.1
n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 0 0 8446.05 11423.8 20512253
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.701 1 8446.05 11423.8 0.636 0.764
sc_temp 0.304 1 10470.72 13123.6 0.254 0.355
sc_rain 0.036 1 16919.56 15417.1 0.000 0.074
sc_rain_48 0.007 1 17069.15 15935.4 -0.031 0.045 

html_summary(read.file = "./data/processed/regression_models/previous_rain6.2.rds",
    n.posterior = 1000)
previous_rain6.2
n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 0 0 11468.2 13337.9 1714781937
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.698 1 11468.2 13337.9 0.633 0.762
sc_temp 0.312 1 13834.2 13973.2 0.262 0.362
sc_rain 0.045 1 23413.4 15796.1 0.008 0.083
sc_rain_24 0.090 1 22867.5 15487.7 0.051 0.129 

html_summary(read.file = "./data/processed/regression_models/night_rain6.3.rds",
    n.posterior = 1000)
night_rain6.3
n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour)
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 4) flat 10000 4 1 5000 2 0 7697.31 11401.3 2086406674
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.699 1 7697.31 11401.3 0.636 0.760
sc_temp 0.572 1 8657.91 12183.9 0.493 0.650
sc_HR 0.318 1 10071.18 12501.7 0.244 0.391
sc_rain_24 0.093 1 14992.10 14450.6 0.055 0.131

Causal model combined effect size plot

Takes the effect sizes (and posteriors) from the right causal models

coef_table <- read.csv("./data/processed/summary_causal_model_table_13-01-23.csv",
    sep = "\t")

coef_table$variable <- coef_table$Label
coef_table$value <- coef_table$Estimate
coef_table$significance <- ifelse(coef_table$CI_low * coef_table$CI_high > 0, "sig",
    "no.sig")

posteriors_l <- lapply(1:nrow(coef_table), function(x) {

    # print(x)
    X <- readRDS(coef_table$Model[x])
    xdrws <- brms::as_draws(X)
    post <- xdrws$`1`[[paste0("b_", coef_table$Variable[x])]]
    out <- data.frame(variable = coef_table$Label[x], value = post, significance = coef_table$significance[x])
    return(out)
})

posteriors <- do.call(rbind, posteriors_l)

coef_table$variable <- factor(coef_table$Label, levels = sort(unique(coef_table$Label),
    FALSE))
posteriors$variable <- factor(posteriors$variable, levels = sort(unique(posteriors$variable),
    TRUE))


fill_values <- c("#FFB9DF", "#FF7598")
fill_values <- adjustcolor(fill_values, alpha.f = 0.5)

# creat plots gg_dists <-
ggplot2::ggplot(data = posteriors, ggplot2::aes(y = variable, x = value, fill = significance)) +
    ggplot2::geom_vline(xintercept = 0, col = "black", lty = 2) + ggdist::stat_halfeye(ggplot2::aes(x = value),
    .width = c(0.95), normalize = "panels", color = "transparent") + ggplot2::scale_fill_manual(values = fill_values,
    guide = "none") + ggplot2::geom_point(data = coef_table) + ggplot2::geom_errorbar(data = coef_table,
    ggplot2::aes(xmin = CI_low, xmax = CI_high), width = 0) + ggplot2::scale_color_manual(values = coef) +
    ggplot2::facet_wrap(~variable, scales = "free_y", ncol = 1) + ggplot2::theme_classic() +
    ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"), plot.margin = ggplot2::margin(0,
        0, 0, 0, "pt"), legend.position = "none", strip.background = ggplot2::element_blank(),
        strip.text = ggplot2::element_blank()) + ggplot2::labs(x = "Effect size",
    y = "Parameter")  #+

# ggplot2::xlim(range(c(posteriors_by_chain$value, 0)) * plot.area.prop)

ggsave(filename = "./figures/summary_effect_sizes_pooled_from_multiple_causal_models.jpeg")

Conditional plots

Measured based on the global model

Rain and temperature

glob_mod_24 <- readRDS("./data/processed/regression_models/global_rain_24.rds")

conditions <- data.frame(sc_temp = c(`Low temperature` = -1, `Mean temperature` = 0,
    `High temperature` = 1))

rain_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_rain", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 520)) +
    ggtitle("Current rain") + labs(x = "Rain", y = "Call activity (calls/hour)")

rain24_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_rain_24", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 520)) +
    ggtitle("Previous 24h rain") + labs(x = "Rain", y = "Call activity (calls/hour)")

glob_mod_48 <- readRDS("./data/processed/regression_models/global_rain_48.rds")

rain48_gg <- plot(conditional_effects(glob_mod_48, effects = "sc_rain_48", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 520)) +
    ggtitle("Previous 48h rain") + labs(x = "Rain", y = "Call activity (calls/hour)")


cowplot::plot_grid(rain_gg, rain24_gg, rain48_gg, nrow = 3)

Relative humidity and temperature

conditions <- data.frame(sc_temp = c(`Low temperature` = -1, `Mean temperature` = 0,
    `High temperature` = 1))

hr_temp_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_HR", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + labs(x = "Relative humidity",
    y = "Call activity (calls/hour)")

hr_temp_gg

conditions <- data.frame(sc_HR = c(`Low humidity` = -1, `Mean humidity` = 0, `High humidity` = 1))

temp_hr_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_temp", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + labs(x = "Temperature",
    y = "Call activity (calls/hour)")

temp_hr_gg

Relative humidity and rain

conditions <- data.frame(sc_rain = c(`Low rain` = -1, `Mean rain` = 0, `High rain` = 1))


rain_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_HR", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 75)) +
    ggtitle("Current rain")

rain24_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_HR", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 75)) +
    ggtitle("Previous 24h rain")

rain48_gg <- plot(conditional_effects(glob_mod_48, effects = "sc_HR", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 75)) +
    ggtitle("Previous 48h rain")

cowplot::plot_grid(rain_gg, rain24_gg, rain48_gg, nrow = 3)

conditions <- data.frame(sc_HR = c(`Low humidity` = -1, `Mean humidity` = 0, `High humidity` = 1))

rain_hr_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_rain", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ggtitle("Current rain") +
    labs(x = "Rain", y = "Call activity (calls/hour)")

rain24_hr_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_rain_24", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ggtitle("Previous 24h rain") +
    labs(x = "Rain previous", y = "Call activity (calls/hour)")

rain48_hr_gg <- plot(conditional_effects(glob_mod_48, effects = "sc_rain_48", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ggtitle("Previous 48h rain") +
    labs(x = "Rain previous", y = "Call activity (calls/hour)")

cowplot::plot_grid(rain_hr_gg, rain24_hr_gg, rain48_hr_gg, nrow = 3)

Moonlight and temperature

conditions <- data.frame(sc_temp = c(`Low temperature` = -1, `Mean temperature` = 0,
    `High temperature` = 1))

moon_temp_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_moonlight", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + labs(x = "Moonlight",
    y = "Call activity (calls/hour)")

moon_temp_gg

conditions <- data.frame(sc_moonlight = c(`Low moonlight` = -1, `Mean moonlight` = 0,
    `High moonlight` = 1))

temp_moon_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_temp", conditions = conditions),
    plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + labs(x = "Temperature",
    y = "Call activity (calls/hour)")

temp_moon_gg

 

Takeaways

 

Sum up results

Next steps

 

Session information

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] brmsish_1.0.0     lunar_0.2-1       ggplot2_3.4.0     knitr_1.42       
##  [5] kableExtra_1.3.4  HDInterval_0.2.2  readxl_1.3.1      posterior_1.3.1  
##  [9] cowplot_1.1.1     ggdist_3.2.0      brms_2.18.0       Rcpp_1.0.9       
## [13] viridis_0.6.2     viridisLite_0.4.1 remotes_2.4.2    
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.4.1      systemfonts_1.0.4    plyr_1.8.7          
##   [4] igraph_1.3.5         splines_4.1.0        crosstalk_1.2.0     
##   [7] TH.data_1.1-0        rstantools_2.2.0     inline_0.3.19       
##  [10] digest_0.6.31        htmltools_0.5.4      fansi_1.0.3         
##  [13] magrittr_2.0.3       checkmate_2.1.0      RcppParallel_5.1.5  
##  [16] matrixStats_0.62.0   xts_0.12.2           sandwich_3.0-1      
##  [19] svglite_2.1.0        prettyunits_1.1.1    colorspace_2.0-3    
##  [22] rvest_1.0.3          textshaping_0.3.5    xfun_0.36           
##  [25] dplyr_1.0.10         callr_3.7.3          crayon_1.5.2        
##  [28] jsonlite_1.8.4       lme4_1.1-27.1        survival_3.2-11     
##  [31] zoo_1.8-11           ape_5.6-2            glue_1.6.2          
##  [34] gtable_0.3.1         emmeans_1.8.1-1      webshot_0.5.4       
##  [37] distributional_0.3.1 pkgbuild_1.4.0       rstan_2.21.7        
##  [40] abind_1.4-5          scales_1.2.1         mvtnorm_1.1-3       
##  [43] DBI_1.1.1            miniUI_0.1.1.1       xtable_1.8-4        
##  [46] diffobj_0.3.4        stats4_4.1.0         StanHeaders_2.21.0-7
##  [49] DT_0.26              htmlwidgets_1.5.4    httr_1.4.4          
##  [52] threejs_0.3.3        ellipsis_0.3.2       pkgconfig_2.0.3     
##  [55] loo_2.4.1.9000       farver_2.1.1         sass_0.4.5          
##  [58] utf8_1.2.2           labeling_0.4.2       tidyselect_1.2.0    
##  [61] rlang_1.0.6          reshape2_1.4.4       later_1.3.0         
##  [64] munsell_0.5.0        cellranger_1.1.0     tools_4.1.0         
##  [67] cachem_1.0.6         cli_3.6.0            generics_0.1.3      
##  [70] ggridges_0.5.4       evaluate_0.20        stringr_1.5.0       
##  [73] fastmap_1.1.0        ragg_1.1.3           oce_1.7-8           
##  [76] yaml_2.3.7           processx_3.8.0       pbapply_1.6-0       
##  [79] nlme_3.1-152         mime_0.12            projpred_2.0.2      
##  [82] formatR_1.11         rstanarm_2.21.3      xml2_1.3.3          
##  [85] compiler_4.1.0       bayesplot_1.9.0      shinythemes_1.2.0   
##  [88] rstudioapi_0.14      gamm4_0.2-6          tibble_3.1.8        
##  [91] bslib_0.4.2          stringi_1.7.12       highr_0.10          
##  [94] ps_1.7.2             Brobdingnag_1.2-9    lattice_0.20-44     
##  [97] Matrix_1.5-1         nloptr_1.2.2.2       markdown_1.3        
## [100] shinyjs_2.1.0        tensorA_0.36.2       vctrs_0.5.2         
## [103] pillar_1.8.1         lifecycle_1.0.3      gsw_1.0-6           
## [106] jquerylib_0.1.4      bridgesampling_1.1-2 estimability_1.4.1  
## [109] httpuv_1.6.6         R6_2.5.1             promises_1.2.0.1    
## [112] gridExtra_2.3        codetools_0.2-18     boot_1.3-28         
## [115] colourpicker_1.2.0   MASS_7.3-54          gtools_3.9.3        
## [118] assertthat_0.2.1     withr_2.5.0          shinystan_2.6.0     
## [121] multcomp_1.4-17      mgcv_1.8-36          parallel_4.1.0      
## [124] grid_4.1.0           coda_0.19-4          minqa_1.2.4         
## [127] rmarkdown_2.20       shiny_1.7.3          base64enc_0.1-3     
## [130] dygraphs_1.1.1.6