Selecting a random subset to match PDF of another or a given distribution

Problem definition: I have a list of patient records from a hospital, who are of a certain distribution of ages. I want to randomly sample patients from this list in such a way that their age distribution resembles that of a known cohort of patients. E.g. patients with heart failure have an age distribution starting at about 40y.o. and going all the way to the max ages in the data set.

Example:

RandomVariate[] does this from a specified numerical distribution, but not from real data. RandomSample[] does it from real data, but only according to the distribution of the data itself. I'm hoping there's some way to use the mixture of the two (perhaps by generating a list of record IDs, sorted/binned by age, and then extracting those parts of the list?) which will solve my problem.

Limitations on functional solutions:

1. No replacement (or available option for no replacement). Each patient's record should be samples only once.
2. One should not be able to sample more data, than available. Which means that if any part of the sample distribution hits the maximum available data for that age/bin then, all other bins should be reduced to never return more data than that (within an allowed margin of error, since this is random sampling and also very small bins, such as "90-year-olds", may not have sufficient data sizes, but inaccuracies in those bins won't significantly impact overall analysis accuracy).
• What's the rationale for the limitations? – mef Apr 29 '18 at 22:36
• Since we're random sampling from a pool of unique individuals (patients), replacement makes no sense. We wouldn't want to sample the same person twice. As for part two, I guess the logic there is that the resulting distribution will not match the desired one of we start getting too many people from plays there they exist and non at all from plays where they don't. This is normal to some extent as the sampling is random. But one generally wants the shape of the distribution to remain similar to the set we're trying to match. – Gregory Klopper Apr 29 '18 at 22:41
• Target distribution isn't always "Normal". I'm trying to achieve similar modality, mean/median/variability/IQR. A small deviation in the resulting PDF is acceptable, but not too much. – Gregory Klopper Apr 29 '18 at 22:42
• My top candidate right now is a bin-wise random sampling with variable width bins as per target values density. – Gregory Klopper Apr 29 '18 at 22:43
• Your limitation "Each patient's record should be samples only once." is not clear to me. If you sample sufficiently long and make this restriction, you will end up with taking all persons into account (provided that the number of patients is finite and countable). That is just trivial. What's the point? Or do you have infinite number of patients? – yarchik May 9 '18 at 11:41

Hopeffully this will work for you:

First I will generate some data that looks like yours

Dist =  SplicedDistribution[{3/5, 2/5}, {0,50, \[Infinity]}, {UniformDistribution[{0,50}],NormalDistribution[60, 15]}];

data1 = Sort[RandomVariate[Dist, 50000]];
data2 = Sort[RandomVariate[WeibullDistribution[7, 75], 50000]];


The Histogram looks like

 Histogram[{data1, data2}, {1}]


Since age is only integers we will group our data and calculate the PMF. The PMF is important because it will give us the weights that we will use in our random sampling. Then we generate a binned list of data1, replace it by ones and multiply by the PMF to obtain the weights.

cases2 = BinCounts[data2, {0, 120, 1}];
pmf2 = cases2/Total[cases2];
list = BinLists[data1, {0, 120, 1}];
func[x_] := 1;
ones = Map[func, list, {2}];
weights = Flatten[ones*pmf2]


Now you sample data1 with the corresponding weights:

sample = RandomSample [weigths -> data1, 3000];
Histogram [sample, {1}]


Its not perfect, but still similar to our cohort distribution.

• Fairly straightforward . – Anton Antonov May 5 '18 at 7:36
• Yes, VERY straightforward! I didn't realize that RandomSample can accept weights. This certainly makes things much simpler. I've the same question as for the other answer (both are really interesting useful answers!). How would I know how much to sample before I run out of reasonable samples and the distribution shape would change too drastically? Ideally, I'd like to get as many samples as available without drastically altering the shape of PMF. – Gregory Klopper May 8 '18 at 22:40
• How would I generalize this for a 2 or more variables (e.g. age, gender, and BMI)? Subselection (select by gender, then by age, then by BMI)? or is there a way to do this simultaneously with a multidimensional weight? – Gregory Klopper May 8 '18 at 22:42
• Run this after your code. This is what I meant above. How do we know when to stop before the sampled distribution is nowhere near the target PMF? s2 = RandomSample[weights -> data1, 2000]; s7 = RandomSample[weights -> data1, 7000]; s15 = RandomSample[weights -> data1, 15000]; s25 = RandomSample[weights -> data1, 25000]; Round[Quantile[#, {.05, .1, .25, .5, .75, .9, .95}], .1] & /@ {s2, s7, s15, s25} // TableForm Histogram[{s2, s7, s15, s25, data1, data2}, {1}, "Probability", ChartLegends -> {"2k", "7k", "15k", "25k", "source", "target"}] – Gregory Klopper May 8 '18 at 22:59
• Just noticed that when I ask for too much data to sample, considering that 0 weights eliminate some data, WM will happily return "RandomSample cannot generate a sample of length 25000, which is greater than the length of the sample set". That's awesome. But it will keep sampling up until that point from ALL available data, returning pretty much a mirror image of the source distribution for very large N. If I knew that my target distribution is always Weibull or any monomodal distribution, I could use PMF value at its median/mode to limit sample count. But what about multimodal target dist? – Gregory Klopper May 8 '18 at 23:59

This problem is known and can even be more generalized than you have described it. Especially, when researching fields like cardiovascular diseases or aortic valve problems, it is often the case that you have a large healthy control group and a small group of patients that presented with conditions like bicuspid aortic valves.

To get any meaningful results from your analysis, you need to have two groups that show similarities in essential properties. Age is one of those properties but by far not the only one. In heart diseases, you additionally want to ensure that your patient and control group show similar distributions regarding e.g., smoking, diabetes, or hyperlipidemia. Each of those properties, so-called confounders, will have a particular distribution in your patient group and you need a control group that matches as close as possible.

A very simple, but widely used approach is Propensity Score Matching. The algorithm works as follows:

1. Your data is in the form {age, 1} for each patient and {age, 0} for each in the control group.
2. You make a logistic regression that gives you for each entry the pscore.
3. You use this pscore to select for each patient the most matching subject from the control group

This approach can lead to duplicates if your control group is not large enough, but there some ways to handle this. First, you can exclude the pair when it leads to duplication. Second, you can loosen your restriction and use the control patient with the second-best pscore.

Let me give a simple example for your specific question.

cohort = {#, 1} & /@ RandomVariate[GammaDistribution[10, 4.], 200];
nonCohort = {#, 0} & /@ RandomReal[{0, 100}, 300];


The two groups have a vastly different distribution

Histogram[{cohort[[All, 1]], nonCohort[[All, 1]]}, {3}]


Now we make a LogitModelFit for all entries combined

ℳ = LogitModelFit[Join[cohort, nonCohort], age, age]


M[age] gives now the propensity-score for each age and we use Nearest to select the best nonCohort entry for each cohort entry:

nf = Nearest[
Association[#[[1]] -> ℳ[#[[1]]] & /@ nonCohort]];
fixedNonCohor = Flatten[nf /@ (ℳ /@ cohort[[All, 1]])];
SmoothHistogram[{cohort[[All, 1]], fixedNonCohor}, 2,
PlotLegends -> {"Original Cohort", "Fixed Non-cohort"}]
`

I understand that this answer is not exactly what you had in mind, but propensity score matching is widely used in this kind of research and I hope it still helps you.

Would there be limitations to selecting with Nearest if there's insufficient data to create the right sample? E.g. in your example if there were no patients over 50 in non-cohort.

Still, there would be a subject in the non-cohort that is the nearest match. However, it is evident that the histograms of age would not match so perfectly after the p-score matching. Additionally, when you have several confounders, the likelihood of finding a perfect match between two subjects from the cohort and non-cohort decreases.

Let me repeat the important points: When you design a research study with experiments on tissue or patients, then you will make sure upfront that you have patients and controls that match in all important properties. What these properties are, depends on your research subject but Age is one property that is used in most cases. So the distribution for age should match. If you are doing research in the field of retinal diseases, another important property would be diabetes as this impacts the retina.

However, in such a design, you are free to define upfront how many patients you need and of what age they need to be. So you select your study subjects to ensure that your statistical analysis can even show the effect you are investigating and that possible effects are not because you compared apples with oranges.

But such an approach is not always possible. If you are dealing with rare conditions, you might have to use what you've got. This seems to be the situation you are in. So you have a small group of subjects that have a disease or condition and you usually have a large group that doesn't have this condition.

P-score matching is now one of the easiest methods to design your study afterward by selecting a sub-group of controls (and probably patients). After the matching, the groups will look like you designed them and make a statistical analysis possible.

Coming back to your question: when you don't have controls (or non-cohort) with an age over 50 and you chose the closed p-score you find anyway, it will make your distributions different and probably introduce effects that are not coming from your medication or experiment, but simply because the groups differ.

The cleanest thing would probably to make a 1:1 matching without repetition. This means, for each patient, you select a control if the p-score matches to a high degree. If this does not exist, you don't use this patient. If several patients match with the same control, you either chose the second best control if the p-scores are close enough, or you throw the patient out as well.

There is no common practise as far as I can tell and different publications used different approaches. You can look for instance over this abstract to see how it is applied in sciences.

(1) seems that your algorithm selects items "with replacement" as it can select the same non-cohort item twice for every cohort item;

Yes, this small example does this. In the research project I was working on, however, I ensured a strict 1:1 match without repetitions. You can easily implement it by throwing out a control you already used. Additionally, don't use the nearest without comparing the difference of the p-score as I explained above.

(2) how would I know how much I can sample before similarities run out? Where's the limitation on the size of the sample?

It depends on the size, the distribution of your non-cohort and of course how good your cohort lies inside the distribution of your non-cohort. If your non-cohort doesn't have patients above 50, while your cohort only contains ages 60-70, this approach will not work at all. The underlying assumption is to have a large non-cohort compared to the cohort group.

Finally, I want to stress the fact that p-score matching is a simple, but not the best way to do it. It seems to me that it is so widely used because of its simplicity and because medical researchers usually don't have a statistician at hand. It's a bit like the infectious widely misused significance level. I'm sure our statistical expert JimB could elaborate on better alternatives.

• VERY interesting approach! Interesting use of logistic regression. As you mentioned, this approach would allow to even match distributions on multiple parameters. Would there be limitations to selecting with Nearest if there's insufficient data to create the right sample? E.g. in your example if there were no patients over 50 in non-cohort. – Gregory Klopper May 8 '18 at 22:09
• another 2 things I'm not sure about: (1) seems that your algorithm selects items "with replacement" as it can select the same non-cohort item twice for every cohort item; (2) how would I know how much I can sample before similarities run out? Where's the limitation on the size of the sample? – Gregory Klopper May 8 '18 at 22:35
• @GregoryKlopper I added some additional information in my answer. – halirutan May 9 '18 at 3:39
• I greatly appreciate your thorough answer! It's very insightful. I'm sure many will find it useful. However, I chose to go with the other answer. Your method is very expensive computationally & suffers from the "replacement" issue. Resulting sample won't conform to target dist once you de-dup items. Reading up on PSM, I see that it's very typical to have just 1 source item match a slew of target items and many source items can match 1 target item. PSM is a GREAT tool for another purpose, but it does not accomplish specifically what I want - reshaping source to match target PDF. – Gregory Klopper May 11 '18 at 20:29
• Those, curious about PSM, might find this interesting: youtu.be/rBv39pK1iEs – Gregory Klopper May 14 '18 at 4:46