Coding an interaction between a nominal and a continuous predictor for logistic regression in MATLAB

Question

So our data is structured as follows:

We have $M$ participants, each participant can be categorized into 3 groups (G $\in {A,B,C}$), and for each participant we have $N$ samples of a continuous variable. And we are trying to predict values that are either 0 or 1.

How would we use matlab to test for an interaction between the continuous variable and the categorical variable in predicting these values?

Answer 1

The easiest way, IMO, is to build the design matrix yourself, as glmfit accepts either a matrix of raw (observed) values or a design matrix. Coding an interaction term isn't that much difficult once you wrote the full model. Let's say we have two predictors, $x$ (continuous) and $g$ (categorical, with three unordered levels, say $g={1,2,3}$). Using Wilkinson's notation, we would write this model as y ~ x + g + x:g, neglecting the left-hand side (for a binomial outcome, we would use a logit link function). We only need two dummy vectors to code the g levels (as present/absent for a particular observation), so we will have 5 regression coefficients, plus an intercept term. This can be summarized as

$$\beta_0 + \beta_1\cdot x +\beta_2\cdot\mathbb{I}_{g=2} +\beta_3\cdot\mathbb{I}_{g=3} + \beta_4\cdot x\times\mathbb{I}_{g=2} + \beta_5\cdot x\times\mathbb{I}_{g=3},$$

where $\mathbb{I}$ stands for an indicator matrix coding the level of $g$.

In Matlab, using the online example, I would do as follows:

x = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]';
g = [1 1 1 1 2 2 2 2 3 3 3 3]';
gcat = dummyvar(g);
gcat = gcat(:,2:3); % remove the first column
X = [x gcat x.*gcat(:,1) x.*gcat(:,2)]; 
n = [48 42 31 34 31 21 23 23 21 16 17 21]';
y = [1 2 0 3 8 8 14 17 19 15 17 21]';
[b, dev, stats] = glmfit(X, [y n], 'binomial', 'link', 'probit');

I didn't include a column of ones for the intercept as it is included by default. The design matrix looks like

    2100           0           0           0           0
    2300           0           0           0           0
    2500           0           0           0           0
    2700           0           0           0           0
    2900           1           0        2900           0
    3100           1           0        3100           0
    3300           1           0        3300           0
    3500           1           0        3500           0
    3700           0           1           0        3700
    3900           0           1           0        3900
    4100           0           1           0        4100
    4300           0           1           0        4300

and you can see that the interaction terms are just coded as the product of x with the corresponding column of g (g=2 and g=3, since we don't need the first level).

The results are given below, as coefficients, standard errors, statistic and p-value (from stats structure):

   int.   -3.8929    2.0251   -1.9223    0.0546
   x       0.0009    0.0008    1.0663    0.2863
   g2     -3.2125    2.7622   -1.1630    0.2448
   g3     -5.7745    7.5542   -0.7644    0.4446
   x:g2    0.0013    0.0010    1.3122    0.1894
   x:g3    0.0021    0.0021    0.9882    0.3230

Now, testing the interaction can be done by computing the difference in deviance from the full model above and a reduced model (omitting the interaction term, that is the last two columns of the design matrix). This can be done manually, or using the lratiotest function which provides Likelihood ratio hypothesis test. The deviance for the full model is 4.3122 (dev), while for the model without interaction it is 6.4200 (I used glmfit(X(:,1:3), [y n], 'binomial', 'link', 'probit');), and the associated LR test has two degrees of freedom (the difference in the number of parameters between the two models). As the scaled deviance is just two times the log-likelihood for GLMs, we can use

[H, pValue, Ratio, CriticalValue] = lratiotest(4.3122/2, 6.4200/2, 2)

where the statistic is distributed as a $\chi^2$ with 2 df (the critical value is then 5.9915, seechi2inv(0.95, 2)). The output indicates a non-significant result: We cannot conclude to the existence of an interaction between x and g in the observed sample.

I guess you can wrap up the above steps in a convenient function of your choice. (Note that the LR test might be done by hand in very few commands!)

---

I checked those results against R output, which is given next.

Here is the R code:

x <- c(2100,2300,2500,2700,2900,3100,3300,3500,3700,3900,4100,4300)
g <- gl(3, 4)
n <- c(48,42,31,34,31,21,23,23,21,16,17,21)
y <- c(1,2,0,3,8,8,14,17,19,15,17,21)
f <- cbind(y, n-y) ~ x*g
model.matrix(f)  # will be model.frame() for glm()
m1 <- glm(f, family=binomial("probit"))
summary(m1)

Here are the results, for the coefficients in the full model,

Call:
glm(formula = f, family = binomial("probit"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7124  -0.1192   0.1494   0.3036   0.5585  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept) -3.892859   2.025096  -1.922   0.0546 .
x            0.000884   0.000829   1.066   0.2863  
g2          -3.212494   2.762155  -1.163   0.2448  
g3          -5.774400   7.553615  -0.764   0.4446  
x:g2         0.001335   0.001017   1.312   0.1894  
x:g3         0.002061   0.002086   0.988   0.3230  

For the comparison of the two nested models, I used the following commands:

m0 <- update(m1, . ~ . -x:g)
anova(m1,m0)

which yields the following "deviance table":

Analysis of Deviance Table

Model 1: cbind(y, n - y) ~ x + g
Model 2: cbind(y, n - y) ~ x * g
  Resid. Df Resid. Dev Df Deviance
1         8     6.4200            
2         6     4.3122  2   2.1078