How to specify an autoregressive correlation structure in a linear model when the data has multiple observations at each time step? I want to fit a multiple regression model that accounts for temporally-correlated errors (i.e., some sort of autoregressive correlation structure, like those provided by ARMA or ARIMA models). However, the problem I have is that my data set isn't a simple continuous time series (i.e., one that has only one data point at every single time step). In my data set, each time step may contain one, several or no data points. 
Because of this, I cannot use standard autoregressive correlation structures, such as those provided by corAR1() or corARMA() in R, or the function auto.arima() with xreg regression.
It might help to be more specific about my data set. I am analysing estimates of biological diversity through geological time (numbers of species, estimated using the fossil record) for roughly equally-sized geographic regions. I am using multiple regression to understand the role of various environmental variables on diversity.
Time points actually correspond to so-called "time bins"—lumped intervals of geological time with roughly equal durations (for example, the Late Cretaceous is one time bin). In some time bins, I have diversity estimates for several different geographic regions. In others, only one estimate is available, or even no estimates, if there are no known fossils. The identities/locations of geographic regions are not constant through time: they vary because the geographic distribution of fossil sites varies through time. 
Is there any way to specify an autoregressive correlation structure for this kind of more complex data set? I'm happy to use whatever modelling approach is necessary.
I think it would be reasonable to treat time steps with multiple data points as representing multiple measurements from the same subject (planet Earth). 
EDIT: The answer posted by user "HCQ" suggested that I a should fit a multilevel model with a random effect for geographic regions. This suggestion doesn't solve the problem, because geographic identities of data points are not constant through time (see revised question, above).
 A: It sounds like you wish to estimate a model of the following form:
$$y_{t,i}=\beta' x_t + \nu_t + \eta_{t,i},$$
where $t$ indexes time, $i$ indexes location and (e.g.):
$$\nu_t=\rho \nu_{t-1} + \varepsilon_t,$$
$$[\eta_{t,1},\dots,\eta_{t,N}]\sim N(0,\Sigma),$$
$$\varepsilon_{t}\sim N(0,\sigma_\varepsilon^2),$$
where:
$$\Sigma_{i,j}=\sigma_\eta^2 \exp{(-\gamma d_{i,j})},$$
and $d_{i,j}$ is (some measure of) the distance between locations $i$ and $j$.
I.e. you have multiple potentially noisy observations of the same unobserved underlying process (with spatially correlated noise).
(This can obviously be generalized to an arbitrary "ARIMA" type persistence structure.)
State space models such as this may be estimated with the Kalman filter. Missing observations (i.e. not observing all locations in each time period) are no problem if estimating with the Kalman Filter.
There is a description of the available tools for state space modelling in R in the paper "State Space Models in R" by Petris and Petrone.
A: I have seen your dataset, I think your data structure has two levels as in the following figure: https://drive.google.com/open?id=1R8cQFGk4Jhw1IHWam6Z2uhcFAZcTCF70.
I believe you should fit a multilevel model which gives geographic regions random effects to explain between-region variance, and I think it is improper to just ignore the variation among different regions which may lead to reduced statistical power. 
If you change your mind to fit a multilevel model, then I think you can easily fit a hierarchical poisson regression model in nlme() with the covariance structure you need.
