Here is a more precise heuristic argument for the probability that $n$ is prime. Let $f(k, n)$ represent our subjective probability that $n$ is $k$-rough, i.e. has no divisors smaller than $k$. We are interested in $f(n, n)$.

We have: $$f(k+1, n) = f(k, n) (1 - \mathbb{P}(\text{k is prime}) \mathbb{P}(\text{$k|n$ given that $k$ is prime and $n$ is $k$-rough}))).$$

From here on out we’ll assume $k$ is prime and write $p$ instead. Your “bad heuristic argument” treats $p|n$ and “$n$ is $p$-rough” as independent. But as you point out they aren’t independent and we can get a better estimate by taking the correlation into account.

To do this we apply Bayes’ rule:
$$\mathbb{P}(\text{$p|n$ given $n$ is $p$-rough}) = \mathbb{P}(p|n)\frac{\mathbb{P}(\text{$n$ is $p$-rough given $p|n$})}{\mathbb{P}(\text{$n$ is $p$-rough})}$$

Note that if $p|n$, then $n$ is $p$-rough if and only if $n/p$ is $p$-rough. So $\mathbb{P}(\text{$n$ is $p$-rough given $p|n$})$ is just $f(p, n/p)$. The denominator $\mathbb{P}(\text{$n$ is $p$-rough})$ is $f(p, n)$. Thus:
$$\mathbb{P}(\text{$p|n$ given $n$ is $p$-rough}) = \frac 1p \frac{f(p, n/p)}{f(p, n)}.$$
This differs from your naive estimate by a factor of $f(p, n/p) / f(p, n)$, reflecting the correlation between different divisibility events. The correction factor is is 1 if $p$ is much smaller than $n$, indicating approximate independence. It’s 0 if $p > \sqrt{n}$, reflecting the common crude approximation where we stop looking for prime factors at $\sqrt{n}$. But you get a softer interpolation between those regimes.

Including this correction factor, we get the recurrence $$f(k+1, n) - f(k, n) = - \mathbb{P}(\text{$k$ is prime}) \frac {f(k, n/k)}{k}.$$

If we don’t assume the prime number theorem is true up to $n$, and just substitute $\mathbb{P}(\text{$k$ is prime}) = f(k, k)$, then it turns out that any solution of this equation satisfies $f(k, k) = (1+o(1))/\log(k)$.

But that doesn’t really answer your question or show that we’ve fixed the bogus independence assumption. The reason it works out is that the resulting differential equation is self-correcting---if the density of primes is more than $1/\log(k)$ for a while then it exerts downwards pressure on $\mathbb{P}(\text{$n$ is prime})$ and vice versa. So even if we make some bogus independence assumptions leading to $O(1)$ multiplicative error, we will still get the correct $\frac {1 + o(1)}{\log n}$ asymptotics.

To answer your question, and show that we’ve actually fixed the problem with the independence assumption in the naive heuristic estimate, let’s instead use the correct distribution for the primes up to $n$ and show that our recurrence gives the correct estimate for $\mathbb{P}(\text{$n$ is prime})$. That is, let’s assume that for $k < n$ we have:

- $\mathbb{P}(\text{$k$ is prime}) = 1/\log(k).$
- $\prod_{p < k}(1 - 1/p) = \exp(-\gamma)/\log(k).$

We can exactly work out $f(k, n)$ for small values of $k$. In this regime, the events $p|n$ are exactly independent, and so $f(k, n) = \prod_{p < k}\left(1 - 1/p\right) = \exp(-\gamma)/\log(k)$.

For large values of $k$ we can approximate $f(k+1, n) - f(k, n) \approx \frac {d}{dk} f(k, n)$.

So now we have a differential equation: $\frac {d}{dk} f(k, n) = -\frac{f(k, n/k)}{k \log k}$, with boundary condition $f(k, n) = \frac{\exp(-\gamma)}{\log k}$ for $ k \ll n$.

I claim the solution is exactly $$f(k, n) = \frac {\omega\left(\frac {\log n}{\log k}\right)}{\log k}$$for $ k \leq n$, and $0$ for $k > n$, where $\omega$ is Buchstab's function. Plugging in $k = n$, and using $\omega(1) = 1$, we get the desired $\mathbb{P}(\text{$n$ is prime}) = \frac 1{\log n}$.

To see that this is a solution:

- If $k \ll n$, then we have $f(k, n) = \frac {\exp(–\gamma)}{\log k}$, since $\omega(u) \rightarrow \exp(-\gamma)$ for large $u$.
- If $k > \sqrt{n}$, then $f(k, n) = \frac {\omega\left(\frac {\log n}{\log k}\right)}{\log k} = \frac {1}{\log n}$, because $\omega(u) = \frac 1u$ for $u < 2$. Thus $\frac{d}{dk} f(k, n) = 0$. This is correct, because $k > k/n$ and hence $f(k, n/k) = 0$.
- If $k$ is large but $k < \sqrt{n}$, then we have $\frac {d}{dk} f(k, n) = \frac {d}{dk} \frac{\omega\left(\frac {\log n}{\log k}\right) \frac {\log n}{\log k}}{\log n}$. Now we can use the fact that the Buchstab function $\omega$ is defined so that $\frac {d}{du} \omega(u) u = \omega(u - 1)$. Applying this with $u = \frac {\log n}{\log k}$, and using the fact that $\frac d{dk} \frac {\log n}{\log k} = -\frac {\log n}{k \log^2 k}$, we have: $$\frac {d}{dk} f(k, n) = -\frac {\omega\left(\frac{\log n}{\log k}-1\right)\frac {\log n}{k \log^2 k}}{\log n} = -\frac {\omega\left(\frac {\log n/k}{\log k}\right)}{k \log^2 k} = - \frac{f(k, n/k)}{k \log k},$$ as desired.

To summarize: if we take into account the correlation of divisibility events using Bayes’ rule in the most naive possible way, we end up with a correction term like $f(k, n/k) / f(k, n)$. Solving the resulting equation gives you the the correct heuristic estimate for the probability that $n$ is $k$-rough, and in particular the probability that $n$ is prime. The missing factor of $\exp(\gamma)$ you are looking for comes from the limit of the Buchstab function, which appears as the solution of the differential equation corresponding to this correction factor.

At first this felt to me like an unexplained coincidence or a bit of magical reverse causality. I've effectively just argued heuristically that for small $k$, $\mathbb{P}(\text{$n$ is prime}) = \frac {\exp(\gamma) \log k}{\log n}\prod_{p < k} (1 - 1/p)$. But why does $\prod_{p < k} (1 - 1/p)$ happen to have an exactly compensating factor of $\exp(-\gamma)$?

This comes directly from the self-correcting nature of the differential equation we mentioned before. As long as $\prod_{p < k} (1 - 1/p) > \exp(-\gamma)/\log(k)$, then $\mathbb{P}(\text{$n$ is prime}) > 1/\log(n)$ and so $\prod_{p < k}(1 - 1/p)$ will shrink faster than $1/\log(k)$. And similarly as long as $\prod_{p < k}(1 - 1/p)< \exp(-\gamma)/\log(k)$ then it will shrink slower than $1/\log(k)$. So I intuitively view the $\exp(\gamma)$ correction from the correlation as being the "logically prior" fact, with the $\exp(-\gamma)$ factor in $\prod_{p < k}(1 - 1/p)$ arising from the control system compensating for that correction.