estimate main effects and interactions

Hi,

I'm trying to replicate Minitab functionality by creating a Pareto Chart of standardised effects.

I'm using data from the example linked here where a DOE was run with 4 inputs and the response was measured. In the link, a full effects table is calculated which has effects of each of the inputs A,B,C & D but also all of the interaction terms.

In Matlab I'm trying to use stepwiselm to fit a model to the table of data but I find that Matlab is removing terms based on their P-value. is there a way to use this function but ask for all linear and interaction terms to be retained?

Secondly, I'm looking at the functions plotEffects() and plotInteraction() to use on the model but is there a function or way to group all main & interaction effects together in a single table, as per table 2 in the link above?

Many thanks, Dan

[dpb]

" I'm trying to use stepwiselm to fit a model to the table of data but I find that Matlab is removing terms based on their P-value. is there a way to use this function but ask for all linear and interaction terms to be retained?"

That's what stepwise does; why using it instead of fitlm if you want a specific model?

I've not investigated the specifics of the rest of the Q?

[Dan Howe]

Thanks for the quick reply,

you're right that if I'm just trying to replicate analysis in the link then I could probably use fitlm. If I use the 'interactions' option under fitlm modelspec I do get interaction effects between say A&B but in the example, terms such as A*B*C are included. Can fitlm include these terms without manually specifying the full formula?

[dpb]

Nope; doesn't appear that that's a named model option....fitlm seems to stop at two-level interactions without writing the model. (That said, I've not actually tried with higher numbers of independent variables to test; you've at least looked at what the 'interactions' value does it seems).

[Dan Howe]

I'll have a look at writing out the model or check out other functions that might have the option of including higher level interactions.

For the second part of my question, the link I mention includes a table of effects that include all individual and interaction terms e.g.:

Term Effect

A -8

B 24

.. ..

A*D 0.75

B*C 4.5

... ..

A*B*C 0.5

A*B*C*D -0.25

Can Matlab produce a similar table of effects once a model has been fit to the data? As in my original post, I can use either plotEffects or plotInteraction but I ideally want to combine them together.

[dpb]

You see the output from fitlm--for the specified model you get the coefficent estimates and statistics for each term in the model; the "effect" is just the magnitude of the coefficient. There's a whole bunch of other anciliary data available from the LinearModel object.

The precise format of a table in the identical form as above you'll undoubtedly have to create yourself from the information; MATLAB doesn't know how any particular user will want the output to look...

[Dan Howe]

Thanks a lot for the support, I didn't appreciate that the coefficients from the model are termed the same as main effects.

Answer 1

[the cyclist]

I was able to replicate the effect estimates, but only via a rather odd normalization scheme of the variables:

data = [1	10	220	10	50	8	70
        2	15	220	10	50	2	60
        3	10	240	10	50	10	89
        4	15	240	10	50	4	81
        5	10	220	12	50	16	60
        6	15	220	12	50	5	49
        7	10	240	12	50	11	88
        8	15	240	12	50	14	82
        9	10	220	10	80	15	69
        10	15	220	10	80	9	62
        11	10	240	10	80	1	88
        12	15	240	10	80	13	81
        13	10	220	12	80	3	60
        14	15	220	12	80	12	52
        15	10	240	12	80	6	86
        16	15	240	12	80	7	79];

A = data(:,2);
B = data(:,3);
C = data(:,4);
D = data(:,5);
Y = data(:,7);

A = normalize(A,'zscore','robust');
B = normalize(B,'zscore','robust');
C = normalize(C,'zscore','robust');
D = normalize(D,'zscore','robust');

Y = 2*Y;

tbl = table(A,B,C,D,Y);

mdl = fitlm(tbl,'Y ~ A + B + C + D + A:B + A:C + A:D + B:C + B:D + C:D + A:B:C + A:B:D + A:C:D + B:C:D + A:B:C:D')

The normalization of the explanatory variables is a standard one, but I have no idea the reason behind needing to multiply the response variable by 2.

A slightly simpler way to specify the model in this case would be

modelMatrix = fullfact([2 2 2 2])-1;
mdl = fitlm([A B C D],Y,modelMatrix)

That's more of a Design of Experiments approach, I guess, but it boils down to the same thing.

[Dan Howe]

I tried this (I couldnt use 'normalize' as I'm on Matlab 2016 but I manually normalised to 0's and 1's using this approach) and after multiplying Y by 2, have this table:

0	0	0	0	140
1	0	0	0	120
0	1	0	0	178
1	1	0	0	162
0	0	1	0	120
1	0	1	0	98
0	1	1	0	176
1	1	1	0	164
0	0	0	1	138
1	0	0	1	124
0	1	0	1	176
1	1	0	1	162
0	0	1	1	120
1	0	1	1	104
0	1	1	1	172
1	1	1	1	158

the output from fitlm with this table, using the full model terms you've declared, gives me the following coefficients

Estimated Coefficients:
                    Estimate      SE    tStat    pValue
                   ___________    __    _____    ______

    (Intercept)            140    0      Inf     NaN   
    A                      -20    0     -Inf     NaN   
    B                       38    0      Inf     NaN   
    C                      -20    0     -Inf     NaN   
    D                       -2    0     -Inf     NaN   
    A:B                      4    0      Inf     NaN   
    A:C                     -2    0     -Inf     NaN   
    B:C                     18    0      Inf     NaN   
    A:D                      6    0      Inf     NaN   
    B:D             1.6117e-13    0      Inf     NaN   
    C:D                      2    0      Inf     NaN   
    A:B:C                    6    0      Inf     NaN   
    A:B:D                   -4    0     -Inf     NaN   
    A:C:D          -7.7716e-14    0     -Inf     NaN   
    B:C:D                   -4    0     -Inf     NaN   
    A:B:C:D                 -4    0     -Inf     NaN   

Can you think of an obvious reason why my coefficients don't match up to the linked example? Thanks!

[the cyclist]

The normalization that worked has -1's and 1's for the explanatory variables. So you need

mine = 2*yours - 1;

[Dan Howe]

that worked perfectly, thank you!

[the cyclist]

Great! Post here if you find out why that normalization scheme makes any sense. :-)

[dpb]

I hadn't actually looked at the link before...makes no sense to me why would have done from what is published there, certainly. Maybe if one looked at the original references could deduce why.