function [sig_table, bias_variance_table, input_params_table] = sim_mixedfx_fpr_power_bias_glmfitmulti(iter, m, n_i, p, b_g, s2, U_g, regressortype, varargin)

% sim_mixedfx_fpr_power_bias_matlablme is a wrapper that generates and fits data for each of a series of iterations, and then collects results (bias, variance, rates of positive findings).  It uses these functions:
%     sim_generate_mixedfx_data1 generates a dataset given parameters, and assuming normally distributed errors
%     matlabLMEstats runs the mixed effects model and collects selected output statistics
%     get_sim_stats collects tables of statistics across iterations
% The false positive rates in this simulation are assessed when the true fixed effects are all zero, and other random effects are unchanged.
% You can:
% 1 - Run this script with different inputs, including sample sizes (between/within), variances, fixed effects, regressor types
% 2 - Copy this script and modify the matlabLMEstats function to return the summary statistics for any mixed effects model you like
% NOTE: optional arguments (varargin) will be passed to glmfit_multilevel, so you can explore options.
%% Collect and display input params

input_params_table = table(iter, m, n_i, p, b_g', s2, {U_g}, {regressortype}, 'VariableNames', {'iter', 'm', 'n_i', 'p', 'b_g', 's2', 'U_g', 'regressortype'});
disp(input_params_table)
disp(U_g)

if ~isempty(varargin)
    disp('Optional inputs')
    disp(varargin{:})
    disp(' ')
end

%% Run simulation


% NOTE:
% iter = 500;

[fixedfx_hat, rfx_hat] = deal(zeros(2, iter));
s2_within = zeros(1, iter);

tic
for i = 1:iter

    % Generate data
    [y, Z_i, x, sim_data] = sim_generate_mixedfx_data1(m, n_i, p, b_g, s2, U_g, regressortype);

    % Fit model - Matlab LME
    [fixedfx_hat(:, i), rfx_hat(:, i), s2_within(:, i), pval(:, i), dfe(:, i), sig(:, i)] = ...
        glmfitstats(y, x, b_g, varargin{:});

    if i == 1
        % run once, print output
        disp('Sample output')
        stats = glmfit_multilevel(y, x, [], varargin{:});
    end

    if mod(i, 100) == 0
        toc
    end

end

toc

[sig_table, bias_variance_table] = get_sim_stats(b_g, U_g, s2, sig, iter, pval, dfe, fixedfx_hat, rfx_hat, s2_within);

end % function




%

function [fixedfx_hat, rfx_hat, s2_within, pval, dfe, sig] = glmfitstats(y, x, b_g, varargin)
% Bundle everything model-specific here
% The rest of the code can run on any model that returns these outputs

% This runs glmfit_multilevel instead of fitlme
stats = glmfit_multilevel(y, x, [], 'noverbose', varargin{:});

%lme = fitlme(sim_data,'Y ~ x01 + (1 + x01|sid)');
% [V_est, s2_est] = covarianceParameters(lme);

% Save estimates to check bias and variance
fixedfx_hat = stats.beta'; % fixedEffects(lme); % out.beta
rfx_hat = stats.var';  % diag(V_est{1}); % out.betastar;
s2_within = NaN; % s2_est; % out.Sigma;

% P-values
pval = stats.p';  % test.pValue;
%test = anova(lme, 'DFMethod', 'Satterthwaite');

dfe = stats.dfe'; % test.DF2;

% if b_g ~= 0, sig reflects power, and we need sign to match.
sig = pval < 0.05 & sign(fixedfx_hat) == sign(b_g);  % P < 0.05 and correct sign
% if b_g = 0, sign reflects FPR, and does not need to match
sig(b_g == 0) = pval(b_g == 0) < 0.05;

end % glmfitstats function

%

function [sig_table, bias_variance_table] = get_sim_stats(b_g, U_g, s2, sig, iter, pval, dfe, fixedfx_hat, rfx_hat, s2_within)

%% Get power/FPR and DFE

percent_pos = sum(sig, 2) ./ iter;  % fpr or power, depending on whether b_g is 0 or not for a parameter.
mean_pval = mean(pval, 2);
mean_dfe = mean(dfe, 2);

sig_table = table(percent_pos, mean_pval, mean_dfe, 'VariableNames', {'percent_pos' 'mean_pval' 'mean_dfe'}, 'RowNames', {'FixedInt' 'FixedSlope'});

%% Get Bias and variance of estimates

% Calculate means, bias (deviations), and variances for estimates across
% iterations

fixedfxmean = mean(fixedfx_hat, 2);
rfxmean = mean(rfx_hat, 2);
s2mean = mean(s2_within, 2);

fixedfxbias = (fixedfx_hat - b_g)';
rfxbias = (rfx_hat - diag(U_g))';
s2bias = (s2_within - s2)';

% Significance test for estimation bias
% Does the model produce the true values on average?
% A significant P-value here indicates systematic bias, which is bad for the model.

[~, ~, ~, fixedfx_bias_stats] = ttest(fixedfxbias);
[~, ~, ~, rfx_bias_stats] = ttest(rfxbias);
[~, ~, ~, s2_bias_stats] = ttest(s2bias);

% Variance of estimates across iterations
% Lower is better

fixedfx_err_var = var(fixedfx_hat, 0, 2);
rfx_err_var = var(rfx_hat, 0, 2);
s2_err_var = var(s2_within, 0, 2);

% Make output table
bias = [mean(fixedfxbias)'; mean(rfxbias)'; mean(s2bias)'];
bias_t = [fixedfx_bias_stats.tstat'; rfx_bias_stats.tstat'; s2_bias_stats.tstat'];
variance = [fixedfx_err_var; rfx_err_var; s2_err_var];

bias_variance_table = table(bias, bias_t, variance, 'VariableNames', {'bias' 'bias_t' 'variance'}, 'RowNames', {'FixedInt' 'FixedSlope' 'RandInt' 'RandSlope' 'ErrVar'});

end % getstats function