R planned comparisons in Zelig negative binomial regression

Question

I have this model:

model <- zelig(dv~(product*intervention), model = "negbin", data = data)

intervention has two levels: neutral(=0), treatment(=1)
product has two levels: product1(=0), product2(=1)

I build f_all to just have one factor with 4 groups for comparison analysis.

Thus I have 4 groups in f_all
1. product1-neutral
2. product1-treatment
3. product2-neutral
4. product2-treament

My interaction hypothesis is that treatment only works for product2.

Zelig gives me my predicted significant interaction.

Yet, I need planned contrasts to test my specific hypothesis: c(-1,1,0,0) and c(0,0,1,-1)

I researched and found a description of doing this with multcomp on this page: post comparisons

The regression output shows my predicted interaction

(Intercept)  1.34223    0.08024  16.728   <2e-16 ***
product      0.08747    0.08025   1.090   0.2757
intervention 0.07437    0.07731   0.962   0.3361
interaction  0.45645    0.22263   2.050   0.0403 * 

However, it said multcomp and the glht function is for linear models, but I am using a negbin model.

3 Questions regarding this problem:
1. Can I do planned comparisons on my negbin model using multcomp?
2. If not what appropriate method is there to do this for my negbin model?
3. Based on R using treatment contrasts per default could I just interpret the interaction coefficient as the contrast comparing product2-neutral versus product2-treatment? Can I then interpret the intervention coefficient as contrast comparing product1-neutral versus product1-treament?

Answer 1

To answer your specific question 1, yes you can do planned comparisons with multcomp even though you are using a generalized linear model. From the package description:

Simultaneous tests and confidence intervals for general linear hypotheses in parametric models, including linear, generalized linear, linear mixed effects, and survival models.

You can easily implement this with the Zelig output (which is an object from the negbin class since Zelig calls the glm.nb function from the MASS package). Here is an example:

library(Zelig)
library(multcomp)
data(sanction)
z.out <- zelig(num ~ target  * coop, model = "negbin", data = sanction)

## construct contrast matrices
hypo.mat <- rbind("coop0:target1 - target0" = c(0, 1, 0, 0),
                  "coop1:target1 - target0" = c(0, 1, 0, 1))
summary(glht(z.out, hypo.mat))

Which gives the following output:

    Simultaneous Tests for General Linear Hypotheses

Fit: zelig(formula = num ~ target * coop, model = "negbin", data = sanction)

Linear Hypotheses:
                             Estimate Std. Error z value Pr(>|z|)
coop0:target1 - target0 == 0  0.04201    0.38908   0.108    0.971
coop1:target1 - target0 == 0  0.09089    0.24811   0.366    0.786
(Adjusted p values reported -- single-step method)

Note that I used different contrasts than you gave. You are putting the contrasts in terms of the vector of groups, but multicomp (and its general form of hypothesis testing) wants contrasts on the model parameters. We can write the model above as

$\log \mu_i = \beta_0 + \beta_1 x_i + \beta_2 z_i + \beta_3 (x_i \times z_i)$

where $E(Y) = \mu_i$ is the expected value of the outcome. Thus, in this model, the hypothesis that the effect of $x_i$ is zero when $z_i$ is 0 is just:

$H_0: \beta_1 = 0$

This leads to the contrast c(0,1,0,0). The hypothesis that the effect of $x_i$ is zero when $z_i$ is 0 is just:

$H_0: \beta_1 + \beta_3 = 0$

This leads to the contrast c(0,1,0,1).

Answer 2

To answer question 3 first, the contrast between product2-neutral and product2-treatment is the sum of the interaction and intervention coefficients, while the contrast between product1-neutral and product1-treatment is just the intervention coefficient. That is, the interaction coefficient is the change in the effect from product1 to product2.

I am not sure about multicomp, but the Zelig package allows for you to easily find what it calls first differences, which are the difference in the expected value of the dependent variable for a given change in an independent variable. Here is some sample code that will help with what you want.

library(Zelig)
data(sanction)
z.out <- zelig(num ~ target  * coop, model = "negbin", data = sanction)

Once you have this model object, you can define the vectors of independent variables you would like to contrast. For instance, say you wanted to find the effect of the target variable when coop is 1. You would simply create vectors using the setx command and then pass it to the sim function:

x.high <- setx(z.out, target = 1, coop = 1)
x.low <- setx(z.out, target = 0, coop = 1)
s.out <- sim(z.out, x = x.low, x1 = x.high)
summary(s.out)

Part of the output will look like this:

First Differences in Expected Values: E(Y|X1)-E(Y|X)
       mean    sd    2.5%  97.5%
1 -0.002548 0.265 -0.7818 0.2627

Thus, in this case, the effect of target when coop is 1 is -0.0025 but is insignificant. You can easily repeat this for when coop is 0 and compare the two "first differences." The sim function estimates various quantities of interest by simulating the parameters of the model (it is akin to a parametric bootstrap). The Zelig documentation has more information on the software and King, Tomz, and Wittenberg (2000) describes the methods behind it. One benefit of this approach is that the answer is on the scale of the dependent variable as opposed to on the scale of the linear predictor.