J-ottings 28: Just say it in J
Vector 17:4, page 31
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 31
An early claim for APL was that the exactness demanded by expressing mathematics in an executable language could lead to significantly greater understanding of underlying algebra and algorithms than was achievable by conventional mathematical notation. Subsequent experience diluted this claim somewhat since the obscurity which even a few lines of code could generate outweighed the inherent discipline of using a programming language as a learning tool. Nevertheless there are some topics in which the expository power of APL, and now even more strongly J, trounces traditional methods. One of the most conspicuous examples of this is Analysis of Variance.
The acronym ANOVA, which is at least as well known by its abbreviation as in its full form, is one of the curiosities of statistical terminology, since it is not a single technique, it is not primarily about variance, and by no stretch of the imagination can it be called a process of analysis. I surmise that many generations of students have tackled ANOVA as an arcane piece of algebra involving mysterious "correction factors", seemingly arbitrary rules concerning denominators, and mantras about "degrees of freedom", "interactions" and "error sums of squares". Yet if ANOVA is seen as series of transformations to numbers which are then later formalised in algebra, then both its rationale and the algebra itself become plain. I am not trying to underestimate the vast scope of the intellectually demanding science of Design of Experiments, but rather to suggest that its foundations might be more clearly presented to beginners through the medium of J. I imagine that if ANOVA were to be reinvented today it would probably be called something like SSP (standing for Sum of Squares Partitioning).
The insight to be obtained from using J does of course assume competence in that language, so that as a starting point the following basic verbs together with the adverb each=.&> can be taken for granted:
mean=.+/ % # NB. Mean of vector
dev=.-mean NB. Deviations from mean
ss=.+/@*: NB. Sum of squares of a vector
ssm=.ss@dev NB. Sum of squares about the mean
First tackle one-way ANOVA which is relevant where there are groups of numbers, not necessarily of equal size, but of sufficient general similarity for it to be interesting to speculate on whether the groups are distinguishably different, or whether any variability present is scattered across all the groups. In order to keep the calculations checkable by mental arithmetic the values used in the following examples are integers. A set of groups of numbers is modelled by a list of lists such as
A common sense way of assessing whether this result shows real differences between groups is to hypothesise no variability within each group, and smooth out any such variability by replacing every value by its group mean. J obliges with the primitive verb copy which, operating as the centre tine of a fork, makes just as many copies of the numbers as are necessary:
(There are redundant parentheses in the above expression in order to make the fork structure stand out.) The sum of squares of these values about the overall mean is totally free of any within-group variation and so, rationally, is a sound measure of the between-group variation:
For variability within groups, a corresponding procedure is to sum the individual sums of squares about the mean within each group:
Given the above results, any subsequent analysis is clear cut, viz. the within group variability overwhelms the between group variability. As a consequence of algebra, which can be investigated separately later, it happens that the sum of bss and wss necessarily equals the sum of squares about the mean which would have arisen with the same set of numbers if there had been no group divisions in the first place. Using raze in, this is:
Hence
aov1=.bss;wss;tss NB. one-way ANOVA
aov1 data1
┌───────┬─────┬───────┐
│3.91667│68.75│72.6667│
└───────┴─────┴───────┘
At this point it is useful to consolidate the verbs which describe this first stage in understanding ANOVA :
bss=.ssm@(#each # mean each) NB. Between sum of sqq
wss=.+/@:(ssm each) NB. Within sum of sqq
tss=.ssm@; NB. Total sum of sqq
aov1=.bss;wss;tss NB. 1-way ANOVA
Try them out again with another example:
In this case there are the same numbers of items within each group, and so items could also be classified by position within group. One way to do this is to define markers for rows and cols (so named because >d2 is a matrix):
In deciding to proceed in this way, the model has subtly changed its character to one where there is on the one hand a single list of values with no inherent defined structure, and on the other a collection of value-free data which defines possible structures within the data proper. (Anyone who has heard a statistician talking about a "design matrix" will see this concept beginning to emerge.)
Using row markers, the result of aov1 is obtained by using the key adverb qualifying the from verb for the purposes of selection, together with a redefinition of the data in a structure-free form :
data2=.;d2
group=.{/.
rows group data2
┌─────┬─────┬─────┐
│6 2 3│2 8 5│5 6 8│
└─────┴─────┴─────┘
cols group data2
┌─────┬─────┬─────┐
│6 2 5│2 8 6│3 5 8│
└─────┴─────┴─────┘
The between groups sums of squares can be found exactly as before:
For convenience the verbs bss and group can be combined as
A further extension to this verb allows the sums of squares to be adjusted simultaneously for any number of design vector(s) as left argument:
The means of these subgroups can be obtained conveniently by a similar small extension to the mean verb:
By developing dev each data-item can be transformed into a deviation from the mean of the group to which it belongs, with ravel applied in order to restore it to single vector form:
deveach=.,@(group(-each)means)
rows deveach data2
2.33333 _1.66667 _0.6666667 _3 3 0 _1.33333 _0.3333333 1.66667
The definition of deveach is not quite complete since it is desirable that the deviations are restored to their positions in the original data vector. Restoring the right order involves using an inverse permutation obtainable through two applications of the grade-up verb :
restore=./:@/:@[ { ]
deve=.[ restore,@(group -&> means)
cols deve data2
1.66667 _3.33333 _2.33333 _2.33333 2.66667 _0.3333333 0.6666667 0.6666667 2.66667
The total sum of squares is
A useful check can be made here since the sums of squares of the row-mean adjusted data:
is just the difference between the total sum of squares and the between-rows sum of squares.
Suppose now that this data had a further "hidden" classification based on a so-called Latin square design whereby a third factor (treatment) was incorporated using the matrix
to overlay the data in a row by column layout.
latin=.0 1 2 2 0 1 1 2 0
latin sq data2
24.6667
latin means data2
7.33333 4 3.66667
(rows;cols;latin) sqq data2
10.6667 2 24.6667
So far 37.333 out of the total of 42 sums of squares have been accounted for. What remains is the sum of squares of the residuals, often called the error sum of squares, which is the sum of squares after the various adjustments have been made in succession:
The form of the above line invites generalisation to an insert adverb:
>(deve&.>)/(latin;cols;rows;data2)
0.6666667 _1 0.3333333 _1 0.3333333 0.6666667 0.3333333 0.6666667 _1
The items in the above vector are the result of successive row, column and treatment smoothing; squaring and summing them leads to a definition in J of the residual sum of squares:
and so the complete partitioning, including the total sum of squares, is given by
The order in which the items in the left argument appear is, of course, immaterial.
The property that the tail of the result is the sum of the curtail (or in plain English the last value is the sum of all the others!) depends on orthogonality in the design. For example a design vector which consists of an arbitrary arrangement of three 0s, three 1s, and three 2s such as
does not possess this property:
nor, more generally, do designs which, although balanced, are not orthogonal. However, such considerations go far beyond the scope of this elementary article which seeks to do no more than break the ground for on the one hand analysis of what such partitioned sums of squares can say about variance, and on the other, extension to increasingly complex structures. This is an appropriate point to consolidate the verbs which are part of the second stage of expounding ANOVA :
group=.{/. NB. .. by design vector
sqq=.(sq each)< NB. Between sum of sqq
sq=.bss@group
means=.mean each@group NB. Group means
rss=.sq@>@(deve&.>/@(,<)) NB. Residual sum of sqq
deve=.[ restore,@(group(-each)means)
restore=./:@/:@[ { ]
aov2=.sqq,rss,sq@] NB. Multi-way ANOVA
The above development involves many J primitives – what is striking is that each of them corresponds to an operation which is a logically necessary step for understanding the underlying process of ANOVA in an abstract way. It could also be argued that for someone who already has a clear understanding of ANOVA, the foregoing discussion might lead towards a fuller appreciation of J!