First of all, I assume you understand that this is meant to be a nonrigorous argument, so there will be a limit to how rigorous I can make my answer.
The intuition here is that $n$ is prime if and only if it is not divisible by any prime $<n$. So we "should" have
$$f(n) \approx \prod_{p < n} \left( 1-1/p \right).$$
Similarly
$$f(n+1) \approx \prod_{p<n+1} \left( 1-1/p \right) = \prod_{p < n} \left( 1-1/p \right) \cdot \left\{ \begin{matrix} \left( 1-1/n \right) \ \mbox{if}\ n\ \mbox{is prime} \\ 1 \ \mbox{if}\ n\ \mbox{is not prime} \end{matrix} \right.$$
$$\approx f(n) \cdot \left\{ \begin{matrix} \left( 1-1/n \right) \ \mbox{if}\ n\ \mbox{is prime} \\ 1 \ \mbox{if}\ n\ \mbox{is not prime} \end{matrix} \right. .$$
Since $n$ is prime with "probability $f(n)$", we interpolate between the two cases next to the brace by writing:
$$f(n+1) \approx f(n) \left( 1-f(n)/n \right).$$
One might argue that it would be better to interpolate with a factor of $(1-1/n)^{f(n)}$, but this will make no difference in the asymptopics as $(1-1/n)^{f(n)} = 1-f(n)/n+O(1/n^2)$.
This argument is famously fishy, because it gives the right answer, but the intermediate point is wrong! The actual asymptopics of $\prod_{p<n} \left( 1-1/p \right)$ do not look like $1/\log n$, but like $e^{-\gamma} /\log n$. I've never seen a good intuitive explanation for why we get the wrong estimate for $\prod_{p<n} \left( 1-1/p \right)$, but the right estimate for the density of the primes.
There are several ideas here, some mentioned in the other answers:
One: When Gauss was a boy (by the dates found on his notes he was approximately 16) he noticed that the primes appear with density $ \frac{1}{\log x}$ around $x$. Then, instead of counting primes and looking at the function $\pi (x)$, lets weight by the natural density and look at $\sum_{p\leq x} \log p$. Since we are weighting by what we think is the density, we expect it to be asymptotic to be $x$.
Two: Differentiation of Dirichlet series. If $$ A(s)=\sum_{n=1}^{\infty} a_{n} n^{-s} $$ then $$ (A(s))'=-\sum_{n=1}^{\infty} a_{n} log(n) n^{-s}$$
The $\log$ term appears naturally in the differentiation of Dirichlet series. Taking the convolution of $\log n$ with the Mobius function (that is multiplying by $\frac{1}{\zeta(s)}$) then gives the $\Lambda(n)$ mentioned above. The $\mu$ function is really the special thing here, not the logarithm.
Expanding on this, there are other weightings besides $\log p$ which arise naturally from taking derivatives. Instead we can look at $\zeta^{''}(s)=\sum_{n=1}^\infty (\log n)^2 n^{-s}$, and then multiply by $\frac{1}{\zeta(s)}$ as before. This leads us to examine the sum $$\sum_{n\leq x} (\mu*\log^2 )(n)$$ (The $*$ is Dirichlet convolution) By looking at the above sum, Selberg was able to prove his famous identity which was at the center of the first elementary proof of the prime number theorem:
$$\sum_{p \leq x} log^2 p +\sum_{pq\leq x}(\log p)(\log q) =2x\log x +O(x).$$
Three: The primes are intimately connected to the zeros of $\zeta(s)$, and contour integrals of $\frac{1}{\zeta(s)}$. (Notice it was featured everywhere here so far) We can actually prove that
$$\sum_{p^k \leq x} \log p= x - \sum_{\rho :\ \zeta(\rho)=0} \frac{x^{\rho}}{\rho} +\frac{\zeta'(0)}{\zeta(0)} $$
Notice that the above is an equality, which is remarkable since the left hand side is a step function. (Somehow, at prime powers all of the zeros of zeta conspire and make the function jump.)
Best Answer
You ask:
As pointed out by the users @wojowu and @PeterHumphries, it is true that the PNT is equivalent to
$$\lim_{x \to \infty} \sum_{n\leq x} \frac{\mu(n)}{n}=0,$$ and it is relatively easy to prove that
$$\lim_{s\rightarrow 1^+} \sum_{n=1}^{\infty} \frac{\mu(n)}{n^s}=0.$$ The real difficulty lies in proving that
$$\lim_{x\rightarrow \infty} \sum_{n\leq x} \frac{\mu(n)}{n}= \lim_{s\rightarrow 1^+} \sum_{n=1}^{\infty} \frac{\mu(n)}{n^s},$$ which is highly nontrivial and requires intricate arguments.
In particular, as pointed out by @TerryTao in the comments: