function L = tdcs_lme(S, cfg)
%TDCS_LME Linear mixed-effects "days x tDCS" model on the Box-B2 vs Box-A2 arm.
%   L = TDCS_LME(S, CFG) subsets one scenario table S to the two anchor groups,
%   codes a binary tDCS factor (Box-B2 = 1, Box-A2 = 0), and fits the LINEAR
%   mixed model
%
%       success ~ day * tDCS + (1|subject)
%
%   replicating the published formulation (a linear mixed-effects model for the
%   number of successful reaches with a days x tDCS interaction, a main effect
%   of days, and a main effect of tDCS). `day` is raw and 0-indexed, and OUR
%   day 0 corresponds to the paper's "Day 1", so the tDCS main effect is the
%   group difference on Day 1 -- directly comparable to the paper's "equal on
%   Day 1" test. (Do not re-index day to 1-based: that would move the main
%   effect's evaluation point off Day 1.)
%
%   L fields:
%     .lme         the LinearMixedModel object
%     .nSubjects   number of subjects (Box-A2 + Box-B2)
%     .nObs        number of sessions
%     .interaction / .day / .tDCS   effect structs, each with .estimate .se
%                  .t .df .p (from Coefficients) and .F .df1 .df2 .Fp (from
%                  ANOVA). For these single-df terms F == t^2 and Fp == p.
%
%   NOTE: this is a LINEAR mixed model (Gaussian on the raw count), matching
%   the paper's method, unlike the Poisson/Binomial GLMMs in tdcs_models.

T = S(ismember(S.group, {cfg.anchorLow, cfg.anchorHigh}), :);
T.group = removecats(T.group);
T.tDCS = double(T.group == cfg.anchorHigh);   % Box-B2 = 1, Box-A2 = 0

lme = fitlme(T, 'success ~ day*tDCS + (1|subject)');
C = lme.Coefficients;
A = anova(lme);

L.lme = lme;
L.nSubjects = numel(unique(T.subject));
L.nObs = height(T);
L.maxDayA = max(T.day(T.group == cfg.anchorLow));    % last day Box-A2 has data
L.maxDayB = max(T.day(T.group == cfg.anchorHigh));   % last day Box-B2 has data
L.interaction = localTerm(C, A, 'day:tDCS');
L.day = localTerm(C, A, 'day');
L.tDCS = localTerm(C, A, 'tDCS');

end

function e = localTerm(C, A, name)
%LOCALTERM Pull one term's coefficient (t/df/p) and ANOVA (F/df/p) stats.
ci = strcmp(C.Name, name);
ai = strcmp(A.Term, name);
e = struct( ...
    'estimate', C.Estimate(ci), 'se', C.SE(ci), ...
    't', C.tStat(ci), 'df', C.DF(ci), 'p', C.pValue(ci), ...
    'F', A.FStat(ai), 'df1', A.DF1(ai), 'df2', A.DF2(ai), 'Fp', A.pValue(ai));
end