Variance of the sum of Bernoulli Random variables? $\newcommand{\var}{\operatorname{var}}$
Let $X_{i}$ be a Bernoulli random variable with paramater $p_{i}$ where $p_i$ itself is a random variable that ranges from $0$ to $1$. The expectation of $p_{i}$ is $\rho$, and $X_{i}$ is independent of $X_{j}$ where $i\neq j$
Let $Y=\sum_{i=1}^{n} X_{i}$
Show that $\var(Y)>n\rho(1-\rho)$
This is what I have done so far:
$$\var(Y)=\var\left(\sum_{i=1}^n X_{i}\right) = \sum_{i=1}^n \var(X_{i})$$
$$\var(X_{i})=E[\var(X_{i}\mid p_{i})]+\var(E[X_i \mid p_i])$$
$$=E[p_i-p_i^2]+\var(p_i)$$
$$=\rho-E[p_i^2]+E[p_i^2]-E[p_i]^2$$
$$=\rho-\rho^2$$
So,
$$\var(Y)=\sum_{i=1}^n (\rho-\rho^{2})=n\rho(1-\rho)$$
What am I doing wrong? This is the variance if $p_i$ were a constant. When $p_i$ varies, surely that should increase the variance of $Y$. Any insight would be appreciated!
 A: Apparently, this insight is not correct - I doubt if what the question asks can be proved at all:

When $p_i$ varies, surely that should increase the variance of $Y$.

Intuitively, for some values of $p_i$, $X_i$ has lower variance than a Bernoulli with parameter $\rho$, and for some other values it has a higher variance. But the effects add up in such a way as if $p_i$ is fixed to $\rho$. 
Here is an alternate proof: 
Since $X_i$ takes only values $0$ and $1$, we have $X_i^2=X_i$ and $E[X_i^2] = E[X_i]$. Therefore, 
$$\operatorname{Var}(X_i) = E[X_i^2]-(E[X_i])^2 = E[X_i](1-E[X_i])$$ 
while by the law of total expectation,
$$E[X_i] = E[E[X_i|p_i]] = E[p_i] = \rho$$
Hence, $$\operatorname{Var}(X_i) =\rho(1-\rho)$$
A: $\newcommand{\var}{\operatorname{var}}$
$\newcommand{\cov}{\operatorname{cov}}$
The original poster has put this link to the question into a comment below the question.  It differs in crucial ways from the question posted above (in its present form).
It says:

Suppose $P(Y_i=1)=1-P(Y_i=0) = \pi$, $i=1,\ldots,n$, where $\{Y_i\}$ are independent.  Let $Y=\sum_i Y_i$.

Then in part (c), which the poster says in comments is what the question is about, it says:

Suppose that heterogeneity exists: $P(Y_i=1\mid\pi)=\pi$ for all $i$, but $\pi$ is a random variable with density function $g(\cdot)$ on $[0,1]$ having mean $\rho$ and positive variance.  Show that $\var(Y)>n\rho(1-\rho)$.

One crucial difference between this and the posted question is that there are not independent $\pi_i$ for $i=1,\ldots,n$, but just one $\pi$ that applies to every $Y_i$.
Thus $Y_i$, $i=1,\ldots,n$ are conditionally independent given $\pi$, but are not marginally independent.  For example, if one observes $Y_1=\cdots=Y_{10}=1$, then given that observation, it is far more probable that $\pi$ is close to $1$ than it was before that observation; hence $P(Y_{11}=1\mid Y_1=\cdots=Y_{10}=1)$ is much bigger than $P(X_{11}=1)$.
Now
\begin{align}
\var(Y) & = \sum_i \var(Y_i) + \sum_{i,j\,:\,i< j} 2\cov(Y_i,Y_j) \\[8pt]
& = n\rho(1-\rho) + \sum_{i,j\,:\,i< j} 2\cov(Y_i,Y_j).
\end{align}
The first term in the last line above follows from what the original poster wrote.
After that it is enough to show that $\cov(Y_i,Y_j)>0$ for $i\ne j$.  Notice that
$$
\begin{align}
\cov(Y_i,Y_j) & = \mathbb E(Y_i Y_j) - (\mathbb E Y_i)(\mathbb E(Y_j) \\
& = P(Y_i = 1 = Y_j) - \pi^2.
\end{align}
$$
So we only need to show that $P(Y_i=1=Y_j)>\pi^2$. We have
$$
P(Y_i=1=Y_j)= P(Y_i=1)P(Y_j=1\mid Y_i=1).
$$
So it's enough to prove that $P(Y_j=1\mid Y_i=1)>\pi=P(Y_j=1)$.
For now I'll leave the rest as an exercise.  Maybe I'll say more later, especially if comments are posted.
