Regarding the premises: I am not sure you can claim with generality that there is no second order phase transition in a model solely from the fact that the Gaussian fixed point is a sink. You would need to know whether there are any other fixed points of the renormalization group transformation with unstable directions that take you from the critical manifold to the Gaussian sink. Keep in mind that 'relevant' and 'irrelevant' interactions are always with reference to a fixed point of the RG transformation, so in principle you could have a Gaussian fixed point at which all interactions are irrelevant but there is another fixed point at which some of the interactions are relevant. If such a critical point exists then you could have a second order transition even if the Gaussian fixed point is a sink. That said, off the top of my head I don't know of any examples in which the Gaussian fixed point is a sink in the first place.
For the other questions:
1) Yes, the RG can identify first order transitions. See, for example, Goldenfeld's "Lectures on Phase transitions and the Renormalization Group," Ch. 9 Sec. 9.6.6 (end) and 9.7. In particular, Goldenfeld states that if a fixed point has an eigenvalues of the RG transformation equal to $\ell^d$, where $\ell$ is the coarse-graining factor in real space and $d$ is the dimension, is indicative of a first order transition. Goldenfeld calls this the "Nienhuis-Nauenberg" criterion, referencing Niehuis and Nauenberg, Phys. Rev. Lett. 35, 477 (1975).
2) You do not need to use RG to calculate statistical moments, so you would calculate them using your microscopic model as given. For example, if you have a nearest-neighbor Ising with Hamiltonian $\mathcal H(\mathbf{s}) = - J_1 \sum_{\langle ij\rangle} s_i s_j + h \sum_i s_i$ you would calculate moments using exactly this form. If you include higher order interactions, like a next-nearest-neighbor term $J_2 \sum_{i,j \in {\rm n.n.n.}} s_i s_j$ then you will get a model with a different numerical set of means, correlations, etc. The next-nearest-neighbor terms being irrelevant at the ferromagnet phase transition is a statement about the interactions under rescaling of the model (the second part of the RG transformation, after coarse-graining), and shows that the ferromagnetic transition of the next-nearest-neighbor model is in the same universality class as the nearest-neighbor model, but this does not come into play in calculating the moments. For the statistical moments you are not rescaling anything, nor do you even need to be anywhere near a fixed point of the rescaled model in general, and so any terms in your Hamiltonian/action contribute to those moments. Thus, adding more interaction terms to your $F(\phi)$ would change the statistical moments of the model, though they may not alter the fixed point structure.
Edit: To elaborate a little bit to answer the follow-up question in the comments:
About the 2nd point: if we compute the statistical moments with the
original nonlinear terms as you mentioned, does this computation
correspond to a specific scale, or its applicability can be extended
to different (presumably larger) scales? I think I am still confused
about the relationship between computing averages in a theory with
some given nonlinearities, with the effect of RG on those nonlinear
terms which, possibly, will generate new nonlinearities.
Let's keep going with the Ising model as an example. Suppose you had an experimental system that were exactly described by a classical nearest neighbor 2d Ising model on a triangular lattice with lattice spacing $a$ and spin-spin coupling $J$ at temperature $T$. If you have an experimental probe with sufficient resolution that you can measure the orientation of individual spins separated only by a lattice spacing (equivalently momenta all the way from $0$ to $\Lambda = \pi/a$), then you have access to every scale of the problem. If you wanted to predict, e.g., the covariance between two spins located at lattice positions $\mathbf{x}_m$ and $\mathbf{x}_n$ then you would compute
$$\langle s_m s_n \rangle = Z^{-1} \sum_{\mathbf{s}} s_i s_j \exp\left(\frac{J}{k_B T} \sum_{\langle ij \rangle} s_i s_j \right).$$
There's no need to renormalize anything, no need to include any higher-order interactions here: the model is exact over all spatial scales (by assumption).
Now, let's suppose that your measurement device cannot resolve the orientation of spins all the way down to separations of $a$, but only down to separations of $R$, such that any spins separated by less than this spacing appear blurred together (essentially, coarse-grained together) into some local magnetization; let the local magnetization at location $\mathbf{x}$ be $M_R(\mathbf{x})$, where I have denoted dependence the resolution scale $R$. The underlying individual spins are still described by the nearest neighbor Ising model, but since the coarse-grained $M_R(\mathbf{x})$ values are the only thing you can measure, you might prefer to have a model for it. In principle, the way you would formally construct such a model would be to calculate
$$\exp(-\beta L[M_R(\mathbf{x}]) = Z^{-1}\sum_{\mathbf{s}} e^{-\beta \mathcal H(\mathbf{s})} \delta\left[\sum_{\mathbf{x}_m \in \mathbf{x} + R} s_m - N_R(\mathbf{x}) M_R(\mathbf{x}) \right],$$
where $\delta$ is a delta function enforcing its argument, $\mathbf{x}_m \in \mathbf{x} + R$ restricts the sum over spins whose position are within the radius $R$ of the position $\mathbf{x}$, $N_R(\mathbf{x})$ counts the number of spins within that radius, and $L[M_R(\mathbf{x})]$ is an effective Hamiltonian for the new coarse degrees of freedom. Why is it labeled $L$? Well, essentially this would be a formal derivation of the Landau free energy; this argument follows Ch. 5 in Goldenfeld's book referenced above. You may know that one does not really construct the Landau free energy this way in practice, but instead appeals to symmetry considerations to work out the terms that would be present in $L$. If we were to derive the Landau free energy this way, however, this is where the effective interactions come into play: those terms that we would include on symmetry considerations would be the terms generated by coarse-graining the Ising model as part of a renormalization group procedure. Now, typically when using the Landau free energy one is just interested in understanding phase transitions, so not all of those higher order interactions generated by this coarse graining procedure are relevant: we know that if we rescale the coarse-grained variables by appropriate powers of the coarse-graining length $R$ to put the model back onto the original scale then many of those higher order terms are unimportant as far as the critical properties of the model are concerned. However, if instead what you wanted to do is accurately calculate things like correlation functions of these coarse-grained spins, then you need those higher order terms if you want to be quantitatively accurate. This is because the following two calculations must agree:
$$\langle M_R(\mathbf{x}) M_R(\mathbf{y}) \rangle = \int \mathcal D M ~M_R(\mathbf{x}) M_R(\mathbf{y}) e^{-\beta L[M]},$$
where I've formally written this as a path integral, but one does not need to assume a continuum limit, and
$$\langle M_R(\mathbf{x}) M_R(\mathbf{y}) \rangle = \frac{1}{N_R(\mathbf{x}) N_R(\mathbf{y})}\sum_{m \in \mathbf{x} + R}\sum_{n \in \mathbf{y} + R} \left[\frac{1}{Z}\sum_{\mathbf{s}} s_m s_n e^{-\beta \mathcal{H}(\mathbf{s})} \right],$$
where I inserted the definition of $M$ in terms of the microscopic spins. If you drop the 'irrelevant' terms from $L$ these two calculations will not agree. As in the answer to another of your questions, the nonlinearities generated are necessary for these two calculations to agree. Dropping them won't change the critical properties of the model (around the critical point corresponding to the ferromagnetic transition), but it would affect the quantitative calculation of the statistical moments/cumulants. If you working with the rescaled variables, $M'_R(\mathbf{x}') \equiv R^{-[M]}M_R(\mathbf{x}/R)$, where $[M] = (2-d)/2$ is the scaling dimension of $M_R$, and you want to calculate, e.g., $\langle M'_R(\mathbf{x}') M'_R(\mathbf{y}')\rangle$, then those higher-order terms in $L$ will be suppressed by powers of $R^{-1}$, and can be neglected without affecting the value of $\langle M'_R(\mathbf{x}') M'_R(\mathbf{y}')\rangle$ too much (because those terms will affect the values of the moments on scales $|\mathbf{x}| \lesssim R$.
I think perhaps where some of the confusion may be coming from the situation in high energy physics, in which there are additional complications arising from the fact that one does not actually know the microscopic model, and so there are additional steps involved in safely taking the continuum limit. In this case, I believe the variables one is working with are more akin to the rescaled variables $M'_R(\mathbf{x}')$, though someone more familiar with the high-energy situation should chime in here.
3) The general rule of thumb is that critical exponents and scaling functions tend to be universal, as well as certain ratios of quantities sometimes, but I don't know of a comprehensive list (to the extent that one could even be compiled). Moreover, even critical exponents are not guaranteed to be universal: Baxter's "Exactly solved models in statistical mechanics" includes an example of the square lattice eight-vertex model (Ch. 10) in which most of the critical exponents vary continuously with the parameters of the model. (However, Baxter notes that others have suggested this might be an artifact of defining the critical exponents as powers of $T-T_c$, and that if one defines exponents as powers of the inverse correlation length $\xi^{-1}$, which essentially divides the normal exponents by the correlation length exponent $\nu$, then these modified exponents do appear to be universal). Another example in which microscopic details show up in critical exponents are networks: in models with phase transitions on networks with degree distributions decaying as $p(n) \sim n^{-\gamma}$ at large degree $n$, with $2 < \gamma < 3$ (such that the mean degree exists but the degree variance is infinite), then the critical exponents are functions of $\gamma$. See, e.g., this paper's Eq. 8 for one example of this result. (Although, in this case $\gamma$ is playing a role similar to the dimension of a lattice, so the exponents are in a sense still 'universal')
Edit: I struck out the last example about networks because thinking about it more the final parenthetical comment interpreting the degree distribution exponent as being similar to dimension is the better way to think about the situation, so it is not really an example of microscopic details affecting the critical exponents. A better example would be this paper by Duclut and Delamotte investigating erosion over intermediate spatial scales. They find that there is a line of fixed points in the renormalization group flow, such that the value of the anomalous exponent is somewhat dependent on the initial conditions of the RG flow (i.e., on details of the microscopic action).
