%% Ashar 2022 LME example
% We'll fit the Ashar dataset with fitlme, Matlab's mixed effects model
% This is Matlab's object-oriented equivalent to lmer in R. It fits a variety of mixed effects models. A list of methods and properties are here:

help LinearMixedModel

% These models, like many other statistics objects implementing linear models, allow you to specify model formulas using Wilkinson notation. There is also a matrix input version.
%% Load the dataset
% This dataset is from:
% Ashar et al. 2022, JAMA Psychiatry
% This study is a randomized clinical trial (RCT) of patients with chronic back pain (>= 4/10 pain for at least 3 months, average duration ~10 years). The trial compares Pain Reprocessing Therapy, Placebo (a single injection into the skin of the back), and Usual Care treatments, with pre-treatment and post-treatment (1 month after treatment initiation), and follow-up to 1 year.

backpain = readtable('Ashar_backpain_simple_long.csv');
backpain

% Remember: You need the dataset on your Matlab path for this to work.
%
%% Fit the model
% FIt a model testing the Group x Time (Post - Pre-treatment) interaction.

lme = fitlme(backpain,'bpi_intensity ~ group*time + (time | id)', 'FitMethod','REML');
disp(lme);

%
%% Diagnose problems!
% There are few things wrong here. What are they?
% It's a great idea to do "sanity checks" to understand what the model is doing and what might be going wrong.
% Hints:
% Look at the dfe. Does it look right?
% Look at the design matrix and the number of fixed effects. Is it as intended?
% Does it makes sense to include healthy individuals scanned only once? Are they being included?
% Are we effectively using the "difference score" method here? Is there another way?
%
%% Issue 1: Fixed effects dfe
% The dfe for the fixed effects is very large - much larger than the number of subjects, and in fact comparable to the total number of observations.
% The fixed effects table does not correct the dfe for effective degrees of freedom, and does not therefore account for reduction in dfe with correlated errors due to random effects!
% To estimate the dfe, we can used the Satterthwaite approximation.

[~, ~, fixedstats] = fixedEffects(lme, 'DFMethod', 'Satterthwaite');
disp(fixedstats)

% Contrasts: Simple (trivial) contrast for the Group x Time interaction
disp('Group x Time')
[PVAL,F,DF1,DF2] = coefTest(lme, [0 0 0 1], 0, 'DFMethod', 'Satterthwaite')

%
%% Issue 2: Mis-coding a predictor
% Look at the design matrix.

X = designMatrix(lme, 'Fixed');

X_table = array2table(X, 'VariableNames', lme.CoefficientNames);
X_table

% There is only 1 column for Group, but we know we need 2 columns to span the space of differences among 3 groups!
% Group is coded with numbers, and is being treated as a numeric variable. It is looking for linear increases (or decreases) in pain as function of the numeric code! If we want lme to treat Group as a categoricall variable, we need to recast it as a categorical variable class.

backpain.group = categorical(backpain.group);

%
% Now we can re-run the model:

lme = fitlme(backpain,'bpi_intensity ~ group*time + (time | id)', 'FitMethod','REML');
disp(lme);

%
% The Group codes are rather awkward. Group 1 is a reference group, so the positive effect of time:group_2 means that there was a greater increase in pain over time in Group 2 (Placebo) than Group 1 (PRT, treatment).
%  How do we do an F-test for Group overall?  The anova() method will do it:

anova(lme, 'DFMethod', 'Satterthwaite')

%
%% Issue 3: The difference score method
% The effect of Time here is the equivalent of the Post vs. Pre difference score method we used earlier, when we created contrasts across the within-person effect of time.
% This kind of subtraction can be super useful, as there may be a lot of individual variability in overall pain scores, and by collecting a baseline score we can "subtract it away", reducing error.  And it is the improvement we are interested in, after all.
% The subtraction method is widely used -- for example, subtracting high vs. low cognitive demand to yield a difference score. The Stroop, Simon, Flanker, Go/No-go, IAT, and many more tests use this method.
% There are a couple of problems with this:
% Interpretability: Treatment is really expected to affect the POST-treatment scores, not the pre-treatment scores. But the contrast is sensitive to variation in both pre- and post-treatment scores.  A negative [post - pre] x [PRT vs. Placebo] score could be driven by lower (better) pain scores post-treatment in the PRT group OR higher baseline scores in the Placebo group.  This makes the contrast hard to interpret. If I find a significant effect, I should check that it's not "driven by" pre-treatment differences.
% Error variance: Random varation in the pre-treatment scores will make the difference scores more noisy too, resulting in higher error variance and lower power. According to the variance sum law, the variance of the difference between two random variables var(A - B) is the sum of var(A) + var(B)!! So difference scores are extra noisy: We've added var(pre-treatment scores) to the variance of our outcome unnecessarily.
% Unwarranted assumptions: The subtraction method assumes a slope of 1 across repeated measures, i.e., pre- and post-scores. This works OK if pre- and post-scores are highly correlated. But they may not be, as there are many factors that could reduce the correlation, including error.  To the extent they are not perfectly correlated, you're subtracting off a random value, increasing noise.
% A solution: Pre-treatment pain as a covariate. Instead, if you control for Pre-treatment pain by including it as a covariate, you estimate a \hat\beta for it. This implicitly subtracts off a fraction of the pre-treatment score (the standardized \hat\beta must be < 1, because it is 1 for a perfect correlation) in direct proportion to the Pre-Post correlation. If they correlation is close to 1, the covariate will be close to the subtraction method. But if it's low (close to 0), the adjustment will be smaller and a fraction of Pre-treatment scores will be subtracted.

% For this the "wide" format is more convenient
backpainwide = readtable('Ashar_backpain_simple_wide.csv');

backpainwide.group = categorical(backpainwide.group);

lme = fitlme(backpainwide,'BPI_Time2 ~ group + BPI_Time1');
disp(lme);

anova(lme, 'DFMethod', 'Satterthwaite')

% Of course, now I have one vector of independent errors (they are not stacked across time points within-person), so there are no random effects. Or, we should say, errors are not grouped by a random factor (participant). So we do not need a mixed effects model to do this: Any OLS GLM will work.
%