function [F1, F2, F3, F4] = pairwise(y, cov, geno, chrid)
%genome scan for joint effects of two markers
%this program computes the four test statistics described in
%Sugiyama et al, 1999
%test 1 compares unconstrained regression model with full interaction 
%	term to a null model
%test 2 is the interaction test
%tests 3 and 4 are for coat-tail effects
%
%INPUT VARIABLES:
%     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
%
%OUTPUT VARIABLES
%All F-matrices are mxm but only the lower triangle (i<j) is used
%to store the F-statistic
%  F1:	the overall F-test
%	F2:	the interaction F-test
%	F3:	the Q1 only test
%  F4:   the Q2 only test
%
%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.
%
%this code cannot handle missing data - use impute.m to fill out
%	a genotype array if needed
%missing data check
if any(any(isnan(geno)))| any(isnan(cov)) | any(isnan(y))
   exit('missing data in geno, yy or cov')
end

%constants and storage
[n,m] = size(geno);
F1 = zeros(m,m);
if nargout > 3
   F2 = zeros(m,m);
   F3 = zeros(m,m);
   F4 = zeros(m,m);
end
if isempty(cov)
   c = 0;
else
   c = size(cov,2);
end

if any(geno(:)>1) %set up intercross indicator
   intercross = 1;
else
   intercross = 0;
end

%null regression
X0 = [ones(n,1),cov];
[rss0, df0] = RSS(y, X0);

%model regressions
if intercross
	for i = 1:m-1
   	for j = i+1:m
         if chrid(i) ~= chrid(j)	%skip markers on same chromosome
            
            %construct the design matrix
         	X = [ones(n,1), cov, geno(:,i)==1, geno(:,i)==2, geno(:,j)==1, geno(:,j)==2,...
               	geno(:,i)==1 & geno(:,j)==1, geno(:,i)==1 & geno(:,j)==2,...
                  geno(:,i)==2 & geno(:,j)==1, geno(:,i)==2 & geno(:,j)==2];
            
            %overall test
				[rss1, df1] = RSS(y, X);
            F1(i,j) = ((rss0-rss1)/(df0-df1))/(rss1/df1);
            
            %secondary tests
            if nargout > 1
				   [rss2, df2] = RSS(y, X(:,[1:c+5]));
               F2(i,j) = ((rss2-rss1)/(df2-df1))/(rss1/df1);
               
				   [rss3, df3] = RSS(y, X(:,[1:c+3]));
               F3(i,j) = ((rss3-rss2)/(df3-df2))/(rss2/df2);
               
				   [rss4, df4] = RSS(y, X(:,[1:c+1,c+4,c+5]));
               F4(i,j) = ((rss4-rss2)/(df4-df2))/(rss2/df2);
            end
      	end
   	end
	end
else	%backcross 
	for i = 1:m-1
   	for j = i+1:m
         if chrid(i) ~= chrid(j)	%skip markers on same chromosome
            
            %create the design matrix
         	X = [ones(n,1), cov, geno(:,i)==1, geno(:,j)==1, geno(:,i)==1 & geno(:,j)==1];
            
            %overall test
				[rss1, df1] = RSS(y, X);
            F1(i,j) = ((rss0-rss1)/(df0-df1))/(rss1/df1);
            
            %secondary tests
            if nargout > 1
				   [rss2, df2] = RSS(y, X(:,[1:c+3]));
               F2(i,j) = ((rss2-rss1)/(df2-df1))/(rss1/df1);
               
				   [rss3, df3] = RSS(y, X(:,[1:c+2]));
               F3(i,j) = ((rss3-rss2)/(df3-df2))/(rss2/df2);
               
				   [rss4, df4] = RSS(y, X(:,[1:c+1,c+3]));
               F4(i,j) = ((rss4-rss2)/(df4-df2))/(rss2/df2);
            end
      	end
   	end
	end
end
 
%%% SUBROUTINE RSS %%%
function [rss,df] = RSS(y,X)
[n,p] = size(X);
df = n - p;
[Q, R]=qr(X,0);
b = R\(Q'*y);
yhat = X*b;		
rss = (y-yhat)'*(y-yhat);