How to calculate StdErr for intercepts and slopes in a multivariate ANCOVA (function ANOVAN() with continuous independent variable)?

I am trying to do multivariate analysis of covariance (mANCOVA) with Matlab. There is a special GUI for it called aoctool(), but it works only for one categorial independent variable (CIV) and one covariate. So if I have two CIVs and/or two covariates, I need something else.

Fortunately, anovan() tool, initially made for n-way ANOVA, could deal with continuos independent variables. So it is possible to perform mANCOVA using this tool and it really works! One can make and tests hypothesis on variance of slopes and treatments, but to provide a comperehensive analysis one needs to compare individual models (slopes or intercepts) to find differenes between them. But stats structure in anovan() outputs contains only values of model coefficients (stats.coeffs) but not their StdErr, SSE etc. So one can not test hypothesis on equality of two slopes (using t-test). Strictly speaking, there are no multcompare() for slopes in mANCOVA.

Can somebody give me an advice (or just any help), how I can borrow data on slopes StdErr or SSE from anovan() with continuous independent variable?

Thanks, Alex.

Answer 1

[Jeff Miller]

If you set up dummy variables to code your categorical independent variables, you can then do analysis of covariance with fitlm. The advance of that is that you can then use regression techniques to get standard errors and confidence intervals for the model coefficients. Maybe there is an easier way, though.

[Alex Sabrekov]

Thanks, Jeff! I tried and it works! Using both anovan() and fitlm() gives possibility to make comperehensive mANCOVA. But there are some issues (probably somebody will be interested):

• fitlm() with categorial variables as a IV by default uses one value of this IV as a reference. So you kind of lose this value from analysis. It could be fixed with dummy variables, and in this case very long model formulations should be made:

mdl1 = fitlm(data1,'y ~ -1 + OO:a +OO:b + OO:c + RO:a + RO:b + RO:c + WP:a + WP:b + WP:c + BW:a + BW:b + BW:c + x:OO + x:RO + x:WP + x:BW')

• when you specify full model including intercept and categorical variables without interaction (in my case OO, RO, a, b and so on) Matlab is complaining on rank deficiency. So only restricted model above with interaction but without intercept and single categorical variables did not report rank deficiency. Results are still ok (I compare with SAS calculation from Milliken and Johnson) but I am nervous about this reports and their possible influence on future results.

But anyway now I can do pairwise comparisons for slopes. Number of degrees of freedom remains unclear, looks like it should not be sizes of subsample for this particular slopes, but total sample size minus parameters - in my case 72 vs 72 and not 18 vs 18. Multiple comparison of slopes is still problematic because it needs customized Bonferroni-type correction but simple correction using Holm method could be made.

[Jeff Miller]

You are welcome, Alex.

Just a comment about rank deficiency that may be relevant: There are various ways to set up dummy variables for coding categorical variables (and there interactions), and some of these produce rank-deficient design matrices even if you leave off the constant term (as you have done with your -1 term in the model)--especially if you have multiple categorical variables. If you are interested, a coding method that avoids rank deficiency problems is described in this online book (with lots of examples).

[Alex Sabrekov]

Thank you very much again, Jeff!

I am wondering, how Matlab would deal with other types of introducing dummy variables. In Matlab in fact this is system of 0 and 1, not +1/0/-1 notation. For me it is more convenient to stay inside Matlab framework. Do you suggest optimal coding method in Chapter 18?

As I see, Milliken uses n-way ANOVA with 3 type of variance to test the overall structure of the model, and then uses DVR modelling to test slopes and estimate models. So possible rank deficiency is not critical for main decisions at least.

[Jeff Miller]

Hi Alex, No, I don't recommend that you change coding methods if you have one that gives you what you need. I just thought the other coding method might be helpful if rank deficiency was blocking you.

[Alex Sabrekov]

As I understand, this rank deficiency is a universal and obligate property of certain linear model formulations, and it doesn't really influence on slope estimates.

Thanks again for the idea of dummy variables! I have heard with half an ear about it, but thought that it is too complicated for me. But now I realize that this is the only way to solve these problems properly. Using fitlm() also gives one an ability to estimate model with a certain confidence interval, that is necesarry for result's reporting. So there is no need to search a new tricky software.