Plot marginals for interaction effects

[AC]

This should probably have a very simple answer, but I am figure it out. Apologies in advance if this is too foolish. I am trying to plot the marginal densities for the net effects from interactive terms. Consider the following:

library(Zelig); data(turnout) #Using a simple binary outcome data from the Zelig package
formula = vote ~ age + race * educate
mod = inla(formula, family="binomial", data=turnout)
summary(mod)

This produces the output:

Fixed effects:
                          mean          sd  0.025quant     0.5quant  0.975quant          kld
(Intercept)       -3.073100135 0.455472486 -3.98322504 -3.067436705 -2.19556728 0.0005631726
age                0.027582979 0.003511253  0.02075451  0.027561629  0.03453359 0.0002608525
racewhite          0.398661635 0.453717767 -0.48119542  0.394968744  1.29955810 0.0001240864
educate            0.225105180 0.035942497  0.15634452  0.224471158  0.29744299 0.0003300071
racewhite:educate -0.001955261 0.040093250 -0.08182132 -0.001537168  0.07549265 0.0001339810

Now, I can plot the approximate marginal for "educate" with:

plot(mod$marginals.fixed$'educate')

But how do I plot "educate + racewhite:educate"? Obviously, if I do "plot(mod$marginals.fixed$'racewhite:educate')", it's just going to plot the density for the additive effect, not the total effect.

[Håvard Rue]

On Sat, 2012-09-15 at 23:36 -0700, AC wrote:

>
> But how do I plot "educate + racewhite:educate"? Obviously, if I do
> "plot(mod$marginals.fixed$'racewhite:educate')", it's just going to
> plot the density for the additive effect, not the total effect.

this is the linear predictor. by default, its not computed, but adding

inla(..., control.predictor = list(compute=TRUE))

adds it. then its there as result$summary.linear.predictor and result
$marginals.linear.predictor



--
Håvard Rue
Department of Mathematical Sciences
Norwegian University of Science and Technology
N-7491 Trondheim, Norway
Voice: +47-7359-3533 URL : http://www.math.ntnu.no/~hrue
Mobile: +47-9260-0021 Email: havar...@math.ntnu.no

R-INLA: www.r-inla.org

[AC]

Thanks for the reply, Håvard. But I'm afraid I'm still lost. So the "control.predictor=list(compute=TRUE)" option yields the marginal linear predictors and summaries for each of the 2000 data points. But how do I calculate/plot the densities for the net effect of "educate" and "racewhite:educate" (i.e. the mean for which would be obviously 0.225 - 0.001 = 0.224)

Thanks, AC

[INLA help]

On Mon, 2012-09-17 at 03:39 -0700, AC wrote:
> Thanks for the reply, Håvard. But I'm afraid I'm still lost. So the
> "control.predictor=list(compute=TRUE)" option yields the marginal
> linear predictors and summaries for each of the 2000 data points. But
> how do I calculate/plot the densities for the net effect of "educate"
> and "racewhite:educate" (i.e. the mean for which would be obviously
> 0.225 - 0.001 = 0.224)

aaah, got it. then its easiest to use 'linear combinations', see the FAQ
for details. here is an artificial example

n = 100
racewhite = 1:n
educate = 1:n
y = rnorm(n)
r = inla(y ~ 1 + racewhite*educate,
data = data.frame(racewhite, educate),
control.predictor = list(compute=TRUE))

## net effect of educate + racewhite:educate
lc = inla.make.lincombs(
educate=educate,
"racewhite:educate" = racewhite * educate)

rr = inla(y ~ 1 + racewhite*educate,
data = data.frame(racewhite, educate),
control.predictor = list(compute=TRUE),
lincomb = lc)

rr$summary.lincomb.derived
rr$marginals.lincomb.derived


--
INLA help <he...@r-inla.org>
R-INLA

[INLA help]

On Mon, 2012-09-17 at 13:44 +0200, INLA help wrote:

## to be pedantic, or in the case where there are factors or contrasts
## involved, here is a revised example. the issue is that the
## model.matrix() sets up the fixed effects, so we need to define the
## linear combinations exactly as it is done there. so we need to run
## it twice, really. (well, that is the easiest option).

n = 100
racewhite = 1:n
educate = 1:n
y = rnorm(n)
r = inla(y ~ 1 + racewhite*educate,
data = data.frame(racewhite, educate),
control.predictor = list(compute=TRUE))

## net effect of educate + racewhite:educate
lc = inla.make.lincombs(

educate= r$model.matrix[, "educate"],
"racewhite:educate" = r$model.matrix[, "racewhite:educate"])

rr = inla(y ~ 1 + racewhite*educate,
data = data.frame(racewhite, educate),
control.predictor = list(compute=TRUE),
lincomb = lc)

[AC]

Could sampling from the posterior densities be an easier alternative (especially for for plotting)?

educate = inla.rmarginal(10000, rr$marginals.fixed$'educate')
educate_racewhite = inla.rmarginal(10000, rr$marginals.fixed$'educate') + inla.rmarginal(10000, rr$marginals.fixed$'racewhite:educate')

plot(density(educate), lty=1)
lines(density(educate_racewhite), lty=2)

[INLA help]

On Wed, 2012-09-19 at 01:51 -0700, AC wrote:
> Could sampling from the posterior densities be an easier alternative
> (especially for for plotting)?
>
> educate = inla.rmarginal(10000, rr$marginals.fixed$'educate')
> educate_racewhite = inla.rmarginal(10000, rr$marginals.fixed
> $'educate') + inla.rmarginal(10000, rr$marginals.fixed
> $'racewhite:educate')
>

don't know, the only thing is that you'll just get a noisy version of
the same marginlal.

if you want a function of two, say, as you do in the above, the sampling
is easier yes. for the above example, then you assume/use that the two
fixed effects are independent. if you want to compute the density of
their sum, you can use lincomb=.., or if you just want their posterior
correlation, you can get that using

control.fixed list(correlation.matrix = TRUE)

see ?control.fixed

[AC]

Got it. This was very helpful. I think I am beginning to get the hang of it now. Many thanks.

AC