A Monte Carlo Approach to the Permutation Test

Enis SINIKSARAN
Assistant Professor of Statistics
Istanbul University,Faculty of Economics,Department of Econometrics
Beyazit-Istanbul,Turkey
saran@istanbul.edu.tr

Created : July, 2001

The permutation test is mainly used to test whether two samples come from the same distribution. Consider the following data from the well known book by Efron [1]. 16 mice were randomly assigned to a treatment group or a control group. The data give their survival times, in days, following a test surgery.  Our aim is to test whether the treatment group prolong the survival.

           ({{Treatment, 94, 197, 16, 38, 99, 141, 23}, {Control, 52, 104, 146, 10, 51, 30, 40, 27, 46}})
        
       A standard approach to the test is to compute a simple t statistic. Assuming both samples come from a normal distribution with identical variances the statistic is 
             
             (Subscript[t, 0])=\(TraditionalForm\`\(\(x\&-\) - \(y\&-\)\)\/\(S \@\(\((1/m)\) + \((1/k)\)\)\)\)
             
where  \(TraditionalForm\`\(x\&-\)\) and \(TraditionalForm\`\(y\&-\)\) are the sample means of the groups, m and k are the sample sizes,  and S is

           S=\(TraditionalForm\`\@\(\(∑\+\(i = 1\)\%m\((x\_i - \(x\&-\))\)\^2 + ∑\+\(j = 1\)\%k\((y\_i - \(y\&-\))\)\^2\)\/\(m + k - 2\)\)\)
           
           In standard approach we know that the p-value of the one sided test  can be computed by the following  formula
           
           p-value = 1 - Φ( (Subscript[t, 0]))
           
where Φ is the cumulative distribution function of the Students' t distribution with  (m + k - 2) degrees of freedom.  For the mouse data, the difference in the means of the samples is 30.63 and p-value is

           p-value = 1 - Φ((30.63/(54.21 (Sqrt[(1/7) + (1/9)])))) = .141
           
which is a not an evidence against the null hypothesis. So we fail to reject  the hypothesis that the treatment group prolong the survival.
           When we do not rely on assuming that the data come from the normal distribution we can appeal to the permutation test. The basic idea behind it is very  simple: If the null hypothesis is true, or equivalently both samples are drawn from the same distribution  any permutation of the elements in sample X and sample Y is equally probable. So we can combine all the m+k elements from both groups together, then draw a sample of size m without replacement to represent the first group, the remaining k elements represents the second group. We compute the difference in the means of both groups for the all permutations. Thus, we will have a large number of differences which constitutes a distribution.  If the original difference (30.63 for the mouse data) is larger than (1-α).100% of the distribution of the differences, then  we can reject the null hypothesis at a α.100% level of significance. Of course, we do not have to restrict ourselves to the one sided test. When  we  apply two-sided test we will observe whether the original difference falls outside the middle (1-α).100% of the distribution of the differences. If so, reject the hypothesis at α.100% level of significance.
           When the sample sizes get larger, listing all permutations will need a very large number of computations. An alternative is to use a Monte Carlo approach. Provided that we choose a large enough number of  replications, we can generate an accurate distribution of the all possible permutations. 
           The following program  is written to perform Monte Carlo replications and to calculate some summary statistics and confidence intervals dealing with the permutation test.  The user function MCPermutOne[data1,data2,n,alpha] gives the distribution of n permutations for the two samples data1 and data2.  Alpha represents the level of significance. Summary statistics in the table are as follows:
MeanDiff: The difference in the means of the original samples data1 and data2
TestStat: The value of the test statistic in the standard approach calculated by the statistics package of Mathematica  
"Hypothesis Testing"        
OneSided P-value: One-sided p-value calculated by the statistics package of Mathematica  
"Hypothesis Testing"
OneSidedMonteCarlo P-value: One sided Monte Carlo p-value .   

            In addition to the summary statistics and the result of the test, MCPermutOne[data1,data2,n,alpha] gives the graphics of  the distribution of all differences. The solid line represents the mean difference of the original samples.
If we count the points above this line and divide n ,we can obtain the one sided Monte Carlo p-value. So we can make the permutation test visually.

One-Sided  Permutation Test

Needs["Statistics`DescriptiveStatistics`"]
Needs["DiscreteMath`Combinatorica`"]
Needs["Statistics`HypothesisTests`"]

\(MCPermutOne[list1_List, list2_List, n_Integer?Positive, 
        alpha_] := 
Module[{jointdata, m, k, montecarlos, 
          montecarloscomp, means, 
          report, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ meandifference, alldifferences, total, pvalue, meandiff, pval, q, 
          teststat}, 

jointdata = Join[list1, list2]; 

        m = Length[list1]; 
k = Length[list2]; 


        montecarlos = Table[RandomKSubset[jointdata, k], {i, 1, n}]; 

        montecarloscomp = 
          Table[Complement[jointdata, montecarlos[\([i]\)]], {i, 1, 
              n}]; 


        alldifferences = 
          Table[Mean[montecarlos[\([i]\)]] - 
              Mean[montecarloscomp[\([i]\)]], {i, 1, n}]; 

        meandifference = 
          N[Max[Mean[list1], 
                Mean[list2]] - \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Min[Mean[list1], Mean[list2]]]; 

        ListPlot[
          alldifferences, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> All, 
          PlotStyle -> {PointSize[
                1/75]}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {"n", \*""\!\(\(x\&-\)\_1\)-\!\(\(x\&-\)\_2\)""}, 
          ImageSize -> {500, 
              300}\ , \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotLabel -> "All differences below and above the original difference", 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Epilog -> {
			\ \ \ Thickness[ .008], 
			\ \ \ Line[{{0, 
                    meandifference}, {1.1*n, 
                    meandifference}}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[
                meandifference, {1.12*n, meandifference}, {\(-1\), 
                  0}]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ }]; 
Do[
          zzz[i] = If[alldifferences[\([i]\)] > meandifference, 0, 1], {i, 1, 
            n}]; 
\ \ \ \ \ total = 
          Sum[zzz[i], {i, 1, n}]; 
\ \ \ \ \ pvalue = 
          N[1 - total/n]; 

report = 
          MeanDifferenceTest[list1, list2, 0, KnownVariance -> None, 
            FullReport -> True, 
            EqualVariances -> 
              True]; 
pval = {OneSidedPValue /. 
              report}; 
teststat = \(Take[
              Flatten[\(\(report[\([1]\)]\)[\([2]\)]\)[\([1]\)]], 
              2]\)[\([2]\)]; 

q = 
          N[Quantile[alldifferences, 
              1 - alpha]]; 


Print[
          TableForm[{{"MeanDiff", 
                meandifference}, 	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
					{"TestStat", 
                teststat}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {"OneSided P-value", 
                pval}, 	
			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {"OneSidedMonteCarlo P-value", 
                pvalue}			\ \ \ 
}]]; 
If[meandifference > q, 
          Print["Reject the null hypothesis at significance level" -> alpha\ ], 
          Print["Fail to reject the null hypothesis at significance level" -> alpha]
];
]\)

Let us apply the program for the mouse data for 100 replications and 0.05 level of significance.

data1 = {94, 197, 16, 38, 99, 141, 23};
data2 = {52, 104, 146, 10, 51, 30, 40, 27, 46};

MCPermutOne[data1, data2, 100,  .05]

1000 replications

MCPermutOne[data1, data2, 1000,  .05]

Two-Sided  Permutation Test

User function MCPermutTwo[data1,data2,n,alpha] gives the similar results with ones of MCPermutOne[]. To calculate the p-value visually, this time, we count the points above the original difference and below the -original difference, and then divide by n. Additionally, it gives another graphics which shows all confidence intervals calculated for each Monte Carlo sample. In this graphics, confidence intervals are represented by a point. The bigger point represents the confidence interval calculated for the original samples. So when we see the most of the points with the bigger point together in the same quadrant, we can fail to reject the null hypothesis.

Needs["Statistics`DescriptiveStatistics`"]
Needs["DiscreteMath`Combinatorica`"]
Needs["Statistics`HypothesisTests`"]
Needs["Statistics`ConfidenceIntervals`"]

MCPermutTwo[list1_List, list2_List, n_Integer?Positive, 
        alpha_]:=
Module[{jointdata, m, k, montecarlos, montecarloscomp, means, 
            report, meandifference, alldifferences, total, pvalue, meandiff, 
            pval, qlower, qupper, teststat},
jointdata = Join[list1, list2];
m = Length[list1];
k = Length[list2];

montecarlos = 
              Table[RandomKSubset[jointdata, k], {i, 1, n}];
montecarloscomp = 
              Table[Complement[jointdata, montecarlos[[i]]], {i, 1, n}];

alldifferences = 
              N[Table[Mean[montecarlos[[i]]] - 
                    Mean[montecarloscomp[[i]]], {i, 1, n}]];

meandifference = 
              N[Max[Mean[list1], 
                    Mean[list2]] -    
                               Min[Mean[list1], Mean[list2]]];

\(ListPlot[
              alldifferences, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> All, 
              PlotStyle -> {PointSize[
                    1/55]}, 
	\ \ \ \ \ \ \ \ \ \ AxesLabel -> {"n", \*""\!\(\(x\&-\)\_1\)-\!\(\(x\&-\)\_2\)""}, \ 
              ImageSize -> {500, 300}\ , 
              PlotLabel -> "All differences below and above the original difference", 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Epilog -> {
			\ \ \ Thickness[ .015], 
			\ \ \ Line[{{0, 
                        meandifference}, {1.1*n, 
                        meandifference}}], 
	\ \ \ \ \ \ \ \ \ \ Line[{{0, \(-meandifference\)}, {1.1*
                          n, \(-meandifference\)}}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[
                    meandifference, {1.12*n, meandifference}, {\(-1\), 
                      0}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[\(-meandifference\), {1.12*n, \(-meandifference\)}, {\(-1\), 
                      0}]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ }]\);

    Do[zzz[i]=If[alldifferences[[i]] < -meandifference ||       
                   alldifferences[[i]] > 
                          meandifference ,0,1],{i, 1, n}];
total = Sum[zzz[i], {i, 1, n}];
     pvalue = N[1 - total/n];

report = 
              MeanDifferenceTest[list1, list2, 0, KnownVariance -> None, 
                FullReport -> True, TwoSided -> True];
pval = {TwoSidedPValue /. report};
teststat = Take[
                  Flatten[report[[1]][[2]][[1]]], 
                  2][[2]];

qlower = 
              N[Quantile[alldifferences,  alpha/2]];
     qupper = N[Quantile[alldifferences,  1 - alpha/2]];

\(Print[
              TableForm[{{"MeanDiff", 
                    meandifference}, 	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
					{"TestStat", 
                    teststat}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {"TwoSided P-value", 
                    
                    pval}	, 	
			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {"TwoSidedMonteCarlo P-value", 
                    pvalue}	
				\ \ \ \ \ \ \ 
}]]\);

If[meandifference < qlower||meandifference > qupper, 
   Print["Reject the null hypothesis at significance \
level" -> alpha ], 
   Print["Fail to reject the null hypothesis at \
significance level" -> alpha]
 ];

allintervals = 
              Table[MeanDifferenceCI[montecarlos[[i]], 
                  montecarloscomp[[i]], ConfidenceLevel ->  .95, 
                  EqualVariances -> True, TwoSided -> True], {i, 1, 
                  n}];

CI = 
              MeanDifferenceCI[list1, list2, 
                ConfidenceLevel -> 1 - alpha, 
                EqualVariances -> True, TwoSided -> True];

\(ListPlot[allintervals, 
              PlotStyle -> {PointSize[1/75]}, 
              Epilog -> {PointSize[1/45], 
                  Point[{CI[\([1]\)], 
                      CI[\([2]\)]}], 
Line[{{CI[\([1]\)], 
                        CI[\([2]\)]}, {1.1*CI[\([1]\)],  .8*
                          CI[\([2]\)]}}], 
Text[{CI[\([1]\)], 
                      CI[\([2]\)]}, {1.1*CI[\([1]\)],  .8*
                        CI[\([2]\)]}, {\(-1\), 0}]
}, 
              PlotRange -> All, AxesOrigin -> {0, 0}, \ 
              ImageSize -> {500, 300}\ , 
              PlotLabel -> "All Monte Carlo confidence intervals with the original confidence interval"
]\);
]

Let us apply the function for the mouse data for 50 replications and 0.05 level of significance.

MCPermutTwo[data1, data2, 
      50,  .05]

1000 replications

MCPermutTwo[data1, data2, 1000,  .05]

References

[1] B. Efron, " An Introduction to the Bootstrap", Chapman & Hall ,1998, page 11.