% This Matlab script has been developped by Jerome G. Prunier.
% 
% Prunier J.G., Kaufmann B., Fenet S., Picard D., Pompanon F.,
% Joly P., Lena, J.P., 2013. Optimizing the trade-off between spatial and
% genetic sampling efforts in patchy populations: towards a better assessment
% of functional connectivity using an individual-based sampling scheme.
% Molecular Ecology (accepted).
% 
% This script can be used to perform spatial Mantel autocorrelation analyses with simple and partial Mantel
% tests with restricted permutations.
%   pop :   the aggregate to which each individual is assigned
%   nperm :   number of permutations
%   DG :      matrix of pairwise genetic distances (dependent  variable)
%   DE :      matrix of pairwise Euclidean distances (explanatory variable I)
%   Dmax :    maximal spatial extent of spatial autocorrelation analysis
%   lag :     extent of distance classes (the number of distant classes is
%               calculated by Dmax/lag)
%   boot :    1 if a confidence interval of Mantel values is to be
%               estimated by boostrap (0 otherwise)
%   nboot :   number of bootstrap iterations (facultative). If boot=1 but nboot is not
%               provided, nboot=1000 (default value);
%   bootprop :   the proportion of sampled plots to randomly discard at each bootstrap iteration

% The script returns a table rablist providing, for each distance class :
%   - the distance class (at mid-point)
%   - a Mantel correlation value r (based on the Pearson coefficient correlation)
%   - a two-tailed p-value p.
%   - the lower and upper confidence limits bouding the 95% confidence
%       interval about r under the null hypothesis of no spatial structure
%       (lower and upper values of Mantel correlation values after permutations)
% If boot=1', the script also returns a vector bootres providing, for each
% distance class, the 95% confidence interval around r as determined by
% boostrap resampling.
% Finally, the script provides a plot of spatial autocorrelation (you will
% be asked for a scale limit for the y-axis).



function [rablist,bootres]=PlotRestrictedMantelAutocor4ok(DG,DE,pop,nperm,dmax,lag,boot,nboot,bootprop)

pop=[pop(:,1),(1:size(pop,1))'];

[rablist,bootres]=RestrictedMantelAutocor(DG,DE,pop,nperm,dmax,lag,boot,nboot,bootprop);

plot(rablist(:,1),rablist(:,2),'k');
hold on;
errorbar(rablist(:,1),rablist(:,2),abs(rablist(:,2)-bootres(:,1)),abs(bootres(:,2)-rablist(:,2)),'Color','k') ;
plot([0,dmax],[0,0],'Color','k')
scatter(rablist(:,1),rablist(:,2),'MarkerEdgeColor','k','MarkerFaceColor','w');
hold on;
for k=1:size(rablist,1)
    if rablist(k,3)<=0.05
        scatter(rablist(k,1),rablist(k,2),'MarkerEdgeColor','k','MarkerFaceColor','k');
    end
end

plot(rablist(:,1),rablist(:,4),'LineStyle',':','Color','k');
plot(rablist(:,1),rablist(:,5),'LineStyle',':','Color','k');
Yup=input('Borne Y sup ? : ');
%Yup=max(abs(rablist(:,2)))+max(abs(bootres(:,2)));
%Yup=1;
ylim([-Yup Yup]);
xlim([0 dmax]);
hold off



function [rablist,bootres,nperlag]=RestrictedMantelAutocor(DG,DE,pop,nperm,dmax,lag,boot,nboot,bootprop)
npop=size(unique(pop(:,1)),1);
nblag=round(dmax/lag);
classe=NaN(1,nblag+1);
classe(1,1)=0;
rablist=zeros(nblag,5);
rablist(:,2)=1;
for i=1:nblag
    classe(1,1+i)=i*lag ;
end
nperlag=NaN(nblag,1);
bootres=NaN(nblag,2);
for k=1:nblag
    disp(['     Class ',num2str(k),'...']);
    DEini=DE;
    DPini=DEini;
    for i=1:size(DEini,2)
        for j=1:size(DEini,2)
            if i>j
                if  DEini(i,j)<classe(k+1) && DEini(i,j)>=classe(k)
                    DPini(i,j)=0;
                    DPini(j,i)=0;
                else DPini(i,j)=1;
                    DPini(j,i)=1;
                end
            end
        end
    end
    if sum(asvect(DPini)==0)==0 || sum(asvect(DPini)==1)==0
        rablist(k,1)=classe(k)+(lag/2);
        rablist(k,2)=0;
        rablist(k,3)=1;
        rablist(k,4)=0;
        rablist(k,5)=0;
        bootres(k,1:2)=zeros(1,2);
        continue
    else
        nperlag(k,1)=(size(DPini,1)*(size(DPini,1)-1)/2)-sum(asvect(DPini));
        DP=DPini;
        [rab,p,UL]=RestrictedMantel4Autocor(DG,DE,DP,pop,nperm);
        rablist(k,1)=classe(k)+(lag/2);
        rablist(k,2)=rab;
        rablist(k,3)=p;
        rablist(k,4)=UL(1);
        rablist(k,5)=UL(2);
        aff=pop;
        
        if boot==1
            disp('          Bootstraping...');
            bootrablist=NaN(size(rablist,1),nboot);
            for bk=1:nboot
                bootpop=(randperm(npop))';
                bootpop=bootpop(1:round((npop-bootprop*npop)));
                
                for i=1:size(bootpop,1)
                    aff2=aff(aff(:,1)==bootpop(i),2);
                    if i==1
                        aff3=aff2;
                    else aff3=[aff3;aff2];
                    end
                end
                [rab]=correlation(DG(aff3,aff3),DP(aff3,aff3),DE(aff3,aff3));
                %[rab]=RestrictedMantel4Autocor(typa,size(bootpop,1),3,DG(aff3,aff3),DE(aff3,aff3),DP(aff3,aff3));
                bootrablist(k,bk)=rab;
            end
            bootrablistk=bootrablist(k,:);
            bootrablistk(isnan(bootrablistk))=[];
            bootres(k,1:2)=quantile(bootrablistk,[0.025,1-0.025]);
            %bootres(k,1)=rablist(k,2)-std(bootrablistk);
            %bootres(k,2)=rablist(k,2)-std(bootrablistk);
        else bootres(k,1:2)=zeros(1,2);
        end
    end
end



function [rr,p,UL]=RestrictedMantel4Autocor(DG,DE,DP,pop,nperm)
npop=size(unique(pop(:,1)),1);
veca=asvect(DG);
vecb=asvect(DP);
vec=asvect(DE);
%Standardization :

veca(vec==0)=[];
vecb(vec==0)=[];
no=length(veca);
nveca=(veca-mean(veca));
nvecb=(vecb-mean(vecb));
sveca=sqrt(sum(nveca.^2)/(no-1));
svecb=sqrt(sum(nvecb.^2)/(no-1));
veca=nveca/sveca;
vecb=nvecb/svecb;
%Pearson correlation coefficient :
nveca=(veca-mean(veca));
nvecb=(vecb-mean(vecb));
sveca=sqrt(sum(nveca.^2)/(no-1));
svecb=sqrt(sum(nvecb.^2)/(no-1));
rr = (nveca*nvecb')/((no-1)*sveca*svecb);
%Permutations :
r = zeros(nperm, 1);
popuni=unique(pop(:,1));

for qsdf=1:nperm
    pop3=randperm(size(unique(pop(:,1)),1))';
    pop3=popuni(pop3);
    tem1=pop(pop(:,1)==pop3(1,1),:);
    for i=2:npop
        tem2=pop(pop(:,1)==pop3(i,1),:);
        tem1=[tem1;tem2];
    end
    DA=DG(tem1(:,2),tem1(:,2));
    avec=asvect(DA);
    avec(vec==0)=[];
    anvec=(avec-mean(avec));
    asvec=sqrt(sum(anvec.^2)/(no-1));
    avec=anvec/asvec;
    anvec=(avec-mean(avec));
    asvec=sqrt(sum(anvec.^2)/(no-1));
    r(qsdf,1)= (anvec*(nvecb)')/((no-1)*asvec*svecb);
end
if rr>=0
    p = (sum(r>=rr))/(nperm+1);
else p = (sum(r<=rr))/(nperm+1);
end

r2(:,2)=sort(r);
r2(:,1)=1:size(r,1);
UL=[quantile(r2(:,2),0.025),quantile(r2(:,2),1-0.025)];

function vect=asvect(mat)
m=size(mat,1);
n=((m-1)*m)/2;
k=1;
vect=zeros(1,n);
for i=(2:m)
    for j = (1:(i-1))
        vect(k)=mat(i,j);
        k=k+1;
    end
end
