I'd like to test for trends in proportions of animals sampled over 12 years at six different beaches so that there is a separate trend test per beach. In the data below 'thisbeach' is the number of animals sampled at that particular beach and 'notthisbeach' is the number of animals sampled at all other beaches.
dat <- data.frame(fBeach = as.factor(rep(c("B6", "B5", "B2", "B1", "B4", "B3"), each=12)),
year = rep(seq(1:12),6),
notthisbeach = c(4990, 1294, 4346, 4082, 4628, 5576, 5939, 5664, 6108, 5195, 5564, 4079, 4694, 1224, 4052,
4019, 4457, 5242, 5259, 5198, 5971, 5208, 5168, 3722, 5499, 1288, 4202, 3988, 4773, 6018,
5952, 6100, 7308, 5821, 6030, 4546, 4698, 1300, 3884, 3943, 4717, 5911, 6110, 6076, 7606,
6138, 6514, 4767, 4830, 1307, 4886, 4327, 5285, 6344, 6627, 5824, 7305, 5991, 6073, 4647,
4584, 1162, 4200, 3956, 4710, 5664, 5533, 4828, 6082, 4697, 4721, 3529),
thisbeach = c(869, 221, 768, 781, 1086, 1375, 1145, 1074, 1968, 1415, 1250, 979, 1165, 291, 1062,
844, 1257, 1709, 1825, 1540, 2105, 1402, 1646, 1336, 360, 227, 912, 875, 941, 933,
1132, 638, 768, 789, 784, 512, 1161, 215, 1230, 920, 997, 1040, 974, 662, 470,
472, 300, 291, 1029, 208, 228, 536, 429, 607, 457, 914, 771, 619, 741, 411,
1275, 353, 914, 907, 1004, 1287, 1551, 1910, 1994, 1913, 2093, 1529))
glmmTMB indicates serial correlation is present;
require(glmmTMB)
require(DHARMa)
require(multcomp)
dat.TMB <- glmmTMB(cbind(notthisbeach,thisbeach) ~ year*fBeach, family = "betabinomial", data=dat)
simres <- simulateResiduals(dat.TMB,plot=T)
res = recalculateResiduals(simres, group = dat$year)
testTemporalAutocorrelation(res, time=unique(dat$year))
Durbin-Watson test
data: simulationOutput$scaledResiduals ~ 1
DW = 0.40903, p-value = 0.0002994
alternative hypothesis: true autocorrelation is not 0
However, I can't seem to find any examples including an autocorrelation structure in a model of this type.
Does anyone have any advice please?
I'm not sure I get the setup of the number of animals at this beach vs that beach, and depending on your research question you may need to do something different. However, basic patterns are easy enough to implement in glmmtmb. The example below shows an ar1.
dat <- data.frame(fBeach = as.factor(rep(c("B6", "B5", "B2", "B1", "B4", "B3"), each=12)),
year = rep(seq(1:12),6),
notthisbeach = c(4990, 1294, 4346, 4082, 4628, 5576, 5939, 5664, 6108, 5195, 5564, 4079, 4694, 1224, 4052,
4019, 4457, 5242, 5259, 5198, 5971, 5208, 5168, 3722, 5499, 1288, 4202, 3988, 4773, 6018,
5952, 6100, 7308, 5821, 6030, 4546, 4698, 1300, 3884, 3943, 4717, 5911, 6110, 6076, 7606,
6138, 6514, 4767, 4830, 1307, 4886, 4327, 5285, 6344, 6627, 5824, 7305, 5991, 6073, 4647,
4584, 1162, 4200, 3956, 4710, 5664, 5533, 4828, 6082, 4697, 4721, 3529),
thisbeach = c(869, 221, 768, 781, 1086, 1375, 1145, 1074, 1968, 1415, 1250, 979, 1165, 291, 1062,
844, 1257, 1709, 1825, 1540, 2105, 1402, 1646, 1336, 360, 227, 912, 875, 941, 933,
1132, 638, 768, 789, 784, 512, 1161, 215, 1230, 920, 997, 1040, 974, 662, 470,
472, 300, 291, 1029, 208, 228, 536, 429, 607, 457, 914, 771, 619, 741, 411,
1275, 353, 914, 907, 1004, 1287, 1551, 1910, 1994, 1913, 2093, 1529))
head(dat)
dim(dat)
require(glmmTMB)
require(DHARMa)
require(multcomp)
# function to test ar
testar <- function(mod, dat) {
simres <- simulateResiduals(mod,plot=T)
res <- recalculateResiduals(simres, group = dat$year)
print(testTemporalAutocorrelation(res, time=unique(dat$year)))
}
mod.TMB <- glmmTMB(cbind(notthisbeach,thisbeach) ~ year*fBeach, family = "betabinomial", data=dat)
testar(mod.TMB, dat)
# results
# Durbin-Watson test
#
# data: simulationOutput$scaledResiduals ~ 1
# DW = 0.40903, p-value = 0.0002994
# alternative hypothesis: true autocorrelation is not 0
mod.TMB.ar <- glmmTMB(cbind(notthisbeach,thisbeach) ~ as.factor(year) + fBeach + ar1(as.factor(year) + 0 | fBeach), family = "betabinomial", data=dat)
testar(mod.TMB.ar, dat)
#
# Durbin-Watson test
#
# data: simulationOutput$scaledResiduals ~ 1
# DW = 1.179, p-value = 0.1242
# alternative hypothesis: true autocorrelation is not 0
VarCorr(mod.TMB.ar)
# Conditional model:
# Groups Name Std.Dev. Corr
# fBeach as.factor(year)1 0.21692 0.464 (ar1)