% function R=phyl_resid(C, x, Y) computes R residuals given C, x, and Y
function R=phyl_resid(C, x, Y)
    
    % find out how many y variables and taxa we have
    [n,m]=size(Y);

    % prepare matrix containing size
    X=ones(n,1); X(:,2)=x;

    % now loop over those variables, each time calculating the regression &
    % residuals

    for i=1:m
        % estimate beta
        beta=(X'*C^-1*X)^-1*(X'*C^-1*Y(:,i));
        % compute residuals and put in ith column of R
        R(:,i)=Y(:,i)-X*beta;
    end
    
% done