Answers to Warren's post:

That is, *any* distance involving a point that has a nan is nan. This seems like a reasonable default behavior.

I agree, it is the way that most functions in other lanaguages/packages handle it.

What should nanpdist(y) be?

Based on your code snippet on StackOverflow and your comment in the github issue, my understanding is this: for any pair, you ignore the coordinates where either has a nan (i.e. compute the distance in a lower dimension). In this case, pdist(y) would be

[nan, 4, 6, 8, 2, 4, 6, 2.83, 5.66, 2.83]

(I'm not sure if you would put nan or something else in that first position.)

Or, if we use the scaling of `n/(n - p)` that you suggested in the github issue, where n is the dimension of the observations and p is the number of "missing" coordinates,

[nan, 8, 12, 16, 4, 8, 12, 2.83, 5.66, 2.83]

Is that correct?

That is what I suggest. The appropriate scaling would have to be checked/discussed in detail as it may differ between distance and similarity measures.

What's the use-case for this behavior? How widely used is it?

I work in bioinformatics and my data set consists of thousands of vectors corresponding to different treatment parameters. Each vector consists of basically the changes in expression levels of a number of genes. I am interested in clustering the treatments, i.e., determine which treatments introduce similar gene expression patterns. Not every treatment leads to significant expression changes, of course, which is why there are missing values. So the vectors have roughly 3000 elements and most of them have about 200 missing values. The data are scaled to follow a normal distribution so I could just replace the missing values with the mean and be done with it but I don't think that's the correct approach. I also don't want the current pdist behavior as it would disregard the majority of my otherwise perfectly valid data. As to the popularity of this use-case: Clustering of gene expression data is very wide-spread, however, usually all gene expression data are considered and thus every treatment consists of a completely filled vector. I can't claim that my current use-case is very popular, it's a slightly new approach. If you think that this behavior has no place in scipy, no problem at all. Best, Moritz

Thank you for your response Nathaniel. I was a bit concerned that by going into the application this would turn into a discussion about the method rather than whether this is a desirable concept for scipy. I suppose it's not possible to fully separate the two issues so I will indulge you. On Fri, Aug 22, 2014 at 4:33 PM, Nathaniel Smith <njs@pobox.com> wrote:

Just as a scientific issue this seems very odd to me and not at all what statisticians usually mean by missing data. Surely if you want to determine "which treatments introduce similar gene expression patterns" then two treatments that both produce no effect on the expression of the same gene should be counted as more similar to each other? If you've measured an expression change to be near 0 then that's a known measured value that happens to be near 0 -- not an unknown value that could be arbitrarily large or small and you have no idea which. (Obviously I don't know any of the details about your setting, but in particular I worry that your reasoning sounds similar to common misconceptions about what "significant" actually means. "Not significantly different from zero" might well be "significantly different from 1000".)

Since I didn't want the discussion to be about the method I tried to describe the situation briefly and did not give you the whole story. My apologies. The real situation is the following: The gene expression data are mapped onto pathways using information on links between proteins and coding genes. The pathway definitions come from a multitude of source databases and were collected in a single database (http://consensuspathdb.org/). Only pathways that have five or more available scores are considered (this is somewhat arbitrary, I suppose). Each pathway is then assigned a mean score. Pathways that have too few scores are not considered. You can read up on more specifics in [1]. So I consider those pathways that did not make the cut-off of 5 scores as "missing values". If all the treatments had missing values at the same pathways, I'd be tempted to just throw those out. We are considering treatments from different studies, however, and the studies report gene expression changes for different genes and consequently different pathways end up having no scores. I still want to be able to compare treatments between different studies. One approach could be to rethink the scoring of pathways and introduce an uncertainty that is larger for pathways with missing scores but since I'm sitting at the end of a pipeline that lands the treatments and pathway response scores in my lap, my preferred way of dealing with this is to simply scale up the distance between treatments where one has a pathway score and it's missing for the other. If this seems unreasonable to you, I'm all ears. It does make sense in my mind. Cheers, Moritz [1] http://toxsci.oxfordjournals.org/content/124/2/278.full in particular in the subsection "pathway response analysis"

Answers to Jaime's post: Warren has already pointed this out, but let me insist: what is nanpdist,

or the nan keyword expected to do? Treat pairs of vectors with NaNs as lower dimensional, removing pairs of entries where either is NaN? Do those results make any real sense? Thinking of euclidean distance for points in 3D space, I have trouble thinking of a practical situation where "if any Z coordinate is missing, just give me the distance of the projections onto the XY plane" would be anything but a misleading result. I presume the case is different for all those other distances I have never needed to use, so I am just curious of the use case.

Please see my answer to Warren about the use-case. In three dimensions this would certainly not make sense but my use-case has over three thousand dimensions. What I have in mind is a scaling factor for distance metrics, as suggested before, and an appropriate consideration of dissimilarity of the missing coordinate in similarity measures.

Looking at your linked post, from an implementation point of view, at the low level function that is actually going to do the heavy lifting, it is probable better to, rather than hardcode a check for NaN-ness, take a 'where' kwarg, as numpy ufuncs already do ( http://docs.scipy.org/doc/numpy/reference/ufuncs.html#optional-keyword-argum...), and build the masking array in a higher level wrapper. This would make it easier to eventually make this functionality work with masked arrays or the like.

I'd be perfectly happy to do so. The hard-coded check is inspired by bottleneck which does exactly that for all its nan* functions. But I agree that a mask is preferable.

As a separate but related issue, I have had this PR open for almost a year now, https://github.com/scipy/scipy/pull/3163, and although me saying I want to complete it is getting old, hopefully whatever you have in mind can fit with the general structure of that.

I haven't fully grasped your code in umath_distance.c.src but that's probably a separate discussion. I also couldn't tell if some of that code is automatically generated or all written by hand.

Lastly, whatever you go for, I don't think you should do anything to pdist that you don't also do for cdist and the individual distance functions.

Noted and agreed. Best, Moritz

So there's quite obviously not a lot of interest in this. I will simply

You might be rushing to a conclusion a bit. The mailing list was down for a good part of the week. And in general, you might want to let people a bit more time to respond --- response times vary a lot, for better or worse. Evgeni

write my own little package in that case. I guess after the weekend I'll close the issue on github unless anyone wants to keep it open.

@Jaime: I've read up on ufuncs and they definitely seem like the way to go. Can you say a bit more on how you generated scipy/spatial/src/umath_distance.c.src? I assume it was generated and not written by hand, so did you do that with Cython or something else not included in the pull request?

_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev