function [p]=PvalMatchResample(samp1,samp2,N)
if nargin<3 N=20000; end
% Matched Pairs
% samp1,samp2 are vectors containing numerical data to be tested corresponding to matched individuals each given one of two randomly assigned treatments.
% The test statistic here is the mean of a numerical variable and the p-values calculated estimate the probabilities of mean differences greater than or equal to those actually observed assuming the null hypothesis, that the mean difference is 0: i.e. the treatment has no overall effect and the actual values seen in the samples are due to chance variation alone.
% List your two samples so that mean(samp1) exceeds mean(samp2) and, if you are running a 1-sided test, your alternative hypothesis must agree that this SHOULD be true. There is no need to run a 1-sided test at all if this cannot be done: the data you collected will not support your alternative hypothesis and the null hypothesis will never be rejected.
% If the probability is very small it is unlikely that the null hypothesis that generated it (namely that the mean difference of the numerical variable is 0) is true, and this counts as evidence in favor of the alternative hypothesis.
% p(1) corresponds to the ``two-sided'' alternative hypothesis, which is the hypothesis that the mean difference is not 0. It is the probability that the magnitude of the mean difference equals or exceeds that which was observed assuming the null hypothesis. p(2) corresponds to the ``one-sided'' alternative hypothesis: that the mean of the first treatment variable exceeds the mean of the second. p(2) estimates the probability that the mean (positive) difference equals or exceeds that which was observed assuming the null hypothesis.
% samp1=[72,65,81,73,53,55,78,43];
% samp2=[45,12,33,24,49,87,10,69];
mean1=mean(samp1)
mean2=mean(samp2)
if mean(samp2)>mean(samp1) error('Examine your data. The mean of the first sample should exceed the mean of the second to run this test. If you are running a 2-sided test the order is irrelevant, so switch the order. But if you are running a 1-sided test we must have the first mean exceeding the second and your alternative hypothesis must agree that this SHOULD be true. There is no need to run a 1-sided test at all if this cannot be done: the data you collected will not support your alternative hypothesis and the null hypothesis will never be rejected.'), end
rng('shuffle') %assign time-based seed to random number generator
n=length(samp1);
diff12=mean(samp1-samp2);
counttwosided=0;
countgreater=0;
for i=1:N
choose=rand(1,n);
for j=1:n
if choose(j)<0.5
randsamp1(j)=samp1(j);
randsamp2(j)=samp2(j);
else
randsamp1(j)=samp2(j);
randsamp2(j)=samp1(j);
end
end
randiff(i)=mean(randsamp1-randsamp2);
if randiff(i)>=diff12
countgreater=countgreater+1;
end
if abs(randiff(i))>=abs(diff12)
counttwosided=counttwosided+1;
end
end
p=[(counttwosided+1)/(N+1),(countgreater+1)/(N+1)];
% (1 is added because the original samples automatically satisfy both conditions.)
figure(77)
histogram(randiff,100);
hold on
plot([-diff12 -diff12 diff12 diff12],[0 N/40 N/40 0])
hold off
end