function [out] = TimeCourseRunRegression(regs,timec,Data)
% runs for each subject, for each time point a separate regression with all
% regressors (regs)

% some lengths extracted here to make code below more compact
ntimep = size(timec.Region(1).v{1},2);   % number of timepoints in each epoch
nregion= size(timec.Region,2);           % number of regions
nregres= size(regs.DesignMatrix{1},2);              % number of regressors

% run the regressions
    betas = nan(length(Data.SubjectIDs),nregion,ntimep,nregres+1); %subs, region, timepoints, regressors
 
for iregion=1:nregion
    for is =  1:length(Data.SubjectIDs) % for each participant
        for ip = 1:ntimep % for each time point
            ydat  = timec.Region(iregion).v{is}(:,ip);   % ydat for the regression, i.e. the intrapolated BOLD values 
            indtr = size(ydat,1); % how many fmri trials we have for each person
            preds=regs.DesignMatrix{is};                        
            betas(is,iregion,ip,1:(nregres+1)) = glmfit(preds,ydat); %running fitglm instead would give more information, but it takes longer to run (and we only need the regression weights
            % Substitute this with the hierarchical fit?            
            
        end        
    end
end



out.betas = betas;

out.labels= regs.Labels;
out.descr = 'subject x region x time point x coefficients';
out.regions= {timec.Region.name};
