How to model multiple multivariate time series for binary classification? The data consist of repeated measurements of multiple individuals, where more variables are measured at each time point. Individuals fall into mutually exclusive groups. Number of measurements may differ with respect to individual, i.e some individuals have longer and some have shorter history. The the grid of time points is equidistant. Each observation has a binary label, typically long history of 0-s and some of them then switch to 1. How to predict the probability of the label being 1 in next step given the history? Ideally the model should be interpretable, i.e. if I say there is 80% chance that the individual A1 will become 1 in next time step and 60% chance it will become 1 in second step given it did not become 1 in first forward step, I would like to have some explanation why there is this chance...
 A: I will use the notation $y_k(t)$ to refer to the observation at time $t$ for the $k$-th individual.
The simplest solution I can see, if we assume little to no measurement error, would be to develop 2 states Markov Chain model with probabilities $p_{ij} = P\big( y(t + 1) = j \, | \, y(t) = i \big)$ such as $p_{00} + p_{01} = 1$ and $p_{10} +p_{11} = 1$.
You can further add individual-dependence by specifying mixed effects on $p_{10}$ and $p_{01}$, for example with $\operatorname{logit}(p_{01}^{(k)}) \sim \mathcal{N}(\mu_{01}, \sigma_{01}^2)$ and  $\operatorname{logit}(p_{10}^{(k)}) \sim \mathcal{N}(\mu_{10}, \sigma_{10}^2)$.
If we are interested in higher order lags (e.g. $t-2$ and $t-1$ to predict $t$) we can simply increase the number of states in the Markov Chain.
This model can be implemented with a standard logistic regression.
In R, if you have a dataframe df with columns y (observation at $t$), y_lag (observation at $t-1$), and k, you could fit this with lme4::glmer(y ~ (y_lag + 1 | k), family = "binomial").
The parametrisation would be different than the one described above, as the intercept coefficient would correspond to $\operatorname{logit}(p_{01})$ and the intercept + slope would correspond to $\operatorname{logit}(p_{11})$.
An alternative would be to develop a bespoke Bayesian model using the probabilistic programming language Stan for example, this is usually what I do in my work.
For example, if we assume measurement errors, we can introduce the latent variable $p_k(t) = P(y_k(t) = 1)$ and assume a mixed effect autoregressive model on $\operatorname{logit}p_k(t)$ (taking the logit ensures that we are in $\mathbb{R}$).
For instance, for a first order autoregressive model: $\operatorname{logit}p_k(t + 1) \sim \mathcal{N}\big(\alpha_k\operatorname{logit}p_k(t) + b_k  , \sigma^2\big)$, where $\alpha_k$ and $b_k$ are the individual-dependent autocorrelation parameter and intercept, respectively. Finally, we can specify a mixed effect on $\alpha_k \sim \mathcal{N}(\mu_\alpha, \sigma_\alpha^2)$ and $b_k \sim \mathcal{N}(\mu_b, \sigma_b^2)$.
With a Bayesian model, we would need set priors but then this is context-specific so I cannot address this here.
