3
$\begingroup$

In my study, children repeat sentences and I count the errors (DV = error rate). Children are divided into two groups (Group factor), and there are two different types of sentences (Condition factor). I am particularly interested in the Group by Condition interaction, which is theoretically important in my field (developmental psycholinguistics).

When I run a by-subjects ANOVA, the "gold standard" in my field, I get a significant Group by Condition interaction. The assumptions of the ANOVA pretty much hold (normal distribution of residuals according to a QQ plot, homogeneity of variance (Levene's), Sphericity (Mauchly's)). The data themselves are not normally distributed, but I gather it's the distribution of the residuals that really matters.

I'm also interested in applying GLM approaches. The raw data (i.e. both groups combined not averaged in any way) are strongly right-ward skewed and are well-approximated by a negative binomial model (I've used STATA nbvargr routine). When a negative binomial regression is run the residuals show a normalish distribution according to a QQ plot (better than the poisson and gaussian models) though not nearly as good as the ANOVA. I also get a significant interaction but it is in the OPPOSITE direction to that of the ANOVA (a negative as opposed to positive coefficient). I think this is because the ANOVA assumes an additive model whereas the negative binomial assumes a multiplicative model (as it employs logs).

My question is, which analysis is "right", i.e. which best describes the underlying processes which result in the observed data? Though it has to be said that the ANOVA fit is excellent, there is a substantial loss of data as they are averaged for each cell of the ANOVA design (participant by group by condition = 50 observations per cell). From this perspective it's no wonder that the fit of the negative binomial model is not quite as good, as the model has to account for 50 times the data. In addition, the negative binomial model attempts to model the data generating process, i.e. a poisson process, whereas I don't think the same thing can be said for the ANOVA.

(NB the issue of how to interpret interactions, and whether to use an additive or multiplicative model is hot issue in Psychiatry (Kendler, K. S., & Gardner, C. O. (2010). Interpretation of interactions: guide for the perplexed. The British Journal of Psychiatry, 197(3), 170–171. doi:10.1192/bjp.bp.110.081331))

Here's some code showing how the interaction effect is negative for the negative binomial model, but positive for the ANOVA.

** ld = error rate, DV
** limp = group (0 = control group, 1 = group of interest)
** hard = experimental condition (0 = "easy" condition, 1 = "hard" condition)
** morphemes = length of sentence in morphemes, reqd as an offset for
***negative binomial

. nbreg ld limp##hard, robust cluster(id) offset(morphemes)

Fitting Poisson model:

[model-fitting procedure omitted for simplicity's sake]

Negative binomial regression                      Number of obs   =       3397
Dispersion           = mean                       Wald chi2(3)    =     116.39
Log pseudolikelihood = -7574.8427                 Prob > chi2     =     0.0000

                                (Std. Err. adjusted for 34 clusters in id_num)
------------------------------------------------------------------------------
             |               Robust
 ld_morph_br |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      1.limp |   2.235774   .2752258     8.12   0.000     1.696341    2.775206
      1.hard |   .7013264   .1065947     6.58   0.000     .4924047    .9102482
             |
   limp#hard |
        1 1  |  -.3988122    .123109    -3.24   0.001    -.6401013   -.1575231
             |
       _cons |  -10.25841   .2499824   -41.04   0.000    -10.74837   -9.768454
   morphemes |   (offset)
-------------+----------------------------------------------------------------
    /lnalpha |    .696452   .0583161                      .5821544    .8107495
-------------+----------------------------------------------------------------
       alpha |   2.006621   .1170184                       1.78989    2.249593
------------------------------------------------------------------------------

Now dataset is collapsed to obtain the mean error rate per participant for each cell of the ANOVA design

. collapse (mean) ld_morph_br, by(id_num limp hard)

. 
. anova ld_morph_br limp / id|limp hard limp#hard, repeated(hard)

                           Number of obs =      68     R-squared     =  0.9819
                           Root MSE      = .323407     Adj R-squared =  0.9621

                  Source |  Partial SS    df       MS           F     Prob > F
             ------------+----------------------------------------------------
                   Model |  181.662443    35  5.19035552      49.62     0.0000
                         |
                    limp |  126.729506     1  126.729506      83.17     0.0000
             id_num|limp |  48.7592906    32  1.52372783   
             ------------+----------------------------------------------------
                    hard |  5.62123672     1  5.62123672      53.74     0.0000
               limp#hard |  .552409685     1  .552409685       5.28     0.0282
                         |
                Residual |   3.3469382    32  .104591819   
             ------------+----------------------------------------------------
                   Total |  185.009381    67  2.76133405   


Between-subjects error term:  id_num|limp
                     Levels:  34        (32 df)
     Lowest b.s.e. variable:  id_num
     Covariance pooled over:  limp      (for repeated variable)

Repeated variable: hard
                                          Huynh-Feldt epsilon        =  1.0323
                                          *Huynh-Feldt epsilon reset to 1.0000
                                          Greenhouse-Geisser epsilon =  1.0000
                                          Box's conservative epsilon =  1.0000

                                            ------------ Prob > F ------------
                  Source |     df      F    Regular    H-F      G-G      Box
             ------------+----------------------------------------------------
                    hard |      1    53.74   0.0000   0.0000   0.0000   0.0000
               limp#hard |      1     5.28   0.0282   0.0282   0.0282   0.0282
                Residual |     32
             -----------------------------------------------------------------

Command "regress" used to obtain coefficients

. regress

      Source |       SS       df       MS              Number of obs =      68
-------------+------------------------------           F( 35,    32) =   49.62
       Model |  181.662443    35  5.19035552           Prob > F      =  0.0000
    Residual |   3.3469382    32  .104591819           R-squared     =  0.9819
-------------+------------------------------           Adj R-squared =  0.9621
       Total |  185.009381    67  2.76133405           Root MSE      =  .32341

------------------------------------------------------------------------------
 ld_morph_br |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      1.limp |   3.139737   .3327826     9.43   0.000     2.461881    3.817593
             |
 id_num#limp |
             | [omitted for simplicity's sake]
             |
      1.hard |   .3947685   .1109275     3.56   0.001     .1688164    .6207205
             |
   limp#hard |
        1 1  |   .3605257   .1568752     2.30   0.028     .0409813    .6800701
             |
       _cons |   .9226158   .2353129     3.92   0.000     .4432992    1.401932
------------------------------------------------------------------------------

Now we have a positive coefficient

Improve this question
$\endgroup$
3
  • $\begingroup$ Given what @Maarten says and what I understand from your question, it could be that the two methods are parameterizing the categorical variables differently - that is, which is the reference category for each independent variable? $\endgroup$
    – Peter Flom
    Commented May 12, 2014 at 10:11
  • $\begingroup$ Both analyses use the same reference category (0 = control group, 1 = population of interest) $\endgroup$
    – Nick Riches
    Commented May 12, 2014 at 12:07
  • $\begingroup$ Also check the categories for the IVs. $\endgroup$
    – Peter Flom
    Commented May 12, 2014 at 12:42

2 Answers 2

Reset to default
4
$\begingroup$

If you only have two categorical variables and you include the interaction between the two, then the nbreg and anova are equivalent in the sense that they produce exactly the same predicted values. As a consequence both models are "equally right" or "equally wrong".

sysuse nlsw88, clear

gen byte black = race == 2 if race < 3

anova wage black##union
predict yhat_anova

nbreg wage black##union
predict yhat_nbreg

tab yhat_anova yhat_nbreg

    Fitted |         Predicted number of events
    values |  5.968046     7.5821   8.614504   8.735929 |     Total
-----------+--------------------------------------------+----------
  5.968046 |       350          0          0          0 |       350 
    7.5821 |         0      1,051          0          0 |     1,051 
  8.614504 |         0          0        151          0 |       151 
  8.735929 |         0          0          0        302 |       302 
-----------+--------------------------------------------+----------
     Total |       350      1,051        151        302 |     1,854

This is not just true for a model with two categorical variables and their interaction; It is true for any fully saturated model, i.e. any model that inlcudes only categorical variables and all interactions including all higher order interactions.

Edit:

In order to see why these models are just different ways of saying the same thing it is easiest to go back at the basics. Consider the example above: Both models model the mean outcome, in this case mean wage. In both models there are only 4 means to model: black non-union, black union, white non-union white union. Both models model these 4 means with 4 parameters (constant, main effect for union, main effect for black, and the interaction term). As a consequence both models impose no constraint on these means.

Lets start with just looking at these means:

table  union black, c(mean wage)

------------------------------
union     |       black       
worker    |        0         1
----------+-------------------
 nonunion |   7.5821  5.968046
    union | 8.735929  8.614504
------------------------------

So, among non-union member there is a noticable difference in mean wage between black and white persons, while among union members this difference is a lot smaller.

Lets look at the regression output:

. reg wage i.union##i.black

      Source |       SS       df       MS              Number of obs =    1854
-------------+------------------------------           F(  3,  1850) =   29.83
       Model |  1472.83336     3  490.944452           Prob > F      =  0.0000
    Residual |  30445.6086  1850  16.4570858           R-squared     =  0.0461
-------------+------------------------------           Adj R-squared =  0.0446
       Total |   31918.442  1853   17.225279           Root MSE      =  4.0567

------------------------------------------------------------------------------
        wage |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval] 
-------------+----------------------------------------------------------------
       union |
      union  |   1.153829   .2648625     4.36   0.000     .6343677    1.673289
     1.black |  -1.614053   .2503572    -6.45   0.000    -2.105066   -1.123041
             |
 union#black |
    union#1  |   1.492629   .4755625     3.14   0.002     .5599329    2.425324
             |
       _cons |     7.5821   .1251339    60.59   0.000     7.336681    7.827518
------------------------------------------------------------------------------

You can see that the mean wage for non-union white persons is the constant. This mean is 1.61 dollars/hour less for black persons. So the wage for black persons is: 

di _b[_cons]+_b[1.black]
5.9680463

Which corresponds with the mean reported in the tabel above. The negative effect of being black is reduced by 1.48 dollars/hour, i.e. black union members earn

di _b[1.black] + _b[1.union#1.black]
-.12142489

dollars/hour less than white union members, a number you can also recover from the table of means above.

Consider a multiplicative model, in this case a Poisson with robust standard errors (for a justification of that choice, see here.

. poisson wage i.union##i.black, irr vce(robust)
note: you are responsible for interpretation of noncount dep. variable

Iteration 0:   log pseudolikelihood = -5239.2159  
Iteration 1:   log pseudolikelihood = -5239.2159  

Poisson regression                                Number of obs   =       1854
                                                  Wald chi2(3)    =     106.56
                                                  Prob > chi2     =     0.0000
Log pseudolikelihood = -5239.2159                 Pseudo R2       =     0.0188

------------------------------------------------------------------------------
             |               Robust
        wage |        IRR   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       union |
      union  |   1.152178   .0360531     4.53   0.000     1.083638    1.225053
     1.black |   .7871232   .0265286    -7.10   0.000     .7368083     .840874
             |
 union#black |
    union#1  |   1.252791   .0759936     3.72   0.000     1.112359    1.410951
             |
       _cons |     7.5821   .1308473   117.39   0.000     7.329932    7.842942
------------------------------------------------------------------------------

Again we can see that the mean wage of white non-union members is 7.5821. Black non-union persons earn $(0.787 - 1) \times 100\% = -21\%$ less then white non-union persons, a number you can recover from the table of means: $( 5.968046 - 7.5821 ) / 7.5821*100 = -21\%$. This negative effect of being black increases, which means becomes less negative, by 25% if someone is a union member. So, the effect of being black for union members is $1.25\times.787=.98$ or black union members earn 2\% less than white union members. A number you could recover from the table of means or the regression output, both of which I will leave as an exercise to the reader.

Improve this answer
$\endgroup$
9
  • $\begingroup$ Thanks for the demonstration. However, my data is slightly different in containing within-subjects factors and I'm also interested in the directionality of effects. I've taken the liberty of adding some code to my questions showing how the effects go in different directions. $\endgroup$
    – Nick Riches
    Commented May 12, 2014 at 12:40
  • 1
    $\begingroup$ While it's true that any number of categorical effects will yield the same predicted values for the outcome in either model, the inference on their corresponding effects will be completely different. This is what OP's original question addressed. $\endgroup$
    – AdamO
    Commented May 12, 2014 at 13:00
  • $\begingroup$ The key thing is that the coefficients mean different things, but you can just compute one set of coefficients from the other. So, the models just say the same thing in different ways. My take on this is that you get most out of these models by treating them as complementary rather then competing represenatations of your data. $\endgroup$
    – Maarten Buis
    Commented May 12, 2014 at 14:06
  • $\begingroup$ @MaartenBuis how exactly would one compute coefficients for one model from the other? Perhaps more importantly, would one be able to also calculate the standard errors? I wouldn't agree that the model coefficients "say the same thing", especially considering that these findings indicate an equal and opposite effect for an interaction parameter. I also definitely would not say they're complementary. All too often I review an analysis where it's clear the statistician has a lack of understanding of the problem, can't justify a single decisive model choice, and presents conflicting findings. $\endgroup$
    – AdamO
    Commented May 12, 2014 at 16:46
  • $\begingroup$ @AdamO I edited my answer to discuss how one can move from one model to the other. As for inference: that is what robust standard errors are for, see e.g. here. I would say that the problem of the papers you reviewed is not so much that the findings are conflicting, but that they haven't been successful at explaining to their readers why they aren't conflicting. $\endgroup$
    – Maarten Buis
    Commented May 13, 2014 at 8:11
2
$\begingroup$

In general, you can't predict what one model might find relative to another. While the effects are additive in the same predictors, each model coefficient should be interpreted completely differently. The correct choice of model, of course, depends upon the nature of the problem. Other model output should be viewed as exploratory and, in this case, spurious.

Negative Binomial Regression has a non-linear link function and a mean-variance relationship. That means that completely different parameters are being fit in the model. The NB GLM estimates a relative rate.

The inference on effects in an ANOVA is completely analogous to the linear regression model. You should be comparing your ANOVA output to that you obtain from regress in Stata, and you'd find that the Wald test p-values for model coefficients are very close the likelihood ratio tests presented in ANOVA. The OLS estimates a mean difference.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks. The main "problem" is whether children in one of the groups exhibit a qualitatively different profile. I guess both analyses point to this. $\endgroup$
    – Nick Riches
    Commented May 12, 2014 at 14:42
  • $\begingroup$ @NickRiches The distribution of the outcome variable is obviously of the main concern. $\endgroup$
    – AdamO
    Commented May 12, 2014 at 16:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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