function [r2] = comp_R2(y,X,indiv)
    n = size(X,1);
    p = size(X,2);

    x = [ones(n,1) X];

    if isempty(indiv)
        indiv =0;
    end

    if indiv == 0
        beta = (x'*x)\(x'*y);
        res  = y-x*beta;
        ssq  = sum(res.^2);
        
        tss  = sum((y-mean(y)).^2);
        
        r2 = 1-ssq./tss;
    else
        r2 = zeros(1,p);
        for i = 1:p
            beta = (x(:,[1 i+1])'*x(:,[1 i+1]))\(x(:,[1 i+1])'*y(:,i));
            res  = y(:,i)-x(:,[1 i+1])*beta;
            ssq  = sum(res.^2);
            
            tss  = sum((y(:,i)-mean(y(:,i))).^2);
            
            r2(:,i) = 1-ssq./tss;
        end

    end
end