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

1. interpret plots of the residuals to look for evidence of autocorrelation and,
2. 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

1. interpret plots of the residuals to look for evidence of autocorrelation and,
2. 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!

• r
• stan
• autocorrelation
• brms

Share Improve this question Follow edited Nov 20, 2024 at 15:00 Jaime Grimm asked Nov 20, 2024 at 3:13 Jaime GrimmJaime Grimm 122 bronze badges 3

• You have samples taken at the same time from different sites. Unless I'm misunderstanding something, that shouldn't be possible. – Roland Commented Nov 20, 2024 at 6:42
• @Roland thanks, I've fixed the dataframe so it better represents the real data. – Jaime Grimm Commented Nov 20, 2024 at 15:00
• OK. 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

Add a comment  | 

1 Answer 1

Sorted by: Reset to default 0

The 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:

1. 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.
2. I ran the above model with autocorrelation specified as:

acformula(~arma(time = seq, gr = species:site, cov=TRUE)