Only prior sampled: Help to setup time series analysis with intrinsic auto-regressive model

  • Operating System: Ubuntu 16.04
  • brms Version: 2.14.0

Dear forum,

I want to use brms (great package!) for an intrinsic auto-regressive model for time series analysis (see e.g. Chen & Ficetola).

With my current setup brms seems to sample from the prior only. Does anyone has an idea what goes wrong here?

To demonstrate the prior-only sampling, there is a sample example below.
With the autoregressive AR(1) model, the autocorrelation is correctly inferred and the true values are correctly predicted by the FitAr1 model.
With the intrinsic auto-regressive model, the prediction interval is the same for all data points.
When shuffling the data (i.e. removing the temporal autocorrelation), the intrinsic auto-regressive model shows the same estimates for the autocorrelation parameter sdcar as before (FitCar2 vs. FitCar1). In contrast, the AR(1) model (FitAr2) correctly infers the lack of autocorrelation.

Thanks a lot for any suggestions,
Torsten

library(brms)

# Helper function similar to 
# https://gitlab.com/wtchen/DNA_CAR_model/-/blob/master/ts_analysis.R
# but faster/less memory hungry
###########################################################
adjacency_matrix <- function(S, Degree = 1) {
  A <- matrix(0, S, S)
  for (i in 1:S) {
    BeforeDiag <- (i - Degree):(i - 1)
    BeforeDiag <- BeforeDiag[BeforeDiag >= (i - Degree) & BeforeDiag > 0 ]
    AfterDiag <- (i + 1):(i + Degree)
    AfterDiag <- AfterDiag[AfterDiag <= S]
    A[i, c(BeforeDiag, AfterDiag)] <- 1
  }
  return(A)
}
# Create some data
################
N <- 100
X <- seq(-5, 5, length.out = N)
Y <- 1 + 2 * X + rnorm(N, sd = 0.1)
Df <- data.frame(X = X, Y = Y, Group = "A")

AdjMat <- adjacency_matrix(N)
rownames(AdjMat) <- colnames(AdjMat) <- Df$Group 

# AR1 Model
FitAr1 <- brm(Y ~ ar(p = 1, gr = Group, cov = FALSE),
             data = Df,
             data2 = list(W = AdjMat),
             chains = 1, cores = 1)
summary(FitAr1)
# True values well predicted
head(Df)
head(predict(FitAr1)) 

# icar Model
FitCar1 <- brm(Y ~ car(W, gr = Group, type = "icar"),
              data = Df,
              data2 = list(W = AdjMat),
              chains = 1, cores = 1)
summary(FitCar1)
# High uncertainty for values
head(Df)
head(predict(FitCar1)) 

# Shuffle data (removes autocorrelation)
################################
Df <- Df[sample(N), ]
# AR1 Model
FitAr2 <- brm(Y ~ ar(p = 1, gr = Group, cov = FALSE),
             data = Df,
             data2 = list(W = AdjMat),
             chains = 1, cores = 1)
# ar[1] with large uncertainty because shuffling removes autocorrelation
summary(FitAr2) 

# icar Model
FitCar2 <- brm(Y ~ car(W, gr = Group, type = "icar"),
              data = Df,
              data2 = list(W = AdjMat),
              chains = 1, cores = 1)
# Same sdcar then for FitCar1
summary(FitCar2)

I think I found the solution: The car model works when the grouping factor gr within car() is unique for each observation.

The only thing to modify in the minimal example would be:

Df <- data.frame(X = X, Y = Y, Group = as.character(1:N))

Is it possible to tag the threat with something like solved?

Best,
Torsten

Good to hear you found a solution.

There should be an option to mark your comment ‘as a solution’ (but I don’t think you can for your original post).