How to plot an interaction term, using a model's coefficients, of three-factorial GEE model with full-order interactions (geepack package)?

Question

Using function geeglm from package geepack (Generalized Estimating Equation), I have modeled counts as being dependent on two nominal (factor) variables, one continuous variable with 3rd-order interactions and one grouping variable. This is the model:

m1<-geeglm(formula = dependent.var  ~  cat.var1 *  cat.var2 * contin.var,  
           family = poisson, id = group, corstr = "exchangeable")

Factor variable cat.var1 has two levels (CD and WL), factor variable cat.var2 has two levels (SAH and SKA). Group is a grouping variable to account for autocorrelation among somehow related subjects. Using anova function, I found that the only significant model’s interaction terms are cat.var1 : contin.var and cat.var2 : contin.var. I would like to plot these two interactions: I have used the codes below, but I am not happy with the fitted lines. Line for CD-level of cat.var1 in the first figure is the same as line for SAH-level of cat.var2 in the second figure (see figs) – therefore I think that they are not correct and I do not know how to solve this. Which model's coefficients I have to use to construct that lines? The model’s coefficients used for producing of the figures are thereinafter.

Codes:

# 1st plot: cat.var1 : contin.var
par(mfrow=c(1,1))
plot(contin.var, dependent.var, type="n")
points(contin.var[cat.var1=="WL"], jitter(dependent.var[cat.var1=="WL"]))
points(contin.var[cat.var1=="CD"], jitter(dependent.var[cat.var1=="CD"]),
       pch=16)
x <- seq(3,7,0.1)
y1 <-exp(-0.1584 + 0.3474*x)
lines(x,y1, lty=2)                                 # CD
y2 <-exp(-0.1584 + 3.7293 +(-0.6685 + 0.3474)*x)
lines(x,y2, lty=1)                                 # WL

# 2nd plot: cat.var2 : contin.var
par(mfrow=c(1,1))
plot(contin.var, dependent.var, type="n")
points(contin.var[cat.var2=="SKA"], jitter(dependent.var[cat.var2=="SKA"]))
points(contin.var[cat.var2=="SAH"], jitter(dependent.var[cat.var2=="SAH"]), 
       pch=16)
x <- seq(3,7,0.1)
y1 <-exp(-0.1584 + 0.3474*x)
lines(x,y1, lty=2)                                 # SAH
y2 <-exp(-0.1584 + -6.9490 + (1.2249 + 0.3474)*x)
lines(x,y2)                                        # SKA

Here are my model's coefficients, used for producing of figures, outputs from summary(m1) (statistically significant terms - outputs from anova(m1) - are denoted with asteriks **):

Coefficients:
                                         Estimate

(Intercept)                                  -0.1584   
cat.var1WL                                    3.7293   
cat.var2SKA                                  -6.9490  
contin.var                                    0.3474   *
cat.var1WL: cat.var2SKA                       3.9970   
cat.var1WL: contin.var                       -0.6685   *
cat.var2SKA: contin.var                       1.2249   **
cat.var1WL: cat.var2SKA: contin.var          -0.6860   

Factor variable *cat.var1* has two levels (CD and WL),
factor variable *cat.var2* has two levels (SAH and SKA).

Here are the data: https://docs.google.com/file/d/0Bz8ojhHeiNclVi1oT0ZwTEtEN2s/edit

Answer 1

First of all, it's great to see the data included with your question!

It sounds like the heart of the question is whether you correctly created curve formulas based on the coefficients from the geeglm package. I don't know the correct way for that model, but using only intercept and contin.var coefficients for those cases can't be correct. Though you only see a coefficient for cat.var1WL in the package output, there is an implied complementary coefficient for cat.var1CD.

In any case, you should be able examine m1.formula to get the full formula.

Given that your model includes categorical interaction, you should show the interaction in the graphs by plotting different curves for each combination (CD&SAH, CD&SKA, WL&SAH, WL&SKA), whether overlaid or in a grid (example below with a smoother).

An alternative is to remove that interaction from the model if you've determined if to be insignificant. Then your provided graphs are more appropriate, once you fix the y formulas.

Answer 2

As the main effects of cat.var1 and cat.var2 were not statistically significant and I was interested about plotting only significant interaction terms (see my question above), I made two separate models to obtain coefficients for construction of plots to show those interaction between explanatory variables, i.e. one model for each significant interaction from the m1. Both models contain two relevant variables and their interaction term. Second factor variable was thus something like "pooled".

# 1st plot for interaction cat.var1:contin.var:

m2<-geeglm(dependent.var ~ cat.var1 * contin.var, family = poisson, id = group, corstr = "exchangeable")

Coefficients:
                                         Estimate

Intercept                                  -3.464
cat.var1WL                                  5.264
contin.var                                  0.912
cat.var1WL: contin.var                     -0.917 

Plot code:

par(mfrow=c(1,1))
plot(contin.var, dependent.var, type="n")
points(contin.var[cat.var1=="WL"], jitter(dependent.var[cat.var1=="WL"]), pch=16)
points(contin.var[cat.var1=="CD"], jitter(dependent.var[cat.var1=="CD"]))
x <- seq(3,7,0.1)
y1 <-exp (-3.464 + 0.912*x)
lines(x,y1, lty=2)                                 # CD
y2 <-exp ((-3.464 + 5.264) + (-0.917 + 0.912)*x)
lines(x,y2, lty=1)                                 # WL
legend(4.82, 15, c("CD", "WL"), pch=c(16, 1), lty=c(2, 1))

Plot:

# 2nd plot for interaction cat.var2:contin.var:

m3<-geeglm(dependent.var ~ cat.var2 * contin.var, family = poisson, id = group, corstr = "exchangeable")

Coefficients:
                                         Estimate

Intercept                                  1.1384
cat.var2SKA                               -4.1958
contin.var                                 0.1182
cat.var2SKA: contin.var                    0.7532

Plot code:

par(mfrow=c(1,1))
plot(contin.var, dependent.var, type="n")
points(contin.var[cat.var2=="SKA"], jitter(dependent.var[cat.var2=="SKA"]))
points(contin.var[cat.var2=="SAH"], jitter(dependent.var[cat.var2=="SAH"]), 
       pch=16)
x <- seq(3,7,0.1)
y1 <-exp (1.1384 + 0.1182*x)
lines(x,y1, lty=2)                                 # SAH
y2 <-exp ((1.1384 + -4.1958) + (0.7532 + 0.1182)*x)
lines(x,y2, lty=1)                                 # SKA
legend(4.82, 15, c("SAH", "SKA"), pch=c(1, 16), lty=c(2, 1))

Plot: