How to calculate $Var(X)$ for the union size? Let's say there is a set of $4$ bags $\{{a,b,c,d\}}$ containing balls of colors $\{{red,blue,green,orange,black\}}$. Balls are assigned to bags in an arrangement that allows a variable number of colors in each bag, and overlaps.
$a: \{{red, green\}}, b: \{{red, blue\}}, c: \{{orange\}}, d: \{{red, black\}}$
When sampling $k$ bags, without replacement, one counts the total union of unique colors. Two examples:

*

*If $k = 1$, then one can sample either $a$, $b$, $c$, or $d$, obtaining $\{{red, green\}}$, $\{{red, blue\}}$, $\{{orange\}}$ or $\{{red, black\}}$ as colors, respectively. Then the number of unique colors is $\{{red, blue\}}$ = $2$ if selecting $a$, or $\{{orange\}}$ = $1$ if selecting $b$.

*If $k = 2$, one can draw $a$ and $b$, obtaining $\{{red, blue\}}$ and $\{{red, blue\}}$ as colors. Then the union size is the number of unique colors, $\{{red, blue, orange\}}$ = $3$.


*

*What's $Var(X)$, when drawing $k$ bags from this set? Checking all options, one can calculate $E(X) = 1.75$ when $k = 1$, $E(X) = 3$ when $k = 2$ $E(X) = 4$ when $k = 3$, and $E(X) = 5$ when $k = 4$.

*Is there a formulation for the underlying distribution that one can use, inferred from $E(X)$ or derived from all provided information? I think this could be a variant of a hypergeometric distribution, with variable weights per member, and duplicates.

 A: The key to this solution is to represent the "bags" with a logical array (of zeros and ones) representing non inclusion of elements in the bags.  Each bag gets its own row and each element of the universal set  $S = \{\text{black}, \ldots, \text{red}\}$ gets its own column.  This matrix will represent the "is not an element of" relation.  That is, for a bag $i$ and element $j$, set $a_{ij}=1$ when $j\notin i$ and otherwise set $a_{ij}=0$ (when $j\in i$).  Here is the matrix $\mathbb{A}=(a_{ij})$ for the example:
$$\mathbb{A} = \begin{array}{l|ccccc}
\text{} & \text{black} & \text{blue} & \text{green} & \text{orange} & \text{red}\\
\hline
\text{a} & 1 & 1 & 0 & 1 & 0 \\
\text{b} & 1 & 0 & 1 & 1 & 0 \\
\text{c} & 1 & 1 & 1 & 0 & 1 \\
\text{d} & 0 & 1 & 1 & 1 & 1
\end{array}$$
The sample space is the collection of all subsets of the bags $\{a,b,c,d\}.$  (It is okay for two or more bags to be identical as sets: everything still works out.)
We are going to define some random variables with values in $\{0,1\}.$  To do this, pick a column $j$ in the matrix and let the value of $X_j$ for any collection of $\mathcal{B}=\{i_1,i_2,\ldots,i_k\}$ be the product of all the entries for those bags:
$$X_j(\mathcal{B}) = \prod_{i\in\mathcal{B}} a_{ij}.$$
For instance, for the collection of bags $\mathcal{B}=\{a,b,c\}$ we will consider the corresponding rows of $\mathcal{A}$ (the first three).  There are five corresponding random variables:
$$\begin{aligned}
X_{\text{black}}(\{a,b,c\}) &= a_{a,\text{black}}\,a_{b,\text{black}}\,a_{c,\text{black}} = (1)(1)(1) = 1; \\
X_{\text{blue}}(\{a,b,c\}) &= a_{a,\text{blue}}\,a_{b,\text{blue}}\,a_{c,\text{blue}} = (1)(0)(1) = 0;
\end{aligned}$$
and so on.
Clearly, the value of $X_j$ is $0$ if and only if $j$ is not an element of at least one bag in the collection.  Consequently, letting $n = |S|$ be the number of columns of $\mathbb{A},$ the cardinality of the union of all bags in a collection $\mathcal{B}$ (its "union count") is given by
$$\text{Union count}(\mathcal{B}) = n - X(\mathcal{B})$$
where
$$X(\mathcal{B}) = \sum_{j\in S} X_j(\mathcal{B}).$$
I will analyze $X,$ because (this should be obvious) it will have the same variance as the union count.
On to the calculations.  We need to find means, variances, and covariances of $X$ when the domain of $X$ is restricted to all subsets of size $k.$  (I will use a superscript $(k)$ to designate this restriction.)  But this is now straightforward.
As a preliminary, for any column $j$ of $\mathbb A$ let $m_j$ be the sum of the column and for any pair of columns $j, j^\prime$ let $m_{jj^\prime}$ be the sum of the (componentwise) product of those columns:

$$m_j = \sum_{i} a_{ij};\quad m_{jj^\prime} = \sum_{i} a_{ij}\,a_{ij^\prime}.\tag{1}$$

(The latter are the entries of $\mathbb{A}^\prime \mathbb{A}.$)
Let $m$ be the number of rows of $\mathbb A$ (the number of bags).

*

*The mean of $X_j^{(k)}$ is the proportion of $k$-subsets of column $j$ of $\mathbb A$ where all the values are $1.$  Since there are $m_j$ ones in that column out of all $m$ entries, this proportion is the number of $k$-subsets of those ones, divided by the total number of possible subsets, $$E\left[X_j^{(k)}\right] = {\binom{m_j}{k}}/{\binom{m}{k}}.$$


*Because the values of every $X_j^{(k)}$ are in $\{0,1\},$ $\left(X_j^{(k)}\right)^2 = X_j^{(k)}.$  Thus,  $$E\left[\left(X_j^{(k)}\right)^2\right] = {\binom{m_j}{k}}/{\binom{m}{k}}.$$


*Similarly, the expectation of $X_j^{(k)}X_{j^\prime}^{(k)}$ is found by counting all bags $i$ where both columns of $\mathbb A$ contain $1$'s.  Thus $$E\left[X_j^{(k)}X_{j^\prime}^{(k)}\right] = \binom{m_{ij}}{k}/{\binom{m}{k}}.$$
We now have all the moments needed to compute the variance-covariance matrix

$$\begin{aligned}
\mathbb{V^{(k)}}_{jj^\prime} &= \operatorname{Cov}\left(X_j^{(k)}, X_{j^\prime}^{(k)}\right) \\
&= \binom{m_{jj^\prime}}{k}/{\binom{m}{k}} - \left(\binom{m_{j}}{k}/{\binom{m}{k}}\right)\left(\binom{m_{j^\prime}}{k}/{\binom{m}{k}}\right)\\
&=\binom{m}{k}^{-2}\left(\binom{m}{k}\,\binom{m_{jj^\prime}}{k} - \binom{m_{j}}{k}\,\binom{m_{j^\prime}}{k}\right).
\end{aligned}\tag{2}$$

Finally, because $X^{(k)}$ is the sum of the $X_j,$ its variance is

$$\operatorname{Var}\left(X^{(k)}\right) = (1,1,\ldots,1)\,\mathbb{V^{(k)}}\,(1,1,\ldots,1)^\prime = \sum_{j, j^\prime\in S}\mathbb{V^{(k)}}_{jj^\prime}.$$

This the sum of all the elements of this matrix, computed using the simple formulas $(1)$ and $(2)$ above.
I doubt there is much, if any, algebraic simplification available, because the variance must depend on the detailed relationship among all the bags: nothing short of making all bag-to-bag comparisons will do the trick.
The computational effort for this approach is proportional to the number of entries in the covariance matrix, on the order of $n^2,$ regardless of how many bags there might be.  Furthermore, $S$ can be restricted to just the union of elements actually appearing in all the bags, potentially decreasing $n.$
The variances in the example are all zero except when $k=1,$ where the variance is $3/16 = 0.1875.$

The following R code implements this algorithm.  After creating $\mathbb A,$ the calculations are noteworthy for their brevity.  Functions mu1, mu2, and variance comprise just five lines of code.  Using them, the calculation of the variances for all $k$ is simmply a matter of computing the variance matrix and summing its entries:
k <- 0:length(Bags)
V <- sapply(variance(k, A), sum)

Here's the full code.  It includes a commented-out section to generate synthetic datasets and at the end it checks the calculation with a direct computation of the random variables $X^{(k)}$ and using the R function var to obtain the answers.  I have run it with examples up to $m=500$ bags, where the calculation still takes less than a second.  The brute-force method used to check the calculation is exponentially slow, however, and so is limited to small values of $m.$
#
# All first moments.
#
mu1 <- function(k, A) {
  k <- pmax(0, pmin(k, nrow(A)))
  exp(lchoose(colSums(A), k) - lchoose(nrow(A), k))
} 
#
# All second moments.
#
mu2 <- function(k, A) {
  lapply(k, function(k) matrix(exp(lchoose(crossprod(A), k) - lchoose(nrow(A), k)), ncol(A)))
}
#
# Variance matrix.
#
variance <- function(k, A) {
  M <- lapply(k, function(k) (function(x) outer(x, x))(mu1(k, A)))
  mapply(function(m2, m) m2 - m, mu2(k, A), M, SIMPLIFY = FALSE)
}
#------------------------------------------------------------------------------#
# 
# EXAMPLES
#
# Represent the element relation for the example.
#
Bags <- list(a = c("red", "green"),
             b = c("red", "blue"),
             c = c("orange"),
             d = c("red", "black"))
S <- sort(unique(unlist(Bags)))
# #
# # Alternatively, create random examples.
# # The algorithms resource requirements scale quadratically in the size of S,
# # which can be as large as `n` (but might be smaller).
# # It scales linearly with `m`, though!
# #
# n <- 5     # Number of elements in the universal set
# m <- 6    # Number of bags
# p <- .4     # Expected proportion in each bag
# S <- paste0("_", seq_len(n))
# # set.seed(17)
# Bags <- lapply(seq_len(m), function(i) S[runif(n) <= p])
# # Supply bag ids for up to 26^2+26=702 bags:
# s <- apply(expand.grid(c(letters), c("", letters)), 1, paste0, collapse="")
# names(Bags) <- s[seq_len(m)] #
#------------------------------------------------------------------------------#
#
# CREATE THE MATRIX A
#
S <- sort(unique(unlist(Bags)))
indices <- seq_along(S)
if (length(intersect(indices, S)) > 0) warning("Integer names of bags may confuse R.")
names(indices) <- S
one <- rep(1, length(S))
A.raw <- sapply(Bags, function(bag) {
  a <- c(one, NA)        # Guarantees a non-empty result
  a[indices[bag]] <- 0
  a
})
A <- matrix(t(A.raw), ncol=length(S)+1)[,-(length(S)+1)] # Strips the NAs
colnames(A) <- S
rownames(A) <- names(Bags)
#------------------------------------------------------------------------------#
#
# The variances of the union counts.
#
k <- 0:length(Bags)
V <- sapply(variance(k, A), sum)
names(V) <- k
#------------------------------------------------------------------------------#
#
# CONFIRMATION
#
# Check with brute force.
#
if (length(Bags) <= 16) {
  UnionCount <- function(a.list) length(unique(unlist(a.list)))
  X <- sapply(k, function(k) {
    i <- combn(rownames(A), k) # Columns are samples of bags w/o replacement
    apply(i, 2, function(i) UnionCount(Bags[i]))
  })
  V.check <- sapply(X, function(x) {
    n <- length(x)
    ifelse(n <= 1, 0, var(x) * (n-1) / n)
  })
  names(V.check) <- k
  #
  # Compare.
  #
  print(results <- rbind(Formula = zapsmall(V), `Brute force` = V.check))
  print(Bags)
  print(A)
  if(!isTRUE(all.equal.numeric(V, V.check))) stop("Discrepant results found.")
} else {
  print(V)
  print(A)
}

