0

我想测试 12 年来在六个不同海滩采样的动物比例趋势，以便每个海滩都有单独的趋势测试。在下面的数据中，“thisbeach”是在该特定海滩采样的动物数量，“notthisbeach”是在所有其他海滩采样的动物数量。

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 表示存在序列相关；

  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

但是，我似乎在这种类型的模型中找不到任何示例，包括自相关结构。

请问有人有什么建议吗？

4

1 回答 1

0

我不确定我是否设置了这个海滩与那个海滩的动物数量，并且根据您的研究问题，您可能需要做一些不同的事情。但是，基本模式很容易在glmmtmb. 下面的示例显示了一个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)
于 2021-12-31T15:47:15.247 回答