function [F1, F2, F3] = covscan(y, cov, geno)
%scan a marker array for associations with a trait
%covscan uses an unconstrained regression model
%INPUT VARIABLES ARE:
%     y: nx1 vector of trait values
%   cov: nx1 vector of covariates,  use [] if none
%  geno: nxm genotype array. 
%OUTPUT VARIABLES
%  F1:	the overall F-test, full model vs null model
%	F2:	the interaction test, full model vs additive model
%	F3:	the coat-tail test, additive model vs covariate only


%constants and storage
[n,m] = size(geno);
ncov = size(cov,2);
F1 = zeros(m,1);
F2 = zeros(m,1);
F3 = zeros(m,1);

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

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

%scan the genome doing regression at each marker				
for i = 1:m
   if intercross
      X1 = [X0, geno(:,i)==1, geno(:,i)==2,];
      X2 = [X1, repmat((geno(:,i)==1),1,ncov).*cov, repmat((geno(:,i)==2),1,ncov).*cov];
   else
      X1 = [X0, geno(:,i)==1];
      X2 = [X1, repmat((geno(:,i)==1),1,ncov).*cov];
   end
   [rss1, df1] = RSS(y, X1);
   [rss2, df2] = RSS(y, X2);

   F1(i) = ((rss0-rss2)/(df0-df2))/(rss2/df2);
   F2(i) = ((rss1-rss2)/(df1-df2))/(rss2/df2);
   F3(i) = ((rss0-rss1)/(df0-df1))/(rss1/df1);
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);