function [cv, pval] = pairpermute(y, cov, geno, chrid, nperm)
%permutation test for joint effects of two markers
%this program computes the permutation test described in
%	Sugiyama et al, 1999
%The test compares unconstrained regression model with full interaction 
%	term to a null model
%
%Input Variable:
%     y: nx1 vector of trait values
%   cov: nx1 vector of covariates,  use [] if none
%  geno: nxm genotyping array, encoded as 0-1 or as 0-1-2
% chrid: mx1 vector = chromosome number associated with markers
% nperm: number of permutations
%OUTPUT VARIABLES
%  F1:	the overall F-test mxm but only the lower triangle (i<j) is used
%  cv: critical values (.8, .9, .95 and .99) from permutation analysis
%  pval: permutation p values for observed maximum F-statistic
%
%Copyright notice
%this code was written by Gary Churchill who retains all rights.
%If this code to generate results in publications please cite
%Sugiyama F, Churchill GA, Higgins DC, Johns C, Makaritsis KP, Gavras H, Paigen B (1999)
%Quantitative trait loci analysis of blood pressure in mice on a high salt regimen:
%Concordance with rat and human loci. Submitted to Nature Genetics.
%
%missing data check
if any(any(isnan(geno)))| any(isnan(cov)) | any(isnan(y))
   exit('missing data in geno, yy or cov')
end

%genome scan

%permutation test
Fmax = zeros(nperm,1);
for j = 1:nperm
   %obtain permutation indices
   ord = randperm(length(y));
   
   %permute the trait values
   ptrait = y(ord);
   
   %permute the covariates
   if ~(isempty(cov))
      pcov = cov(ord,:);
   else
      pcov = [];
   end
   
   %compute overall F and record the maximum value
   F1 = pairwise(ptrait, pcov, geno, chrid);
   Fmax(j) = max(F1(:));
end

%sort maximal F's
Fmax = sort(Fmax);

%compute critical values 
cv=Fmax([ceil(.8*nperm),ceil(.9*nperm),ceil(.95*nperm),ceil(.99*nperm)]);

%compute the p-value
if nargout>1
	F1 = pairwise(y, cov, geno, chrid);
	pval = sum(Fmax(:) > max(F1(:)))/nperm;
end