When (in the context of 'non-independent tests') does Bonferroni correction *not* yield the appropriate FWER? I have a question about the oft-heard suggestion that a Bonferroni correction is appropriate in the context of posing multiple independent hypotheses and inappropriate elsewhere, in contrast to statements like:

One of the most widespread of these misunderstandings is that the method would be based on an assumption of independence between p-values. This misunderstanding comes from a frequently used, but deficient motivation for Bonferroni, saying that the probability of making a false rejection if all $m_0$ p-values of true hypotheses are independent, and we perform each test at level... By this reasoning, Bonferroni seems like a method that only provides approximate FWER control and that requires an assumption of independence for its validity. In fact, the method of Bonferroni provides exact FWER control under any dependence structure of the p-values.

To explore this question numerically, I wrote a small R script where I simulated multivariate normal data according to some correlation matrix (sampled from a flat LKJ), estimated sample correlations among that data, and computed the corresponding p-values. I then estimated regression coefficients between these samples $x$ and independently generated data (i.e. from $y = a + 0 \cdot x + e$), and then computed the corresponding p-values. After Bonferroni correction, the FWER for both was at the desired value, $\alpha$. This should also apply if the correlation structure / non-independence were in $y$, given the symmetry between $y \sim x$ and $x \ sim y$.
So is this whole 'non-independence between tests' thing actually a thing w/ the Bonferroni? I'm aware there are other motivations for disfavoring the correction, e.g. it only bounding the FWER from above, and with sufficiently many tests over finite data you get ever higher false-negative rates as the significance threshold goes to 0. If the dependence structure between tests is irrelevant to Bonferroni correction, where did this advised requirement of "multiple independent hypotheses" come from?
(sorry if this question is ill-posed -- maybe I'm misunderstanding the term 'independence' here. I usually just fit hierarchical models, cite Gelman, Hill, & Yajima '12, and call it a day, but am working more in a NHST world lately and am trying to wrap my head around different multiple comparisons strategies)

Neither of the positions in the quote is correct: the Bonferroni method does not assume independence of the p-values, but it is also not an exact method.  The Bonferroni correction is actually a highly conservative (non-exact) correction and so it is usually the case that the true family-wise error rate is substantially below the desired value $\alpha$ under this procedure.  In fact, the only time the Bonferroni method will give you the exact family-wise error rate you stipulate is when the events that the individual hypothesis are rejected are mutually exclusive, which never happens in actual statistical practice.
You say that you did a simulation and it gave you a family-wise error rate at the desired level $\alpha$.  Presumably what you mean is that the actual family-wise error rate was bounded by the desired level $\alpha$ --- it would be highly unusual (and indicative of a mistake) if it was at that exact value.

Let me illustrate how the Bonferroni method works, and what it relies on mathematically.  To do this, suppose we are conducting $m$ comparisons.  Let $p_i$ denote the p-value for comparison $i$ and let $\mathbf{H}_0$ denote the null hypothesis for all the comparisons.  Under the Bonferroni method we reject the null hypothesis if $p_i \leqslant \alpha/m$.  Consequently, under this method, the error rate for the $i$th comparison is the probability of Type I error, which is:
$$\text{ER}_i \equiv \mathbb{P} \Big( p_i \leqslant \frac{\alpha}{m} \Big| \mathbf{H}_0 \Big) = \frac{\alpha}{m},$$
and the family-wise error rate for the set of all $m$ comparisons is the probability of Type I error for at least one of those comparisons, which is:
$$\text{FWER}_m \equiv \mathbb{P} \Big( \min (p_1,...,p_m) \leqslant \frac{\alpha}{m} \Big| \mathbf{H}_0 \Big).$$
Without any assumption of independence, you can apply Boole's inequality to get:
$$\begin{align}
\text{FWER}_m 
&= \mathbb{P} \Big( \min (p_1,...,p_m) \leqslant \frac{\alpha}{m} \Big| \mathbf{H}_0 \Big) \\[6pt]
&= \mathbb{P} \Big( p_1 \leqslant \frac{\alpha}{m}, \cdots , p_m \leqslant \frac{\alpha}{m} \Big| \mathbf{H}_0 \Big) \\[6pt]
&\leqslant \sum_{i=1}^m \mathbb{P} \Big( p_i \leqslant \frac{\alpha}{m} \Big| \mathbf{H}_0 \Big) \\[6pt]
&= \sum_{i=1}^m \frac{\alpha}{m} \\[6pt]
&= \alpha. \\[6pt]
\end{align}$$
As you can see, the use of Boole's inequality gives this result a highly conservative character.  The only time you will get $\text{FWER}_m = \alpha$ exactly is when the events for rejection of the null hypothesis are mutually exclusive, which never really happens in practice.  Depending on the dependency (or independence) between the p-values for the comparisons, the exact family-wise error rate will usually be substantially less than the desired level $\alpha$.
