Implement hierarchical API calls

[jmarshallnz]

We need to implement a facility for random effects/hierarchical models with heterogeneous random effects, e.g. for the AR1 model. This needs to be done via a method that is familiar for the user (e.g. via corAR1()) but that might be tricky. Instead we might use some modification of the formula as is being done in flexLmer? Or just come up with something ourselves I guess...

(If we use corAR1() we can probably pull the group + time variable(s) out without using the rest of it).