%% Assessing bias and variance in mixed effects model estimates using simulation
% This script simulates power and false positive rates given a set of input parameters.
% It strings together several scripts we've used so far. it is a "wrapper" for
% sim_mixedfx_fpr_power_bias_matlablme
% with different values for sample size.
% About sim_mixedfx_fpr_power_bias_matlablme
% 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.
% We've built up these functions over the last several scripts so you can see how they work and modify them.  They are designed to be modular.
% 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
%
% Table of Contents
% Specify input parameters
% Test run
% Power, FPR, bias, and variance as a function of number of subjects
% Activities
%
%
%% Specify input parameters

n_i = 20;           % len
m = 5;              % sub
p = 2;              % num params
b_g = [1 3]';       % fixed effects: Intercept and slope    
s2 = 25;            % residual var
U_g = [3 0; 0 9];   % Cov matrix of random intercept and slope
regressortype = 'continuous';

% These are params whose estimates we will compare later in the report
fixedint = b_g(1);
fixedslope = b_g(2); 
 
withinerr_std = sqrt(s2);

randintvar = U_g(1, 1);
randslopevar = U_g(2, 2);
 

%
%% Test run
% Generate a dataset, fit the model, and plot the results

iter = 50;
[sig_table, bias_variance_table] = sim_mixedfx_fpr_power_bias_matlablme(iter, m, n_i, p, b_g, s2, U_g, regressortype);

% [sig_table, bias_variance_table] = sim_mixedfx_fpr_power_bias_glmfitmulti(iter, m, n_i, p, b_g, s2, U_g, regressortype);

%
%% Power, FPR, bias, and variance as a function of number of subjects
% This cell block runs the simulation:

iter = 1000;
m_vals = [5 10 20 40 100 200];

[sig_tables, bv_tables] = deal({});
eval_time_per_model = zeros(1, length(m_vals));
b_g_null = zeros(size(b_g)); % global null w/no fixed effects

for i = 1:length(m_vals)

    timekeeper = clock;

    [sig_tables{i}, bv_tables{i}] = sim_mixedfx_fpr_power_bias_matlablme(iter, m_vals(i), n_i, p, b_g, s2, U_g, regressortype);

    eval_time_per_model(i) = etime(clock, timekeeper) ./ iter;

    fpr_tables{i} = sim_mixedfx_fpr_power_bias(iter, m_vals(i), n_i, p, b_g_null, s2, U_g, regressortype);

end

%
% This cell block aggregates results into a table and plots figures

%% Collect the results by statistic of interest

[power, fpr, dfe, bias_t, variance] = deal([]);

for i = 1:length(m_vals)

    power(:, i) = sig_tables{i}.percent_pos;
    fpr(:, i) = fpr_tables{i}.percent_pos;

    dfe(:, i) = sig_tables{i}.mean_dfe;
    
    bias_t(:, i) = bv_tables{i}.bias_t;

    variance(:, i) = bv_tables{i}.variance;

end

% Make a summary table
tablenames = {'N' 'power' 'fpr' 'dfe' 'bias_t' 'variance' 'eval_time'};
% rownames = mat2cell(m_vals, 1, ones(length(m_vals), 1));
% rownames = cellfun(@num2str, rownames, 'UniformOutput', false); 
summary_table = table(m_vals', power', fpr', dfe', bias_t', variance', eval_time_per_model', 'VariableNames', tablenames);
disp(summary_table)

% Plot the results 

figure; 
subplot(1, 2, 1); hold on; set(gca, 'FontSize', 16)
plot(m_vals, power, 'LineWidth', 3);
xlabel('Sample size (m)'); ylabel('Power');
legend({'Intercept' 'Slope'})
title('Power')

subplot(1, 2, 2); hold on; set(gca, 'FontSize', 16)
plot(m_vals, fpr, 'LineWidth', 3);
xlabel('Sample size (m)'); ylabel('False Pos Rate');
plot([min(m_vals) max(m_vals)], [0.05 0.05], 'k--');
title('False positive rate')

% Plot bias and var   

figure; 
subplot(1, 2, 1); hold on; set(gca, 'FontSize', 16)
plot(m_vals, bias_t, 'LineWidth', 2);
xlabel('Sample size (m)'); ylabel('Bias (t-score)');
legend(bv_tables{1}.Properties.RowNames')
plot([min(m_vals) max(m_vals)], [0 0], 'k--');
title('Bias in estimates')

subplot(1, 2, 2); hold on; set(gca, 'FontSize', 16)
plot(m_vals, variance, 'LineWidth', 2);
xlabel('Sample size (m)'); ylabel('Variance');
legend(bv_tables{1}.Properties.RowNames')
title('Variance in estimates')

subplot(1,2,1)
legend("Position", [0.25717,0.23647,0.19821,0.21548])

% Plot DFE and eval time   

figure; 
subplot(1, 2, 1); hold on; set(gca, 'FontSize', 16)
plot(m_vals, dfe, 'LineWidth', 2);
xlabel('Sample size (m)'); ylabel('Error DF');
plot([0 max(m_vals)], [0 max(m_vals)], 'k--');
title('Degrees of freedom')

subplot(1, 2, 2); hold on; set(gca, 'FontSize', 16)
plot(m_vals, eval_time_per_model, 'LineWidth', 2);
xlabel('Sample size (m)'); ylabel('Eval time');
title('Evaluation time per model')

%

% Save like this:

% save mixedfx_sim3_glmfitmultilevel summary_table sig_tables bv_tables



%
%% Activities
% Choose something else to vary or explore, and modify the simulation accordingly.
% Regressor type: Do categorical predictor variables (i.e., those from an experimental design) yield different power or FPR levels than continuous predictors?
% Number of trials within-person: Which is more important -- the number of participants or the number of trials within-person? Are there diminishing marginal returns of adding more of one or both?
% Error variance: Which has a bigger impact -- changing noise variance or variance of random effects? Why? What is the relationship between the two?
% Missing data: What is the impact on power and bias when some participants have unequal amounts of data (e.g., missing data)? What is the relative power of the full LME vs. the two-stage summary statistics approach as variation in missing data across participants increases?
% Heteroscedasticity: What is the impact on power and bias when some subjects have a larger error variance than others?
% Mismodeling: What is the impact of leaving a grouping variable (random effect) out of the model?  e.g., running an random intercept-only model if there is a random slope for participant?
% Compare performance of different models: Compare the power, bias, variance, and DFE above with those from a different model:
% Here is a summary of some different models you could select from:
% Two-stage summary stats:  Unweighted glmfit_multilevel
% Two-stage precision-weighted: weighted glmfit_multilevel
% Two-stage precision-weighted bootstrap: weighted glmfit_multilevel boot
% Matab MFX with REML
% LMER in R with ML
% LMER in R with REML
% statsmodels.lm in python
% bambi in python
% igls
% rigls
%
%