8
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ @Thislstheld What code are you using for fitting your logistic regression? $\endgroup$ – chl Jun 7 '11 at 12:14
  • $\begingroup$ @chl - Either glmfit with a 'link', 'logit' or mnrfit would work, no particular preference for either. $\endgroup$ – mpacer Jun 8 '11 at 0:04
  • $\begingroup$ I don't think this is a statistic question, but a programming question, which is placed more properly at StackOverflow... $\endgroup$ – Manoel Galdino Jun 22 '11 at 23:50
  • $\begingroup$ @Manoel as the reply by @chl shows, this question is fundamentally statistical in nature. If questions like this were off-topic, then so would be about half of the questions on this site! $\endgroup$ – whuber Jun 23 '11 at 1:57
9
+50
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Btw, for others who wish to use this, for your reference you need to divide by -2 not 2, otherwise it will return an error. See en.wikipedia.org/wiki/Deviance_%28statistics%29 $\endgroup$ – mpacer Sep 6 '11 at 12:18
  • $\begingroup$ Also, many thanks - this was incredibly helpful both in my understanding of the problem and the general class of problems that I was dealing with :) . $\endgroup$ – mpacer Sep 6 '11 at 12:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.