Bernoulli distribution with random means Let $S = \frac{X_1 + \cdots + X_n}{n}$ where the $X_i$ are IID Bernoulli distributed with mean $p$, then $E[S] = p$ and $Var(S) = \frac{p(1-p)}{n}$.
Now consider the slightly more complex setup where $S'=\frac{1}{nk}(\sum_{1 \leq i \leq k} \sum_{1\leq j \leq n} X_{ij})$ where $\{X_{ij} | 1 \leq j \leq n \}$ are IID Bernoulli samples with a mean of $p_i \forall 1 \leq i \leq k$.
The $p_i$ are IID sampled from the uniform distribution on $[0,1]$.
How do I go about calculating mean and Var of $S'$? How do I formalize this?

Since all variables are IID, this is quite straightforward.  First we compute the conditional moments:
$$\begin{align}
\mathbb{E}(S' | \mathbf{p}) 
&= \mathbb{E} \bigg( \frac{1}{nk} \sum_{i=1}^k \sum_{j=1}^n X_{ij} \bigg|\mathbf{p}\bigg) \\[6pt]
&= \frac{1}{nk} \sum_{i=1}^k \sum_{j=1}^n \mathbb{E}(X_{ij}|p_i) \\[6pt]
&= \frac{1}{nk} \sum_{i=1}^k \sum_{j=1}^n p_i \\[6pt]
&= \frac{1}{k} \sum_{i=1}^k p_i, \\[12pt]
\mathbb{V}(S' | \mathbf{p}) 
&= \mathbb{V} \bigg( \frac{1}{nk} \sum_{i=1}^k \sum_{j=1}^n X_{ij} \bigg| \mathbf{p}\bigg) \quad \quad \quad \quad \\[6pt]
&= \frac{1}{n^2 k^2} \sum_{i=1}^k \sum_{j=1}^n \mathbb{V}(X_{ij}|p_i) \\[6pt]
&= \frac{1}{n^2 k^2} \sum_{i=1}^k \sum_{j=1}^n p_i (1-p_i) \\[6pt]
&= \frac{1}{n k^2} \sum_{i=1}^k p_i (1-p_i). \\[12pt]
\end{align}$$
Then we use the law of iterated expectation and the law of iterated variance to compute the unconditional moments:
$$\begin{align}
\mathbb{E}(S') 
&= \mathbb{E}(\mathbb{E}(S' | \mathbf{p})) \\[6pt]
&= \mathbb{E} \bigg( \frac{1}{k} \sum_{i=1}^k p_i \bigg) \\[6pt]
&= \frac{1}{k} \sum_{i=1}^k \mathbb{E}(p_i)) \\[6pt]
&= \frac{1}{k} \sum_{i=1}^k \frac{1}{2} \\[6pt]
&= \frac{1}{2}, \\[12pt]
\quad \quad \quad \quad \quad
\mathbb{V}(S') 
&= \mathbb{E}(\mathbb{V}(S' | \mathbf{p})) + \mathbb{V}(\mathbb{E}(S' | \mathbf{p})) \\[6pt]
&= \mathbb{E} \bigg( \frac{1}{n k^2} \sum_{i=1}^k p_i (1-p_i) \bigg) + \mathbb{V} \bigg( \frac{1}{k} \sum_{i=1}^k p_i \bigg) \\[6pt]
&= \frac{1}{n k^2} \sum_{i=1}^k \mathbb{E}(p_i (1-p_i)) + \frac{1}{k^2} \sum_{i=1}^k \mathbb{V}(p_i) \\[6pt]
&= \frac{1}{n k^2} \sum_{i=1}^k \bigg( \frac{1}{2} - \frac{1}{3} \bigg) + \frac{1}{k^2} \sum_{i=1}^k \frac{1}{12} \\[6pt]
&= \frac{1}{n k^2} \sum_{i=1}^k \frac{1}{6} + \frac{1}{k^2} \sum_{i=1}^k \frac{1}{12} \\[6pt]
&= \frac{1}{6 n k} + \frac{1}{12 k} \\[6pt]
&= \frac{1}{12 k} \cdot \frac{n+2}{n}. \\[6pt]
\end{align}$$
As a sanity check, we observe that $\mathbb{V}(S') \rightarrow 0$ as $\min(k,n) \rightarrow \infty$, which is what we would expect from the law of large numbers.
