hw6_2.m

%% 2. Randomization test whether significant difference exists between 
%  the density of Oncaea Venusta and Canthocalanus pauper. (Assume all 
%  data are independent and use all 34 stations.) Report your p-value.

% sp1 is Oncaea venusta, sp2 is canthocalanus pauper
sp1=copepod_comp(169,:);
sp2=copepod_comp(15,:);

mu_sp12=mean(sp1-sp2);

% mix two species
sp12=[sp1,sp2]';

% randomization with 5000 times
sp12_perm=[];
for i=1:r
    [a,b]=sort(rand(length(sp12),1));
    resample_obj=sp12(b);
    sp12_perm=[sp12_perm,resample_obj];
end

rsp1_perm=sp12_perm(1:34,:);
rsp2_perm=sp12_perm(35:68,:);

% calculate the permutation means
%m_rsp1=mean(rsp1_perm);
%m_rsp2=mean(rsp2_perm);
%mrsp12_mean=(m_rsp1-m_rsp2);
mrsp12_mean=mean(rsp1_perm-rsp2_perm);
 
% percentile method to calculate p-value
p_pvalue=2*(1-sum(mrsp12_mean<=mu_sp12)/r
 
% t-test
%se_sp12=sumsq(mrsp12_mean-mu_sp12)/(r-1);
%t_sp12=(mean(mrsp12_mean)-mu_sp12)/se_sp12/(r-1);
%sp12_pv=tcdf(t_sp12,r-1)