# Repeated measures ANOVA vs. factorial ANOVA with subject factor: understanding "error strata" and Error() term in aov

Consider repeated measures ANOVA (RM-ANOVA) with one within-subject factor A and several measurements per subject for each level of A.

It is closely related to two-way ANOVA with two factors: A and subject. They use identical decomposition of the sum of squares into four parts: A, subject, A⋅subject, and residual. However, two-way ANOVA tests the effect of A by comparing SS of A with the residual SS, while RM-ANOVA tests the effect of A by comparing SS of A with the A$\cdot$subject interaction SS.

Why the difference?

1. Does this difference automatically follow from the repeated-measures structure of the data, or is it some convention?
2. Does this difference between two-way ANOVA and RM-ANOVA correspond to testing two different nulls? If so, what exactly are they and why would we use different nulls in these two cases?
3. Two-way ANOVA's test can be understood as an F-test between two nested models: the full model, and the model without A. Can RM-ANOVA be understood in a similar way?

(If there is only one measurement per subject for each level of A, then the distinction sort of disappears because A$\cdot$subject and residual variation cannot be disentangled: Is one-way repeated measures ANOVA equivalent to a two-way ANOVA?)

Demonstration

I will use toy data d2 generated in http://dwoll.de/rexrepos/posts/anovaMixed.html. The same webpage shows correct syntax for RM-ANOVA.

# Discarding between-subject factors and leaving only one within-subject factor
d = d2[d2$Xb1=='CG' & d2$Xb2 == 'f', c(1,4,6)]


(See reproducible version here on pastebin.) Data look like that:

     id Xw1     Y
1    s1   A  28.6
2    s1   A  96.6
3    s1   A  64.8
4    s1   B 107.5
5    s1   B  77.3
6    s1   B 120.9
7    s1   C 141.2
8    s1   C 124.1
9    s1   C  88.0
10   s2   A  86.7
...


Here is two-way ANOVA: summary(aov(Y ~ Xw1*id, d))

             Df Sum Sq Mean Sq F value   Pr(>F)
Xw1           2  95274   47637  16.789 3.73e-07 ***
id           19  31359    1650   0.582    0.913
Xw1:id       38  71151    1872   0.660    0.929
Residuals   120 340490    2837


Here is RM-ANOVA: summary(aov(Y ~ Xw1 + Error(id/Xw1), d))

Error: id
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 19  31359    1650

Error: id:Xw1
Df Sum Sq Mean Sq F value   Pr(>F)
Xw1        2  95274   47637   25.44 9.73e-08 ***
Residuals 38  71151    1872

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 120 340490    2837


Note the identical SS decomposition, but two-way ANOVA tests Xw1 against the residual, while RM-ANOVA tests Xw1 against the Xw1:id interaction.

Why?

This question is related to How to write the error term in repeated measures ANOVA in R: Error(subject) vs Error(subject/time). If we try using Error(id) instead of Error(id/Xw1) in the example above, then Xw1 will get tested against Xw1:id interaction lumped together with the residual variation.

(The same issue arises in factorial RM-ANOVA with multiple within-subject factors, where each factor or interaction is tested against its own "error term" aka "error stratum". These error strata are always given by the corresponding interaction with the block/plot/subject variable id.)

• Relevant thread: r.789695.n4.nabble.com/AOV-and-Error-td865845.html - but no real answer there. Jun 20, 2017 at 12:30
• Okay, I re-read @JakeWestfall's paper jakewestfall.org/publications/JWK.pdf and realized that the whole issue boils down to RM-ANOVA treating subject effect (and all its interactions!) as random, whereas 2-way ANOVA treats it as fixed. I have to think more about it to figure out all the details. Jun 21, 2017 at 8:05
• For point (2), the null hypothesis is exactly what makes the ratio of the expected mean squares of the corresponding two sum of squares equal to one and the noncentrality parameter corresponding to both sum of squares equal to 0. This it so that the $p$ value for the $F$ statistic is computable. It isn't currently clear to me why we can accomplish all 3 of these goals in the nulls we're clasiccaly used to seeing in ANOVA, but it seems like we only need to focus on the ratio of the EMS when the effects are random and the numerator SS noncentrality parameter when the (numerator) effect is fixed. Jun 23, 2017 at 17:15
• These comments relate to cochran's theorem (en.wikipedia.org/wiki/Cochran%27s_theorem). (The book I use as an ANOVA reference calls this "Bhat's Lemma", by the way.) Jun 23, 2017 at 17:19
• Similar question here, Understanding the split plot, but no terrific answer there yet either Jun 28, 2017 at 2:52

... two-way ANOVA tests the effect of A by comparing SS of A with the residual SS, while RM-ANOVA tests the effect of A by comparing SS of A with the A⋅subject interaction SS.

1) Does this difference automatically follow from the repeated-measures structure of the data, or is it some convention?

It follows from the repeated-measures structure of the data. The basic principle of analysis of variance is that we compare the variation between levels of a treatment to the variation between the units that received that treatment. What makes the repeated measure case somewhat tricky is estimating this second variation.

In this simplest case, the thing we're interested in are the differences between the levels of A. So how many units have we measured that difference on? It's the number of subjects, not the number of observations. That is, each subject gives us an additional independent piece of information about the difference, not each observation. Adding more repeated measures increases the accuracy of our information about each subject, but doesn't give us more subjects.

What the RM-Anova does when using the A--subject interaction as the error term is to correctly use the variation in differences between levels of A between subjects as the variation to test the A level effect. Using the observational error instead uses the variation in the repeated measures on each individual, which is not correct.

Consider a case where you take more and more data on just a couple individuals. If using the observation level error, you would eventually reach statistical significance, even though you only have a couple individuals. You need more individuals, not more data on them, to really increase the power.

2) Does this difference between two-way ANOVA and RM-ANOVA correspond to testing two different nulls? If so, what exactly are they and why would we use different nulls in these two cases?

Nope, same null hypothesis. What's different is how we estimate the test statistic and its null distribution.

3) Two-way ANOVA's test can be understood as an F-test between two nested models: the full model, and the model without A. Can RM-ANOVA be understood in a similar way?

Yes, but not perhaps in the way you're hoping for. As you see in the output from aov, one way of thinking about these kinds of models is that they're really several models in one, with one model for each level.

One can fit the models for higher levels individually by averaging the data over the lower levels. That is, an RM-Anova test for A is equivalent to a standard Anova on the averaged data. Then one can compare models in the usual way.

> library(plyr)
> d2 <- ddply(d, ~Xw1 + id, summarize, Y=mean(Y))
> a1 <- aov(Y ~ id, d2)
> a2 <- aov(Y ~ Xw1+id, d2)
> anova(a1, a2)
Analysis of Variance Table

Model 1: Y ~ id
Model 2: Y ~ Xw1 + id
Res.Df   RSS Df Sum of Sq      F    Pr(>F)
1     40 55475
2     38 23717  2     31758 25.442 9.734e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Alternatively, one can fit the full aov with all the data but without the term of interest, and then compare the fit with the full aov with the term of interest, but then to compare models you need to pick out the level of the model you've changed (here the id:Xw1 level) and then you can compare those two models.

> summary(aov(Y ~ 1 + Error(id/Xw1), d))

Error: id
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 19  31359    1650

Error: id:Xw1
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 40 166426    4161

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 120 340490    2837
> (F <- ((166426 - 71151)/2) / (71151/38))
 25.44202
> pf(F, 2, 38, lower=FALSE)
 9.732778e-08

• (+1) Thanks for taking the time to write this up! This is an interesting perspective that allows us to gain some intuition about why it's natural to compare against the interaction sum of squares in the repeated measures case. However, it seems to fail in elucidating the details of the testing, since you wrongly (according to the arguments in my response) claim that the null hypotheses are the same. The last paragraph of my answer writes what I derived the null hypotheses to be. Please let me know if you think I'm mistaken! Jun 28, 2017 at 15:27
• I think we need to distinguish between what's being tested and what's an assumption of the null hypothesis (which is part of what I mean when I say the null distribution is different). The σ^2_{id∗Xw1}=0 you have is not actually being tested, you can have data where that's not at all true but if X_{w1j} is exactly equal to 0 for all j then you won't reject the null. Jun 28, 2017 at 17:07
• The question is, what do you conclude when you reject the null? In either case, you're concluding that you have evidence that the group means are different. You don't conclude that either the group means are different OR the variance is large. That is, the null hypothesis in both cases is simply that all the group means are the same. What changes is the test statistic we use to test that and the distribution of that test statistic. Jun 28, 2017 at 19:29
• I've realized I'm confused by your whole line of reasoning. A null hypothesis isn't derived, it's simply stated apriori, and then one chooses a test statistics and determines its distribution under the null. In both these cases, the null hypothesis is simply that all the group means are equal. Jun 29, 2017 at 1:59
• @Aaron In chat, amoeba kindly pointed out that I seem to have misunderstood your response to question 2. I interpreted you as saying that in the repeated measures case, the null hypotheses corresponding to the test statistics with MSE in denom or MS_inter in the denom are the same. (Indeed, my final paragraph that I pointed you to is in the setting of repeated measures.) However, it seems now like that was not what you were saying. My mistake! amoeba and I deleted our comments to keep this from misleading future readers. Jun 29, 2017 at 13:11

This note depends on results contained in Moser's Linear Models: A Mean Model Approach. I'll cite some results from this book in what follows. When I saw your question, I started looking through the book: this note is just the way my thoughts were organized after.

Let $y \sim \mathcal{N}_n(\mu, \Sigma)$ be the response, with $\mu$ containing the fixed effects and $\Sigma$ containing the random effects.

Take $y^T A_i y$ to be the sums of squares corresponding to each term (covariates and interactions) in the model. Note that these sums of squares are invariant to whether terms are fixed or random. Assume that each $A_i$ is symmetric and idempotent, which will be true in most models of interest.

When it holds that $$I = \sum_i A_i,$$ which amounts to the sums of squares corresponding to a decomposition into orthogonal subspaces since we've assumed the $A_i$ are projectors, and $$\Sigma = \sum_i c_i A_i,$$ by Cochran's theorem (lemma 3.4.1), \begin{equation} y^T A_i y \sim c_i \chi^2_{d_i} (\mu^T A_i \mu/c_i), \end{equation} for $d_i = tr(A_i)$, and $y^T A_j y$ is independent of $y^T A_k y$ for $j \ne k$.

The term $$\tilde{F} = \frac{y^T A_j y / d_j}{y^T A_k y / d_k} \sim \frac{c_j \chi_{d_j}^2(\mu^T A_j \mu / c_j)/d_j}{c_k \chi_{d_k}^2(\mu^T A_k \mu / c_k)/d_k}$$ is indeed a (central) $F$ statistic if and only if \begin{align*} \frac{c_j}{c_k} & = 1 , \tag{1} \\ \mu^T A_j \mu & = 0 , \tag{2} \\ \mu^T A_k \mu & = 0 , \textrm{ and }\tag{3} \end{align*} When these three conditions are satisfied, we can compute $p$-values corresponding to the statistic $\tilde{F}$. These terms basically just aid in computability since the $c_i$'s depend on variance components and the the noncentrality parameters depend on the mean $\mu$. The second condition ensures that $\tilde{F}$ will have (at least) a noncentral $F$ distribution. Under the second condition, the third condition gives that $\tilde{F}$ has a central $F$ distribution.

The expected mean squares ($\mathrm{EMS}$) corresponding to the $i^\mathrm{th}$ sum of squares $y^T A_i y$ is $$\mathrm{EMS}_i := \frac{1}{tr(A_i)} \mathbb{E} [y^T A_i y] = \frac{tr(A_i \Sigma) + \mu^T A_i \mu}{tr(A_i)} = c_i + \frac{\mu^T A_i \mu}{tr(A_i)},$$ where $tr(A_i\Sigma) = c_i \, tr(A_i)$ due to cor 3.1.2. The ratio $$\frac{\mathrm{EMS}_j}{\mathrm{EMS}_k} = \frac{c_j + \frac{\mu^T A_j \mu}{tr(A_j)}}{c_k + \frac{\mu^T A_k \mu}{tr(A_k)}} = 1$$ if conditions $(1)$, $(2)$, and $(3)$ hold. This is why people inspect the ratio of $\mathrm{EMS}$ when determining which sums of squares to divide to form a $F$ statistic to test a particular null hypothesis.

We use conditions $(1), (2)$, and $(3)$ to specify the null hypothesis. In my experience, when the term (corresponding to $j$) that we're interested in testing is random, we make the null hypothesis be $c_j/c_k = 1$, and, when it's fixed, we make the null hypothesis be $y^T A_j y = 0$. In particular, these amount to us being able to choose $k$ so that the rest of conditions $(1), (2)$ and $(3)$ are satisfied. Such a choice of $k$ isn't always possible, which leads to Behrens-Fisher-like difficulties.

This doesn't explain anything particularly related to the problem at hand, but that just amounts to computing $\mu$ and $\Sigma$. I hope this is seen to be a useful way of thinking about the problem. Note that example 4.4.1 works out what all of the quantities above are in the two-way ANOVA example.

The difference is due to the problem structure and not due to convention. These different approaches (two-way vs repeated measure) change $\mu$ and $\Sigma$, which changes the EMS, which changes which $k$ we choose to construct the test.

Let's consider the model \begin{equation} y_{ijk} = \mu_0 + \mathrm{id}_i + \mathrm{Xw1}_{j} + \mathrm{id * Xw1}_{ij} + \mathrm{R(id * Xw1)}_{k(ij)}, \end{equation} where $i$ denotes the level of $\mathrm{id}$, etc. Here $k$ denotes which of the 3 replicates are being considered.

We now introduce some helpful vector notation: write $y = (y_{111}, y_{112}, y_{113}, y_{121}, \dots y_{20,3,3})$. Since this data is balanced, we can make us of kronecker product notation. (As an aside, I was told that Charlie Van Loan once called the kronecker product "the operation of the 2000s!") Define $\bar{J} \in \mathbb{R}^{m \times m}$ to be the matrix with all entries equal to $\frac{1}{m}$ and $C=I-\bar{J}$ to be the centering matrix. (The centering matrix is so named since, for instance, $\|C x \|_2^2 = \sum_i (x_i - \bar{x})^2$ for a vector $x$.)

With this kronecker product notation under out belt, we can find the matrices $A_i$ mentioned above. The sum of squares correspoding to $\mu_0$ is \begin{equation} SS(\mu_0) = n (\bar{y}_{\cdot\cdot\cdot})^2 = \|(\bar{J} \otimes \bar{J} \otimes \bar{J}) y\|_2^2 = y^T (\bar{J} \otimes \bar{J} \otimes \bar{J}) y, \end{equation} where the first component $\bar{J} \in \mathbb{R}^{20 \times 20}$, the second is in $\mathbb{R}^{3 \times 3}$, and the third is in $\mathbb{R}^{3 \times 3}$. Generally speaking, the matrices in those components will always be of that size. Also, the sum of squares due to $\mathrm{id}$ is \begin{equation} SS(\mathrm{id}) = \sum_{ijk} (\bar{y}_{i\cdot\cdot} - \bar{y}_{\cdot \cdot \cdot})^2 = \|(C \otimes \bar{J} \otimes \bar{J}) y\|_2^2 = y^T (C \otimes \bar{J} \otimes \bar{J}) y. \end{equation} Notice that $SS(\mathrm{id})$ does indeed measure the variation among levels of $\mathrm{id}$. Similarly, the other matrices are $A_{Xw1} = \bar{J} \otimes C \otimes \bar{J}$, $A_{id * Xw1} = C \otimes C \otimes \bar{J}$, and $A_{R()} = I \otimes I \otimes C$.

This is shown to be consistent with aov by running code to give, for instance, the residual sum of squares $SS(\mathrm{R(id * Xw1)}) = y^T A_{R()} y$:

mY <- c()
for(j in 1:(nrow(d)/3)) {
mY <- c(mY, rep(mean(d$Y[3*(j-1)+(1:3)]), 3)) } sum((d$Y - mY)^2) #this is the residual sum of squares


At this point, we have to make some modeling choices. In particular, we have to decide whether $\mathrm{id}$ is a random effect. Let's first suppose that it isn't a random effect, so that all effects besides the replication are fixed. Then \begin{equation} \mathbb{E} [y_{ijk}] = \mu_{ij} = \mu_0 + \mathrm{id}_i + \mathrm{Xw1}_{jk} + \mathrm{id * Xw1}_{ij} \end{equation} and $R(\mathrm{id * Xw1})_{k(ij)} \sim_{iid} \mathcal{N}(0, \sigma^2)$. Notice that there's no dependence between distinct observations. In vector notation, we can write $$y \sim \mathcal{N} (\mu, \Sigma)$$ for $\mu = \mathbb{E} [y] = (\mu_{11}, \mu_{12}, \dots, \mu_{20,3}) \otimes \mathbf{1}_{3}$ and $\Sigma = \sigma^2 (I \otimes I \otimes I)$.

Noticing that the sum of all $5$ of the $A$'s defined above is the identity, we know by cochran's theorem that, among other things, $$SS(\mathrm{Xw1}) = y^T A_{Xw1} y \sim \sigma^2 \chi^2_{(19)(1)(1)} (\mu^T A_{Xw1} \mu / \sigma^2)$$ and $$SS(\mathrm{R(id * Xw1)}) = y^T A_{R()} y \sim \sigma^2 \chi^2_{(20)(3)(2)} (\mu^T A_{R()} \mu / \sigma^2)$$ and these sums of squares are independent.

Now, in line with what we discussed above, we want conditions $(1), (2),$ and $(3)$ to hold. Notice that condition $(1)$ holds (because there's no other variance components to complicate things.) What's really cool to notice now is that $\mu^T A_{R()} \mu = 0$, since $\mu$ is constant along this third "component" that is being centered by $A_{R()}$. This means that $(3)$ is behind us. Therefore we only have to fret about condition $(2)$: if we assume it (as a null hypothesis) then we're assuming that $0 = \mu^T A_{Xw1} \mu = \sum_{ijk} (\mu_{ij} - \bar{\mu}_{i \cdot})^2$, which is the same as $\mu_{ij} = \bar{\mu}_{i \cdot}$ for all $i,j$, which is the same as $\mathrm{Xw1}_j = 0$ and $\mathrm{id * Xw1}_{ij} = 0$ for all $i,j$ (since the mean level is in the other terms.)

In summary, the null hypothesis can be seen to just be testing whether a noncentrality parameter is zero, which is equivalent to effects concerning the covariate being zero. The repeated measures case follows a similar line of reasoning, where we instead make the modeling choice that the $\mathrm{id}$ effect is random. There, condition $(1)$ will become the null hypothesis.

Related to the R command, like you mention in the comments to the original post, this error term just specifies which terms are to be considered as random effects. (Note that all terms that are to be included in the model should be plainly input or input inside the Error() term. This is why there's a difference between id/Xw1 = id + id:Xw1 and id being in the Error term. Non-included terms are lumped in with the error in the sense that $A_{R()} + A_{id * Xw1}$ is relabeled as $A_{R()}$.)

Here's the explicit details related to the repeated measures case where the terms related to $\mathrm{id}$ (which are $\mathrm{id}$ and $\mathrm{id * Xw1}$) are random. We'll see that this is the more interesting case.

There we have the same sum of squares matrices (since they don't depend on whether a factor is fixed or random.) The covariance matrix there is \begin{align*} \Sigma & \stackrel{(a)}{=} \sigma^2_{id} (I \otimes J \otimes J) + \sigma^2_{id * Xw1} (I \otimes C \otimes J) + \sigma^2_{R()} (I \otimes I \otimes I) \\ & = \sigma^2_{id} (3)(3) (A_{\mu_0} + A_{id}) + \sigma^2_{id * Xw1} (3) (A_{Xw1} + A_{id * Xw1}) + \sigma^2_{R()} (A_{\mu_0} + A_{id} + A_{Xw1} + A_{id * Xw1} + A_{R()}) \\ & = ((3)(3)\sigma^2_{id} + \sigma^2_{R()})A_{\mu_0} + ((3)(3)\sigma^2_{id} + \sigma^2_{R()}) A_{id} + ((3)\sigma^2_{id * Xw1} + \sigma^2_{R()}) A_{Xw1} + ((3)\sigma^2_{id * Xw1} + \sigma^2_{R()}) A_{id * Xw1} + \sigma^2_{R()} A_{R()}, \end{align*} where $J$ is the matrix of all ones. The first and last summand on the right hand side of equality (a) offer intuitive explanations: the first summand shows that there's an additional source of correlation among observations with the same $\mathrm{id}$, and the third summand shows, as in the two-way example, the base source of variation. This second summand is less intuitive, but among observations with the same \mathrm{id}, it can be seen as increasing variation between observations with same $\mathrm{Xw1}$ while decreasing variation between observations with different $\mathrm{Xw1}$, due to the shape of $I \otimes C \otimes J$.

Also, since all of the terms related to $\mathrm{id}$ are random, the mean is just due to $\mathrm{Xw1}$, so that $\mathbb{E} [y_{ijk}] = \mu_{j} = \mu_0 + \mathrm{Xw1}_j$, or $\mu = \mathbf{1} \otimes (\mu_1, \mu_2, \mu_3) \otimes \mathbf{1}$.

Notice that, related to condition $(1)$: we have $$\frac{c_{Xw1} }{c_{id * Xw1}} = \frac{(3)\sigma^2_{id*Xw1} + \sigma^2_{R()}}{(3)\sigma^2_{id*Xw1} + \sigma^2_{R()}} = 1,$$ while $$\frac{c_{Xw1}}{c_{R()}} = \frac{(3)\sigma^2_{id*Xw1} + \sigma^2_{R()}}{\sigma^2_{R()}} \neq 1.$$ Further, related to condition $(3)$ both $\mu^T A_{Xw1*id} \mu = 0$ and $\mu^T A_{R()} \mu = 0$. Also, related to condition $(2)$: we see that \begin{align*} \mu^T A_{Xw1} \mu & = \|A_{Xw1} \mu\|_2^2 \\ & = \|(\bar{J} \otimes C \otimes \bar{J}) (\mathbf{1} \otimes (\mu_1, \mu_2 \mu_3)' \otimes \mathbf{1}) \|_2^2 \\ & = (20)(3) \|C (\mu_1, \mu_2 \mu_3)'\|_2^2 \\ & = (20)(3) \sum_j (Xw1_j)^2. \end{align*}

Therefore, if the denominator sum of squares was the residual $\mathrm{R(id * Xw1)}$ like before, there would be both conditions $(1)$ and $(2)$ in the null hypothesis---since those are the two conditions that aren't satisfied without assumptions. However, if we were to use denominator sum of squares as the interaction, since condition $(1)$ is already satisfied, the null hypothesis would just be condition $(2)$. So, as you mention in your question, these different denominators just amount to different null hypotheses.

This analysis technique we use allows the choice of which null hypothesis is being tested to be transparent. Indeed, we can see this by writing out the conditions mentioned in the previous paragraph more explicitly. Using the denominator as the residual sum of squares forces us to test $Xw1_j = 0$ for all $j$ and $\sigma^2_{id * Xw1} = 0$, while using the denominator as the interaction sum of squares allows us to simply test $Xw1_j = 0$ for all $j$.

• +1. Wow, thanks a lot. It will take me some time to digest this answer. I am not very familiar with the mathematical theory of hypothesis testing in linear models, so this is a bit hard to understand. I might come back to you with some questions in the following days. I was more expecting to get an answer in the style of the example on pp. 2-3 of this paper jakewestfall.org/publications/JWK.pdf, where expected mean squares are computed in several fixed-vs-random situations and everything follows from there. It looks like you are talking about the same thing, but more formal. Jun 23, 2017 at 21:59
• I've included an example. (They can get pretty long to write out!) I think it takes some time to get comfortable with kronecker product manipulations, but, after that, this is more easily understandable. Also, I keep finding typos in the answer. Please let me know if you think there's any! Jun 27, 2017 at 20:50
• Whew, that's a lot of math! The question seems much more conceptual to me, I'll see if I can find the time add an answer in words. Jun 28, 2017 at 2:53
• @Aaron since amoeba asked for a comprehensive answer and about extending this problem to other scenarios, I figured it would be worthwhile to provide a full explanation of $F$ tests in ANOVA. The answer got notationally heavy just because there's a lot of computation involved when doing it in a fully generalizable way. (Although, to be clear, the most math involved is evaluating the norm of a projected vector.) I'd be very interested to see a more conceptual answer that fully explains the intricacies that I introduced (more than a little) notation to explain. Please do post if you have time! Jun 28, 2017 at 3:58