Bill van Zwet’s contributions to his profession are, in their own way, comparable to his very substantial research, research supervision and teaching contributions. He had several leadership roles in professional societies, edited leading theoretical journals, started up or revived important conference series, founded a research center and played a huge role in providing colleagues in Central and Eastern Europe access to the western world of probability and statistics. And to each of these several activities, he brought farsightedness, energy and wit, qualities that shone through when he was interviewed about what he’d done (Statist. Sci. 24 (2009) 87–115). In fact, because Bill provided a wealth of detail in “the interview,” this article focuses more on the recollections of some of the people with whom he interacted.

Willem van Zwet made deep and influential contributions to probability and statistics, which we review in this paper. Bickel and Götze collaborated with him on his major contributions to higher order asymptotics of nonlinear statistics and on resampling and the bootstrap. We relate this work to his remarkable development of the properties of the Hoeffding expansion for symmetric statistics as well as Fourier analytic tools. Fiocco and De Gunst were his students. We describe how in their theses and subsequent papers with him they developed statistical inference in two subtle stochastic models, the contact process and cell population development under a plausible regime. We also touch on his solutions of intriguing problems not related directly to his main interests.

Willem van Zwet was supervisor of sixteen PhD students. All of them pursued academic careers and most of them became full professor. Below are some stories of PhD students Wim Albers, Cees Diks, Ronald Does, Marta Fiocco, Sara van de Geer, Mathisca de Gunst, Chris Klaassen, Hein Putter, Aad van der Vaart, Marten Wegkamp and Martien van Zuijlen with in addition a contribution by Nelly Litvak who was guided by Willem after her PhD.

We study the problem of independence testing given independent and identically distributed pairs taking values in a σ-finite, separable measure space. Defining a natural measure of dependence $\mathit{D}(\mathit{f})$ as the squared ${\mathit{L}}^2$-distance between a joint density f and the product of its marginals, we first show that there is no valid test of independence that is uniformly consistent against alternatives of the form $\{\mathit{f}:\mathit{D}(\mathit{f})\ge {\mathit{\rho }}^2\}$. We therefore restrict attention to alternatives that impose additional Sobolev-type smoothness constraints, and define a permutation test based on a basis expansion and a U-statistic estimator of $\mathit{D}(\mathit{f})$ that we prove is minimax optimal in terms of its separation rates in many instances. Finally, for the case of a Fourier basis on ${[0,1]}^2$, we provide an approximation to the power function that offers several additional insights. Our methodology is implemented in the $\mathtt{R}$ package $\mathtt{USP}$.

Bayesian nonparametric regression under a rescaled Gaussian process prior offers smoothness-adaptive function estimation with near minimax-optimal error rates. Hierarchical extensions of this approach, equipped with stochastic variable selection, are known to also adapt to the unknown intrinsic dimension of a sparse true regression function. But it remains unclear if such extensions offer variable selection consistency, that is, if the true subset of important variables could be consistently learned from the data. It is shown here that variable consistency may indeed be achieved with such models at least when the true regression function has finite smoothness to induce a polynomially larger penalty on inclusion of false positive predictors. Our result covers the high-dimensional asymptotic setting where the predictor dimension is allowed to grow with the sample size. The proof utilizes Schwartz theory to establish that the posterior probability of wrong selection vanishes asymptotically. A necessary and challenging technical development involves providing sharp upper and lower bounds to small ball probabilities at all rescaling levels of the Gaussian process prior, a result that could be of independent interest.

Spectral clustering is one of the most popular algorithms to group high- dimensional data. It is easy to implement and computationally efficient. Despite its popularity and successful applications, its theoretical properties have not been fully understood. In this paper, we show that spectral clustering is minimax optimal in the Gaussian mixture model with isotropic covariance matrix, when the number of clusters is fixed and the signal-to-noise ratio is large enough. Spectral gap conditions are widely assumed in the literature to analyze spectral clustering. On the contrary, these conditions are not needed to establish optimality of spectral clustering in this paper.

Proposed by Donoho (Ann. Statist. 25 (1997) 1870–1911), Dyadic CART is a nonparametric regression method which computes a globally optimal dyadic decision tree and fits piecewise constant functions in two dimensions. In this article, we define and study Dyadic CART and a closely related estimator, namely Optimal Regression Tree (ORT), in the context of estimating piecewise smooth functions in general dimensions in the fixed design setup. More precisely, these optimal decision tree estimators fit piecewise polynomials of any given degree. Like Dyadic CART in two dimensions, we reason that these estimators can also be computed in polynomial time in the sample size N via dynamic programming. We prove oracle inequalities for the finite sample risk of Dyadic CART and ORT, which imply tight risk bounds for several function classes of interest. First, they imply that the finite sample risk of ORT of order $\mathit{r}\ge 0$ is always bounded by $\mathit{C}\mathit{k}\frac{log\mathit{N}}{\mathit{N}}$ whenever the regression function is piecewise polynomial of degree r on some reasonably regular axis aligned rectangular partition of the domain with at most k rectangles. Beyond the univariate case, such guarantees are scarcely available in the literature for computationally efficient estimators. Second, our oracle inequalities uncover minimax rate optimality and adaptivity of the Dyadic CART estimator for function spaces with bounded variation. We consider two function spaces of recent interest where multivariate total variation denoising and univariate trend filtering are the state of the art methods. We show that Dyadic CART enjoys certain advantages over these estimators while still maintaining all their known guarantees.

Multivariate extreme value theory is concerned with modeling the joint tail behavior of several random variables. Existing work mostly focuses on asymptotic dependence, where the probability of observing a large value in one of the variables is of the same order as observing a large value in all variables simultaneously. However, there is growing evidence that asymptotic independence is equally important in real world applications. Available statistical methodology in the latter setting is scarce and not well understood theoretically. We revisit nonparametric estimation and introduce rank-based M-estimators for parametric models that simultaneously work under asymptotic dependence and asymptotic independence, without requiring prior knowledge on which of the two regimes applies. Asymptotic normality of the proposed estimators is established under weak regularity conditions. We further show how bivariate estimators can be leveraged to obtain parametric estimators in spatial tail models, and again provide a thorough theoretical justification for our approach.

Let ${\mathit{X}}_1,\dots ,{\mathit{X}}_{\mathit{n}}$ be i.i.d. random variables sampled from a normal distribution $\mathit{N}(\mathit{\mu },\mathrm{\Sigma })$ in ${\mathbb{R}}^{\mathit{d}}$ with unknown parameter $\mathit{\theta }=(\mathit{\mu },\mathrm{\Sigma })\in \mathrm{\Theta }:={\mathbb{R}}^{\mathit{d}}\times {\mathcal{C}}_{+}^{\mathit{d}}$, where ${\mathcal{C}}_{+}^{\mathit{d}}$ is the cone of positively definite covariance operators in ${\mathbb{R}}^{\mathit{d}}$. Given a smooth functional $\mathit{f}:\mathrm{\Theta }\mapsto {\mathbb{R}}^1$, the goal is to estimate $\mathit{f}(\mathit{\theta })$ based on ${\mathit{X}}_1,\dots ,{\mathit{X}}_{\mathit{n}}$. Let

$$\mathrm{\Theta }(\mathit{a};\mathit{d}):={\mathbb{R}}^{\mathit{d}}\times \{\mathrm{\Sigma }\in {\mathcal{C}}_{+}^{\mathit{d}}:\mathit{\sigma }(\mathrm{\Sigma })\subset [1/\mathit{a},\mathit{a}]\},\mathit{a}\ge 1,$$

where $\mathit{\sigma }(\mathrm{\Sigma })$ is the spectrum of covariance Σ. Let $\hat{\mathit{\theta }}:=(\hat{\mathit{\mu }},\hat{\mathrm{\Sigma }})$, where $\hat{\mathit{\mu }}$ is the sample mean and $\hat{\mathrm{\Sigma }}$ is the sample covariance, based on the observations ${\mathit{X}}_1,\dots ,{\mathit{X}}_{\mathit{n}}$. For an arbitrary functional $\mathit{f}\in {\mathit{C}}^{\mathit{s}}(\mathrm{\Theta })$, $\mathit{s}=\mathit{k}+1+\mathit{\rho },\mathit{k}\ge 0,\mathit{\rho }\in (0,1]$, we define a functional ${\mathit{f}}_{\mathit{k}}:\mathrm{\Theta }\mapsto \mathbb{R}$ such that

$$\underset{\mathit{\theta }\in \mathrm{\Theta }(\mathit{a};\mathit{d})}{sup}{\Vert {\mathit{f}}_{\mathit{k}}(\hat{\mathit{\theta }})-\mathit{f}(\mathit{\theta })\Vert }_{{\mathit{L}}_2({\mathbb{P}}_{\mathit{\theta }})}{\lesssim }_{\mathit{s},\mathit{\beta }}\Vert \mathit{f}{\Vert }_{{\mathit{C}}^{\mathit{s}}(\mathrm{\Theta })}\left[\right(\frac{\mathit{a}}{\sqrt{\mathit{n}}}\vee {\mathit{a}}^{\mathit{\beta }\mathit{s}}{\left(\sqrt{\frac{\mathit{d}}{\mathit{n}}}\right)}^{\mathit{s}})\wedge 1],$$

where $\mathit{\beta }=1$ for $\mathit{k}=0$ and $\mathit{\beta }>\mathit{s}-1$ is arbitrary for $\mathit{k}\ge 1$. This error rate is minimax optimal and similar bounds hold for more general loss functions. If $\mathit{d}={\mathit{d}}_{\mathit{n}}\le {\mathit{n}}^{\mathit{\alpha }}$ for some $\mathit{\alpha }\in (0,1)$ and $\mathit{s}\ge \frac{1}{1-\mathit{\alpha }}$, the rate becomes $\mathit{O}({\mathit{n}}^{-1/2})$. Moreover, for $\mathit{s}>\frac{1}{1-\mathit{\alpha }}$, the estimator ${\mathit{f}}_{\mathit{k}}(\hat{\mathit{\theta }})$ is shown to be asymptotically efficient. The crucial part of the construction of estimator ${\mathit{f}}_{\mathit{k}}(\hat{\mathit{\theta }})$ is a bias reduction method studied in the paper for more general statistical models than normal.

Additive regression is studied in a very general setting where both the response and predictors are allowed to be non-Euclidean. The response takes values in a general separable Hilbert space, whereas the predictors take values in general semimetric spaces, which covers a very wide range of nonstandard response variables and predictors. A general framework of estimating additive models is presented for semimetric space-valued predictors. In particular, full details of implementation and the corresponding theory are given for predictors taking values in Hilbert spaces and/or Riemannian manifolds. The existence of the estimators, convergence of a backfitting algorithm, rates of convergence and asymptotic distributions of the estimators are discussed. The finite sample performance of the estimators is investigated by means of two simulation studies. Finally, three data sets covering several types of non-Euclidean data are analyzed to illustrate the usefulness of the proposed general approach.

Entropy integrals are widely used as a powerful empirical process tool to obtain upper bounds for the rates of convergence of global empirical risk minimizers (ERMs), in standard settings such as density estimation and regression. The upper bound for the convergence rates thus obtained typically matches the minimax lower bound when the entropy integral converges, but admits a strict gap compared to the lower bound when it diverges. Birgé and Massart (Probab. Theory Related Fields 97 (1993) 113–150) provided a striking example showing that such a gap is real with the entropy structure alone: for a variant of the natural Hölder class with low regularity, the global ERM actually converges at the rate predicted by the entropy integral that substantially deviates from the lower bound. The counter-example has spawned a long-standing negative position on the use of global ERMs in the regime where the entropy integral diverges, as they are heuristically believed to converge at a suboptimal rate in a variety of models.

The present paper demonstrates that this gap can be closed if the models admit certain degree of “set structures” in addition to the entropy structure. In other words, the global ERMs in such set structured models will indeed be rate-optimal, matching the lower bound even when the entropy integral diverges. The models with set structures we investigate include (i) image and edge estimation, (ii) binary classification, (iii) multiple isotonic regression, (iv) s-concave density estimation, all in general dimensions when the entropy integral diverges. Here, set structures are interpreted broadly in the sense that the complexity of the underlying models can be essentially captured by the size of the empirical process over certain class of measurable sets, for which matching upper and lower bounds are obtained to facilitate the derivation of sharp convergence rates for the associated global ERMs.

For multivariate distributions in the domain of attraction of a max-stable distribution, the tail copula and the stable tail dependence function are equivalent ways to capture the dependence in the upper tail. The empirical versions of these functions are rank-based estimators whose inflated estimation errors are known to converge weakly to a Gaussian process that is similar in structure to the weak limit of the empirical copula process. We extend this multivariate result to continuous functional data by establishing the asymptotic normality of the estimators of the tail copula, uniformly over all finite subsets of at most D points (D fixed). An application for testing tail copula stationarity is presented. The main tool for deriving the result is the uniform asymptotic normality of all the D-variate tail empirical processes. The proof of the main result is nonstandard.

Two-stage enrichment designs can be used to target the benefiting population in clinical trials based on patients’ biomarkers. In the case of continuous biomarkers, we show that using a bivariate model that treats biomarkers as random variables more accurately identifies a treatment-benefiting enriched population than assuming biomarkers are fixed. Additionally, we show that under the bivariate model, the maximum likelihood estimators (MLEs) follow a randomly scaled mixture of normal distributions. Using random normings, we obtain asymptotically standard normal MLEs and construct hypothesis tests. Finally, in a simulation study, we demonstrate that our proposed design is more powerful than a single stage design when outcomes and biomarkers are correlated; the model-based estimators have smaller bias and mean square error (MSE) than weighted average estimators.

In matrix-valued datasets the sampled matrices often exhibit correlations among both their rows and their columns. A useful and parsimonious model of such dependence is the matrix normal model, in which the covariances among the elements of a random matrix are parameterized in terms of the Kronecker product of two covariance matrices, one representing row covariances and one representing column covariance. An appealing feature of such a matrix normal model is that the Kronecker covariance structure allows for standard likelihood inference even when only a very small number of data matrices is available. For instance, in some cases a likelihood ratio test of dependence may be performed with a sample size of one. However, more generally the sample size required to ensure boundedness of the matrix normal likelihood or the existence of a unique maximizer depends in a complicated way on the matrix dimensions. This motivates the study of how large a sample size is needed to ensure that maximum likelihood estimators exist, and exist uniquely with probability one. Our main result gives precise sample size thresholds in the paradigm where the number of rows and the number of columns of the data matrices differ by at most a factor of two. Our proof uses invariance properties that allow us to consider data matrices in canonical form, as obtained from the Kronecker canonical form for matrix pencils.

We establish adaptive results for trend filtering: least squares estimation with a penalty on the total variation of $(\mathit{k}-1)$th order differences. Our approach is based on combining a general oracle inequality for the ${\ell }_1$-penalized least squares estimator with “interpolating vectors” to upper bound the “effective sparsity.” This allows one to show that the ${\ell }_1$-penalty on the kth order differences leads to an estimator that can adapt to the number of jumps in the $(\mathit{k}-1)$th order differences of the underlying signal or an approximation thereof. We show the result for $\mathit{k}\in \{1,2,3,4\}$ and indicate how it could be derived for general $\mathit{k}\in \mathbb{N}$.

In time series analysis there is an apparent dichotomy between time and frequency domain methods. The aim of this paper is to draw connections between frequency and time domain methods. Our focus will be on reconciling the Gaussian likelihood and the Whittle likelihood. We derive an exact, interpretable, bound between the Gaussian and Whittle likelihood of a second order stationary time series. The derivation is based on obtaining the transformation which is biorthogonal to the discrete Fourier transform of the time series. Such a transformation yields a new decomposition for the inverse of a Toeplitz matrix and enables the representation of the Gaussian likelihood within the frequency domain. We show that the difference between the Gaussian and Whittle likelihood is due to the omission of the best linear predictions outside the domain of observation in the periodogram associated with the Whittle likelihood. Based on this result, we obtain an approximation for the difference between the Gaussian and Whittle likelihoods in terms of the best fitting, finite order autoregressive parameters. These approximations are used to define two new frequency domain quasi-likelihood criteria. We show that these new criteria can yield a better approximation of the spectral divergence criterion, as compared to both the Gaussian and Whittle likelihoods. In simulations, we show that the proposed estimators have satisfactory finite sample properties.

Inference after model selection has been an active research topic in the past few years, with numerous works offering different approaches to addressing the perils of the reuse of data. In particular, major progress has been made recently on large and useful classes of problems by harnessing general theory of hypothesis testing in exponential families, but these methods have their limitations. Perhaps most immediate is the gap between theory and practice: implementing the exact theoretical prescription in realistic situations—for example, when new data arrives and inference needs to be adjusted accordingly—may be a prohibitive task. In this paper, we propose a Bayesian framework for carrying out inference after variable selection. Our framework is very flexible in the sense that it naturally accommodates different models for the data instead of requiring a case-by-case treatment. This flexibility is achieved by considering the full selective likelihood function where, crucially, we propose a novel and nontrivial approximation to the exact but intractable expression. The advantages of our methods in practical data analysis are demonstrated in an application to HIV drug-resistance data.

Privacy-preserving data analysis is a rising challenge in contemporary statistics, as the privacy guarantees of statistical methods are often achieved at the expense of accuracy. In this paper, we investigate the tradeoff between statistical accuracy and privacy in mean estimation and linear regression, under both the classical low-dimensional and modern high-dimensional settings. A primary focus is to establish minimax optimality for statistical estimation with the $(\mathit{\epsilon },\mathit{\delta })$-differential privacy constraint. By refining the “tracing adversary” technique for lower bounds in the theoretical computer science literature, we improve existing minimax lower bound for low-dimensional mean estimation and establish new lower bounds for high-dimensional mean estimation and linear regression problems. We also design differentially private algorithms that attain the minimax lower bounds up to logarithmic factors. In particular, for high-dimensional linear regression, a novel private iterative hard thresholding algorithm is proposed. The numerical performance of differentially private algorithms is demonstrated by simulation studies and applications to real data sets.

This paper considers distributed statistical inference for general symmetric statistics in the context of massive data with efficient computation. Estimation efficiency and asymptotic distributions of the distributed statistics are provided, which reveal different results between the nondegenerate and degenerate cases, and show the number of the data subsets plays an important role. Two distributed bootstrap methods are proposed and analyzed to approximation the underlying distribution of the distributed statistics with improved computation efficiency over existing methods. The accuracy of the distributional approximation by the bootstrap are studied theoretically. One of the methods, the pseudo-distributed bootstrap, is particularly attractive if the number of datasets is large as it directly resamples the subset-based statistics, assumes less stringent conditions and its performance can be improved by studentization.

A considerable portion of the work on mixed orthogonal arrays applies specifically to arrays of strength 2. Although strength $\mathit{t}=2$ is arguably the most important case for statistical applications, there is an urgent need for better methods for $\mathit{t}\ge 3$. However, the knowledge on the existence of arrays for $\mathit{t}\ge 3$ is rather limited. In this paper, new construction methods for symmetric and asymmetric orthogonal arrays (OAs) with high strength are proposed by using lower strength orthogonal partitions of spaces and OAs. A positive answer is provided to the open problem in Hedayat, Sloane and Stufken (Orthogonal Arrays: Theory and Applications (1999) Springer) on developing better methods and tools for the construction of mixed orthogonal arrays with strength $\mathit{t}\ge 3$. Not only are the methods straightforward, but also they are useful for constructing symmetric or asymmetric OAs of arbitrary strengths, numbers of levels and various sizes. The constructed OAs can be utilized to generate more OAs. The resulting OAs have a high degree of flexibility and many other desirable properties. Some selective OAs are tabulated for practical uses.

Structural causal models (SCMs), also known as (nonparametric) structural equation models (SEMs), are widely used for causal modeling purposes. In particular, acyclic SCMs, also known as recursive SEMs, form a well-studied subclass of SCMs that generalize causal Bayesian networks to allow for latent confounders. In this paper, we investigate SCMs in a more general setting, allowing for the presence of both latent confounders and cycles. We show that in the presence of cycles, many of the convenient properties of acyclic SCMs do not hold in general: they do not always have a solution; they do not always induce unique observational, interventional and counterfactual distributions; a marginalization does not always exist, and if it exists the marginal model does not always respect the latent projection; they do not always satisfy a Markov property; and their graphs are not always consistent with their causal semantics. We prove that for SCMs in general each of these properties does hold under certain solvability conditions. Our work generalizes results for SCMs with cycles that were only known for certain special cases so far. We introduce the class of simple SCMs that extends the class of acyclic SCMs to the cyclic setting, while preserving many of the convenient properties of acyclic SCMs. With this paper, we aim to provide the foundations for a general theory of statistical causal modeling with SCMs.

The proposal and study of dependent Bayesian nonparametric models has been one of the most active research lines in the last two decades, with random vectors of measures representing a natural and popular tool to define them. Nonetheless, a principled approach to understand and quantify the associated dependence structure is still missing. We devise a general, and not model-specific, framework to achieve this task for random measure based models, which consists in: (a) quantify dependence of a random vector of probabilities in terms of closeness to exchangeability, which corresponds to the maximally dependent coupling with the same marginal distributions, that is, the comonotonic vector; (b) recast the problem in terms of the underlying random measures (in the same Fréchet class) and quantify the closeness to comonotonicity; (c) define a distance based on the Wasserstein metric, which is ideally suited for spaces of measures, to measure the dependence in a principled way. Several results, which represent the very first in the area, are obtained. In particular, useful bounds in terms of the underlying Lévy intensities are derived relying on compound Poisson approximations. These are then specialized to popular models in the Bayesian literature leading to interesting insights.

This paper delivers improved theoretical guarantees for the convex programming approach in low-rank matrix estimation, in the presence of (1) random noise, (2) gross sparse outliers and (3) missing data. This problem, often dubbed as robust principal component analysis (robust PCA), finds applications in various domains. Despite the wide applicability of convex relaxation, the available statistical support (particularly the stability analysis in the presence of random noise) remains highly suboptimal, which we strengthen in this paper. When the unknown matrix is well conditioned, incoherent and of constant rank, we demonstrate that a principled convex program achieves near-optimal statistical accuracy, in terms of both the Euclidean loss and the ${\ell }_{\infty }$ loss. All of this happens even when nearly a constant fraction of observations are corrupted by outliers with arbitrary magnitudes. The key analysis idea lies in bridging the convex program in use and an auxiliary nonconvex optimization algorithm, and hence the title of this paper.

Delayed-acceptance Metropolis–Hastings and delayed-acceptance pseudo-marginal Metropolis–Hastings algorithms can be applied when it is computationally expensive to calculate the true posterior or an unbiased stochastic approximation thereof, but a computationally cheap deterministic approximation is available. An initial accept–reject stage uses the cheap approximation for computing the Metropolis–Hastings ratio; proposals which are accepted at this stage are subjected to a further accept–reject step, which corrects for the error in the approximation. Since the expensive posterior, or the approximation thereof, is only evaluated for proposals which are accepted at the first stage, the cost of the algorithm is reduced and larger scalings may be used.

We focus on the random walk Metropolis (RWM) and consider the delayed-acceptance RWM and the delayed-acceptance pseudo-marginal RWM. We provide a framework for incorporating relatively general deterministic approximations into the theoretical analysis of high-dimensional targets. Justified by diffusion-approximation arguments, we derive expressions for the limiting efficiency and acceptance rates in high-dimensional settings. Finally, these theoretical insights are leveraged to formulate practical guidelines for the efficient tuning of the algorithms. The robustness of these guidelines and predicted properties are verified against simulation studies, all of which are strictly outside of the domain of validity of our limit results.

When the response mechanism is believed to be not missing at random (NMAR), a valid analysis requires stronger assumptions on the response mechanism than standard statistical methods would otherwise require. Semiparametric estimators have been developed under the parametric model assumptions on the response mechanism. In this paper, a new statistical test is proposed to guarantee model identifiability without using instrumental variable assumption. Furthermore, we develop optimal semiparametric estimation for parameters such as the population mean. Specifically, we propose two semiparametric optimal estimators that do not require any model assumptions other than the response mechanism. Asymptotic properties of the proposed estimators are discussed. An extensive simulation study is presented to compare with some existing methods. We present an application of our method using Korean labor and income panel survey data.

Regular $2^{\mathit{n}-\mathit{p}}$ designs are also known as single flat designs. Parallel flats designs (PFDs) consisting of three parallel flats (3-PFDs) are the most frequently utilized PFDs, due to their simple structure. Generalizing to f-PFD with $\mathit{f}>3$ is more challenging. This paper aims to study the general theory for the f-PFD for any $\mathit{f}\ge 3$. We propose a method for obtaining the confounding frequency vectors for all nonequivalent f-PFDs, and to find the least G-aberration (or highest D-efficiency) f-PFD constructed from any single flat. PFDs are particularly useful for constructing nonregular fraction, split-plot or randomized block designs. We also characterize the quaternary code design series as PFDs. Finally, we show how designs constructed by concatenating regular fractions from different families may also have a parallel flats structure. Examples are given throughout to illustrate the results.