# Mixed effects model with an ARIMA correlation structure

I'm fitting a mixed effects model where each cluster is a time series. For this purpose, I'm using the R package nlme, whose main module (also nlme) may get a correlation structure argument of class corStruct. The package also comes with some implementations, one of them is corARMA which designs an ARMA(p,q) structure.

The problem is, my clusters are typically not stationary, and usually a single differencing $y_t'=y_t-y_{t-1}$ is required to get them to a stationary state, meaning I'm actually in need of an ARIMA(p,1,q) model.

The easiest solution would obviously be to manually transform my data and just use the supported ARMA(p,q) class. Another, more difficult solution, may be to dig in the code and program my own ARIMA subclass - I'd prefer to avoid such hacking as this seems like a deep rabbit hole.

However, both proposed solutions are mere patches that may cause some estimation errors down the road. My model contains some fixed effects, covariates that essentially aim to model and explain the trend $T_t$ component, the same trend that causes non-stationarity, which I'm removing with differencing.

So my main question is, how may I estimate the ARMA and covariates correctly?

The strategy I've gathered so far is as follows:

1. Estimate ARMA parameters on differenced $y_t'$ clusters.
2. Plug the ARMA parameters estimates into a corARMA object, then pass it as an argument to nlme and compute a model over the original $y_t$, so as to estimate the fixed and random effects.

Will this yield correct estimations? If not, why, and what should I do to fix my algorithm?

Thank you