function [F1, F2, F3, F4] = covpairwise(y, cov, geno, chrid)
%genome scan for joint effects of two markers with a covariate
%this program computes the four test statistics 
%used to compare models of uniform order
%test 1 compares the covariate-only model to the full three-way interaction model
%test 2 compare covariate-only model to an addtive model with the two qtl
%tests 3 compares the additive model to all pairwise interactions models
%test 4 compares the all-pairwise model to the three-way model
%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 zero order vs first order test
%	F3:	the first order vs second order test
%  F4:   the the second order vs three-way test

%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)
   exit('requires a non-empty cov');
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
            X1 = [X0, geno(:,i)==1, geno(:,i)==2, geno(:,j)==1, geno(:,j)==2];
            X2 = [X1,geno(:,i)==1 & geno(:,j)==1, geno(:,i)==1 & geno(:,j)==2,...
                  geno(:,i)==2 & geno(:,j)==1, geno(:,i)==2 & geno(:,j)==2,...
                  (geno(:,i)==1).*cov, (geno(:,i)==2).*cov];
            X3 = [X2, (geno(:,i)==1 & geno(:,j)==1) .* cov,...
               	(geno(:,i)==1 & geno(:,j)==2) .* cov,...
                  (geno(:,i)==2 & geno(:,j)==1) .* cov,...
                  (geno(:,i)==2 & geno(:,j)==2) .* cov];
            
            %overall test
				[rss0, df0] = RSS(y, X0);
            [rss3, df3] = RSS(y, X3);
            F1(i,j) = ((rss3-rss0)/(df3-df0))/(rss3/df3);

           
            %secondary tests
            if nargout > 1
               [rss1, df1] = RSS(y, X1);
               F2(i,j) = ((rss1-rss0)/(df1-df0))/(rss1/df1);
               
				   [rss2, df2] = RSS(y, X2);
               F3(i,j) = ((rss2-rss1)/(df2-df1))/(rss2/df2);
               
               F4(i,j) = ((rss3-rss2)/(df3-df2))/(rss3/df3);
            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
            
            %construct the design matrix
            X1 = [X0, geno(:,i)==1, geno(:,j)==1];
            X2 = [X1,geno(:,i)==1 & geno(:,j)==1, (geno(:,i)==1).*cov];
            X3 = [X2, (geno(:,i)==1 & geno(:,j)==1) .* cov];
            
            %overall test
				[rss0, df0] = RSS(y, X0);
            [rss3, df3] = RSS(y, X3);
            F1(i,j) = ((rss3-rss0)/(df3-df0))/(rss3/df3);

            %secondary tests
            if nargout > 1
               [rss1, df1] = RSS(y, X1);
               F2(i,j) = ((rss1-rss0)/(df1-df0))/(rss1/df1);
               
				   [rss2, df2] = RSS(y, X2);
               F3(i,j) = ((rss2-rss1)/(df2-df1))/(rss2/df2);
               
               F4(i,j) = ((rss3-rss2)/(df3-df2))/(rss3/df3);
            end
      	end
   	end
	end
end

%[df0 df1 df2 df3 df4]

%%% 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);