¿Cómo trazar un término de interacción, utilizando los coeficientes de un modelo, de un modelo GEE trifactorial con interacciones de orden completo (paquete geepack)?

Question

Uso de la función geeglm del paquete geepack (Ecuación de estimación generalizada), he modelado los recuentos como dependientes de dos variables nominales (factoriales), una variable continua con interacciones de tercer orden y una variable de agrupación. Este es el modelo:

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

Variable factorial cat.var1 tiene dos niveles (CD y WL), la variable factorial cat.var2 tiene dos niveles (SAH y SKA). Grupo es una variable de agrupación para tener en cuenta la autocorrelación entre sujetos relacionados de alguna manera. Utilizando anova encontré que los únicos términos de interacción significativos del modelo son cat.var1 : contin.var y cat.var2 : contin.var . Me gustaría trazar estas dos interacciones: He utilizado los códigos de abajo, pero no estoy contento con las líneas ajustadas. Línea para el nivel CD de cat.var1 en la primera figura es la misma que la línea para el nivel SAH de cat.var2 en la segunda figura (ver figuras) - por lo que creo que no son correctas y no sé cómo resolverlo. ¿Qué coeficientes del modelo tengo que utilizar para construir esas líneas? Los coeficientes del modelo utilizados para la elaboración de las cifras figuran a continuación.

Códigos:

# 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

Aquí están los coeficientes de mi modelo, utilizados para producir las cifras, salidas de summary(m1) (términos estadísticamente significativos - resultados de anova(m1) - se indican con asteriscos **):

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).

Aquí están los datos: https://docs.google.com/file/d/0Bz8ojhHeiNclVi1oT0ZwTEtEN2s/edit

Answer 1

En primer lugar, ¡es estupendo ver los datos incluidos en tu pregunta!

Parece que el meollo de la cuestión es si has creado correctamente las fórmulas de las curvas basadas en los coeficientes del paquete geeglm. No conozco la forma correcta para ese modelo, pero usando sólo intercept y contin.var los coeficientes para esos casos no pueden ser correctos. Aunque sólo se ve un coeficiente para cat.var1WL en la salida del paquete, hay un coeficiente complementario implícito para cat.var1CD .

En cualquier caso, debería poder examinar m1.formula para obtener la fórmula completa.

Dado que su modelo incluye una interacción categórica, debe mostrar la interacción en los gráficos trazando diferentes curvas para cada combinación (CD&SAH, CD&SKA, WL&SAH, WL&SKA), ya sea superpuestas o en una cuadrícula (ejemplo de abajo con un suavizador).

Una alternativa es eliminar esa interacción del modelo si se ha determinado que es insignificante. En ese caso, los gráficos que ha proporcionado son más apropiados, una vez que haya arreglado el y fórmulas.

Answer 2

Como los principales efectos de cat.var1 y cat.var2 no eran estadísticamente significativas y me interesaba trazar sólo los términos de interacción significativos (véase mi pregunta anterior), hice dos modelos separados para obtener los coeficientes para la construcción de gráficos que mostraran esas interacciones entre las variables explicativas, es decir, un modelo para cada interacción significativa del m1 . Ambos modelos contienen dos variables relevantes y su término de interacción. La segunda variable factorial era, por tanto, algo así como "agrupada".

# Primer gráfico para la interacción cat.var1:contin.var:

m2<-geeglm(var.dependiente ~ var.cat1 * var.contin, familia = poisson, id = grupo, corstr = "intercambiable")

Coefficients:
                                         Estimate

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

Código de la parcela:

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))

Trama:

# Segundo gráfico para la interacción cat.var2:contin.var:

m3<-geeglm(var.dependiente ~ var.cat2 * var.contin, familia = poisson, id = grupo, corstr = "intercambiable")

Coefficients:
                                         Estimate

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

Código de la parcela:

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))

Trama: