Comparing observed data to simulated data?

Question

I have a set of observed values (shown with black dots in the figure) that I would like to compare to some simulated data (100 simulated datasets shown as box plots with quartiles, extremes (excluding outliers) as whiskers, and outliers as white dots).

The observed values are outside of the 95% confidence intervals of the simulated dataset so it seems pretty obvious that there is a significant difference between the observed and the simulated; but in some cases, the difference might not be so obvious. Which statistical test (preferably in R) can I do on this data to get a P-value for the significant difference between the observed and the simulated sets?

Answer 1

Building on the comment by @whuber, a straightforward test that would deliver reasonable results would be a Chi square test that compares your observed counts in each bin of "number of shared mutations" with the expected counts in that bin. This is particularly the obvious way to do it if you know how the data were simulated and hence you can directly use the expected counts in the bin; but even if you don't it is a reasonable pragmatic approach.

In the current case the expected number with zero shared mutations is zero, which will make a Chi square statistic infinite (because its calculation involves dividing by the expected number in each bin), and give a p-value of zero, which is to be expected. Basically this reflects that it is literally impossible for the observed data to have been generated from a distribution that gives zero probability to zero mutations.

> horses <- data.frame(
+ expected = c(0,270,410,230,80,10),
+ observed = c(800,230,40,10,5,1),
+ mutations =c(0,1,2,3,4,"lots")
+ )
> 
> horses$expected.scaled <- horses$expected * 
+ sum(horses$observed) / sum(horses$expected)
> 
> horses
  expected observed mutations expected.scaled
1        0      800         0            0.00
2      270      230         1          293.22
3      410       40         2          445.26
4      230       10         3          249.78
5       80        5         4           86.88
6       10        1      lots           10.86
> 
> X <- with(horses,
+ sum((observed-expected.scaled)^2/expected.scaled)
+ )
> X
[1] Inf
> 1-pchisq(X,5) 
[1] 0

Answer 2

Consider the frequency distribution from each simulation run (and the observed data) as a 9-dimensional vector with $j$th element being the number of observations with $j$ shared mutations. Now you have an observed vector $z$, and a set of simulated vectors $x_i$, $i=1,\ldots,100$. Your goal is to determine if the vector $z$ could have come from the distribution defined by $x_i$'s.

A possible approach is to look at deviations from the mean of the simulations. You could use Euclidean distance, or Euclidean distance on the square-roots (which has variance-stabilizing properties), or Mahalanobis distance, or anything else you come up with. Now your simulated data will give you a random sample of 100 of such distances $||x_i-\bar{x}||$, which you can compare to the observed distance $||z-\bar{x}||$ to get a p-value.

If you are braver, you can even make multivariate joint normality assumptions, fit a multivariate normal to the simulations (i.e. get its mean and variance matrix) and calculate a p-value based on that. This would be the parametric equivalent of using the Mahalanobis distance.

Answer 3

One way to test and get a p-value by comparing actual data to simulated data (though in a way very different from what you show above) is discussed in:

 Buja, A., Cook, D. Hofmann, H., Lawrence, M. Lee, E.-K., Swayne,
 D.F and Wickham, H. (2009) Statistical Inference for exploratory
 data analysis and model diagnostics Phil. Trans. R. Soc. A 2009
 367, 4361-4383 doi: 10.1098/rsta.2009.0120

And is implemented in the vis.test function in the TeachingDemos package for R.

One approach more along the lines that you are showing would be to compute a p-value for each x-value above as the proportion of simulated points that are as extreme or more extreme than the observed point. Then use meta-analysis techniques to combine the p-values (9 in the case of your plot above), one option is that the negative log of a p-value (under the null) follows a chi-square distribution with 2 degrees of freedom and you can add chi-square statistics (though the degrees of freedom should probably be discounted since I am guessing that the parts above are not independent). Or if you can define "more extreme" in terms of the whole set of 9 points from each simulation then you could just do the simple p-value instead of combining 9.