# Comparing observed data to simulated data?

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?

-
It is unusual for boxplots to show confidence intervals. Don't these show quartiles, fences, and any outliers? It also seems worthwhile asking why you would want a p-value when the difference is so blatantly obvious: clearly the simulation does not reproduce the observed values. –  whuber Jul 20 '12 at 15:34
Hi whuber, Sorry I am mistaken, the whiskers represent the maximum and minimum (excluding outliers). The reason I want to do this, is that in some of my datasets it will be a much closer fit, so then a p-value become very useful. –  blJOg Jul 20 '12 at 16:09
Could you comment on the nature of the simulation? In principle, it is merely drawing data from an underlying distribution. We could probably bypass the simulation altogether by knowing the distribution. This leads to very simple tests, such as the Chi-squared--and obviates the need to run the simulation, too! –  whuber Jul 20 '12 at 18:13
So I have observed a number of mutations in 48 horses and determined how many mutations are shared between pairs of horses. 800 pairs of horses have no mutations in common, approx 200 pairs of horses share 1 mutation etc... For each observed mutation, the mutation is randomly assigned to the same number of horses it was found in. I then determine based on this random assignment, how many pairs of horses share 0...8 mutations if mutations were spread randomly in the horse population. –  blJOg Jul 22 '12 at 9:51

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

-

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.

-

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.

-