One heuristic is to replace the $n^{th}$ roots of unity by $n$ iid elements $\zeta_1,\dots,\zeta_n$ of the unit circle, drawn uniformly at random. For any sum $\zeta_{i_1} + \dots + \zeta_{i_k}$ of $k$ of these with $i_1 < \dots < i_k$, a standard Fourier analytic calculation (Esseen concentration inequality) shows that if $k \geq 5$ (to make the $k$-fold convolution power of the Fourier transform of the unit circle have absolutely integrable Fourier transform), then $$ {\bf P}( |\zeta_{i_1} + \dots + \zeta_{i_k}| < r ) \ll_k r^2 $$ for $0 < r < 1$. Taking the union bound over all $\binom{n}{k}$ different sums, we see that the probability that $k$ distinct elements of the $\zeta_1,\dots,\zeta_n$ sum to something of magnitude less than $r$ is $O_k( \binom{n}{k} r^2 )$. This shows that with probability $\gg 1$, we should have the lower bound $$ |\zeta_{i_1} + \dots + \zeta_{i_k}| \gg_k \binom{n}{k}^{-1/2} \asymp_k n^{-k/2}$$ for all $i_1 < \dots < i_k$. This is for a fixed $n$ (which after averaging corresponds to predicting a lower bound $\gg_k n^{-k/2}$ for a positive density fraction of $n$, rather than for all $n$). To get a bound that is true almost surely for all but finitely many $n$ (in the spirit of the "strong" law of large numbers, as opposed to the "weak" law), the Borel-Cantelli lemma suggests that we need to get the failure probability for a fixed $n$ down to something like $1/n^{1+\varepsilon}$, and so the lower bound now worsens from $\gg_k n^{-k/2}$ to $\gg_{k} n^{-(k+1+\varepsilon)/2}$.
These bounds are fairly reversible in this model due to all the independence, and with a little more effort one should be able to show that the bounds of $\gg_k n^{-k/2}$ (for a positive density set of $n$) and $\gg_k n^{-(k+1+\varepsilon)/2}$ (for all but finitely many $n$) cannot be significantly improved except possibly for the $\varepsilon$ factor, but I have not attempted to work this out rigorously.
To complete the heuristic analysis one should also treat the cases where there are collisions amongst the $i_1,\dots,i_k$, which can be handled by the same methods as long as there are at least $5$ distinct values of the $i_1,\dots,i_k$; the cases of only four or fewer values need to be treated by more direct computations but I think they should give lower order contributions in the $k \geq 5$ regime.
EDIT: while the analysis of the above model remains accurate within the scope of that model, I have just realised that the rotation invariance symmetry of the $n^{th}$ roots of unity leads to a larger lower bound prediction than the above naive model. The key point is that for the original $n^{th}$ roots of unity problem, one can always rotate one of the roots to equal $1$. So one should really be looking at sums of the form $1 + \zeta_{i_2} + \dots + \zeta_{i_k}$ with $1 < i_2 < \dots < i_k$. This effectively lowers $\binom{n}{k}$ to $\binom{n-1}{k-1}$ and raises all the exponents in the previous analysis by $1/2$, thus the revised prediction would be a lower bound of $\gg_k n^{-(k-1)/2}$ for a positive fraction of $n$ and $\gg_k n^{-(k+\varepsilon)/2}$ for all but finitely many $n$. There are some additional symmetries coming from cyclic permutation of the roots and conjugation symmetry but these only affect the $k$-dependent implied constants and not the exponents. In principle, the further Galois group symmetries of the problem could lead to additional refined corrections to the model, but I tentatively am of the opinion that this will not change the predicted exponents further.