Binomial-binomial is binomial? Let $X\mid Y\thicksim\operatorname{Binomial}(Y,p)$, and let $Y\thicksim\operatorname{Binomial}(n,q)$. I have seen it written that:

Claim: $X$ is marginally $\operatorname{Binomial}(n,pq)$.

There is a simple justification: first, $Y$ members of the population are chosen at random with probability $q$, and then each of these is selected with probability $p$, so that each member is randomly selected in the end with probability $pq$. That is satisfying.
However, my computation of the marginal variance of $X$ disagrees with the answer that should be true if the claim above is true. If the claim above is true, we should have $\operatorname{var} X = npq(1-pq)$. However, I'm getting $$\begin{aligned}
\operatorname{var} X
&= E[\operatorname{var}(X\mid Y)] + \operatorname{var}(E[X\mid Y])\\
&=E[Yp(1-p)] + \operatorname{var}(pY)\\
&=np(1-p) + p^2 n q (1-q).
\end{aligned}$$
This does not equal $npq(1-pq)$. To see this, plug in $p=0.5,q=0.5$, and in that case $npq(1-pq)=0.1875n$ and $np(1-p) + p^2 n q (1-q)=0.3125n$.
I'm wondering if someone can check my math, or point to any resource showing that $X$ is really marginally Bin(n,pq). Thanks.

As Ben points out, you've made an algebraic error and the result is correct. This process is called binomial thinning and, if you search for that expression, you'll find many mentions of it in the published literature. The process applies not just to binomial random variables, but also to multinomial, Poisson and negative binomial.
Suppose that we have binomial, Poisson or negative binomial random variables:

*

*$Y_1\sim {\rm Binomial}(n,q)$

*$Y_2\sim {\rm Poisson}(\lambda)$

*$Y_3\sim {\rm Negative\ Binomial}(\mu,\phi)$, i.e., with mean $\mu$ and variance $\mu+\phi\mu^2$
We can view each of these random variables as counting events from a random process.
Suppose now that the individual events are not all observed but are randomly intercepted so that on average $p$ of them get through and are observed while the others are lost.
In other words, we "thin out" the random processes by keeping each of the original events with probability $p$:

*

*$X_1|Y_1 \sim {\rm Binomial}(Y_1,p)$

*$X_2|Y_2 \sim {\rm Binomial}(Y_2,p)$

*$X_3|Y_3 \sim {\rm Binomial}(Y_3,p)$
The resulting "thinned" distributions have the following marginal distributions:

*

*$X_1 \sim {\rm Binomial}(n,pq)$

*$X_2 \sim {\rm Poisson}(p\lambda)$

*$X_3 \sim {\rm Negative\ Binomial}(p\mu,\phi)$
The effect is to scale down the expected value of the distribution by factor $p$ without otherwise changing the distributional form.
An example of binomial thinning from my own use is the thinCounts function of the edgeR package ( https://rdrr.io/bioc/edgeR/man/thinCounts.html ) which can be used to generate RNA-seq read counts for reduced sequencing depths.
