Aims. Stellar flares serve as crucial indicators stellar magnetic activity. Radio emissions were detectable across all stages of stellar evolution. We define radio stars in our paper as stars with radio continuum emission in the frequency region of about 1.4 - 375 GHz from the catalog published by Wendker (1995, 2015). We also include detected radio stars from the LOFAR Two-metre Sky Survey (LoTSS) and the Australian Square Kilometre Array Pathfinder (ASKAP), and VLA Sky Survey (VLASS) surveys, and present comparative discussions. We used the light curves from the Transiting Exoplanet Survey Satellite (TESS) survey coupled with low- and medium-resolution spectra from the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) survey to investigate the magnetic activity and statistical properties of radio stars. Methods. We cross-matched the radio stars from the radio star catalog and the LoTSS, ASKAP, and VLASS surveys with the TESS survey. For the stars from the radio star catalog, we obtained a matched sample of 1,537 stars (hereafter Sample 1) and downloaded their TESS light curves, which include 4,001 light curves at a 120s cadence. For the stars from the LoTSS, ASKAP, and VLASS surveys, we obtain a matched sample of 98 stars (hereafter Sample 2) and downloaded their 120s TESS light curves, totaling 260. Subsequently, we employed repeated fitting techniques to distinguish stellar background light curves from flare events. For Sample 1, we successfully identify 12,155 flare events occurring on 856 stars. For Sample 2, we identify a total of 3,992 flare events on 86 stars. Furthermore, by cross-referencing our samples with the Gaia survey, the TESS Input Catalog, and the LAMOST survey, we obtain additional stellar parameters, facilitating the determination of relationships between stellar and flare parameters. Results. For stars in Sample 1, within the 12,155 flare events observed on the 856 flare-active radio stars, a majority of more than 97$<!PCT!>$ have durations of less than 2 hours, while for stars in Sample 2, all 3,992 flare events have a duration of less than 2 hours. We calculated the flare occurrence percentage for each flare-active radio star, observing a decrease as effective temperature increased for both Sample 1 and Sample 2. We derived values of the power-law index (for Sample 1, the alpha value is approximately 1.50pm 0.11 for single stars and 1.38pm 0.09 for binary stars, for Sample 2, the alpha value is about 1.47pm 0.11 for single stars and 1.42pm 0.09 for binary stars) for the cumulative flare frequency distribution. In both Sample 1 and Sample 2, stars with lower effective temperatures tend to exhibit increased activity. Using LAMOST spectra and the Gaia DR3 chromospheric activity index, we note that the Halpha equivalent width and Gaia Ca II IRT activity index of flare-active radio stars is significantly larger than that of nonflaring stars. An intriguing finding is our potential identification of a coronal rain candidate through asymmetry observed in the Halpha line.