Testing independence of n binary variables I have 6 binary (dummy) variables corresponding to different possible states that could apply to the subjects, and each subject has two of the six possible states, ie the sum of all of their dummy variables = 2. I am interested in how well each of the pairs are related, ie which pairs tend to appear together, but I'm not sure the best way to approach this problem. I suppose that if each state was independent of the others, there would be the same number of each possible combination, which would be my null hypothesis. Would I be doing a chisq test comparing the observed vs the expected values? I have individual subject-level data as well, but summarized as below:
   A   B   C   D   E   F
A  -  44  99 129  71   3
B  -   -  88  52   5   0
C  -   -   -  41   2   0
D  -   -   -   -  10   1 
E  -   -   -   -   -   0
F  -   -   -   -   -   -

It appears that A-D are more commonly paired than E-F, but I'm not sure how to formally test or express this statistically. I'd really appreciate any help, or even just pointing me to learn about a topic/test. Thank you.
 A: This is not a standard situation because it concerns unordered pairs.  It needs a model and some analysis.

The model describes the state pairs when there is no association between states.  One plausible and flexible model supposes each state $s$ is associated with a constant but unknown probability $\pi_s$.  (We introduce these probabilities, and allow them to differ, for otherwise our analysis would be unable to distinguish lack of independence of states from a non-uniform distribution of state frequencies.)
A subject's states are presumed to arise by selecting one state according to this probability distribution and then, independently, to select a second state from the truncated probability distribution of all the remaining states.  Thus, there are two ways any unordered pair $\{s,t\}$ may arise: either $s$ was selected first or $t$ was.  The total probability of this pair is
$$\pi_{\{s,t\}} = \Pr(s)\Pr(t\mid s) + \Pr(t)\Pr(s\mid t) = \pi_s \frac{\pi_t}{1-\pi_s} + \pi_t \frac{\pi_s}{1-\pi_t}.$$
When there are many subjects and most unordered pairs are well-represented in the data (which is the case in the question), a maximum likelihood estimate will work well.  To obtain it we need the log likelihood, which (up to a constant) is
$$\Lambda = \sum_{\{s,t\}} n_{\{s,t\}} \log \pi_{\{s,t\}}$$
where $n_{\{s,t\}}$ designates the count of subjects exhibiting states $s$ and $t.$
There is no evident closed form solution for the minimum, but it can be obtained numerically by searching over all but one value of the parameters $\log\pi_{s}$ (the last value is determined by all the rest since their exponentials must sum to unity by the Law of Total Probability).  A good starting value for a search of the minimum would be to estimate the $\pi_s$ from the relative total frequencies of all pairs $\{s,t\}.$  Values of $\pi_s$ that achieve a minimum are written $\hat\pi_s$ and the pair-probabilities computed from them (using the foregoing formula) are written $\hat\pi_{\{s,t\}}.$
Comparing this to the "saturated model" that includes one parameter for each possible pair gives a standard Chi-squared test.  Specifically, when the total number of subjects is $n,$ the estimated count of the pair $\{s,t\}$ is $\hat\pi_{\{s,t\}}n$ which is compared to the observed count $n_{\{s,t\}}$ by means of the chi-squared residual
$$r_{\{s,t\}} = \frac{n_{\{s,t\}} - \hat\pi_{\{s,t\}}n}{\sqrt{\hat\pi_{\{s,t\}}n}}.$$
Residuals much larger than $1$ in size are unusual when the model is correct: thus, inspecting the array of residuals for large values is informative.
The chi-squared statistic is the sum of their squares,
$$\chi^2 = \sum_{\{s,t\}} r_{\{s,t\}}^2.$$
In terms of the number $m$ of possible states, the degrees of freedom $\nu$ are the number of pairs $\binom{m}{2}$ minus the number of parameters $m,$ equal to $9$ when $m=6$ in the question.  The p-value of this test is obtained by computing the upper tail area for $\chi^2$ using the Chi-squared distribution with $\nu$ degrees of freedom.
Some simulations bear out the efficacy of this technique.  For instance, 500 datasets of 545 subjects with six states were generated by the model and the resulting p-values appeared to be uniformly distributed, as one would hope.  The code is appended at the end.

In the question, $n=545$ and $m=6.$  The Maximum Likelihood estimate is
$$(\hat\pi_A, \hat\pi_B, \ldots, \hat\pi_F) = (0.356, 0.162, 0.203, 0.205, 0.071, 0.003).$$
The chi-squared statistic of $135.6$ with $9$ degrees of freedom has an astronomically low p-value (around $10^{-24}$), indicating this model is somehow incorrect.  Given that its basic assumptions are plausible, the natural conclusion is that the assumption of independence among the states is wrong.
It may be of interest to see how this independence assumption fails.  The residuals of largest size are $r_{\{A,B\}}=-4.5,$ $r_{\{A,E\}}=5.8,$ $r_{\{B,C\}}=6.7,$ and $r_{\{C,E\}} = -3.8.$  This means there are too few $\{A,B\}$ and $\{C,E\}$ pairs and too many $\{A,E\}$ and $\{B,C\}$ pairs to be consistent both with the observed frequencies of states found among all the subjects and the assumption that states are independent.

R Code
#
# MLE for a lower triangular count matrix of unordered pairs.
#
mle.pairs <- function(X, ...) {
  p.hat <- log((colSums(X)+1)/(sum(X) + ncol(X)))
  p.hat <- p.hat[-length(p.hat)]

  iterations <- 0
  while(iterations < 1e4) {
    fit <- nlm(lambda, p.hat, X=X, ...)
    iterations <- iterations + fit$iterations
    if (fit$code < 4) break
    p.hat <- fit$estimate
  }
  fit$iterations <- iterations
  fit$p.hat <- exp(fit$estimate)
  fit$p.hat <- c(fit$p.hat, 1-sum(fit$p.hat))
  fit$pi.hat <- outer(fit$p.hat, fit$p.hat, function(x,y) x * y * (1/(1-x) + 1/(1-y)))
  fit
}
#
# Chi-squared test for a count matrix of unordered pairs.
#
chisq.test.unordered <- function(X) {
  N <- sum(X)
  fit <- mle.pairs(X)
  e <- N * fit$pi.hat 
  r <- ((X-e) / sqrt(e))[lower.tri(X)]
  chi.sq <- sum(r^2)
  df <- length(r) - length(p)
  p.value <- ifelse(fit$code < 4, pchisq(chi.sq, df, lower.tail=FALSE), NA)
  c(list(stat=chi.sq, df=df, p.value=p.value, residuals=r), fit)
}
#
# Negative log likelihood as a function of the log probabilities `p`.
#
lambda <- function(p, X) {
  i <- lower.tri(X)
  p <- c(p, log(1-sum(exp(p))))
  pi <- outer(p, p, function(x,y) x + y + log(1/(1-exp(x)) + 1/(1-exp(y))))
  -sum(X[i] * pi[i])
}
#
# Simulation check.
# The p-values ought to be close to uniformly distributed under the null.
#
n <- 545
p <- rgamma(6, 2)
p <- p/sum(p)
k <- length(p)

pp <- c(k, 1) %*% (combn(length(p), 2) - 1) + 1
pi <- outer(p, p, function(x,y) x * y * (1/(1-x) + 1/(1-y)))
pi <- pi[lower.tri(pi)]

sim <- replicate(5e2, {
  Y <- sample(pp, n, replace=TRUE, prob=pi)
  X <- matrix(tabulate(Y, nbins=k^2), k)
  chisq.test.unordered(X)$p.value
})
hist(sim)
plot(ecdf(sim))
ks.test(sim, punif)
#
# The test applied to the data in the question.
#
X <- matrix(0, 6, 6)
X[upper.tri(X)] <- c(44,99,88,129,52,41,71,5,2,10,3,0,0,1,0)
X <- t(X)

fit <- chisq.test.unordered(X) # Inspect this for the results.
#
# Example: display the residuals.
#
residuals <- X
residuals[lower.tri(residuals)] <- fit$residuals
residuals

A: Your intuition is right, you test observed/expected, and you can do using a chi-squared test, but you have to derive the right expected/observed values. 
What you have is a contingency table. There any many way to test association between variable (e.g. https://en.wikipedia.org/wiki/Contingency_table#Measures_of_association). The easiest and intuitive method is probably using the odds ratio (see wiki link). There is a very good online compendium explaining how to do this in R: https://xiangzhu.github.io/stanford-stats110/twoway_tab.
Hope this helps!
Ciao,
Emilio
