# Fitting a multilevel AR1 in R

I have some short grouped time series data. I would like to fit a dynamic multilevel regression model in R, with random coefficients for the mean and first order auto-correlation in each group, and with no cross correlation between the two variance parameters; i.e. this model:

$y_{i,t} - \mu_{i} = \rho_{i} (y_{i,t-1} - \mu_{i}) + \epsilon_{i,t}$

$\mu_{i} = \mu + v_{\mu}$

$\rho_{i} = \rho + v_{\rho}$

$\epsilon_{i,t} \sim N(0,\sigma_\epsilon^2), v_{\mu} \sim N(0,\sigma_\mu^2), v_{\rho} \sim N(0,\sigma_\rho^2)$

The best I can do so far in R is:

library(nlme)
library(dplyr)

#create toy data set
df0 <- Orthodont %>%
group_by(Subject) %>%
mutate(lag1=lag(distance)) %>%
filter(!is.na(lag1))

#multilevel model, mean (not Subject specific mean) centered
m1 <- lme(fixed = distance ~ I(lag1 - mean(distance)), data=df0,
random= list(Subject = pdDiag(~ + I(lag1 - mean(distance)))) )
m1
# Linear mixed-effects model fit by REML
#   Data: df0
#   Log-restricted-likelihood: -164.8976
#   Fixed: distance ~ I(lag1 - mean(distance))
#              (Intercept) I(lag1 - mean(distance))
#               25.7907622                0.8289975
#
# Random effects:
#  Formula: ~+I(lag1 - mean(distance)) | Subject
#  Structure: Diagonal
#          (Intercept) I(lag1 - mean(distance)) Residual
# StdDev: 7.892818e-05                0.3031237 1.675277
#
# Number of Observations: 81
# Number of Groups: 27


This 1) does not estimate ($\mu$) and 2) centers the distance lag ($y_{i,t-1}$) on the population mean ($\mu$) rather than the subject ($i$) specific mean ($\mu_i$) that I desire.

Does anyone know how to fit a multilevel AR(1) model that uses the subject (or group) means rather than the overall mean?