We present two new connections between the inhomogeneous stochastic higher spin six vertex model in a quadrant and integrable stochastic systems from the Macdonald processes hierarchy. First, we show how Macdonald $q$-difference operators with $t=0$ (an algebraic tool crucial for studying the corresponding Macdonald processes) can be utilized to get $q$-moments of the height function $\mathfrak{h}$ in the higher spin six vertex model first computed in arXiv:1601.05770 using Bethe ansatz. This result in particular implies that for the vertex model with the step Bernoulli boundary condition, the value of $\mathfrak{h}$ at an arbitrary point $(N+1,T)\in\mathbb{Z}_{\ge2}\times\mathbb{Z}_{\ge1}$ has the same distribution as the last component $λ_N$ of a random partition under a specific $t=0$ Macdonald measure. On the other hand, it is known that $\mathbf{x}_N:=λ_N-N$ can be identified with the location of the $N$th particle in a certain discrete time $q$-TASEP started from the step initial configuration. The second construction we present is a coupling of this $q$-TASEP and the higher spin six vertex model (with the step Bernoulli boundary condition) along time-like paths providing an independent probabilistic explanation of the equality of $\mathfrak{h}(N+1,T)$ and $\mathbf{x}_N+N$ in distribution. Combined with the identification of averages of observables between the stochastic higher spin six vertex model and Schur measures (which are $t=q$ Macdonald measures) obtained recently in arXiv:1608.01553, this produces GUE Tracy--Widom asymptotics for a discrete time $q$-TASEP with the step initial configuration and special jump parameters.