Unless specifying something like a Causal Forests, which is quite inappropriate here, it is somewhat questionable using a panel format with Random Forests. I would think that you could just treat the bands of each image as an independent measurement and include them as separate parameters (eg., y1, b1_t1, b1_t2, b2_t1, b2_t2, ...). This would require assigning the raster values of all images to each observation and not worry about specifying the data in a temporal model structure. Keep in mind that you are working with a scaled hypervolume and not pairwise additive functions so, you do not need to think in terms of each parameter having a one-to-one linear relationship with y. I would however test the model with the inclusion of a factorial parameter indicating the time period that the plot was collected (eg, year).
An interesting aspect of this approach would be to look at the estimated probabilities for each class and identify the observations that do not fit well into any class. You could then look at the distribution(s) of the temporal measures to identify stochasticity in the temporal data and how it relates to class transition and emergent "fuzzy" classes.
You could also use some application of feature engineering (following work such as Fanaee-T & Gama 2013) to derive metrics representing the temporal process thus, collapsing your time series into a design matrix. One R library that comes to mind in supporting this type of analysis is
timetk (Dancho & Vaughan 2018), using
recipes for the data processing pipeline. This workflow fits into the new tidymodels framework and is aimed at supporting time-series data in machine learning models. This however, may be a challenge to apply across a raster stack to provide metrics supporting a spatial estimate.
Fanaee-T, H, & J. Gama (2013) Event labeling combining ensemble detectors and background knowledge. Progress in Artificial Intelligence. pp. 1-15, Springer Berlin Heidelberg.