library(ggplot2)
nitt <- 10000
bs = 30
library(MCMCglmm)
## Loading required package: Matrix
## Loading required package: coda
## Loading required package: ape
deg_dat <- readRDS("../data/processed/degradation_parameters.RDS")
deg_dat$species.height <- as.numeric(deg_dat$species.stratum) - 0.5
pca <- prcomp(deg_dat[, c("blur.ratio", "envelope.correlation", "SNR", "excess.attenuation")], scale. = TRUE)
deg_dat$PC1 <- pca$x[, 1]
summary(pca)
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.4470 1.0243 0.7095 0.59459
## Proportion of Variance 0.5235 0.2623 0.1258 0.08839
## Cumulative Proportion 0.5235 0.7858 0.9116 1.00000
pca$rotation
## PC1 PC2 PC3 PC4
## blur.ratio 0.4897179 -0.5138268 0.5505624 0.4393626
## envelope.correlation -0.5207786 0.4589055 0.4319763 0.5758402
## SNR -0.5298089 -0.4128413 0.5127233 -0.5347705
## excess.attenuation -0.4563645 -0.5957767 -0.4973856 0.4351886
agg_deg <- aggregate(cbind(blur.ratio, envelope.correlation, SNR, excess.attenuation, PC1) ~ distance + species.stratum + playback.stratum + species, data = deg_dat, mean)
ggplot(agg_deg, aes(x = distance, y = blur.ratio, group = species, col = species.stratum)) +
geom_point() +
geom_line() +
scale_color_viridis_d() +
labs(x = "Distance (m)", y = "Blur ratio", col = "Species stratum") +
facet_wrap(~ playback.stratum) +
theme_classic(base_size = bs)
ggplot(agg_deg, aes(x = distance, y = envelope.correlation, group = species, col = species.stratum)) +
geom_point() +
geom_line() +
scale_color_viridis_d() +
labs(x = "Distance (m)", y = "Envelope correlation", col = "Species stratum") +
facet_wrap(~ playback.stratum) +
theme_classic(base_size = bs)
ggplot(agg_deg, aes(x = distance, y = SNR, group = species, col = species.stratum)) +
geom_point() +
geom_line() +
scale_color_viridis_d() +
labs(x = "Distance (m)", y = "Signal-to-noise ratio", col = "Species stratum") +
facet_wrap(~ playback.stratum) +
theme_classic(base_size = bs)
ggplot(agg_deg, aes(x = distance, y = excess.attenuation, group = species, col = species.stratum)) +
geom_point() +
geom_line() +
scale_color_viridis_d() +
labs(x = "Distance (m)", y = "Excess attenuation", col = "Species stratum") +
facet_wrap(~ playback.stratum) +
theme_classic(base_size = bs)
ggplot(agg_deg, aes(x = distance, y = PC1, group = species, col = species.stratum)) +
geom_point() +
geom_line() +
scale_color_viridis_d() +
labs(x = "Distance (m)", y = "PC1 (overall degradation)", col = "Species stratum") +
facet_wrap(~ playback.stratum) +
theme_classic(base_size = bs)
blur ratio
mod1 <- MCMCglmm(blur.ratio ~ height + species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod2 <- MCMCglmm(blur.ratio ~ height * species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod3 <- MCMCglmm(blur.ratio ~ species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod4 <- MCMCglmm(blur.ratio ~ height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
br.mods <- list(mod1, mod2, mod3, mod4)
saveRDS(br.mods, "blur.ratio.models.RDS")
br.mods <- readRDS("../blur.ratio.models.RDS")
dics <- sapply(br.mods, function(x) x$DIC)
dics <- dics - min(dics)
names(dics) <- paste0("mod", 1:4)
sort(dics)
## mod4 mod1 mod2 mod3
## 0.0000000 0.1085204 2.3235323 23.3519324
summary(br.mods[[which.min(dics)]])
##
## Iterations = 3001:9991
## Thinning interval = 10
## Sample size = 700
##
## DIC: -32380.2
##
## G-structure: ~transect
##
## post.mean l-95% CI u-95% CI eff.samp
## transect 0.0001261 1.194e-05 0.0004238 700
##
## ~signal.type
##
## post.mean l-95% CI u-95% CI eff.samp
## signal.type 0.0002393 2e-04 0.0002821 700
##
## ~distance
##
## post.mean l-95% CI u-95% CI eff.samp
## distance 0.003712 5.606e-05 0.008039 700
##
## ~species
##
## post.mean l-95% CI u-95% CI eff.samp
## species 0.003096 0.001217 0.005879 700
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.002461 0.002399 0.002536 700
##
## Location effects: blur.ratio ~ height
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.143471 0.077015 0.206211 700.0 0.0143 *
## height -0.007437 -0.009025 -0.005472 825.5 <0.001 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
excess.attenuation
mod1 <- MCMCglmm(excess.attenuation ~ height + species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod2 <- MCMCglmm(excess.attenuation ~ height * species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod3 <- MCMCglmm(excess.attenuation ~ species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod4 <- MCMCglmm(excess.attenuation ~ height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
ea.mods <- list(mod1, mod2, mod3, mod4)
saveRDS(ea.mods, "excess.attenuation.models.RDS")
ea.mods <- readRDS("../excess.attenuation.models.RDS")
dics <- sapply(ea.mods, function(x) x$DIC)
dics <- dics - min(dics)
names(dics) <- paste0("mod", 1:4)
sort(dics)
## mod2 mod1 mod4 mod3
## 0.000000 7.261564 7.516969 262.118659
summary(ea.mods[[which.min(dics)]])
##
## Iterations = 3001:9991
## Thinning interval = 10
## Sample size = 700
##
## DIC: 13224.8
##
## G-structure: ~transect
##
## post.mean l-95% CI u-95% CI eff.samp
## transect 0.1192 0.01779 0.3618 700
##
## ~signal.type
##
## post.mean l-95% CI u-95% CI eff.samp
## signal.type 0.008934 0.006482 0.01128 490.4
##
## ~distance
##
## post.mean l-95% CI u-95% CI eff.samp
## distance 39.61 0.8046 89.19 700
##
## ~species
##
## post.mean l-95% CI u-95% CI eff.samp
## species 0.06383 0.02233 0.1183 755.6
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.2043 0.1988 0.2097 423.3
##
## Location effects: excess.attenuation ~ height * species.height
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 5.30505 0.59968 11.23076 700 0.0429 *
## height -0.21207 -0.23316 -0.18897 700 <0.001 **
## species.height -0.04132 -0.16778 0.09561 700 0.5029
## height:species.height -0.03078 -0.04655 -0.01637 1044 <0.001 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
envelope.correlation
mod1 <- MCMCglmm(envelope.correlation ~ height + species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod2 <- MCMCglmm(envelope.correlation ~ height * species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod3 <- MCMCglmm(envelope.correlation ~ species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod4 <- MCMCglmm(envelope.correlation ~ height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
ec.mods <- list(mod1, mod2, mod3, mod4)
saveRDS(ec.mods, "envelope.correlation.models.RDS")
ec.mods <- readRDS("../envelope.correlation.models.RDS")
dics <- sapply(ec.mods, function(x) x$DIC)
dics <- dics - min(dics)
names(dics) <- paste0("mod", 1:4)
sort(dics)
## mod1 mod2 mod4 mod3
## 0.00 45602.23 45611.01 45866.70
summary(ec.mods[[which.min(dics)]])
##
## Iterations = 3001:9991
## Thinning interval = 10
## Sample size = 700
##
## DIC: -32378.41
##
## G-structure: ~transect
##
## post.mean l-95% CI u-95% CI eff.samp
## transect 0.0001216 1.105e-05 0.0003533 700
##
## ~signal.type
##
## post.mean l-95% CI u-95% CI eff.samp
## signal.type 0.0002408 0.0002007 0.0002885 700
##
## ~distance
##
## post.mean l-95% CI u-95% CI eff.samp
## distance 0.002425 4.159e-05 0.006577 700
##
## ~species
##
## post.mean l-95% CI u-95% CI eff.samp
## species 0.002917 0.001171 0.005227 700
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.00246 0.002389 0.002527 750
##
## Location effects: blur.ratio ~ height + species.height
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.120336 0.062820 0.195388 782.7 0.0114 *
## height -0.007497 -0.009357 -0.005756 700.0 <0.001 **
## species.height 0.016972 -0.012072 0.045922 700.0 0.2257
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Signal-to-noise ratio
mod1 <- MCMCglmm(SNR ~ height + species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod2 <- MCMCglmm(SNR ~ height * species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod3 <- MCMCglmm(SNR ~ species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod4 <- MCMCglmm(SNR ~ height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
snr.mods <- list(mod1, mod2, mod3, mod4)
saveRDS(snr.mods, "snr.models.RDS")
snr.mods <- readRDS("../snr.models.RDS")
dics <- sapply(snr.mods, function(x) x$DIC)
dics <- dics - min(dics)
names(dics) <- paste0("mod", 1:4)
sort(dics)
## mod1 mod2 mod4 mod3
## 0.000000 1.162794 1.557138 7.924739
summary(ec.mods[[which.min(dics)]])
##
## Iterations = 3001:9991
## Thinning interval = 10
## Sample size = 700
##
## DIC: -32378.41
##
## G-structure: ~transect
##
## post.mean l-95% CI u-95% CI eff.samp
## transect 0.0001216 1.105e-05 0.0003533 700
##
## ~signal.type
##
## post.mean l-95% CI u-95% CI eff.samp
## signal.type 0.0002408 0.0002007 0.0002885 700
##
## ~distance
##
## post.mean l-95% CI u-95% CI eff.samp
## distance 0.002425 4.159e-05 0.006577 700
##
## ~species
##
## post.mean l-95% CI u-95% CI eff.samp
## species 0.002917 0.001171 0.005227 700
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.00246 0.002389 0.002527 750
##
## Location effects: blur.ratio ~ height + species.height
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.120336 0.062820 0.195388 782.7 0.0114 *
## height -0.007497 -0.009357 -0.005756 700.0 <0.001 **
## species.height 0.016972 -0.012072 0.045922 700.0 0.2257
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
PC1
mod1 <- MCMCglmm(PC1 ~ height + species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod2 <- MCMCglmm(PC1 ~ height * species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod3 <- MCMCglmm(PC1 ~ species.height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
mod4 <- MCMCglmm(PC1 ~ height, random = ~ transect + signal.type + distance + species, data = deg_dat, verbose = FALSE, nitt = nitt)
pc.mods <- list(mod1, mod2, mod3, mod4)
saveRDS(pc.mods, "pc1.models.RDS")
pc.mods <- readRDS("../pc1.models.RDS")
dics <- sapply(pc.mods, function(x) x$DIC)
dics <- dics - min(dics)
names(dics) <- paste0("mod", 1:4)
sort(dics)
## mod4 mod1 mod2 mod3
## 0.0000000 0.1372839 1.1420333 7.0963914
summary(pc.mods[[which.min(dics)]])
##
## Iterations = 3001:9991
## Thinning interval = 10
## Sample size = 700
##
## DIC: 24206.85
##
## G-structure: ~transect
##
## post.mean l-95% CI u-95% CI eff.samp
## transect 0.06689 0.005389 0.1973 700
##
## ~signal.type
##
## post.mean l-95% CI u-95% CI eff.samp
## signal.type 0.0464 0.03762 0.05618 549.8
##
## ~distance
##
## post.mean l-95% CI u-95% CI eff.samp
## distance 5.622 0.1368 22.69 700
##
## ~species
##
## post.mean l-95% CI u-95% CI eff.samp
## species 0.5279 0.2169 0.9457 789.5
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.5849 0.5683 0.6026 700
##
## Location effects: PC1 ~ height
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.17256 -2.44449 2.51412 700 0.811
## height -0.05787 -0.08423 -0.03259 700 <0.001 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1