A statistical approach to determine if data are missing at random I have a large set of feature vectors which I will use to attack a binary classification problem (using scikit learn in Python). Before I start to think about imputation, I am interested in trying to determine from the remaining parts of the data if the missing data are 'missing at random' or missing not at random.
What is a sensible way to approach this question?

It turns out a better question is to ask if the data is 'missing completely at random' or not.  What is a sensible way to do that?
 A: A method I use is a shadow matrix, in which the dataset consists of indicator variables where a 1 is given if a value is present, and 0 if it isn't. Correlating these with each other and the original data can help determine if variables tend to be missing together (MAR) or not (MCAR). Using R for an example (borrowing from the book "R in action" by Robert Kabacoff):
#Load dataset
data(sleep, package = "VIM")

x <- as.data.frame(abs(is.na(sleep)))

#Elements of x are 1 if a value in the sleep data is missing and 0 if non-missing.
head(sleep)
head(x)

#Extracting variables that have some missing values.
y <- x[which(sapply(x, sd) > 0)]
cor(y)

#We see that variables Dream and NonD tend to be missing together. To a lesser extent, this is also true with Sleep and NonD, as well as Sleep and Dream.

#Now, looking at the relationship between the presence of missing values in each variable and the observed values in other variables:
cor(sleep, y, use="pairwise.complete.obs")

#NonD is more likely to be missing as Exp, BodyWgt, and Gest increases, suggesting that the missingness for NonD is likely MAR rather than MCAR.

A: This sounds quite doable from a classification standpoint.
You want to classify missing versus non-missing data using all other features. If you get significantly better than random results, then your data aren't missing at random.
A: This is not possible, unless you managed to retrieve missing data. You cannot determine from the observed data whether the missing data is missing at random (MAR) or not at random (MNAR). You can only tell whether the data is clearly not missing completely at random (MCAR). Beyond that only appeal to plausibility of MCAR or MAR as opposed to MNAR based on what you know (e.g. reported reasons for why data is missing). Alternatively, you might be able to argue that it does not matter too much, because the proportion of missing data is small and under MNAR very extreme scenarios would have to happen for your results to be overturned (see "tipping point analysis").
A: You want to know whether there is some correlation of a value being missed in feature and the value of any other of the features.
For each of the features, create a new feature indicating whether the value is missing or not (let's call them "is_missing" feature).
Compute your favourite correlation measure (I suggest using here mutual information) of the is_missing features and the rest of the features.
Note the if you don't find any correlation between two features, it is still possible to have a correlation due to group of features (a value is missing as a function of XOR of ten other features).
It you have a large set of features and a large number of values, you will get false correlations due to randomness. 
Other than the regular ways of coping with that (validation set, high enough threshold) You can check if the correlations are symmetric and transitive.
If they are, it is likely that they are true and you should further check them. 
A: I found the information I was talking about in my comment.
From van Buurens book, page 31, he writes
"Several tests have been proposed to test MCAR versus MAR. These tests
are not widely used, and their practical value is unclear. See Enders (2010, pp. 17–21) for an evaluation of two procedures. It is not possible to test MAR versus MNAR since the information that is needed for such a test is missing."
A: There is a useful package called finalfit check here it has a missing_pairs(outcome VAR, explanatory Vars) were you can explore patterns of missingness and decide whether data is MCAR or MAR. It produces pairs plots to show relationships between missing values and observed values in all variables.
