1 - Yes, in a linear circuit model the performance of the components is entirely characterized by their S matrices. You can extract the power imbalance, insertion loss, splitting ratio etc. from the S matrix. For example, in a 1x2 MMI, the transmission into output 2 would be |S_{21}|^2, transmission into port 3 is |S_{31}|^2, reflection back into port 1 would be |S_{11}|^2, insertion loss would be 10\textrm{log}_{10}(|S_{21}|^2 + |S_{31}|^2), etc.

2 - I would still recommend using the N port S-parameter element. You can use the experimental results to obtain the S parameters. This is a bit tricky if your experimental results are based on measuring power (like insertion loss and power imbalance), which lack phase information (unless the results are from some sort of interference experiment, I suppose). In other words, power measurements can give the magnitudes of the S parameters, but not the real and imaginary parts.

You could assume that the S parameters are purely real, but you should keep in mind that this will not include the phase information, which could play a role in any interference or resonance effects. However, if you are using the MMI as a power splitter, for example, this method would be sufficient.

3, 4 - You could use the method I mentioned above and calculate the S parameters based on experimental measurements. You could also obtain the S parameters using simulation and then verify that the simulation results are reasonable by comparing them to the experimental results. The third option you mentioned, using a scripted element and an analytical model, is also possible. It is up to you how exactly you want to model the devices in your circuit.

Let me know if you have any questions.