PSPP Utilities

About this Page

This page shows how to develop utilities for PSPP using its matrix sub-language. It currently includes programs for Mahalanobis Distance, Trimmed Mean, Cronbach's Alpha, and Median Absolute Distance.

PSPP includes a sub-language that can do matrix calculations. It is available in psppire's Syntax Editor window and in the command line pspp program. It is a compact environment for numerical computation and for creating matrix files, making it useful for deriving statistics and utilities directly.

This sub-language is also useful for learning matrix math and linear algebra. It provides enough functionality to compute many college-level assignments and for experimenting with algorithms step-by-step.

The examples use PSPP's matrix sub-language, but the programs are short and self-contained. You don't need prior matrix programming experience to follow them.

An Appendix with more detail is provided at the end of this page.

This is the process we'll use:

matrix process

Here, I'll do four short PSPP utility examples: a trimmed mean program, a Mahalanobis Distance program, Cronbach's Alpha probram, and a median absolute deviation (MAD) program. They are not very long, but make good examples and allow hand computation to see how things work. They are intended to introduce you to the PSPP matrix sub-language and show some useful computation.

In the trimmed-mean and MAD examples below, the matrix sub-language is used mainly for scalar computations: determining indices, computing medians, evaluating conditions, and performing loops. The only true matrix operations are reading a column vector with GET, slicing it with indexing, and summing it with CSUM(). This illustrates how Matrix can be used effectively even when the underlying computation is not inherently matrix-based.

How to Compute a Trimmed Mean

The trimmed mean is a robust measure of central tendency. Instead of using all values, a trimmed mean removes a fixed percentage of the smallest and largest observations and computes the mean of the remaining data. This reduces the influence of outliers and produces a more stable estimate of the "typical" value when the dataset contains extreme observations.

In this example, 20% of the values are trimmed–10% from the low end and 10% from the high end. The Matrix program sorts the data, determines how many values to trim, slices the vector accordingly, and computes the mean of the trimmed subset.

Results: 
Reading free-form data from INLINE.
+--------+------+
|Variable|Format|
+--------+------+
|x       |F8.0  |
+--------+------+

Value of k
  1

Value of start
  2

Value of end 
  8

Sorted data
   10
   12
   13
   14
   15
   16
   17
   19
  100

Trimmed data
  12
  13
  14
  15
  16
  17
  19

Trimmed mean
 15.1428571429

Here, k determines the 20% number of values to remove, start gives the starting index in the vector, and end gives the ending index. Thus we remove the 10 and the 100 from the data and then compute the mean.

Matrix Syntax: 
DATA LIST LIST /x .
BEGIN DATA
10
12
13
14
14
16
17
19
100
END DATA.

SORT CASES BY x.

MATRIX.
    GET x.

    COMPUTE k = TRUNC(0.20 * NROW(x)).
    COMPUTE start = k + {1}.
    COMPUTE end = NROW(x) - k.
    PRINT k /TITLE="Value of k".
    PRINT start /TITLE="Value of start".
    PRINT end /TITLE="Value of end".

    COMPUTE x_trim = x(start:end) .

    COMPUTE tmean = CSUM(x_trim) / NROW(x_trim).

    PRINT x /TITLE="Sorted data".
    PRINT x_trim /TITLE="Trimmed data".
    PRINT tmean /TITLE="Trimmed mean".
END MATRIX.

This program shows several things:

  1. DATA LIST LIST was used to read the data and make it available to the matrix language.
  2. Selecting values from a matrix uses parentheses for indexing instead of square brackets.
  3. To sum data in pspp matrix, you use RSUM, CSUM, or MSUM. Since we are reading a vertical vector of values, CSUM was used.
  4. To create "variables" use COMPUTE. You can't just assign to variable names directly.
  5. The NROW() matrix function is used to count the rows.
  6. 0.20 is used as the percentage of values to trim from the data, split equally between beginning and end of the data.

Mahalanobis Distance

Mahalanobis distance measures how far an observation is from the center of a multivariate dataset, taking into account the scale and correlation of the variables. It’s useful when your variables move together — for example, petal length and petal width in the Iris data (Fisher, 1936).

People use it to:

Unlike ordinary Euclidean distance, Mahalanobis distance adjusts for covariance, so it gives a more realistic measure of “how unusual” a point is in the context of all variables at once.

This example uses the Fisher Iris data set. String variables cannot be used in Matrix, so the long Species variable is ignored. Also the ID variable has not been used since it doesn't represent iris data. To more easily get the data into the Matrix subsystem, the iris.csv file has been converted to a pspp system data file.

Matrix Syntax 
GET FILE="Iris.sav".

MATRIX.

  GET x /VARIABLES=SepalLengthCm SepalWidthCm PetalLengthCm PetalWidthCm .

  COMPUTE n = NROW(x).
  COMPUTE k = NCOL(x).

  * Column means.
  COMPUTE mu = CSUM(x) / n.

  * Center.
  COMPUTE one = MAKE(n,1,1).
  COMPUTE Xc = x - one * mu.

  * Column variances.
  COMPUTE var = MAKE(1,k,0).

  LOOP j = 1 TO k.
    COMPUTE col = Xc(:,j).
    COMPUTE var(1,j) = (T(col)*col) / (n - 1).
  END LOOP.

  COMPUTE sd = SQRT(var).

  * Standardize.
  COMPUTE denom = one * sd.
  COMPUTE z = MAKE(n,k,0).

  LOOP j = 1 TO k.
    COMPUTE z(:,j) = Xc(:,j) / sd(1,j).
  END LOOP.

  * Covariance of standardized data.
  COMPUTE S = T(z) * z / (n - 1).

  COMPUTE Sinv = INV(S).

  COMPUTE M = z * Sinv * T(z).
  COMPUTE D2all = DIAG(M).

  PRINT D2all /FORMAT=F8.3.

END MATRIX.

This example shows several things about pspp's Matrix sub-language. First, it doesn't have a MEAN function, so that had to be done as CSUM(x)/NROW(x). PSPP Matrix also lacks a COV() function for the the covariance. That was also done using matrix commands. Finally, the data were standardized. (Note: SPSS's Matrix implementation doesn't have MEAN or COV functions, either, according to the manual.)

By performing all computations inside Matrix, the same missing‑data rules apply consistently to both the mean vector and the covariance matrix. Matrix uses strict listwise deletion: any case with a missing value in any variable is excluded entirely from the analysis. This produces a covariance matrix based on a single, coherent set of complete cases, which is numerically more stable and reliably invertible than matrix files constructed using pairwise deletion or AGGREGATE‑based macros.

To use the the Mahalanobis distance matrix program yourself, add your own data set where Iris.sav is, before the Matrix command. Put your (numeric) variable names in place of the ones there in the example (SepalLengthcm, etc.). Then run the program with pspp.

Output

Mahalanobis distances for the 150 iris cases, wrapped into columns for readability.

D2All
   2.152          2.810          2.078          2.447          2.509
   3.957          2.907          1.846          3.359          2.338
   3.314          2.816          2.566          3.588          8.769
   9.836          5.806          2.339          4.530          3.518
   2.622          2.980          3.674          2.244          5.373
   2.399          1.784          2.210          2.496          2.473
   1.966          4.590          8.537          7.346          2.338
   3.398          5.326          2.338          3.143          1.897
   2.570          11.394         3.337          3.737          4.595
   2.956          4.535          2.368          3.049          1.921
   4.446          .630           3.014          3.608          2.043
   3.324          1.290          4.286          2.725          3.168
   7.647          .434           7.512          1.548          1.040
   2.784          3.496          3.296          7.441          2.154
   3.269          1.358          2.379          4.495          1.599
   2.334          4.186          1.423          .319           2.086
   2.786          3.093          .845           2.517          6.254
   2.817          1.795          6.557          1.682          1.996
   4.951          1.188          1.353          4.545          1.485
   2.705          1.190          .648           4.831          .629
   8.888          2.967          2.447          3.568          2.342
   6.066          10.141         7.588          3.327          7.285
   2.249          1.241          2.432          4.769          11.387
   5.843          1.981          12.860         7.211          4.381
   4.018          5.350          8.799          1.811          2.880
   5.899          1.455          1.134          1.667          7.066
   5.698          13.140         2.542          2.404          12.887
   9.696          8.146          3.277          1.467          3.492
   5.947          12.466         2.967          3.083          7.401
   9.079          4.078          1.764          7.637          3.491

Interpreting Mahalanobis distance:

For k variables, compare D² to a chi‑square distribution with k degrees of freedom. For 4 variables (as in the Iris example), values above 9.5 are unusual, above 13 are strong outliers, and above 18 are extreme. The iris data set has several values that are considered unusual, bordering on strong outliers, based on this.

Validation was performed using R and the results were the same.

Note: SAS's mahalanobis function uses an internal scaling that does not match the standard textbook definition used in R or this program.

See Johnson & Wichern, Applied Multivariate Statistical Analysis for the formula used here.

Cronbach's Alpha

This statistic is used with survey or scale variables to assess the internal consistency of the items and whether they represent a single concept. The statistic usually is between zero and one. First, create some test data.

Test Data

This illustrates using INPUT PROGRAM and a LOOP to create random data for testing the Cronbach's Alpha Matrix code. This data will change with each new session.

INPUT PROGRAM.
  LOOP id = 1 TO 200.
    COMPUTE base = RV.NORMAL(50,10).
    COMPUTE item1 = base + RV.NORMAL(0,3).
    COMPUTE item2 = base + RV.NORMAL(0,3).
    COMPUTE item3 = base + RV.NORMAL(0,3).
    COMPUTE item4 = base + RV.NORMAL(0,3).
    END CASE.
  END LOOP.
END FILE.
END INPUT PROGRAM.
EXECUTE.
Matrix Syntax 
MATRIX.
  GET x /VARIABLES=item1 item2 item3 item4.

  * Number of cases and variables.
  COMPUTE n = NROW(x).
  COMPUTE k = NCOL(x).

  * Column means.
  COMPUTE ones_n = MAKE(n,1,1).
  COMPUTE mu = (T(ones_n) * x) / n.

  * Expand mean row to n rows. 
  COMPUTE mu_expanded = ones_n * mu. 

  * Center the data.
  COMPUTE xc = x - mu_expanded.

  * Covariance matrix.
  COMPUTE S = (T(xc) * xc) / (n - 1). 

  * Item variances (diagonal).
  COMPUTE itemvars = DIAG(S).

  * Sum of item variances.
  COMPUTE ones_k = MAKE(k,1,1).
  COMPUTE sumvars = T(ones_k) * itemvars.

  * Total variance (sum of all elements in S). 
  COMPUTE totalvar = T(ones_k) * S * ones_k.

  * Cronbach alpha.
  COMPUTE alpha = (k/(k-1)) * (1 - (sumvars/totalvar)).
  PRINT alpha /FORMAT=F8.3.
END MATRIX.
Results 
alpha
     .978

For Cronbach's Alpha values on 0.90 and higher are considered evidence of strong internal consistency among the items, though extremely high values sometimes reflect item redundancy.

Results were verified with other programming languages, for example a SAS program:

proc corr data=dataset_name alpha;
    var item1 item2 item3 item4;
run;

This pspp Matrix program matched the SAS output.

Also note that the RELIABILITY pspp procedure provides Cronbach's Alpha:

PSPP Output:

    Reliability Statistics
╭────────────────┬──────────╮
│Cronbach's Alpha│N of Items│
├────────────────┼──────────┤
│             .98│         4│
╰────────────────┴──────────╯

This Matrix example is included as a demonstration of a complete, working Matrix program. Even though PSPP provides Cronbach's Alpha, the example show how to compute it manually and illustrates how Matrix handles the underlying statistical operations.

Next, a Matrix program for computing Median Absolute Deviation is shown.

How to Compute Median Absolute Deviation (MAD)

The median absolute deviation (MAD) is a robust measure of spread that uses absolute deviations from the median, making it resistant to outliers.

The data represent wait times at a customer service desk. Most issues appear to be solved in a few minutes, but one person spent over half an hour at the desk. Thus our data has a large outlier which would influence the computation of mean.

This example computes the strict definition of MAD and not the SD-scaled version, which multiplies the MAD value by 1.4826 to make it consistent with the standard deviation under normality. The program as shown outputs the unscaled MAD value; modifying it for the SD-scaled version should be easy if you need it.

This uses DATA LIST FREE to read the data. After the data are read, the values that we don't need (like -1) are removed. Then the cases are sorted for computing the median value. This program checks to see if we have and even number of row or odd numbers of row values and then computes the median based on that with basic math. Function used include MOD(), NROW(), and ABS(). There is no matrix sort that we can use so a bubble sort is added to the program.

Results: 
Median Wait Time (6+4/2 = 5 for 9 values once the -1 is removed)
  5

Median Absolute Deviation (MAD)
  1
Matrix Syntax: 
DATA LIST FREE /wait.
BEGIN DATA
4 5 6 5 4 6 5 5 -1 32
END DATA.

* Remove walk-offs (-1).
SELECT IF wait <> -1.

* Sort wait times.
SORT CASES BY wait.

MATRIX.
    GET wait /VARIABLES=wait.

    COMPUTE n = NROW(wait).

    * Determine odd/even sample size.
    COMPUTE is_even = (MOD(n,2) = 0).

    * Median calculation.
    DO IF is_even = 1.
        COMPUTE med = (wait(n/2) + wait(n/2 + 1)) / 2.
    ELSE.
        COMPUTE med = wait((n+1)/2).
    END IF.

    * Absolute deviations.
    COMPUTE dev = ABS(wait - med).

    * Bubble sort deviations.
    LOOP i = 1 TO n-1.
        LOOP j = 1 TO n-i.
            DO IF dev(j) > dev(j+1).
                COMPUTE temp = dev(j).
                COMPUTE dev(j) = dev(j+1).
                COMPUTE dev(j+1) = temp.
            END IF.
        END LOOP.
    END LOOP.

    * MAD calculation with odd/even branching.
    COMPUTE is_even_dev = (MOD(n,2) = 0).

    DO IF is_even_dev = 1.
        COMPUTE mad = (dev(n/2) + dev(n/2 + 1)) / 2.
    ELSE.
        COMPUTE mad = dev((n+1)/2).
    END IF.

    PRINT med /TITLE="Median Wait Time".
    PRINT mad /TITLE="Median Absolute Deviation (MAD)".
END MATRIX.

Verifying PSPP MAD Results

Independent numberical tools or external calculations were used only to validate that the PSPP Matrix implementation produces the correct strict MAD value. They are included here for completeness and transparency. PSPP users do not need these external tools—the Matrix program shown above computes the strict MAD directly.

I. Independent verification code:

mad(c(4,4,5,5,5,5,6,6,32))
x <- c(4, 4, 5, 5, 5, 5, 6, 6, 32)
median(abs(x - median(x)))

Output:

> median(abs(x - median(x)))
[1] 1

II. Independent verification code:

import numpy as np

x = np.array([4, 4, 5, 5, 5, 5, 6, 6, 32])
mad_strict = np.median(np.abs(x - np.median(x)))
print(mad_strict)

Output:

1.0

III. Independent verification code:

data x;
  input value;
  datalines;
  4
  4
  5
  5
  5
  5
  6
  6
  32
  ;
run;

/* Compute median */
proc sql noprint;
  select median(value) into :med from x;
  quit;

  /* Compute absolute deviations */
data dev;
  set x;
  absdev = abs(value - &med);
run;

/* Compute strict MAD (median of abs deviations) */
proc sql noprint;
  select median(absdev) into :mad from dev;
  quit;

  %put Strict MAD = &mad;

Output:

100        %put Strict MAD = &mad;
 Strict MAD =        1

These examples demonstrated how PSPP's Matrix facility can be used to compute custom statistics using scalar arithmetic, branching, and matrix operations. Users can adapt these patterns to implement trimmed means, percentiles, quartiles, and other robust statistics like median absolute deviation (MAD).

Appendix: Matrix Sub-Language Reference

How to Use Matrix Inside PSPP

Although Matrix is designed for matrix algebra, some programs use it primarily for scalar arithmetic. PSPP treats scalars as 1x1 matrix files, so functions such as NROW(), MOD(), ABS(), arithmetic operators, relational operators, and branching constructs all work naturally on scalar values. In practice, Matrix often serves as a compact procedural language for computing custom statistics, even when the calculations involve only simple scalars rather than full matrix operations.

PSPP Matrix Capabilities

The PSPP Matrix sub-language is a general matrix computation and programming environment. It supports control flow (DO IF, ELSE IF, ELSE, END IF, LOOP, BREAK) and numerical operations such as modular arithmetic (MOD). It provides matrix algebra (addition, subtraction, multiplication, transposition via TRANSPOS, inversion via INV and GINV) and advanced linear-algebra routines including singular value decomposition (CALL SVD), Cholesky factorization (CHOL), solving linear systems (SOLVE), and computing eigenvectors and eigenvalues (CALL EIGEN).

It operates independently of the active dataset. Instead of working with variables and cases, it works directly with matrix files, including row vectors, column vectors, and scalars. Only short strings (8 bytes) are supported. A "variable" is simply an object created with COMPUTE, and indexing uses parentheses rather than square brackets.

Note: PSPP has the FLIP command that transposes variables and cases but that works on the active data set. It is not a Matrix command like TRANSPOS or T, which switch rows and columns within the Matrix environment.

How to Create Matrix Files

Matrix files can be created in several ways. The most common is using GET to pull variables from the active dataset into Matrix. Matrix files can also be created directly inside the Matrix block using COMPUTE, or read from external files using READ. PSPP also provides a separate Matrix DATA procedure that defines a matrix using inline numeric data between BEGIN DATA and END DATA. Unlike COMPUTE or READ, MATRIX DATA is a top-level command and cannot be placed inside a MATRIX block. Matrix files are a PSPP file format.

How to Use Full Matrix Files -- Example

Up to this point the examples have been run from vector data. Now we try a small full matrix example with a matrix operation that you probably expected (inverse). Note that to create B below the COMPUTE command is needed. The following was run in the psppire GUI on Windows 11.

Matrix Syntax:

DATA LIST LIST /x y.
BEGIN DATA
1 2
3 4
END DATA.

MATRIX.
    GET A /VARIABLES=x y.
    COMPUTE B = INV(A).
    PRINT B.
END MATRIX.

Results:

DATA LIST LIST /x y.
Reading free-form data from INLINE.
╭────────┬──────╮
│Variable│Format│
├────────┼──────┤
│x       │F8.0  │
│y       │F8.0  │
╰────────┴──────╯
BEGIN DATA
1 2
3 4
END DATA.

MATRIX.
    GET A /VARIABLES=x y.
    COMPUTE B = INV(A).
    PRINT B.
    COMPUTE C = INV(B).
    PRINT C.
B
 -2.0000000000  1.0000000000
  1.5000000000  -.5000000000

C
  1.0000000000  2.0000000000
  3.0000000000  4.0000000000
END MATRIX.

Note that we get back the original matrix if the first inversion is inverted (invert the inversion).

For additional reference, the PSPP User Manual includes a full section on the Matrix sub-language, covering matrix, matrix data, matrix files, and mconvert.

[PSPP Syntax] [PSPP Analysis-Exploratory] [PSPP Analysis-Inferential & Modeling] [PSPP Merging] [PSPP Multi-Record] [[PSPP Utils]] [PSPP Graphics]

If you have suggestions, comments, or corrections, you can open an issue on the Github repository Issues list