glmmTMB: AR1 models fail to converge I am trying to utilize the first-order autocorrelation [AR(1)] covariance structure abilities of the glmmTMB package (described here by Kasper Kristensen) to model experimental time series data collected from multiple locations. However, the models consistently fail to converge, with the warning message: In fitTMB(TMBStruc) :  Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting').
To diagnose the problem, I created some simulated AR1 data (see code below). However, the failure-to-converge problem persists even when varying the time series length, data distribution and link function, dispersion, zero inflation, and model random effects structure.
Is this a problem with my syntax, a data issue, or a glmmTMB bug? Many thanks for the assistance.
# Simulation test of AR(1) time series using glmmTMB

# ------------------------------------------------------------------------------------
# Install or load packages

#install.packages('bbmle')
#install.packages('dplyr')
#install.packages('ggplot2')
#install.packages('patchwork')
#install.packages('Matrix')
#install.packages('devtools')
#library('devtools')
#devtools::install_github("glmmTMB/glmmTMB/glmmTMB")

library(bbmle)
library(dplyr)
library(ggplot2)
library(glmmTMB)
library(patchwork)

# ------------------------------------------------------------------------------------
# Set parameters
t = 40       # Length of time series
phi = 0.8    # Amount of autocorrelation
stdev = 0.1  # Standard deviation
eff <- 3    # Effect of a treatment relative to a control

# Simulate data
sim.dat <- expand.grid(time = 1:t,
                      treatment = c("control", "manipulate"),
                      location = c("site1", "site2", "site3", "site4", "site5"))
sim.dat$plot <- paste(sim.dat$location, sim.dat$treatment, sep = "_")

# List description for AR(1) model
ar.sim <- list(order = c(1, 0, 0), ar = phi, sd = stdev)

# Simulate AR(1) data for each plot (= unique combination of location and treatment)
set.seed(1);  ar1.sim_site1.control      <- arima.sim(n = t, model = ar.sim)
set.seed(2);  ar1.sim_site2.control      <- arima.sim(n = t, model = ar.sim)
set.seed(3);  ar1.sim_site3.control      <- arima.sim(n = t, model = ar.sim)
set.seed(4);  ar1.sim_site4.control      <- arima.sim(n = t, model = ar.sim)
set.seed(5);  ar1.sim_site5.control      <- arima.sim(n = t, model = ar.sim)
set.seed(6);  ar1.sim_site1.manipulate   <- arima.sim(n = t, model = ar.sim) + eff
set.seed(7);  ar1.sim_site2.manipulate   <- arima.sim(n = t, model = ar.sim) + eff
set.seed(8);  ar1.sim_site3.manipulate   <- arima.sim(n = t, model = ar.sim) + eff
set.seed(9);  ar1.sim_site4.manipulate   <- arima.sim(n = t, model = ar.sim) + eff
set.seed(10); ar1.sim_site5.manipulate   <- arima.sim(n = t, model = ar.sim) + eff

sim.dat$response <- c(as.vector(ar1.sim_site1.control),
                      as.vector(ar1.sim_site1.manipulate),
                      as.vector(ar1.sim_site2.control),
                      as.vector(ar1.sim_site2.manipulate),
                      as.vector(ar1.sim_site3.control),
                      as.vector(ar1.sim_site3.manipulate),
                      as.vector(ar1.sim_site4.control),
                      as.vector(ar1.sim_site4.manipulate),
                      as.vector(ar1.sim_site5.control),
                      as.vector(ar1.sim_site5.manipulate)
                      )

# Add random noise
set.seed(100)
sim.dat$response <- sim.dat$response + runif(nrow(sim.dat), min = -0.5, max = 0.5)

# Convert simulated data to positive integers for poisson model
sim.dat$response <- round(sim.dat$response * 100) 
sim.dat$response <- sim.dat$response - min(sim.dat$response)
hist(sim.dat$response)

# Plot simulated data
ggplot(data = sim.dat, aes(x = time, y = response, color = treatment)) +
  geom_abline(intercept = 0, slope = 0, linetype = "dashed") +
  geom_line() +
  facet_wrap(~ location, ncol = 1, scales = 'fixed')

# ------------------------------------------------------------------------------------
# Fit models using glmmTMB

# Fixed-effects model, no correlation structure
sim.mod1 <- glmmTMB(response ~ treatment, family = poisson, data = sim.dat)

# Mixed-effect model (random intercept), no correlation structure
sim.mod2 <- glmmTMB(response ~ treatment + (1 | location), family = poisson, data = sim.dat)

# Fixed-effects model, AR(1) correlation structure
sim.mod3 <- glmmTMB(response ~ treatment + ar1(time + 0 | location), family = poisson, data = sim.dat)
# WARNING: Model convergence problem; non-positive-definite Hessian matrix.

# Mixed-effect model (random intercept), AR(1) correlation structure
sim.mod4 <- glmmTMB(response ~ treatment + (1 | location) + ar1(time + 0 | location), family = poisson, data = sim.dat)
# WARNING: Model convergence problem; non-positive-definite Hessian matrix.

AICtab(sim.mod1, sim.mod2, sim.mod3, sim.mod4)

# ------------------------------------------------------------------------------------
# Calculate and visualize ACF

# Define function
plot.acf.fun <- function(sim.mod){
  sim.dat$resid <- resid(sim.mod, type = "pearson")

  acf.dat <- sapply(unique(sim.dat$location), function(x){
acf(sim.dat$resid[sim.dat$location == x], lag.max = length(unique(sim.dat$time)) / 3, plot = FALSE)$acf
  })

  pacf.dat <- sapply(unique(sim.dat$location), function(x){
pacf(sim.dat$resid[sim.dat$location == x], lag.max = length(unique(sim.dat$time)) / 3, plot = FALSE)$acf
  }
  )

  acf.dat <- data.frame(acf.dat)
  pacf.dat <- data.frame(pacf.dat)

  colnames(acf.dat) <- (unique(sim.dat$location))
  colnames(pacf.dat) <- (unique(sim.dat$location))

  acf.dat <- acf.dat %>%
    dplyr::mutate(lag = 1:nrow(acf.dat) - 1) %>%
    tidyr::gather(key = "location", value = "acf", -lag)

  pacf.dat <- pacf.dat %>%
    dplyr::mutate(lag = 1:nrow(pacf.dat)) %>%
    tidyr::gather(key = "location", value = "pacf", -lag)

  acf.dat <- dplyr::left_join(acf.dat, pacf.dat, by = c("lag", "location"))

  # Calculate critical value (based on the lowest length of time series available)
  acf.dat$crit <- qnorm((1 + 0.95)/2) / sqrt(length(unique(sim.dat[sim.dat$location == "a1", ]$time)))

  # Plot ACF by location
  p1 <- ggplot(data = acf.dat, aes(x = lag, y = acf)) +
    ggtitle("Autocorrelation by location") +
    facet_wrap(~ location) +
    geom_bar(stat = "identity", width = 0.1, color = "black", fill = "black") +
    geom_hline(yintercept = 0) +
    geom_line(aes(y = crit), linetype = "dashed") +
    geom_line(aes(y = -crit), linetype = "dashed") +
    scale_y_continuous(breaks = seq(-10, 10, by = 2)/10, name = "ACF") +
    scale_x_continuous(breaks = 0:max(acf.dat$lag), name = "Lag") +
    theme_classic() +
    theme(aspect.ratio = 1)

  # Plot average PACF
  p2 <- ggplot(data = acf.dat[!is.na(acf.dat$pacf), ], aes(x = lag, y = pacf)) +
ggtitle("Average partial autocorrelation across locations") +
stat_summary(fun.data = mean_cl_boot) +
geom_hline(yintercept = 0) +
geom_line(aes(y = crit), linetype = "dashed") +
geom_line(aes(y = -crit), linetype = "dashed") +
coord_cartesian(ylim = c(-0.4, 1.0)) +
scale_y_continuous(breaks = seq(-1, 1, by = 0.2), name = "PACF") +
scale_x_continuous(limits = c(0.95, max(acf.dat$lag)), breaks = 1:max(acf.dat$lag), name = "Lag") +
    theme_classic() +
    theme(aspect.ratio = 1)

  p1 + p2
}

plot.acf.fun(sim.mod1) # Fixed-effects model, no correlation structure
plot.acf.fun(sim.mod2) # Mixed-effect model (random intercept), no correlation structure
plot.acf.fun(sim.mod3) # Fixed-effects model, AR(1) correlation structure
plot.acf.fun(sim.mod4) # Mixed-effect model (random intercept), AR(1) correlation structure



 A: I was able to get your example to run by turning time into a factor variable (disappointing, I know) :)
Here's a working example based loosely on the Ben Bolker's post here
library(tidyverse)
library(glmmTMB)
library(gsarima)

# experimental design
t <- 20
locs <- 5
treats <- c("control","manipulate")
eff <- 3
N <- t * locs * length(treats)

# gsarima parameters (see gsarima vignette)
ar <- 0.8
intercept <- 3
frequency <- 1
X=matrix(c(rep(intercept, N+length(ar))), ncol=1)

# Simulate poisson AR(1)
y.sim <- garsim(n=(N+length(ar)), phi=ar, beta=c(1), link= "identity",
                family= "poisson", minimum = 0, X=X)
y<-y.sim[(1+length(ar)):(N+length(ar))]
tsy<-ts(y, freq=frequency)

# Create experimental data
exp_data <- expand_grid(location = 1:locs,
                        time = factor(1:t),
                        treatment = factor(treats)) %>% 
  mutate(response = tsy,
         response = ifelse(treatment == "manipulate",
                           response + eff,
                           response))

mod_ar1 <- glmmTMB(response ~ treatment + (1|location) + 
                     ar1(time + 0|location),
                   data=exp_data,family=poisson)

