APL and Nested Arrays — a Dream for Statistical Computation
Vector 14:2, page 58
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 58
Introduction¶
Alan Sykes
Many papers have been produced in the last 10 or more years in the APL community extolling the virtues of APL for statistical computing so why another one? Such papers have rightly stressed APL’s array-handling operations, the use of user-defined operators, the often transparent flow from mathematical notation to APL’s equivalent computational notation, together with specific tools with a statistical bent, like the dyadic domino.
Nested arrays however figure relatively little in this. Of course the current Windows-oriented computing medium means that a function’s argument could be, for example, a list of properties of an object. For example drawing a polygonal line in Dyalog APL requires an incantation such as:
Nested arrays play a crucial role here. Similarly, our use of nested arrays in APL programming tends only to exploit nested array structures for passing complex arrays in and out of user-defined functions.
However, most of statistics involves number crunching, calculating means and standard deviations. This paper discusses a particular statistical application for which nested numeric vectors and matrices are a natural data structure for which means and variances need to be calculated. This application is Multiple Imputation.
Because we will be dealing with numeric vectors or matrices, which have at least one scalar, which when disclosed, is itself a numeric vector, the principle of scalar extension will play a fundamental role in what follows. The following code illustrates this simple idea.
Introducing Multiple Imputation¶
Among the greatest users of statistical techniques are the Market Researchers who often collect information about customers through a survey. Social scientists have used surveys too, many of them huge, to investigate specific social scientific questions relevant to modern day life. The bugbear of all such survey work, no matter how much work is put into preparing the right kind of questions and the best survey design, is non-response. Multiple Imputation [1] is, in the words of its originator (D. Rubin), the technique that replaces each missing value with two or more acceptable values representing a distribution of possibilities.
The need for it stems from the fact that:
- to ignore item non-response is potentially disastrous for subsequent inferences made from the database;
- if single imputation is used, then complete-data methods of analysis may be used and the need for specialist computer programs to handle non-response is avoided;
- single imputation, however well chosen, means that inferences based on the imputed data will be too sharp since the extra variability due to the unknown missing values is not being taken into account.
Multiple imputation proceeds by:
- imputing a modest number m (between 2 and 10) missing values by some appropriate means;
- producing the required statistical analysis using all m different databases;
- using the variation in the m different results to tease out the extra variability due to the unknown missing values.
APL to the Rescue¶
With this scenario in mind, one can begin to see the role APL can play in multiple imputation calculations. Rubin suggests for example that the imputed values ‘are stored in an auxiliary matrix with one row for each missing value and m columns’. This implies straightaway that specialist programs will have to be constructed. For example, calculating a simple mean of a variable that has been imputed, would presumably require a program that loops round each imputation, which it would have to know where to find, in order to produce the m means.
In APL, we do not need any auxiliary matrix, nor do we need a specialist program - an appropriately ordinary MEAN program should suffice! Take a simple example. We have a two-variable database and the Y-variable (first column) has two missing values which have been imputed m=2 times:
| Y | X | |
|---|---|---|
| 1 | 10 | 8 |
| 2 | ? | 9 |
| 3 | 14 | 11 |
| 4 | ? | 13 |
| 5 | 16 | 16 |
| 6 | 15 | 18 |
| 7 | 20 | 6 |
| 8 | 4 | 4 |
| 9 | 18 | 20 |
| 10 | 22 | 25 |
| Repetition One | Repetition Two | |
|---|---|---|
| Unit 2 | 10 | 14 |
| Unit 4 | 16 | 14 |
The important point to realise is that we have two databases here - one with each unknown value replaced by 10 and 16 respectively, and one with them replaced by 14 and 14. Thus if we require the mean of the Y column we must do one calculation substituting the values 10 and 16 for the two queries, and a further calculation substituting the values 14 and 14. In APL, the original database plus the auxiliary matrix of imputed values may be stored as one matrix and means can be calculated directly:
data←10 2⍴10 8 0 9 14 11 0 13 16 16 15 18 20 6 4 4 18 20
22 25
data[2;1]←⊂10 14
data[4;1]←⊂16 14
data
10 8
10 14 9
14 11
16 14 13
16 16
15 18
20 6
4 4
18 20
22 25
R←MEAN X
[1] R←(+⌿X)÷⊖⍴⍴X
MEAN data[;1]
14.5 14.7
MEAN data
14.5 14.7 13
DISP MEAN data[;1]
┌─────────┐
│14.5 14.7│
└─────────┘
DISP MEAN data
┌─────────┬──┐
│14.5 14.7│13│
└─────────┴──┘
As you can see, this works superbly, irrespective of whether we require just the mean of the y-data or require means of both columns. The only proviso, from an APL point of view, is that the number of imputations is the same for each missing value. But this is also required from a statistical point of view so there is no restriction in its application. In a similar fashion, many routine statistical calculations in APL translate to the nested-array scenario painlessly. For example, take the functions SSQ (Sum of squares), SPR (sum of products), SDEV (Standard deviation) given by:
R←SSQ X;T
[1] R←T+.×T←X-MEAN X
R←X SPR Y
[1] R←(X-MEAN X)+.×Y-MEAN Y
R←SDEV X;T
[1] R←((+⌿T×T←X-(⍴X)⍴MEAN X)÷¯1+⊖⍴⍴X)*0.5
It is easy to see that these simple functions (that our students would write as part of their first assignment in APL Statistical Computing) do what we want, so we can use them to calculate the slope and intercept of a regression line of Y on X:
Regression Line Slope:
Regression Line Intercept:
In each case the returning result is an enclosed vector of length 2 because the number of imputations is 2.
Of course, advanced functions of the APL interpreter (e.g. ⌹) will not operate on nested arrays. Nevertheless, as we show in the next section, it is possible to write simple first-generation APL code in such a way that the dyadic use of ⌹ can be achieved, and in addition, that same code will work quite naturally for nested arrays!
APL Code for Dyadic Domino¶
The most common statistical use for the Dyadic Domino is to calculate the least-squares regression coefficients in a multiple linear regression problem. If Y1 is the vector of ‘y-values’ taken from the first copy of Y and D is the design matrix of linear predictors, i.e.
then we need the following important calculations:
- the least-squares estimates of the coefficients of the linear predictors;
- the residual Sum of Squares;
- the covariance matrix of the regression estimates;
B
7.417910448 0.5447761194
RSS
135.1940299
COV
8.794338011 ¯0.5464932799
¯0.5464932799 0.04203794461
(Note that since Y1 is the first copy of Y then the values of B equal the first values of the slope and intercept previously calculated.)
To have this facility for a general design matrix D, which together with Y may contain enclosed vectors of the same length, we appeal to the efficient and elegantly designed Beaton’s algorithm (see [2]) which is similar to Gaussian elimination. It can be coded easily in APL and the coding works without alteration for nested arrays.
R←Y LS X;FL;X;Y;Z;K;TV;T;XTX;B;SS;D;A;B0
[1] K←1 ⋄ Z←X,Y ⋄ Z←(⍉Z)+.×Z ⋄ D←⊖⍴⍴Z
[2] ⍝ do Beaton recipe on sums of squares matrix
[3] ⍝ NB - at the moment it does nothing about checking whether to
[4] ⍝ cope for ill-conditioning if the first column in design
[5] ⍝ matrix is all 1's - it pivots on all columns!
[6] LP:Z←Z-(TV∘.×TV←Z[K;])÷T←Z[K;K]
[7] Z[;K]←Z[K;]←TV÷|T
[8] Z[K;K]←-÷T
[9] →((¯1+1 ⍴Z)≥K←K+1)/LP
[10] ⍝ now allocate separate results for clarity
[11] XTX←-¯1 0↓0 ¯1↓Z ⋄ SS←Z[D;D] ⋄ B←¯1↓Z[;D]
[12] R←B SS XTX
We use nested arrays only to return the three arrays as one, following the specification above - i.e. the first element consists of the regression coefficients, the second equals the residual sum of squares, and the third equals the estimated covariance matrix of the regression coefficients. Again, we illustrate it using Y1 equal to the vector of first y-values in Y and also with Y2 as the vector of the second y-values in Y, suppressing all but two decimal places in the output.
DISP 2 DP Y1 LS D
┌─────────┬──────┬───────────┐
│7.42 0.54│135.19│ 8.79 ¯0.55│
│ │ │¯0.55 0.04│
└─────────┴──────┴───────────┘
Note that the first component gives the intercept and slope values, the next is the regression sum of squares and finally we have the 2 by 2 covariance matrix. Similarly for Y2:
DISP 2 DP Y2 LS D
┌────────┬──────┬───────────┐
│8.14 0.5│129.59│ 8.43 ¯0.52│
│ │ │¯0.52 0.04│
└────────┴──────┴───────────┘
The same code, containing no second generation idioms apart from the output of the three arrays now works for the nested Y vector:
DISP 2 DP Y LS D
┌────────────────────┬───────────────┬─────────────────────────┐
│┌─────────┬────────┐│┌─────────────┐│┌───────────┬───────────┐│
││7.42 8.14│0.54 0.5│││135.19 129.59│││ 8.79 8.43 │¯0.55 ¯0.52││
│└─────────┴────────┘│└─────────────┘│├───────────┼───────────┤│
│ │ ││¯0.55 ¯0.52│ 0.04 0.04 ││
│ │ │└───────────┴───────────┘│
└────────────────────┴───────────────┴─────────────────────────┘
Note that the 2 by 2 covariance matrix is nested because there are two residual sum of squares and not because the design matrix is nested. However, there is no reason, either theoretical or practical why the design matrix should not be nested. Consider replacing the value 25 in the [10;2] position of D by the enclosed vector 25 30:
D[10;2]←⊂25 30
DISP 2 DP Y LS D
┌─────────────────────┬───────────────┬─────────────────────────┐
│┌─────────┬─────────┐│┌─────────────┐│┌───────────┬───────────┐│
││7.42 8.76│0.54 0.44│││135.19 126.76│││ 8.79 6.89 │¯0.55 ¯0.39││
│└─────────┴─────────┘│└─────────────┘│├───────────┼───────────┤│
│ │ ││¯0.55 ¯0.39│ 0.04 0.03 ││
│ │ │└───────────┴───────────┘│
└─────────────────────┴───────────────┴─────────────────────────┘
As is to be expected, all first values of each element are the same as before, whereas second values are different, as for those results the 10th x-value is now 30 and not 25.
Generalized Linear Models¶
As is well documented in previous APL papers, (see e.g.[3]) it is possible using an iteratively re-weighted least-squares procedure to extend the range of validity of regression-type problems to situations where the ‘y-data’ is not normally distributed. A typical example might be where each y-value represents the number of successes out of a number of trials and where the aim is to explain how the probability of success varies with the predictor values. By a direct emulation of the software of ASLGREG (see [5]) for this particular example of a generalized linear model, we can achieve estimates for a ‘logistic regression’ problem.
RES←Y LOGREG X;T;MU;TP;P;N;ETA;DETA;C;WT;WV;BETA;V;D
[1] ⍎(1=⍴⍴Y)/'N←(⍴Y)⍴1'
[2] ⍎(1≠⍴⍴Y)/'N←Y[;2]⋄Y←Y[;1]'
[3] ⍝start
[4] MU←Y ⋄ C←1 ⋄ DEV←⌊/⍳0 ⋄ P←¯1 ↑⍴X
[5] ETA←⍟TP÷1-TP←(TP+0.0001×TP=0)÷N+N=TP←MU
[6] ⍝ now start iteration
[7] LP:DETA←((1+D)*2)÷N×D←*ETA
[8] WV←ETA+(Y-MU)×DETA
[9] WT←N×TP×1-TP
[10] WT←WT*0.5
[11] BETA←1⊃BT←(WT×WV)LS X×(⍴X)⍴P/WT
[12] ETA←X+.×BETA ⋄ MU←N×TP←D÷1+D←*ETA
[13] DEVNEW←+/2×(Y×⍟(Y+Y=0)÷MU+0.001×MU=0)+(N-Y)×⍟(N-Y
-Y=N)÷(N-MU)+0.001×0=N-MU
[14] ⍎(0.00001<⌈/|DEV-⊃DEVNEW)/'C←C+1⋄DEV←⊃DEVNEW⋄→LP'
[15] RES←BETA DEVNEW((3⊃BT)×(-/⍴X)÷2⊃BT)
The left argument is either a vector of 0s and 1s (successes out of 1 trial) or a matrix whose first column contains the number of successes and the second column the corresponding total number of trials. The right argument is a matrix of predictor variables. A simple example shows how both can be nested.
DISP Y
┌───┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
│0 1│0│1│0│1│1│0│1│1│1│
└───┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
X←1,[1.5]⍳10
X[1;2]←⊂1 5
DISP X
┌─┬───┐
│1│1 5│
├─┼───┤
│1│ 2 │
├─┼───┤
│1│ 3 │
├─┼───┤
│1│ 4 │
├─┼───┤
│1│ 5 │
├─┼───┤
│1│ 6 │
├─┼───┤
│1│ 7 │
├─┼───┤
│1│ 8 │
├─┼───┤
│1│ 9 │
├─┼───┤
│1│10 │
└─┴───┘
DISP 3 DP Y LOGREG X
┌──────────────────────────┬─────────────┬─────────────────────────────┐
│┌────────────┬───────────┐│┌───────────┐│┌─────────────┬─────────────┐│
││¯2.29 ¯1.548│0.528 0.444│││9.803 10.32│││ 3.233 3.903 │¯0.545 ¯0.671││
│└────────────┴───────────┘│└───────────┘│├─────────────┼─────────────┤│
│ │ ││¯0.545 ¯0.671│ 0.114 0.136 ││
│ │ │└─────────────┴─────────────┘│
└──────────────────────────┴─────────────┴─────────────────────────────┘
As before, the first element of the result is the nested vector of ‘regression coefficients’. The second vector is now the ‘deviance’ (equivalent to the residual sum of squares) and the final element is the covariance matrix of the regression estimates. These results may be verified by repeated use of ASLGREG. (Note however that accuracy of the results is not as high as it should be - see the comments in the function LS.)
Other Utilities¶
To finish with we look at some other programming issues associated with the use of nested arrays for multiple imputation. (It is hoped to write up the details of an actual application of multiple imputation elsewhere as they are beyond the scope of this paper.)
Identifying rows of a nested matrix¶
Suppose X is a nested matrix (in this case with two columns). We wish to identify each row in terms of its position in a set of unique rows (given by the columns of T). Standard first-generation look-up code works wonderfully once we have found T. For this, we first have to convert X to a non-nested version which includes all possible rows:
DISP X
┌───┬───┐
│1 2│2 │
├───┼───┤
│ 1│1 │
├───┼───┤
│ 2│2 1│
├───┼───┤
│ 1│2 │
└───┴───┘
⊃,/(,¨¨)/X
1 2
2 2
1 1
2 2
2 1
1 2
Then a simple UNIQUE function which identifies unique elements of a vector or unique rows of a matrix will do the job:
Now we can identify which rows of the original (nested) X matrix are which:
The interpretation of this is clear - the first element which is nested tells us that the first row of X (which is 1 2 or 2 2 depending on which copy we take) consists of the 2nd or 4th elements in the unique listing. Similarly for all the other elements.
So indexing the rows of a nested matrix requires converting the nested matrix to a non-nested matrix. Calculating the median of each column of a matrix requires a similar conversion. For example, the function NVTOMAT (NMTOMAT) converts a nested vector (matrix) to a matrix.
∇ r←NVTOMAT v;m
[1] ⍝ converts nested vector to matrix
[2] m←⌈/,⊃¨⍴¨v
[3] r←↑m⍴¨v
∇
∇NMTOMAT∇
∇ r←NMTOMAT mat;m;d
[1] ⍝ converts nested matrix to fatter matrix
[2] d←⍴mat
[3] m←⊃¨⍴¨mat ⋄ m←d⍴⌈⌿m+m=0
[4] r←(d[1],+/m[1;])⍴⊃,/,m⍴¨mat
∇
To illustrate these functions consider the nested vector v and the nested matrix m below:
DISP v
┌─────┬─┬──┐
│1 2 3│9│99│
└─────┴─┴──┘
NVTOMAT v
1 2 3
9 9 9
99 99 99
DISP m
┌─────┬─┐
│1 2 3│2│
├─────┼─┤
│ 3│4│
├─────┼─┤
│ 5│6│
└─────┴─┘
NMTOMAT m
1 2 3 2
3 3 3 4
5 5 5 6
Calculating Medians¶
Norman Thomson in his St Petersburg paper ([4]) refers to the median function as a function which is ‘not dimensionally extendible’ by which he is suggesting that (unlike the MEAN program) it is not possible to program MEDIAN so that it works on an array of any dimension without looping. This is not true as the following program illustrates.
R←MEDIAN X;ORD;ND;IND;D
[1] D←⍴X
[2] ⍎(1=⍴D)/'X←((⍴X),1)⍴X' ⋄ ND←⍴X
[3] IND←ND⍴⍳ND[2] ⋄ X←,⍉X ⋄ IND←,⍉IND
[4] X←X[ORD←⍋X] ⋄ IND←IND[ORD]
[5] X←X[⍋IND]
[6] X←⍉(⌽ND)⍴X
[7] R←0.5+.×X[⌈0.5×0 1+D[1];]
[8] ⍎(1=⍴D)/'R←⊖⍴R'
D
1 8
1 9
1 11
1 13
1 16
1 18
1 6
1 4
1 20
1 25
MEDIAN D
1 12
If we wish to extend the usefulness of the program MEDIAN to nested arrays of the type encountered in multiple imputation, then there is no real way of avoiding converting the nested vector to a matrix, use MEDIAN to calculate the median value for each column, and then convert this back to have the appropriate nested structure.
DISP Y
┌──┬─────┬──┬─────┬──┬──┬──┬─┬──┬──┐
│10│10 14│14│16 14│16│15│20│4│18│22│
└──┴─────┴──┴─────┴──┴──┴──┴─┴──┴──┘
IMEDIAN Y
15.5 14.5
DISP IMEDIAN D,Y
┌─┬──┬─────────┐
│1│12│15.5 14.5│
└─┴──┴─────────┘
r←IMEDIAN arr;i;m;k
[1] ⍝computes median of nested vector/matrix
[2] ⍎(1≥|≡arr)/'r←MEDIAN arr⋄→0'
[3] →(vec,mat)[⍴⍴arr]
[4] vec:r←MEDIAN NVTOMAT arr ⋄ →0
[5] mat:m←⌈⌿m←⊃¨⍴¨arr ⋄
[6] m←m+m=0
[7] m←⊃,/m ¨1
[8] r←m⊂MEDIAN NMTOMAT arr
Conclusion¶
Nested numeric arrays in APL form a natural tool for use in database work involving multiply imputed values. We have been pleasantly surprised how much basic APL code works with the data structures discussed in this paper and how quickly other facilities can be constructed. The simple idea of scalar extension at the heart of both first- and second-generation APL ensures that the statistician can largely operate with nested arrays in the sure knowledge that the results will be what they obviously should be with nested elements containing the results corresponding to each choice of elements in the nested elements of the original array. In short, the illustration of the MEAN program operating on a nested vector or matrix says it all!
References¶
- Rubin, D.B. (1987), Multiple Imputation for Nonresponse in Surveys. Wiley, New York.
- Sykes A.M. and Ansell J.I., Beaton’s Recipe for Least-Squares. Vector Vol. 8, No. 1, July 1991, pp36–44.
- Sykes A.M., Generalized Linear Models in APL, Proceedings of APL86 (Manchester) pp 48–70.
- Thomson N., Some Proposals for APL2 Specification of Statistical Algorithms, APL92 Conference Proceedings (St Petersburg), APL Quote Quad Volume 23 Number 1, July 1992. pp 243–255.
- Morgan E., ASLGREG: Proposed Additional Features, Predictions with Confidence Limits, Vector Vol. 10 No. 3, January 1994, pp 81–93.