Complementing Terry's nice response, let me try to explain from a more elementary point of view why Selberg's formula (1) is natural, and why it is true.
It is natural to formulate the PNT in terms of the von Mangoldt function $\Lambda(n)$, because it is supported on the prime powers, where it is given by the simple analytic formula $\Lambda(p^r)=\log p$. In addition, $\Lambda(n)$ is intimately related to the Möbius function $\mu(n)$ and the Riemann-zeta function via the equivalent formulae $$ \Lambda=\mu\ast\log,\qquad 1\ast\Lambda=\log,\qquad \sum_{n=1}^\infty\frac{\Lambda(n)}{n^s}=-\frac{\zeta'}{\zeta}(s).$$ Selberg's insight was that a natural and slightly more accessible quantity arises if we convolve $\mu$ not with $\log$, but a power of $\log$: $$ \Lambda_k=\mu\ast\log^k,\qquad 1\ast\Lambda_k=\log^k,\qquad \sum_{n=1}^\infty\frac{\Lambda_k(n)}{n^s}=(-1)^k\frac{\zeta^{(k)}}{\zeta}(s).$$ It is a nice exercise that $\Lambda_k(n)$ is nonnegative, hence at most $\log^k(n)$, and it is supported on numbers with at most $k$ distinct prime factors. Now (1) is essentially the same as $$ \sum_{n\leq x}\Lambda_2(n)=2x\log x+O(x), $$ i.e. an analogue of the PNT for numbers with at most two distinct prime factors (which makes the sum smoother). This explains why (1) is so natural and useful to look at. It is another matter, and the ingenuity of Selberg and Erdős, that one can derive information from the distribution of numbers with at most two distinct prime factors for the distribution of prime powers.
In general, we have the following generalization of Selberg's formula (1): $$ \psi_k(x):=\sum_{n\leq x}\Lambda_k(n)=x P_k(\log x)+O(x),\qquad k\geq 2,$$ where $P_k(t)\in\mathbb{R}[t]$ is a polynomial of degree $k-1$ with leading coefficient $k$. Let me sketch a proof of this more general formula. First, using the definitions, $$ \psi_k(x)=\sum_{n\leq x}\mu(n)L_k(x/n),$$ where $L_k$ is the summatory function of $\log^k$: $$ L_k(x):=\sum_{n\leq x}\log^k(n).$$ A key point is that $L_k(x)$ is a smooth sum, in fact by simple calculus $$ L_k(x)=xR_k(\log x)+O(\log^k(x)), $$ where $R_k(t)\in\mathbb{R}[t]$ is a polynomial of degree $k$ with leading coefficient $1$. Now the idea is to compare $L_k$ with the summatory functions of the generalized divisor functions: $$ T_j(x):=\sum_{n\leq x}\tau_j(n),\qquad \tau_j:=1\ast\dots\ast 1. $$ These are easier to analyze than $\Psi_k(x)$$\psi_k(x)$, in fact one can prove in an elementary way that $$ T_j(x)=xS_j(\log x)+O(x^{1-1/j}), \qquad j\geq 1,$$ where $S_j(t)\in\mathbb{R}[t]$ is a polynomial of degree $j-1$ with leading coefficient $1/(j-1)!$. Here we used the convention that $\tau_1:=1$, so that $T_1(x)=\lfloor x\rfloor=x+O(1)$ as claimed. It follows that $L_k(x)$ can be well-approximated with a linear combination of the $T_j(x)$'s: $$ L_k(x)=k!T_{k+1}(x)+a_kT_k(x)+\dots+a_1T_1(x)+O(x^{1-1/k}).$$ Plugging this back into our formula for $\psi_k$, we get $$ \psi_k(x)=\sum_{n\leq x}\mu(n)\Bigl\{k!T_{k+1}(x/n)+a_kT_k(x/n)+\dots+a_1T_1(x/n)\Bigr\}+O(x).$$ However, for each $j\geq 2$, $$ \sum_{n\leq x}\mu(n)T_j(x/n)=\sum_{n\leq x}(\mu\ast\tau_j)(n)=\sum_{n\leq x}\tau_{j-1}(n)=T_{j-1}(x), $$ and also $$ \sum_{n\leq x}\mu(n)T_1(x/n)=\sum_{n\leq x}(\mu\ast 1)(n)=1, $$ hence in fact $$\begin{align*} \psi_k(x)&=k!T_k(x)+a_kT_{k-1}(x)+\dots +a_2 T_1(x)+a_1+O(x)\\ &=x\Bigl\{k!S_k(\log x)+a_kS_{k-1}(\log x)+\dots + a_2 S_1(\log x)\Bigr\} + O(x).\end{align*} $$ This proves the claimed approximation for $\psi_k(x)$ with the polynomial $$P_k(t):=k!S_k(t)+a_kS_{k-1}(t)+\dots+a_2 S_1(t).$$ This polynomial has real coefficients, degree $k-1$ and leading coefficient $k!/(k-1)!=k$, where we used that $k\geq 2$ throughout. The proof is complete.