# P-values differ in tests for normality between Mathematica and SAS

I am testing the following data for normality, in Mathematica and SAS:

data = {91.6667, 87.5, 54.1667, 75., 83.3333, 87.5, 83.3333, 62.5, 70.8333,70.8333, 75., 87.5, 66.6667, 87.5, 70.8333, 54.1667, 79.1667,95.8333, 70.8333, 37.5, 75., 79.1667, 75., 66.6667, 87.5, 75., 58.3333, 87.5, 91.6667, 91.6667, 54.1667, 87.5, 75., 79.1667, 95.8333, 66.6667, 91.6667, 62.5, 79.1667, 70.8333, 87.5, 91.6667, 54.1667, 87.5, 58.3333, 66.6667, 75., 83.3333, 70.8333, 58.3333, 75., 62.5, 41.6667, 37.5, 91.6667, 54.1667, 91.6667, 83.3333, 62.5, 87.5, 29.1667, 87.5, 91.6667, 75., 58.3333, 75.};


In Mathematica, my code is:

m = N[Mean[data]];
std = N[StandardDeviation[data]];
mo = Commonest[data][[1]];
Grid[{{Labeled[m, "Mean: ", Left]}, {Labeled[mo, "Mode: ", Left]}, {Labeled[std, "Standard Dev.: ", Left]}}, Alignment -> Left]
d = NormalDistribution[m, std];
H = DistributionFitTest[data, d, "HypothesisTestData"];
H["TestDataTable", All]


The output is:

In SAS, I'm using the Univariate procedure on the same data:

proc univariate data=knowledgescorescsv normal plots mu0=73.9268;
var Post;
run;


The output is:

As you can see, the mean, mode, standard deviation and Anderson-Darlin, Cramer-von Mises and Shapiro-Wilk statistics are the same in Mathematica and SAS. However, their corresponding p-values differ substantially, enough that Mathematica's p-values would lead me to fail to reject the null hypothesis while SAS p-values for these three statistics would lead me to reject the null hypothesis (i.e., that the data is normally distributed).

Why would the p-values for these statistics be so vastly different between Mathematica and SAS?

• …but, where's data? Feb 17 '16 at 2:02
• @J.M. Hi, in Mathematica, I'm importing the data from an Excel file on my Desktop, while in SAS, I'm uploading the data in .csv format. Feb 17 '16 at 2:05
• All well and good, but I was nudging you to post your data here or some other convenient location so potential helpers can run their own tests. Feb 17 '16 at 2:07
• What is a good way to post the data here? Could I attach it? Feb 17 '16 at 2:16
• If data is not too big, you can paste the contents of the *.csv file here; otherwise, you can use Pastebin. Feb 17 '16 at 2:17

## 3 Answers

I should have looked at your code but I didn't, because I saw that the test-statistics for the tests in question were the same.

What I did then is to look up how p-values are computed for Anderson-Darling, which seems to be not as standardized as one might think. I came across this site where a formula is given to compute the p-value from the test-statistic. In your case, this formula was

p[ad_] := Exp[1.2937 - 5.709 ad + 0.0186 ad^2];
p[1.275759]
(* 0.00258163 *)


which supported the result that is given by SAS. At this point, I DID look on your code.

So what you do is to use the test wrong. By giving the (estimated) normal distribution as second parameter, Mathematica does not make the usual single sample normality test but it tests against your given distribution.

This results in the same statistic-value, but obviously the p-values are computed differently. Therefore, the simple solution is

DistributionFitTest[data, Automatic, {"TestDataTable", All}]


• Thank you for looking into this @haliruta. I initially used "Automatic" for the second parameter until I saw the small p-values and attempted using the estimated parameter. The reason I did so was an example in the Mathematica documentation for ref/DistributionFitTest were it does the same thing in order to "analyze whether a dataset is drawn from a normal distribution." The function used is: "H = DistributionFitTest[data, d = NormalDistribution[82, 9], "HypothesisTestData"];". A follow-up is, when would one use the NormalDistribution[] parameter and why is it valid in this example? Feb 17 '16 at 3:53
• I'm not a statistician, but purely from the standpoint of using the p-value as indicator whether your sample was drawn from a known distribution it makes a difference for me if I (a) KNOW that I test for normality or (b) I test against another distribution (which happens to be a normal distribution). In case (b), Mathematica obviously accounts for that by adjusting the p-value accordingly. If you want ta detailed answer, I would give support@wolfram.com a shot. Feb 17 '16 at 4:20
• @Lucas Another, even better suggestion is to ask the question, why p-values might be different for tests against a specific distribution on Cross Validated. Refer in your question to this Q&A here and mention Mathematica, but don't make the question specific to Mathematica so that it isn't off topic there! The reason is, over there the user whuber is mod and very knowledgeable AND he is a high-rep user on this site as well. Therefore, he knows Mathematica very good and probably won't hesitate to answer. Feb 17 '16 at 4:24
• I will attempt posting the question on the Mathematica stackexchange as I could refer directly to the example given in the Mathematica documentation. Otherwise, I would really need to have clarity on how to make the question more generic as you suggest. Feb 17 '16 at 4:55
• @Lucas Just post the question on Cross Validated and comment here giving me the link. I could look over it and improve it if it is really necessary. Another thing: Can you make a distribution test against a specific distribution in SAS as well or is this what SAS is already doing? Otherwise, you might want to see whether there are different p-values as well. Feb 17 '16 at 5:12

Yes, the P-values are being computed differently but there's a good reason for doing so and no bug report to Support needs to be made.

First generate some data:

SeedRandom[12345];
data = RandomVariate[NormalDistribution[], 2500];


Now test to see if the sample is from a normal distribution (and by that I mean any normal distribution: an unknown mean and unknown standard deviation):

h = DistributionFitTest[data, NormalDistribution[μ, σ],
"HypothesisTestData"];
h["TestDataTable", All]


$$\begin{array}{l|ll} \text{} & \text{Statistic} & \text{P-Value} \\ \hline \text{Anderson-Darling} & 1.28718 & 0.00211436 \\ \text{Cramer-von Mises} & 0.184634 & 0.00785495 \\ \text{Pearson }\chi ^2 & 29.6667 & 0.000242055 \\ \text{Shapiro-Wilk} & 0.930615 & 0.0011702 \\ \end{array}$$

The same output is obtained with NormalDistribution[μ,σ] is replaced with Automatic.

But now suppose we wanted to test against a specific normal distribution. In fact suppose we test against a normal distribution with the same mean and standard deviation as the sample mean and sample standard deviation.

h = DistributionFitTest[data,
NormalDistribution[Mean[data], StandardDeviation[data]],
"HypothesisTestData"];
h["TestDataTable", All]


$$\begin{array}{l|ll} \text{} & \text{Statistic} & \text{P-Value} \\ \hline \text{Anderson-Darling} & 1.27575 & 0.239873 \\ \text{Cramer-von Mises} & 0.183053 & 0.302856 \\ \text{Pearson }\chi ^2 & 29.6667 & 0.000971023 \\ \text{Shapiro-Wilk} & 0.930615 & 0.0011702 \\ \end{array}$$

So we see that all of the P-values (except the Shapiro-Wilk) have increased because we've stacked the deck by testing against a normal distribution with the same mean and standard deviation. In other words we're much less likely to reject the null hypothesis.

So the P-values differ but that's because we're testing two different hypotheses. The first is if the sample might have been drawn from an unspecified normal distribution and the other hypothesis deals with a specific normal distribution.

(I've only kept the tests that matched what SAS produced and why the Shapiro-Wilk P-value does not change is left as an exercise for the reader.)

Minor point: DistributionFitTest[myData, Automatic, {"TestDataTable", All}] blows up if myData includes negative numbers, but the fix is trivial:

DistributionFitTest[N[myData], Automatic, {"TestDataTable", All}]

• You should probably report this to Support. Apr 7 '18 at 0:25
• The original command DistributionFitTest[-data, Automatic, {"TestDataTable", All}] using the data above (but made all negative with the minus sign) works fine on Mathematica 11.1, Windows 10. What version of Mathematica are you using? Would you supply a specific dataset where DistributionFitTest blows up?
– JimB
Apr 7 '18 at 5:21
• @JimB I'm using V11.3;myData = {-4, -3, -3, -2, -2, -2, -1, -1, -1, -1, 0, 0, 0, 1, 1, 2}; DistributionFitTest[myData, Automatic, {"TestDataTable", All}] I see now that it's more a matter of MMCA leaving some statistics in expanded form, and N[] just collapses that into a single real. Apr 13 '18 at 15:59
• Yes. When you use integers (or rather rational numbers) without decimal points, some of these statistics give you the "exact" number representation. Also, most if not all of these measures assume a continuous distribution and when there are lots of ties (especially all integer values), they may not be appropriate and therefore misleading.
– JimB
Apr 13 '18 at 16:12