The answer by @ZeroTheHero is perfectly good and gives an intuition to the why. Namely $a_-a_+$ has a factor of $n$ because it is a piece of the Hamiltonian and the energy is linear in the number of excitations.
I will show an alternative approach which is purely algebraic as it depends on the commutator. It was proposed in a comment, I will just spell it out in detail.
We have
$$
[a_-, a_+] = 1\,.
$$
By induction we can show
$$
[a_-, a_+^m] = m a_+^{m-1}
$$
Let's work through the steps
- Base case: $m=1$ is the commutation relation
- Inductive case:
$$
[a_-,a_+^{m+1}] = [a_-,a_+^m] a_+ + a_+^m [a_-,a_+] = m\, a_+^m + a_+^m = (m+1)\,a_+^m\,.
$$
Now that we proved this identity we can build the Hilbert space of the harmonic oscillator by applying powers of $a_+$ to the ground state
$$
\psi_n = \frac{a_+^n}{\sqrt{n!}}\psi_0\,.
$$
From this we get
$$
\begin{aligned}
a_+ a_- \,\psi_n &= \frac{a_+ a_- a_+^n}{\sqrt{n!}}\psi_0
\\ &= \frac{a_+ (a_- a_+^n - a_+^n a_-)}{\sqrt{n!}}\psi_0
\\ &= \frac{a_+ [a_-, a_+^n]}{\sqrt{n!}}\psi_0
\\ &= \frac{n\, a_+^n}{\sqrt{n!}}\psi_0
\\ &= n \,\psi_n
\end{aligned}
$$
The second line is due to $a_-\psi_0 = 0$.
You might also wonder why I defined the state with a $\sqrt{n!}$. I did not have to, but it's a convention to make them unit normalized $|\psi_n|^2 = 1$. This can also be proved using the commutator rules
$$
\begin{aligned}
|\psi_n|^2 &= \langle 0 | \frac{a_-^n a_+^n}{n!} | 0\rangle
\\ &= \langle 0 | \frac{a_-^{n-1} [a_-, a_+^n]}{n!} | 0\rangle
\\ &= n\, \langle 0 | \frac{a_-^{n-1} a_+^{n-1}}{n!} | 0\rangle
\\ &= \langle 0 | \frac{a_-^{n-1} a_+^{n-1}}{(n-1)!} | 0\rangle
\\ &= |\psi_{n-1}|^2
\end{aligned}
$$
By induction this can be taken all the way down to $|\psi_0|^2$ which is unit normalized by definition.