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
LinearMixedModel Fitted linear mixed effects model.
LME = FITLME(...) fits a linear mixed effects model to data. The fitted
model LME is a LinearMixedModel modeling a response variable as a
linear function of fixed effect predictors and random effect
predictors.
LinearMixedModel methods:
coefCI - Coefficient confidence intervals
coefTest - Linear hypothesis test on coefficients
predict - Compute predicted values given predictor values
random - Generate random response values given predictor values
partialDependence - Partial Dependence values.
plotPartialDependence - Plot Partial Dependence.
plotResiduals - Plot of residuals
designMatrix - Fixed and random effects design matrices
fixedEffects - Stats on fixed effects
randomEffects - Stats on random effects
covarianceParameters - Stats on covariance parameters
fitted - Fitted response
residuals - Various types of residuals
response - Response used to fit the model
compare - Compare fitted models
anova - Marginal tests for fixed effect terms
disp - Display a fitted model
LinearMixedModel properties:
FitMethod - Method used for fitting (either 'ML' or 'REML')
MSE - Mean squared error (estimate of residual variance)
Formula - Representation of the model used in this fit
LogLikelihood - Log of likelihood function at coefficient estimates
DFE - Degrees of freedom for residuals
SSE - Error sum of squares
SST - Total sum of squares
SSR - Regression sum of squares
CoefficientCovariance - Covariance matrix for coefficient estimates
CoefficientNames - Coefficient names
NumCoefficients - Number of coefficients
NumEstimatedCoefficients - Number of estimated coefficients
Coefficients - Coefficients and related statistics
Rsquared - R-squared and adjusted R-squared
ModelCriterion - AIC and other model criteria
VariableInfo - Information about variables used in the fit
ObservationInfo - Information about observations used in the fit
Variables - Table of variables used in fit
NumVariables - Number of variables used in fit
VariableNames - Names of variables used in fit
NumPredictors - Number of predictors
PredictorNames - Names of predictors
ResponseName - Name of response
NumObservations - Number of observations in the fit
ObservationNames - Names of observations in the fit
See also fitlme, fitlmematrix, LinearModel, GeneralizedLinearModel, NonLinearModel.
Documentation for 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
backpain = 341×6 table
| | id | time | group | bpi_intensity | gender | is_patient |
|---|
| 1 | 12 | 1 | 3 | 2.5000 | 1 | 1 |
|---|
| 2 | 12 | 2 | 3 | 2.7500 | 1 | 1 |
|---|
| 3 | 14 | 1 | 3 | 2.5000 | 2 | 1 |
|---|
| 4 | 14 | 2 | 3 | 2 | 2 | 1 |
|---|
| 5 | 15 | 1 | 3 | 2.2500 | 1 | 1 |
|---|
| 6 | 15 | 2 | 3 | 2.7500 | 1 | 1 |
|---|
| 7 | 18 | 1 | 2 | 4.5000 | 1 | 1 |
|---|
| 8 | 18 | 2 | 2 | 2.2500 | 1 | 1 |
|---|
| 9 | 23 | 1 | 2 | 2.5000 | 2 | 1 |
|---|
| 10 | 23 | 2 | 2 | 2.2500 | 2 | 1 |
|---|
| 11 | 24 | 1 | 1 | 4.5000 | 1 | 1 |
|---|
| 12 | 24 | 2 | 1 | 0.5000 | 1 | 1 |
|---|
| 13 | 29 | 1 | 1 | 3.7500 | 2 | 1 |
|---|
| 14 | 29 | 2 | 1 | 0.2500 | 2 | 1 |
|---|
| 15 | 32 | 1 | 1 | 5.7500 | 1 | 1 |
|---|
| 16 | 32 | 2 | 1 | 0.2500 | 1 | 1 |
|---|
| 17 | 38 | 1 | 1 | 5.5000 | 2 | 1 |
|---|
| 18 | 38 | 2 | 1 | 1.7500 | 2 | 1 |
|---|
| 19 | 51 | 1 | 2 | 4 | 1 | 1 |
|---|
| 20 | 51 | 2 | 2 | 3.7500 | 1 | 1 |
|---|
| 21 | 61 | 1 | 2 | 4.5000 | 2 | 1 |
|---|
| 22 | 61 | 2 | 2 | 1.2500 | 2 | 1 |
|---|
| 23 | 67 | 1 | 1 | 5.2500 | 2 | 1 |
|---|
| 24 | 69 | 1 | 2 | 1.2500 | 1 | 1 |
|---|
| 25 | 69 | 2 | 2 | 0.5000 | 1 | 1 |
|---|
| 26 | 72 | 1 | 3 | 4 | 1 | 1 |
|---|
| 27 | 72 | 2 | 3 | 4.5000 | 1 | 1 |
|---|
| 28 | 79 | 1 | 3 | 4.5000 | 2 | 1 |
|---|
| 29 | 79 | 2 | 3 | 3.2500 | 2 | 1 |
|---|
| 30 | 85 | 1 | 2 | 1 | 1 | 1 |
|---|
| 31 | 87 | 1 | 2 | 4.5000 | 2 | 1 |
|---|
| 32 | 87 | 2 | 2 | 2.2500 | 2 | 1 |
|---|
| 33 | 96 | 1 | 2 | 6.2500 | 2 | 1 |
|---|
| 34 | 96 | 2 | 2 | 6.7500 | 2 | 1 |
|---|
| 35 | 102 | 1 | 1 | 5.7500 | 1 | 1 |
|---|
| 36 | 102 | 2 | 1 | 3 | 1 | 1 |
|---|
| 37 | 103 | 1 | 2 | 3.7500 | 2 | 1 |
|---|
| 38 | 103 | 2 | 2 | 2.7500 | 2 | 1 |
|---|
| 39 | 104 | 1 | 3 | 2 | 2 | 1 |
|---|
| 40 | 104 | 2 | 3 | 1.7500 | 2 | 1 |
|---|
| 41 | 108 | 1 | 3 | 2.2500 | 2 | 1 |
|---|
| 42 | 108 | 2 | 3 | 2 | 2 | 1 |
|---|
| 43 | 112 | 1 | 1 | 2 | 1 | 1 |
|---|
| 44 | 112 | 2 | 1 | 1 | 1 | 1 |
|---|
| 45 | 113 | 1 | 1 | 5.2500 | 2 | 1 |
|---|
| 46 | 113 | 2 | 1 | 3.5000 | 2 | 1 |
|---|
| 47 | 115 | 1 | 1 | 3.7500 | 2 | 1 |
|---|
| 48 | 115 | 2 | 1 | 0.7500 | 2 | 1 |
|---|
| 49 | 117 | 1 | 1 | 5.2500 | 2 | 1 |
|---|
| 50 | 117 | 2 | 1 | 4 | 2 | 1 |
|---|
| 51 | 125 | 1 | 2 | 3 | 2 | 1 |
|---|
| 52 | 125 | 2 | 2 | 2.5000 | 2 | 1 |
|---|
| 53 | 142 | 1 | 1 | 2 | 1 | 1 |
|---|
| 54 | 142 | 2 | 1 | 0.5000 | 1 | 1 |
|---|
| 55 | 146 | 1 | 2 | 3 | 2 | 1 |
|---|
| 56 | 146 | 2 | 2 | 3.7500 | 2 | 1 |
|---|
| 57 | 147 | 1 | 3 | 4.7500 | 2 | 1 |
|---|
| 58 | 147 | 2 | 3 | 2 | 2 | 1 |
|---|
| 59 | 148 | 1 | 2 | 3.5000 | 1 | 1 |
|---|
| 60 | 148 | 2 | 2 | 2 | 1 | 1 |
|---|
| 61 | 165 | 1 | 2 | 3.7500 | 1 | 1 |
|---|
| 62 | 172 | 1 | 2 | 3.7500 | 1 | 1 |
|---|
| 63 | 172 | 2 | 2 | 5.5000 | 1 | 1 |
|---|
| 64 | 192 | 1 | 3 | 4.7500 | 2 | 1 |
|---|
| 65 | 192 | 2 | 3 | 5 | 2 | 1 |
|---|
| 66 | 194 | 1 | 3 | 3.2500 | 2 | 1 |
|---|
| 67 | 194 | 2 | 3 | 3 | 2 | 1 |
|---|
| 68 | 195 | 1 | 1 | 2.7500 | 2 | 1 |
|---|
| 69 | 195 | 2 | 1 | 0 | 2 | 1 |
|---|
| 70 | 219 | 1 | 3 | 4.5000 | 1 | 1 |
|---|
| 71 | 221 | 1 | 2 | 1 | 1 | 1 |
|---|
| 72 | 221 | 2 | 2 | 1 | 1 | 1 |
|---|
| 73 | 228 | 1 | 3 | 2.5000 | 1 | 1 |
|---|
| 74 | 228 | 2 | 3 | 1.2500 | 1 | 1 |
|---|
| 75 | 248 | 1 | 2 | 5.2500 | 1 | 1 |
|---|
| 76 | 248 | 2 | 2 | 4.5000 | 1 | 1 |
|---|
| 77 | 262 | 1 | 1 | 3.2500 | 2 | 1 |
|---|
| 78 | 262 | 2 | 1 | 0.7500 | 2 | 1 |
|---|
| 79 | 295 | 1 | 3 | 3.2500 | 2 | 1 |
|---|
| 80 | 295 | 2 | 3 | 4.2500 | 2 | 1 |
|---|
| 81 | 312 | 1 | 1 | 3.2500 | 1 | 1 |
|---|
| 82 | 312 | 2 | 1 | 0.7500 | 1 | 1 |
|---|
| 83 | 319 | 1 | 1 | 5 | 1 | 1 |
|---|
| 84 | 319 | 2 | 1 | 0.2500 | 1 | 1 |
|---|
| 85 | 324 | 1 | 1 | 5.5000 | 1 | 1 |
|---|
| 86 | 324 | 2 | 1 | 1.7500 | 1 | 1 |
|---|
| 87 | 330 | 1 | 3 | 4.5000 | 1 | 1 |
|---|
| 88 | 330 | 2 | 3 | 2 | 1 | 1 |
|---|
| 89 | 339 | 1 | 3 | 4.5000 | 2 | 1 |
|---|
| 90 | 341 | 1 | 2 | 4 | 1 | 1 |
|---|
| 91 | 341 | 2 | 2 | 2.5000 | 1 | 1 |
|---|
| 92 | 348 | 1 | 3 | 3.5000 | 2 | 1 |
|---|
| 93 | 348 | 2 | 3 | 2 | 2 | 1 |
|---|
| 94 | 352 | 1 | 3 | 5.2500 | 2 | 1 |
|---|
| 95 | 364 | 1 | 3 | 5.5000 | 2 | 1 |
|---|
| 96 | 364 | 2 | 3 | 1.7500 | 2 | 1 |
|---|
| 97 | 378 | 1 | 1 | 3 | 2 | 1 |
|---|
| 98 | 378 | 2 | 1 | 1 | 2 | 1 |
|---|
| 99 | 389 | 1 | 3 | 3 | 2 | 1 |
|---|
| 100 | 389 | 2 | 3 | 4 | 2 | 1 |
|---|
| ⋮ |
|---|
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);
Linear mixed-effects model fit by REML
Model information:
Number of observations 286
Fixed effects coefficients 4
Random effects coefficients 302
Covariance parameters 4
Formula:
bpi_intensity ~ 1 + time*group + (1 + time | id)
Model fit statistics:
AIC BIC LogLikelihood Deviance
958.91 988.05 -471.46 942.91
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)'} 7.4115 0.48139 15.396 282 3.0697e-39 6.4639 8.3591
{'time' } -3.4186 0.29129 -11.736 282 3.6424e-26 -3.992 -2.8452
{'group' } -1.1237 0.22222 -5.0569 282 7.6861e-07 -1.5611 -0.68631
{'time:group' } 1.0362 0.1337 7.7504 282 1.6691e-13 0.77305 1.2994
Random effects covariance parameters (95% CIs):
Group: id (151 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std' } 0.97441 NaN NaN
{'time' } {'(Intercept)'} {'corr'} -0.14395 NaN NaN
{'time' } {'time' } {'std' } 0.33244 NaN NaN
Group: Error
Name Estimate Lower Upper
{'Res Std'} 0.87779 NaN NaN
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)
Fixed effect coefficients: DFMethod = 'Satterthwaite', Alpha = 0.05
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)'} 7.4115 0.48139 15.396 145.67 2.331e-32 6.4601 8.3629
{'time' } -3.4186 0.29129 -11.736 136.4 2.0222e-22 -3.9946 -2.8426
{'group' } -1.1237 0.22222 -5.0569 144.41 1.273e-06 -1.5629 -0.6845
{'time:group' } 1.0362 0.1337 7.7504 135.71 1.9285e-12 0.77182 1.3006
% Contrasts: Simple (trivial) contrast for the Group x Time interaction
[PVAL,F,DF1,DF2] = coefTest(lme, [0 0 0 1], 0, 'DFMethod', 'Satterthwaite')
PVAL = 1.9285e-12
F = 60.0683
DF1 = 1
DF2 = 135.7109
Issue 2: Mis-coding a predictor
Look at the design matrix.
X = designMatrix(lme, 'Fixed');
X_table = array2table(X, 'VariableNames', lme.CoefficientNames);
X_table
X_table = 341×4 table
| | (Intercept) | time | group | time:group |
|---|
| 1 | 1 | 1 | 3 | 3 |
|---|
| 2 | 1 | 2 | 3 | 6 |
|---|
| 3 | 1 | 1 | 3 | 3 |
|---|
| 4 | 1 | 2 | 3 | 6 |
|---|
| 5 | 1 | 1 | 3 | 3 |
|---|
| 6 | 1 | 2 | 3 | 6 |
|---|
| 7 | 1 | 1 | 2 | 2 |
|---|
| 8 | 1 | 2 | 2 | 4 |
|---|
| 9 | 1 | 1 | 2 | 2 |
|---|
| 10 | 1 | 2 | 2 | 4 |
|---|
| 11 | 1 | 1 | 1 | 1 |
|---|
| 12 | 1 | 2 | 1 | 2 |
|---|
| 13 | 1 | 1 | 1 | 1 |
|---|
| 14 | 1 | 2 | 1 | 2 |
|---|
| 15 | 1 | 1 | 1 | 1 |
|---|
| 16 | 1 | 2 | 1 | 2 |
|---|
| 17 | 1 | 1 | 1 | 1 |
|---|
| 18 | 1 | 2 | 1 | 2 |
|---|
| 19 | 1 | 1 | 2 | 2 |
|---|
| 20 | 1 | 2 | 2 | 4 |
|---|
| 21 | 1 | 1 | 2 | 2 |
|---|
| 22 | 1 | 2 | 2 | 4 |
|---|
| 23 | 1 | 1 | 1 | 1 |
|---|
| 24 | 1 | 1 | 2 | 2 |
|---|
| 25 | 1 | 2 | 2 | 4 |
|---|
| 26 | 1 | 1 | 3 | 3 |
|---|
| 27 | 1 | 2 | 3 | 6 |
|---|
| 28 | 1 | 1 | 3 | 3 |
|---|
| 29 | 1 | 2 | 3 | 6 |
|---|
| 30 | 1 | 1 | 2 | 2 |
|---|
| 31 | 1 | 1 | 2 | 2 |
|---|
| 32 | 1 | 2 | 2 | 4 |
|---|
| 33 | 1 | 1 | 2 | 2 |
|---|
| 34 | 1 | 2 | 2 | 4 |
|---|
| 35 | 1 | 1 | 1 | 1 |
|---|
| 36 | 1 | 2 | 1 | 2 |
|---|
| 37 | 1 | 1 | 2 | 2 |
|---|
| 38 | 1 | 2 | 2 | 4 |
|---|
| 39 | 1 | 1 | 3 | 3 |
|---|
| 40 | 1 | 2 | 3 | 6 |
|---|
| 41 | 1 | 1 | 3 | 3 |
|---|
| 42 | 1 | 2 | 3 | 6 |
|---|
| 43 | 1 | 1 | 1 | 1 |
|---|
| 44 | 1 | 2 | 1 | 2 |
|---|
| 45 | 1 | 1 | 1 | 1 |
|---|
| 46 | 1 | 2 | 1 | 2 |
|---|
| 47 | 1 | 1 | 1 | 1 |
|---|
| 48 | 1 | 2 | 1 | 2 |
|---|
| 49 | 1 | 1 | 1 | 1 |
|---|
| 50 | 1 | 2 | 1 | 2 |
|---|
| 51 | 1 | 1 | 2 | 2 |
|---|
| 52 | 1 | 2 | 2 | 4 |
|---|
| 53 | 1 | 1 | 1 | 1 |
|---|
| 54 | 1 | 2 | 1 | 2 |
|---|
| 55 | 1 | 1 | 2 | 2 |
|---|
| 56 | 1 | 2 | 2 | 4 |
|---|
| 57 | 1 | 1 | 3 | 3 |
|---|
| 58 | 1 | 2 | 3 | 6 |
|---|
| 59 | 1 | 1 | 2 | 2 |
|---|
| 60 | 1 | 2 | 2 | 4 |
|---|
| 61 | 1 | 1 | 2 | 2 |
|---|
| 62 | 1 | 1 | 2 | 2 |
|---|
| 63 | 1 | 2 | 2 | 4 |
|---|
| 64 | 1 | 1 | 3 | 3 |
|---|
| 65 | 1 | 2 | 3 | 6 |
|---|
| 66 | 1 | 1 | 3 | 3 |
|---|
| 67 | 1 | 2 | 3 | 6 |
|---|
| 68 | 1 | 1 | 1 | 1 |
|---|
| 69 | 1 | 2 | 1 | 2 |
|---|
| 70 | 1 | 1 | 3 | 3 |
|---|
| 71 | 1 | 1 | 2 | 2 |
|---|
| 72 | 1 | 2 | 2 | 4 |
|---|
| 73 | 1 | 1 | 3 | 3 |
|---|
| 74 | 1 | 2 | 3 | 6 |
|---|
| 75 | 1 | 1 | 2 | 2 |
|---|
| 76 | 1 | 2 | 2 | 4 |
|---|
| 77 | 1 | 1 | 1 | 1 |
|---|
| 78 | 1 | 2 | 1 | 2 |
|---|
| 79 | 1 | 1 | 3 | 3 |
|---|
| 80 | 1 | 2 | 3 | 6 |
|---|
| 81 | 1 | 1 | 1 | 1 |
|---|
| 82 | 1 | 2 | 1 | 2 |
|---|
| 83 | 1 | 1 | 1 | 1 |
|---|
| 84 | 1 | 2 | 1 | 2 |
|---|
| 85 | 1 | 1 | 1 | 1 |
|---|
| 86 | 1 | 2 | 1 | 2 |
|---|
| 87 | 1 | 1 | 3 | 3 |
|---|
| 88 | 1 | 2 | 3 | 6 |
|---|
| 89 | 1 | 1 | 3 | 3 |
|---|
| 90 | 1 | 1 | 2 | 2 |
|---|
| 91 | 1 | 2 | 2 | 4 |
|---|
| 92 | 1 | 1 | 3 | 3 |
|---|
| 93 | 1 | 2 | 3 | 6 |
|---|
| 94 | 1 | 1 | 3 | 3 |
|---|
| 95 | 1 | 1 | 3 | 3 |
|---|
| 96 | 1 | 2 | 3 | 6 |
|---|
| 97 | 1 | 1 | 1 | 1 |
|---|
| 98 | 1 | 2 | 1 | 2 |
|---|
| 99 | 1 | 1 | 3 | 3 |
|---|
| 100 | 1 | 2 | 3 | 6 |
|---|
| ⋮ |
|---|
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);
Linear mixed-effects model fit by REML
Model information:
Number of observations 286
Fixed effects coefficients 6
Random effects coefficients 302
Covariance parameters 4
Formula:
bpi_intensity ~ 1 + time*group + (1 + time | id)
Model fit statistics:
AIC BIC LogLikelihood Deviance
952.38 988.73 -466.19 932.38
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } 6.598 0.31156 21.177 280 4.434e-60 5.9847 7.2113
{'time' } -2.628 0.18597 -14.131 280 1.3627e-34 -2.9941 -2.2619
{'group_2' } -2.0483 0.4391 -4.6647 280 4.7925e-06 -2.9127 -1.1839
{'group_3' } -2.2628 0.4385 -5.1603 280 4.6771e-07 -3.126 -1.3996
{'time:group_2'} 1.7695 0.2628 6.7331 280 9.3944e-11 1.2521 2.2868
{'time:group_3'} 2.0878 0.25945 8.047 280 2.4295e-14 1.5771 2.5986
Random effects covariance parameters (95% CIs):
Group: id (151 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std' } 0.94359 3.3893e-16 2.627e+15
{'time' } {'(Intercept)'} {'corr'} 0.14217 -1 1
{'time' } {'time' } {'std' } 0.19832 2.7884e-141 1.4105e+139
Group: Error
Name Estimate Lower Upper
{'Res Std'} 0.87048 0.00020409 3712.8
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')
ans =
ANOVA marginal tests: DFMethod = 'Satterthwaite'
Term FStat DF1 DF2 pValue
{'(Intercept)'} 448.47 1 145.69 2.5654e-46
{'time' } 199.69 1 135.92 1.8529e-28
{'group' } 16.168 2 144.82 4.5931e-07
{'time:group' } 37.109 2 135.44 1.4098e-13
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
for it. This implicitly subtracts off a fraction of the pre-treatment score (the standardized 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);
Linear mixed-effects model fit by ML
Model information:
Number of observations 135
Fixed effects coefficients 4
Random effects coefficients 0
Covariance parameters 1
Formula:
BPI_Time2 ~ 1 + group + BPI_Time1
Model fit statistics:
AIC BIC LogLikelihood Deviance
423.39 437.92 -206.7 413.39
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)'} -1.0502 0.3473 -3.0238 131 0.0030043 -1.7372 -0.36312
{'group_2' } 1.6587 0.23937 6.9292 131 1.7252e-10 1.1851 2.1322
{'group_3' } 2.0183 0.23594 8.554 131 2.7305e-14 1.5515 2.485
{'BPI_Time1' } 0.60256 0.074734 8.0627 131 4.1045e-13 0.45472 0.7504
Random effects covariance parameters (95% CIs):
Group: Error
Name Estimate Lower Upper
{'Res Std'} 1.1187 0.99288 1.2604
anova(lme, 'DFMethod', 'Satterthwaite')
ans =
ANOVA marginal tests: DFMethod = 'Satterthwaite'
Term FStat DF1 DF2 pValue
{'(Intercept)'} 9.1433 1 135 0.0029886
{'group' } 40.98 2 135 1.2354e-14
{'BPI_Time1' } 65.007 1 135 3.5533e-13
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.