I aggree, it would be fascinating to add such a feature. I would probably iterate the simulated RNA molecues, once expressed, and carry out Bernoulli trials for each intron to be retained, with empirically estimated values. There are a couple of technical issues though, currently the most important that the annotation is read only twice, once in the beginning and then at the sequencing step. I put it as a request for enhancement.
A currently possible workaround:
Read in the annotation and perform the gene expression. At this stage, the .PRO file already contains expression values (i.e., number of molecules) for each transcript. Iterate over these numbers and determine the location and frequency of retained introns. For each combination of retained introns, add a corresponding description to the initial GTF file, e.g. exons with transcript_ID "myID_unspliced1", and a line containing the fully spliced length plus the retained introns, and the frequency of observations to the .PRO file. Continue the simulation with these files.
It is a bit of work, but it should give a realistic model. Modelling intron retention at the stage of sequencing in my opinion can produce misleading results, as the fragmentation and RT have already been carried out on a molecule with different physical attributes.