  Modeling of longitudinal data often requires diffusion models that
incorporate overall time-dependent, nonlinear dynamics of multiple components
and provide sufficient flexibility for subject-specific modeling. This
complexity challenges parameter inference and approximations are inevitable. We
propose a method for approximate maximum-likelihood parameter estimation in
multivariate time-inhomogeneous diffusions, where subject-specific flexibility
is accounted for by incorporation of multidimensional mixed effects and
covariates. We consider $N$ multidimensional independent diffusions $X^i =
(X^i_t)_{0\leq t\leq T^i}, 1\leq i\leq N$, with common overall model structure
and unknown fixed-effects parameter $\mu$. Their dynamics differ by the
subject-specific random effect $\phi^i$ in the drift and possibly by (known)
covariate information, different initial conditions and observation times and
duration. The distribution of $\phi^i$ is parametrized by an unknown
$\vartheta$ and $\theta = (\mu, \vartheta)$ is the target of statistical
inference. Its maximum likelihood estimator is derived from the continuous-time
likelihood. We prove consistency and asymptotic normality of $\hat{\theta}_N$
when the number $N$ of subjects goes to infinity using standard techniques and
consider the more general concept of local asymptotic normality for less
regular models. The bias induced by time-discretization of sufficient
statistics is investigated. We discuss verification of conditions and
investigate parameter estimation and hypothesis testing in simulations.
