I will list 5 issues with F12/R12 methods, and then try to explain them the best I can:
Need for an auxiliary basis set
Most (if not all) F12/R12 methods require more than just a standard single-particle Gaussian basis set. For one of the most common elements (carbon) in electron structure calculations, you can see some of these auxiliary basis sets on Basis Set Exchange, for example by typing "F12" into the search field:
Notice that they only go up to QZ. This is a problem if you want to get close to the basis set limit: advocates for F12/R12 usually say that using these methods gains you 2-3 levels of Z (meaning that QZ-F12 might be as good as 6Z or 7Z without F12), but you will unlikely get results at the quality of 8Z using QZ-F12. A few days ago I posted on arXiv this paper in which we went up to aug-cc-pCV8Z for the carbon atom, and still basis set incompleteness error was the largest source of error in our final ionization energy: specifically, extrapolations to the CBS limit using FCI level calculations at aug-cc-pCV7Z and aug-cc-pCV8Z gave extrapolated values that were more than 12 cm$^{-1}$ away from the aug-cc-pCV8Z energy. That is 12x bigger than our targeted accuracy of $\pm$1 cm$^{-1}$, so ideally we would be doing calculations at 9Z or 10Z or even even higher. Since F12 methods are claimed to gain you only 2-3 levels of Z, this would mean we would need auxiliary basis sets for 6Z-F12 (at the very least) to beat our 8Z results. I have still never seen an cc-pCV5Z-F12 basis set although I do know that aug-cc-pV5Z-F12 have existed (albeit only since 2015). So to compete with my aug-cc-pCV8Z calculation, F12 methods are stuck with aug-cc-pCVQZ-F12 or aug-cc-pV5Z-F12, which is not as good enough. The auxiliary basis set also makes F12 method non-variational!
Loss of accuracy due to density fitting and other approximations
Let's say that an aug-cc-pCV5Z-12 basis set did exist, and did gain us 3Z, would it be able to compete with my aug-cc-pCV8Z calculation? This is where I would say that F12 might be good for gaining you 3Z in terms of giving 5Z quality results with a DZ basis set (where your goal is to get milli-hartree or kcal/mol precision) but not for getting 8Z quality results with a 5Z or even 6Z basis set when your goal is micro-hartree or cm$^{-1}$ precision. There's so many approximations (such as density fitting) that go into F12 methods, that can be ignored at the milli-hartree or kcal/mol level but not at the micro-hartree or cm$^{-1}$ level. Apart from density fitting, let me show this nice table from a 2018 paper which shows how many terms are ignored from the "full" F12 theory, in various variants of F12 which are implemented to make their computational cost more reasonable:
Empirical (or chosen in an ad hoc way) parameter $\gamma$
Many F12/R12 methods rely on choosing a value for $\gamma$ in the correlation factor:
\begin{equation}
\tag{1}
f_{12} = -\frac{1}{\gamma}e^{-\gamma r_{12}}
\end{equation}
This isn't such a big problem when aiming for kcal/mol precision (although some might say it's a bit inelegant), but if you really want to reach the basis set limit to reach "spectroscopic" accuracy, you may have to optimize this parameter, which then becomes very expensive.
Increased cost due to needing the density matrix
Some electronic structure methods that are implemented within the F12/R12 scheme may rely on a density matrix rather than just the wavefunction, for example in the case of FCIQMC-F12 where the cost of using the F12/R12 scheme is at least 2x the CPU cost and 2x the RAM cost.
Not all post-HF methods are implemented with F12/R12
Many of the most common post-HF methods have been implemented within some variant of an F12/R12/CT scheme. These post-HF methods include MP2, GF2, CCSD, CCSD(2), CCSD(3), CCSD(T), CCSDT, CCSDT(2), CCSDTQ, CASSCF, CASPT2, MRCI, MRCI+Q, DMRG and FCIQMC.
But many post-HF methods have not yet been implemented within an F12/R12 schemes. These include (to the best of my knowledge), CCSDT(Q) and any level of CC beyond CCSDTQ; CC(2), CC(3), CC(4), etc.; MR-ACPF, MR-AQCC, RASSCF, RASPT2, GASSCF, GASPT2, etc.; HCI, iFCI, MBE-FCI, etc. If you're aiming for reaching the CBS limit with high-accuracy, you probably also want to use a high-accuracy method to treat electron correlation, and the state-of-the-art post-HF methods are not implemented within the F12/R12 scheme (and when some of them are, such as FCIQMC, it is not very "black-box" as you may have to use several codes and really know what you're doing: the paper I mentioned needed to use a tweaked version of GAMESS plus a tweaked version of MPQC or a tweaked version of a code written by Takeshi Yanai which is not packaged for public use, in addition to NECI). Recently a new CT (canonical transcorrelation) method was developed for FCIQMC, but has the issue that you pointed out, about needing 3-electron integrals (but not 4-electron integrals!).
Conclusion: F12 methods certainly have their place and are of great value to people aiming for kcal/mol (or less) accuracy. If aiming for high accuracy (for example cm$^{-1}$) you will suffer from various limitations such as: not having the required auxiliary basis set to reach your desired level of basis set convergence, losing accuracy due to approximations in F12/R12 implementations, not being able to use state-of-the-art post-HF methods to treat the electron correlation very accurately, and/or the added cost will not be worth it and you might be better off just using a much larger basis set with no attempt to introduce explicit correlation, such as the aug-cc-pCV8Z that I used in the above mentioned paper.
