Using INLA for predicting a binary variable

[Samuel Molano]

Hello,

I'm trying to use INLA to predict a binary variable (called empty, that can be either 0 or 1) over a grid of points, using a spatial approach.

This is the code that I run:

formula <- empty ~ -1 + Intercept + f(field, model = spde)

mod <- inla(formula, data=inla.stack.data(stack, spde=spde),
            family="binomial",
            control.predictor=list(A=inla.stack.A(stack), compute=TRUE, link = 1))

index.pred = inla.stack.index(stack,"pred")$data
lp.mean.pred = mod$summary.linear.predictor[index.pred, "mean"]
lp.sd.pred = mod$summary.linear.predictor[index.pred, "sd"]

However, the lp.mean.pred vector takes negative values sometimes, which is of course absurd, since I'm trying to predict a binary variable. I then tried to transform the empty variable into a factor but I get this error:

Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts, : The response for family[1] is not of type 'numeric|list|matrix'; don't know what to do. The inla program failed and the maximum number of tries has been reached.

I want, if possible that the model predicts P(empty=1) if possible, as a logistic regression does for example. What am I doing wrong? Am I using the wrong family?

Best regards,

Samuel

[INLA help]

The linear predictor is in a transformed scale, check ‘fitted.values'

—

Haavard Rue

Helpdesk

he...@r-inla.org

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/474d1869-fd1a-4edd-85a1-dfcaf56fa574n%40googlegroups.com.

[Samuel Molano]

Hello,

Thank you for your answer. Before, I had tried to fix this by doing it manually by myself, by taking the inverse of the logit link function (link = 1) at -lp.mean.pred. We have then P(empty=1) = 1/(1+exp(-lp.mean.pred)).  As I understand, this should be the same as mod$summary.fitted.values[index.pred, "mean"]. However, I obtain different values but that are highly correlated (Pearson's correlation coefficient of 0.9918 for my example).

Do you see a difference between the two methods?

Thank you.

Best regards,

Samuel Molano

[Finn Lindgren]

From regular probability theory, E(g(X)) is not equal to g(E(X)), which is why you don't get exactly the same thing.
The posterior mean of the linear predictor cannot be fed into a nonlinear (inverse) link function if what you want is the posterior mean of the nonlinear transformation of the linear predictor.

Finn

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/a6c9e803-c802-43cd-a0cb-af18db370548n%40googlegroups.com.

--

Finn Lindgren
email: finn.l...@gmail.com

[Helpdesk (Haavard Rue)]

doing a plug-in is almost correct, but we do the proper integration, like

p.mean = E ( 1/ (1 + exp(-eta)))

over the posterior for 'eta'

> https://groups.google.com/d/msgid/r-inla-discussion-group/a6c9e803-c802-43cd-a0cb-af18db370548n%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org