function [BETA_COEF,dof,E,PValue] = bnpm_glm(data,confound,col_conf, brain_mask, type,no_subj,nr,BETA_COEF,E,dof,PValue)
%-----------------------------------------------------------------------%
%                         BnPM GLM                                      %                         
%                                                                       %
%  It fits a GLM voxel-wise including imaging covariates for one slice  %
%  Non-imaging covariates  can be included in the analysis.             %
%                                                                       %
%-----------------------------------------------------------------------%
% nonparametric method, haven't consider confound covariat

[n,ng] = size(data); % ng - number of groups
[n,nc] = size(confound); % nc - number of confounds
Tsubj  = sum(no_subj);
[M,N,s] = size( data{1} );
Data = [];

if strcmp(type,'ANCOVA')
    X = wfu_bpm_design_mat(no_subj); 
end

if strcmp(type,'REGRESSION')
    X = ones(no_subj,1);
end

%----- Append non-imaging confounds -----------%
if ~isempty(col_conf)      
    X = [X col_conf];
end

% ----------- The analysis is performed on voxel by voxel basis --------- %

[I,J]          = find(brain_mask > 0);
for kmv = 1:length(I)
    
    m = I(kmv);
    n = J(kmv);    
    if ~isempty(confound)
        X1 = [];   % for confounds
        for c = 1:nc
            x1 = [];
            for g = 1:ng
                x11 = confound{c}{g}(m,n,:);
                x1 = [x1;x11(:)];
            end 
            X1 = [X1 x1];
        end
        X2 = [X X1];     % complete Design Matrix
    else
        X2 = X;
    end
    % ------- creating the data vector ---------- %
    y = [];
    for g  = 1:ng
        y1 = data{g}(m,n,:);
        y  = [y;y1(:)];
    end
    
    X3 = X2;
    X3(:,1) = [];
    [beta stats] = nonpreg(X3,y);
    BETA_COEF(m,n,:) = beta     ;
    e = sqrt(stats.res);
    Rscan            = size(y)       ;
    Tsubj            = Rscan(1)      ;
    p                = rank(X2)      ;
    E(m,n,:)         = e        ;
    dof(m,n) = Tsubj - p;
    PValue(m,n,:) = stats.p;
end
end

