See this good mathematical treatment of undersampling by Prof. Fowler at Binghamton University, using the comb function as you have done in the time domain, taking the z-transform of the resulting function, and then evaluating the z-transform on the unit circle to view the frequency spectrum.
IV-03 Transform View_2007.pdf
As for an intuitive understanding, I like to explain this by first understanding aliasing when going from the continuous time domain to the discrete time domain (analog-digital conversion), and then showing that decimation is just a re-sampling of the digital signal. The two converge to be the same when the former case is considered to be a discrete signal with an infinite sampling rate.
To review sampling, we note that sampling is multiplication in time of a signal with an impulse train.

And making note that the sampling process in time as a stream of impulses, results in a stream of impulses in frequency:

Multiplication in time is convolution in frequency. Convolution is simplified when dealing with impulses as the result is a simple shift.
I have personally preferred to make note that each impulse at location $\omega_n$ in the complex frequency domain is the transform of $e^{j\omega_nt}$, and knowing that $e^{j\phi}$ is just a phasor with magnitude 1 and angle $\phi$, therefore each impulse in the frequency domain is a spinning phasor in time at rate $\omega_n$. Knowing this, the concept of positive and negative frequencies now make complete sense (direction of phasor rotation) as is clear with Euler's identify decomposing a cosine wave into two complementary spinning phasors ($2cos\omega_n t = e^{+j\omega_n t}+e^{-j\omega_n t}$) as shown by our input spectrum in the plot below. With this view I make note that exponents add when you multiply, so instead of the possible mental complexity of shifting through convolution, you can alternatively simply add the frequencies for each pair of impulses considered - for example, consider the 20Hz impulse specifically in the sampling spectrum in the figure below, this in the time domain is a spinning phasor at 20Hz rate, when multiplied in time with the input spectrum (two spinning phasors at +3 Hz and -3 Hz) produces 20+3 = 23 Hz and 20-3 = 17 Hz in the output spectrum. (the result is the same: a shift, just the latter offers an alternate explanation). This visual view, and thinking, has explained for myself so many operations in DSP in a more intuitive way, so I want to emphasize that point.
Note the resulting digital spectrum in the example shown in the figure below, sampling a 3 Hz sine wave with a 20 Hz sampling clock $F_s$. Due to the shift property mentioned, the input spectrum is replicated at all the impulse locations in frequency of the sampling clock. Since the digital spectrum is repeated over the interval from $-F_s/2$ to $+F_s/2$ (or equivalently from 0 to $F_s$), we only need to define it over that interval. However with regards to multi-rate signal processing with decimation and interpolation operations, I find it very helpful in getting an intuitive feel to extend it in the digital domain as I have done here (basically converting it to an analog equivalent spectrum).

Now consider the same scenario, but this time sampling a 23 Hz Sine Wave with a 20 Hz clock; the same digital spectrum will result! It should also be clear that this will also happen at 43 Hz, 63 Hz, 83 Hz ... Specifically the following frequency bands in the input spectrum will all alias into our baseband bandwidth B: $nF_s±B$, where B is a frequency band of interest at baseband, and n = ...-3,-2,1,0,1,2,3...

Carrying this explanation to decimation, we make note that decimation is simply resampling a signal that is already digital. Just as we noted with the original digital spectrum that it uniquely extends from 0 to $F_s$, the new digital spectrum from a decimation by M will extend from 0 to $F_s/M$ while the spectrum itself remains unchanged on a unit Hz domain, giving the property that the spectrum stretches on a normalized digital frequency domain. Just as we identified with the original analog input spectrum frequency bands that would alias into the digital spectrum at $nF_s±B$, the new digital spectrum from a decimation by M will have frequency bands in the original digital spectrum that will alias into the final digital spectrum at $nF_s/M±B$ where n = 0 to M-1 (since $nF_s = \textrm{modulo}(n,M)F_s$.