Run a repeated measures anova as a mixed effects model using fitlme

Hi, I am using a repeated measures anova to run a model in which I am looking at the effects of 3 factors, 2 levels each (instruction, picture, time) on the DV which is activation with 22 subjects. Instruction is either "look" or "negative", picture is either "alcohol" or "food", and time is either "pre" or post"

So I start with a table with notations described as above:

t =

22×9 table

subject look_alc_pre look_alc_post negative_alc_pre negative_alc_post look_food_pre look_food_post negative_food_pre negative_food_post

_______ ____________ _____________ ________________ _________________ _____________ ______________ _________________ __________________

1 12.133 -0.58848 1.7663 -1.1068 1.1367 -12.784 -7.5119 -7.8238

2 11.246 -1.3831 11.598 -11.673 6.3135 -10.428 -0.54937 2.0609

3 -3.5529 8.2517 -12.154 18.181 -12.304 5.343 -12.638 15.727

4 -29.929 -17.964 -18.838 -23.523 -27.619 -20.193 -25.449 -18.107

5 -8.6783 5.1554 -10.748 4.1162 -12.658 -0.76942 -7.6787 11.932

6 42.469 60.953 32.816 71.571 14.558 45.931 20.959 38.89

7 11.843 5.0678 6.8864 9.7936 -1.2638 -3.7865 -3.7561 6.1002

8 42.138 -1.6591 33.319 2.8538 20.252 3.1166 21.965 6.3843

9 21.812 7.4502 16.274 -4.9925 23.128 -11.01 12.453 -3.9199

10 -0.46039 12.549 12.595 8.6604 8.4653 -6.2358 -1.1716 21.91

11 32.655 -6.289 23.065 -6.1028 16.403 -4.4612 15.165 -13.452

12 14.303 -8.3735 25.575 -0.040629 16.255 4.1393 19.625 9.6402

13 -5.7241 7.7997 6.1318 14.951 -7.9766 -3.5728 -14.098 4.4836

14 16.11 -4.0946 13.814 0.54614 0.24486 -0.64993 6.3281 -5.7614

15 -13.502 -4.5322 -4.9492 -7.0545 -4.4366 -6.1812 -17.884 -6.3731

16 10.407 -2.5284 -3.9457 7.8309 -0.78716 -6.3991 2.2326 -0.62589

17 13.739 3.0792 5.6006 -2.0368 6.5132 1.4708 3.9475 3.1808

18 15.826 -4.8494 7.9326 -16.506 2.8565 3.8328 7.6348 -4.6485

19 -8.8854 -0.71365 2.9053 -6.4694 2.485 3.3196 0.74676 4.0554

20 0.8475 -2.5931 -1.7847 -6.2106 3.6092 0.33887 -0.42818 -2.1065

21 2.3565 9.5125 -8.8134 3.2384 8.0217 8.9698 17.694 5.6276

22 -6.819 7.9659 -19.151 8.3581 -4.2582 15.056 -11.426 -6.0232

Then the following code successfully runs the model:

     WithinStructure=table([1 1 2 2 1 1 2 2]',[1 1 1 1 2 2 2 2]',[1 2 1 2 1 2 1 2]','VariableNames',{'instruction','picture','time'});
    WithinStructure.instruction = categorical(WithinStructure.instruction);
    WithinStructure.picture = categorical(WithinStructure.picture);
    WithinStructure.time = categorical(WithinStructure.time);
    rm=fitrm(t,'look_alc_pre-negative_food_post~1','WithinDesign',WithinStructure);
    ranovatbl=ranova(rm,'WithinModel','instruction*picture*time');

Howevever I will need to re-run this as a mixed effects model (including fitlme) as I will be including more sujects in the "pre" time. But I am confused as to how to set up the model (just using the original data above):

lmeModel = fitlme(t,'look_alc_pre-negative_food_post ~ 1+instruction*picture*time+(1|subject)');

But I keep getting the following error

Error using classreg.regr.LinearMixedFormula>parseStr
Unable to understand the character vector or string scalar 'look_alc_pre-negative_food_post ~ 1+instruction*picture*time+(1|subject)'.

Error in classreg.regr.LinearMixedFormula (line 166)
                [f.str,f.tree] = parseStr(modelSpec);

Error in LinearMixedModel.fit (line 2275)
        model.Formula = classreg.regr.LinearMixedFormula(formula,ds.Properties.VariableNames);

Error in fitlme (line 233)
    lme = LinearMixedModel.fit(ds,formula,varargin{:});

Not sure what I am doing incorrectly! I have attached the data. THank you!!

[the cyclist]

It would be helpful to everyone here if you provided complete code that creates the table(s) you used in the two analyses (or at least for the LME).

Although each of us can probably correctly infer some of what you did (e.g. look_alc_pre --> look_alcohol_pre_stack), there are still some puzzles to me, such as what

look_alcohol_pre_stack(:,a)

means. What is "a"? (Maybe I missed it?)

Please make it as easy as possible for us to help you.

[bsriv]

Hi apologies I just added the table directly which has all the relevant information

[the cyclist]

Can you upload the table in a MAT file instead? (Otherwise, we would need to cut & paste, define all the columns, etc.) You can use the paper clip icon in the INSERT section of the toolbar.

[bsriv]

Just uploaded. Thank you!

[bsriv]

Hi any follow up on this? Thanks!

Answer 1

[the cyclist]

Here is a fairly massive re-write of your analysis. Almost all of this code is restructuring the data into the tidy format. I used some fairly sophisticated table juju here. But the bottom line is to convert the information currently encoded in column names into variables. I display the top few rows of the restructured table, which makes it obvious what I did.

The model fit is trivial, once you format the data properly.

% Load the data. (I re-saved this MAT file with the variable named "tbl", because "table" is a keyword, and it was causing a conflict.)
load("exampledata.mat","tbl")

% Stack the variables into the "tidy" format, so that variable values are not encoded into column names.
tbl = stack(tbl,["look_alc_pre","look_alc_post","negative_alc_pre","negative_alc_post","look_food_pre","look_food_post","negative_food_pre","negative_food_post"],"IndexVariableName","index","NewDataVariableName","value");

% Use the dummy encoding that was created for the repeated measures analysis,
% and replicate it for each subject (but make them categorical)
WithinStructure=table(categorical([1 1 2 2 1 1 2 2])',categorical([1 1 1 1 2 2 2 2])',categorical([1 2 1 2 1 2 1 2])','VariableNames',{'instruction','picture','time'});
WithinStructure = repmat(WithinStructure,[22,1]);

% Concatenate the table with the dummy encoding
tbl = [tbl,WithinStructure];

% Removing the index variable, because that information is encoded in the dummy encoding
tbl = removevars(tbl,"index");

% Display what the first few rows of the table look like now
head(tbl)

    subject     value      instruction    picture    time
    _______    ________    ___________    _______    ____

       1         12.133         1            1        1  
       1       -0.58848         1            1        2  
       1         1.7663         2            1        1  
       1        -1.1068         2            1        2  
       1         1.1367         1            2        1  
       1        -12.784         1            2        2  
       1        -7.5119         2            2        1  
       1        -7.8238         2            2        2  

% Fit a model -- It might not be the one you want
lmeModel = fitlme(tbl,'value ~ 1 + instruction*picture*time + (1|subject)')

lmeModel = 
Linear mixed-effects model fit by ML

Model information:
    Number of observations             176
    Fixed effects coefficients           8
    Random effects coefficients         22
    Covariance parameters                2

Formula:
    value ~ 1 + instruction*picture + instruction*time + picture*time + instruction:picture:time + (1 | subject)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    1375.7    1407.4    -677.84          1355.7  

Fixed effects coefficients (95% CIs):
    Name                                      Estimate    SE        tStat       DF     pValue      Lower      Upper  
    {'(Intercept)'                   }         7.7425     3.0555      2.5339    168    0.012194     1.7103     13.775
    {'instruction_2'                 }        -2.2927     2.9754    -0.77054    168     0.44206    -8.1667     3.5814
    {'picture_2'                     }        -5.0634     2.9754     -1.7017    168    0.090653    -10.937    0.81063
    {'time_2'                        }        -4.4599     2.9754     -1.4989    168     0.13577    -10.334     1.4141
    {'instruction_2:picture_2'       }        0.80274     4.2079     0.19077    168     0.84894    -7.5044     9.1099
    {'instruction_2:time_2'          }         1.9367     4.2079     0.46026    168     0.64593    -6.3705     10.244
    {'picture_2:time_2'              }         2.0103     4.2079     0.47774    168     0.63346    -6.2969     10.317
    {'instruction_2:picture_2:time_2'}         2.1034     5.9509     0.35346    168     0.72418    -9.6447     13.852

Random effects covariance parameters (95% CIs):
Group: subject (22 Levels)
    Name1                  Name2                  Type           Estimate    Lower     Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        10.393      7.4791    14.442

Group: Error
    Name               Estimate    Lower     Upper 
    {'Res Std'}        9.8684      8.8256    11.034

[bsriv]

Thank you!!!!!

Answer 2

[the cyclist]

Please see my comment about helping us help you.

That being said, I don't think you can specify the response variable as

look_alc_pre-negative_food_post

because that is not valid Wilkinson notation. I'm surprised that that worked in fitrm. I'm frankly not sure what is going on there. Maybe there is an interpretation of that that I don't understand.

If you intended to specify a difference (with a minus sign), that won't work. You'd need to define a new variable for the difference. (I have not tried to figure out if that is just a typo.)

[bsriv]

Hi it looks like this is standard for fitrm for ModelSpec https://www.mathworks.com/help/stats/fitrm.html#bt9ic32-modelspec

But for fitlme you're right the notation is different, but i am still unsure

[the cyclist]

Oh, that's interesting! I've never used fitrm, and made assumptions. Sorry for the misinformation there.