To start, note that the generation of new interactions has everything to do with coarse-graining, it only appears in RG calculations because coarse-graining is the first step of an RG procedure. The essential issue is that if one has a system comprising multiple interacting random variables and you want to describe the behavior of just a subset of those random variables, then that description of the subset must contain effective interactions between random variables that do not necessarily have any direct interaction. This is true whether you have a field theory, multivariate distribution, or even a bivariate distribution. I'll leave aside the issue of relevant and irrelevant interactions here, which come up in the second step of the RG procedure, rescaling.
I always like to go back to simpler models in statistical mechanics to think about these things. As a simple example, let's consider the classical Ising model and a one-step coarse graining in real space by integrating out every other spin. In the full Hamiltonian for the model in which all of the spins are observed, suppose each spin only interacts with its nearest neighbors. e.g., in $d=1$
$$\mathcal H(\mathbf{s}) = -J \sum_{i} s_i s_{i+1} - h s_i$$
with $P(\mathbf{s}) = \exp(-\beta \mathcal H(\mathbf{s}))/Z$ and $Z$ the normalization factor. The coarse-grained Hamiltonian for the 1d example would be defined, say, by summing over all of the even sites:
$$\exp(-\beta\mathcal H'(\mathbf{s})) = \sum_{s_{2n}} \exp(-\beta \mathcal H(\mathbf{s})).$$
Note that if we know the full Hamiltonian, we can in principle calculate any statistical moment of the distribution: mean orientation of each spin, spin-spin correlations, etc. Our coarse-grained Hamiltonian must give the exact same values for these statistics; e.g., $\langle s_1 s_5 \rangle$ must be the same value whether we compute it using the full Hamiltonian $\mathcal H(s)$ or the coarse-grained Hamiltonian $\mathcal H'(s)$.
We know it cannot possibly be the case that our coarse-grained Hamiltonian contains only independent spins, since then all of the higher order statistics like correlations would be $0$. It must be the case that the coarse-grained Hamiltonian has developed effective interactions between spins that did not directly interact in the fully resolved model. One might hope that these effective interactions only couple what were formerly next-nearest neighbor spins, but often that is not exactly the case: any interaction allowed by symmetry has developed, although the relative magnitudes of these interactions are not all the same---far away spins would have weaker interactions than nearby spins in this one-step coarse-grained model.
As an even simpler example, let's consider a zero-dimensional model containing just two random variables $x$ and $y$, both taking values on $\mathbb{R}$. Suppose the probability density describing the interactions between these two random variables is
$$\rho(x,y) \propto \exp\left(-x^2 - x^4 - y^4 - x^2y^2\right).$$
Some things to note about this distribution: it is even in both $x$ and $y$, and it does converge. Suppose we only care about the variable $x$ (say because it is easier to observe), so we would be just as happy with a model for only $x$. The proper way to obtain this model is to marginalize out $y$. In this case one can perform the quartic integral in terms of an infinite series:
$$\int_{-\infty}^\infty dy~e^{-x^2y^2 - y^4} = \frac{1}{2}\sum_{p=0}^\infty \frac{(-x^2)^{p}}{p!} \Gamma\left(\frac{2p+1}{4}\right) = e^{\frac{x^4}{8}} |x| K_{1/4}\left(\frac{x^4}{8} \right);$$
for completeness I've added the representation of the series in terms of the modified Bessel function $K_\alpha(x)$, but the series is more useful for our general purpose. Now let us write the marginalized distribution in the form
$$\rho(x) \propto \exp\left(-x^2 - x^4 + \ln\left[\frac{1}{2}\sum_{p=0}^\infty \frac{(-x^2)^{p}}{p!} \Gamma\left(\frac{2p+1}{4}\right) \right] \right).$$
Thus, we see that even in this $0d$ example marginalizing out one of the random variables has resulted in a marginal distribution for $x$ that has not only modified the coefficients of the existing terms (if we extract the $x^2$ and $x^4$ terms by expanding the logarithm in a formal series), but has also "generated" even powers of all higher orders---note that there are no odd powers, consistent with the symmetry of the full distribution. These higher order interactions must be there to ensure that the statistical moments $\langle x^m \rangle = \int dx~\rho(x) x^m $ match the moments obtained by calculating $\langle x^m \rangle = \int dx dy~\rho(x,y) x^m$. (You can check numerically that, e.g., the second moment evaluates to $\sim 0.21$ using both the full and marginalized distributions. For the latter, the representation of the series in terms of the Bessel function is useful in, e.g., Mathematica).