AR(1) correlation structure on errors with MixedModels.jl

Specific Domains Statistics
1

Hi all, I’m looking at the MixedModels.jl package and I’m curious if it (or something similar in Julia) can fit the following type of model I used nlme to fit in R. The part I can’t figure out is if there is a simple way to specify an AR(1) correlation structure on the errors like is available via corAR1: AR(1) Correlation Structure in nlme: Linear and Nonlinear Mixed Effects Models in R.

The model is from Simon Wood’s GAM book, for the curious.

library(nlme)
data(Loblolly)
Loblolly$age <- Loblolly$age - mean(Loblolly$age)
lmc <- lmeControl(niterEM=500,msMaxIter=100)
m0 <- lme(
  height ~ age + I(age^2) + I(age^3),
  Loblolly,
  random = list(Seed = ~ age + I(age^2) + I(age^3)), 
  correlation = corAR1(form = ~ age|Seed),
  control=lmc
)
2

You can’t currently do this with MixedModels.jl (see my notes here)

It looks like Metida.jl supports this though: Examples · Metida.jl

1 Like
3

Thanks, is this a fundamental constraint? I’m curious how the formulation of the numerical problem makes it difficult to use different covariance structures, since they are necessary for my case. Could you point me to where I can learn more about the implementation in the package?

4

The package documentation has details on the parameter estimation.

The biggest issue for AR1 is that the formulation depends rather heavily on the residual error term being a constant multiple of the identity matrix. In theory, it’s possible to implement this via iterative reweighting, like is done for GLMMs, but I have a very long backlog of FOSS improvements to do and likely won’t get around to this anytime soon…

1 Like
5

Thanks! I’ll read that page soon.

I see a PIRLS loop here for something MixedModels.jl/src/generalizedlinearmixedmodel.jl at main · JuliaStats/MixedModels.jl · GitHub, would that be possible to adapt for estimation of correlation terms as well? I implemented something like this as a one-off piece of code in the past for something similar in Julia, perhaps I could see how to work it out.

1 Like
6

MixedModel.jl and Metida.jl works with different implementation of optimizing REML/ML function. MixedModels using the penalized, iteratively reweighted, least squares (PIRLS) algorithm to determine parameters and it very fast (but works only with random effects and not support repeated[just now, maybe in theory it can be implemented]), Metida.jl - direct optimize REML function and it’s slowly, but can fit any covariance structure - random and repeated effects. If you want AR(1) - it is fully supported for random and repeated effects in Metida.jl and now it works even with rank deficient fixed effects.

1 Like