I have a set of ten samples, comprised of five pairs of twins, and I have calculated the pairwise distances between all ten samples. The question I would like to answer with this data is:

Are twins more similar to each other than to a random sample

or, more formally

Are the distances between paired samples significantly lower than the distance between unpaired samples

On the face of it, it seems like an easy comparison: all paired distances go in one bin, all unpaired go in the other, and compare. However, I recognize that this breaks all kinds of rules regarding independence. A quick simulation in R suggests that including all pairwise distances significantly under-estimates the significance of the p-values for the comparison. Unfortunately, I have no idea how to solve this problem.

I can think of a number of options, such as: only including one non-twin comparison for each twin pair, averaging all non-twin distances for each sample, or comparing each point to a centroid, but none of these seem entirely reasonable. Is there a standard way for doing such a comparison? Or are pairwise distances simply not the right way to answer our hypothesis?

EDIT: The overall study design is a bit complicated, which is why I didn't post it here originally. However, if working with distances is not possible, then it will be necessary to consider exactly what data I'm working with. I'm including a quick description of the study design below:

We're working on something very similar to a RNA-seq study, where RNA is collected from each of our samples, amplified, and then sequenced (note that only a small portion of the total RNA pool is sequenced). After sequencing, you have a random number of molecules from each sample which we bin according to which gene the RNA molecule was produced from. The output is a sample-by-gene count matrix, which we want to use as an estimate for the relative likelihood that each gene will be transcribed into RNA. Generally, these counts are modeled by the combination of two distributions: a Gamma distribution for the biological noise, and a Poisson distribution for the number of molecules actually sequenced. The resulting distribution is usually modeled as a negative binomial. Whether these assumptions hold for our own data is debatable, but for now I'm going to assume they do.

Making the above assumptions, there are some well-established methods for group comparisons of individual genes. However, quantifying the overall similarity of samples is apparently not common. The major problem with any standard measures of similarity is that variance increases in samples with fewer total molecules, and the samples with the highest total number of molecules are always most similar. So currently, to estimate the differences between samples, we are subsampling down to the smallest total read count and calculating distances on the resulting count matrix. The subsampling is then repeated N times, and an average distance matrix is calculated.

There's probably better ways to quantify the overall similarity, but they need to take these distributions into account.


1 Answer 1


I would use a simple multilevel model. You have an intercept plus a random intercept for each twin pair. The model will estimate the variance associated with twin pairs (that is between pairs), and the residual variance (which is the variance within pairs). You can get confidence intervals for the variances, if you like.

  • $\begingroup$ This may be a workable solution, but I'm struggling with exactly how to implement it. The distances I calculate are actually the average of the distances from multiple repetitions (which is why I was looking for a solution using just the distances). I can imagine creating the multilevel model and averaging the intercepts from each run, but I don't know if there's a way to get significance from this average. Additionally, I'm unclear on how to use the intercepts to show the magnitude of the difference between twin-pair distances and non-twin distances. $\endgroup$
    – cbolen
    Nov 12, 2014 at 1:52
  • $\begingroup$ You don't need to average. Just use (all of the) original scores. Then fit a multilevel model. The variance of the intercepts is between pairs, the error variance is within pairs. Can you say more about your study, your data & your goals? $\endgroup$ Nov 12, 2014 at 3:07
  • $\begingroup$ I don't think using all the data is the correct solution. I'll try and give a bit more information about the setup we're working with, but it's a bit too long to explain in a comment down here. I'll add a description of the experimental design to my original post, above. $\endgroup$
    – cbolen
    Nov 12, 2014 at 19:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.