I have a dataset of qPCR reads (relative density) from eDNA for seven species at two sites. The samples were taken over a period of two days. I anticipate some autocorrelation between samples due to some DNA persisting in the sampling equipment between sites.
I am modeling DNA quantity between sites with species ID as a random effect using the brm() function from the brms package in R. In this framework, I'm unsure how to
- interpret plots of the residuals to look for evidence of autocorrelation and,
- account for autocorrelation in the model prediction
Reproducible data
data <- data.frame(site = c(rep("siteA", times=35), rep("siteB", times=35)),
seq = c(rep(1:5, each=7), rep(6:10, each=7)),
species = rep(c("a", "b", "c", "d","e", "f", "g")),
dnaquant = rgamma(70, 5))
Where seq is the order the samples were taken in and the variable I'm hoping to use to account for autocorrelation. I suspect my trouble might be coming from how my dataframe is specified (see error below) - because this is eDNA, I have information for multiple species coming from one sample (seq). Therefore, the seq value is not unique.
I am running a gamma hurdle model to account for zero inflation in my real dataset. So far my model runs without accounting for autocorrelation, specified like this:
m4 <- brm(bf(dnaquant ~ site + (1 + site|species), hu ~ site),
data = data, family = hurdle_gamma(), chain = 2, cores = 2)
I am trying to add an autocorrelation component like this:
m4 <- brm(bf(dnaquant ~ site + (1 + site|species), hu ~ site) + acformula(~arma(time = seq, gr = species, cov=TRUE), data = data, family = hurdle_gamma(), chain = 2, cores = 2)
And am getting, Error: Time points within groups must be unique.
I'd appreciate any thoughts!
I have a dataset of qPCR reads (relative density) from eDNA for seven species at two sites. The samples were taken over a period of two days. I anticipate some autocorrelation between samples due to some DNA persisting in the sampling equipment between sites.
I am modeling DNA quantity between sites with species ID as a random effect using the brm() function from the brms package in R. In this framework, I'm unsure how to
- interpret plots of the residuals to look for evidence of autocorrelation and,
- account for autocorrelation in the model prediction
Reproducible data
data <- data.frame(site = c(rep("siteA", times=35), rep("siteB", times=35)),
seq = c(rep(1:5, each=7), rep(6:10, each=7)),
species = rep(c("a", "b", "c", "d","e", "f", "g")),
dnaquant = rgamma(70, 5))
Where seq is the order the samples were taken in and the variable I'm hoping to use to account for autocorrelation. I suspect my trouble might be coming from how my dataframe is specified (see error below) - because this is eDNA, I have information for multiple species coming from one sample (seq). Therefore, the seq value is not unique.
I am running a gamma hurdle model to account for zero inflation in my real dataset. So far my model runs without accounting for autocorrelation, specified like this:
m4 <- brm(bf(dnaquant ~ site + (1 + site|species), hu ~ site),
data = data, family = hurdle_gamma(), chain = 2, cores = 2)
I am trying to add an autocorrelation component like this:
m4 <- brm(bf(dnaquant ~ site + (1 + site|species), hu ~ site) + acformula(~arma(time = seq, gr = species, cov=TRUE), data = data, family = hurdle_gamma(), chain = 2, cores = 2)
And am getting, Error: Time points within groups must be unique.
I'd appreciate any thoughts!
Share Improve this question edited Nov 20, 2024 at 15:00 Jaime Grimm asked Nov 20, 2024 at 3:13 Jaime GrimmJaime Grimm 12 bronze badges 3 |1 Answer
Reset to default 0The dummy dataset I provided didn't represent the real issue in my data well (my bad). I have resolved the issue with my real data taking the following steps:
- My qPCR was run in triplicate for each sample so I had 3 species estimates for each sample. I took the mean of those three estimates.
- I ran the above model with autocorrelation specified as:
acformula(~arma(time = seq, gr = species:site, cov=TRUE)
m4
with ac fits without errors. I get a number of warnings but that is to be expected with this dataset where the DV doesn't depend on the IVs. – Roland Commented Nov 20, 2024 at 15:09