This is a famously difficult problem in combinatorics. The following is the outline of an algebraic proof by Richard Stanley, using the theory from his textbook Algebraic Combinatorics.
The proof requires the theory of graded posets. A poset is a set equipped with a partial order. If $x,y$ are poset elements, we say $x$ covers $y$ if $x> y$, and there is no poset element $z$ with $x>z>y$. Then, a poset is called graded if there exists a rank function defined on the poset, such whenever $x$ covers $y$, the rank of $x$ is one more than the rank of $y$. We also usually normalize things so all minimal elements of the poset have rank zero. Finally, a poset is called rank-unimodal if the sequence of sizes of the ranks is unimodal, increasing up to a widest point in the poset's middle and decreasing thereafter.
A common example of a poset is $B_n$, the Boolean poset of the set of subsets of an $n$ element set. Stanley's big theorem concerns a certain quotient poset of $B_n$. Given a subgroup $G$ of $S_n$, the symmetric group on the same $n$-element set, we get an induced action on $B_n$ in a natural way. We then turn the set of orbits of this action into a poset, by saying for all $S,T\in B_n$ that $\text{orbit}(S)\le \text{orbit}(T)$ whenever there exists permutations $g_1,g_2\in G$ such that $g_1(S)\subseteq g_2(T)$. This poset on orbits is denoted $B_n/G$.
Theorem: (Stanley) Let $G$ be a subgroup of $S_n$. Then the quotient poset $B_n/G$ is rank unimodal.
I omit the proof, but you can read it in Stanley's book hosted by Stanley for free online, Theorem 5.8.
Therefore, in order to prove that the coefficient sequence of $\binom{n}{k}_q$ is rank unimodal, it suffices to find a poset where the size of the $i^{th}$ rank is the coefficnet of $q^i$ in $\binom{n}k_q$, for each $0\le i\le k(n-k)$, and prove that this poset can be realized as the quotient of a Boolean poset by a group action.
For positive integers $m,n$, let $L(m,n)$ be the poset of integer partitions with at most $n$ parts, where every part is at most $m$, where $(\lambda_1,\dots,\lambda_n)\le (\mu_1,\dots,\mu_n)$ if $\lambda_i\le \mu_i$ for all $i\in \{1,\dots,n\}$. Equivalently, these are lattice paths in an $m\times n$ rectangle, ordered by subset inclusion of the regions beneath the path. It is well known that the rank of path is the area of the region beneath, and that the numbers of paths with each rank exactly match the coefficients of $\binom{m+n}{m}_q$.
To realize this as a quotient of a Boolean poset, we will use $B_{m\times n}$, where we think of this as subsets of the $m\times n$ boxes forming the lattice rectangle spanned by $(0,0)$ and $(m,n)$. The permutation group $G$ acting on this array of boxes is as follows: $G$ is the set of permutations given by arbitrarily permuting within each row, and then permuting the rows themselves. It can be shown for any subset $S$ of the grid, the orbit of $S$ under this action will contain exactly one subset where every row is left-justified, and every columns is bottom-justified, which exactly corresponds to the region below some lattice paths. In this way, we can see that $L(m,n)$ is just the quotient of $B_{m\times n}$ by this action. Applying Stanley's theorem, we conclude the coefficents of $\binom{m+n}{m}_q$ are unimodal.
Over the years, other proofs have been given.
Igor Pak and Greta Panova gave a proof using Schur functions: https://arxiv.org/abs/1306.5085
Kathy O'hara gave a combinatorial proof. Here is an exposition of her proof by Doron Zeilberger: https://www.jstor.org/stable/2325177