Multinomial data with covariates using INLA

I’ve recently started using the new R-INLA package developed by Havard Rue, Janine Illian Finn Lindgren et al. It’s quite easy to start using, has a lot of depth and is extremely useful for spatial data. For another project, I was trying to figure out how to fit a multinomial distribution in INLA for predicting age structures as a function of covariates. While a good example exists on the website, it wasn’t immediately clear to me how to add covariates to the model. After a discussion with Havard Rue, I came up with this example which both shows how to add covariates and a comparison of analysis of simple data using INLA versus another existing R package.

I apologize in advance for any crummy table or equation formatting. I haven’t quite figured out LaTeX withouth mathjax on WordPress.com.

The analysis is based on this example describing the use of the multinom() function from the nnet package.

require(foreign)
require(nnet)
require(xtable)
ml <- read.dta("http://www.ats.ucla.edu/stat/data/hsbdemo.dta")
ml$prog2 <- relevel(ml$prog, ref = "academic")

This is what the data look like:

id	cid	write
45.00	1	35.00
108.00	1	33.00
15.00	1	39.00
67.00	1	37.00
153.00	1	31.00
51.00	1	36.00

Here, “id” indicates the student-level indicator and “prog” indicates the possible choices or outcomes that we are interested in estimating probabilities for. We will use both of these to index the random effects appropriately. The continuous covariate “write” will be used in the analysis.

Because it makes the assumption that everything is based on a Gausian random field, INLA likes is when everything is centered and normalized. To meet those assumptions, we will centre the data. As well, we need to explicitly state the number of individuals in each category for each observation. For these data, we have one choice per person, so we need to fill in some zeros:

ml$number = 1
m1 = ml
m1$number = 0
m2 = m1
m1$prog2 = ifelse(ml$prog2 == "academic", "general", ifelse(m1$prog2 == "general", 
    "vocation", "academic"))
m2$prog2 = ifelse(m1$prog2 == "academic", "general", ifelse(m1$prog2 == "general", 
    "vocation", "academic"))
mI = rbind(ml, m1, m2)

The centering function and transformation:

my.centre = function(x) x - mean(x)
mI$C.write = my.centre(m2$write)
ml$C.write = my.centre(ml$write)

With the data massaging out of the way, we can run some models. First, with multinom(),

test2 <- multinom(prog2 ~ 1 + C.write, data = ml)

the results look like this:

	(Intercept)	C.write
general	-0.77	-0.07
vocation	-0.86	-0.12

Now we can compare results from INLA, providing an example of INLA code. The basic multinomial model is provided by the example in the FAQ section of the INLA site , which allows us to estimate an intercepts-only model. In order to add covariates, we need to recognize that the Poisson approximation to the Multinomial distributoin assumes that the observed counts \( Y_{ij} \) in each of the \( J \) age categories are assumed to be Poisson random variables with means satisfying the equation:

\[ \mbox{log}(\mu_{ij})=\nu + \theta_i + \alpha_j +x_i \beta_j. \]

where \( \nu \) is a common intercept, assumed to be 0, \( \theta_i \) relates to the number of samples drawn for each Multinomial variable, \( \alpha_j \) specifies the relative strength in each category and \( \beta_j \) is a vector of \( J-1 \) covariates specific for each of the covariates in matrix \( x_i \). Thus, to estimate the influence of a covariate, we need to specify that we want \( J-1 \) covariates relating to the response in the relevant categories, assuming that the $J$th category will be treated as the reference (0).

We therefore need some more data manipulation to get the appropriate dummy variables. Here, our reference is “academic”, and we assume the coefficient for the reference group is 0.

mI$C.write_pg = with(mI, C.write * (prog2 == "general"))
mI$C.write_pv = with(mI, C.write * (prog2 == "vocation"))

Next, a formula. INLA extends the standard formula() structured used by linear models in R using the f() notation to define random effects. In this model, we define random effects for both the Observation and category levels, which allows us to do the approximation of the Multinomial as a set of poisson distributions.

formula = number ~ -1 + f(id, model = "iid", hyper = list(prec = list(initial = log(1e-06), 
    fixed = TRUE))) + 
f(prog2, model = "iid", constr = FALSE, hyper = list(prec = list(initial = log(0.001), 
    fixed = TRUE))) + C.write_pg + C.write_pv

Note that we have added two covariates C.write_pg and C.write_pv to estimate coefficients specifically for each category level (except the reference). Now we can pass the formula and data to INLA:

require(INLA)
r <- inla(formula, family = "poisson", data = mI, control.compute = list(config = TRUE))

Some brief notes on notation: the formula, family and data are self-explanatory, but the control.compute allows us to generate model predictions if we need them. Now we can compare the intercepts and coefficients for the multinom() and INLA models:

## compare intercepts Need to make sure they both sum to 1
INLA_Mn = r$summary.random$prog2$mean
INLA_res2 = exp(INLA_Mn)/sum(exp(INLA_Mn))

Multi_Mn = coef(test2)[, 1]
test2_res = exp(Multi_Mn)/(1 + sum(exp(Multi_Mn)))
MN_res2 = c(1 - sum(test2_res), test2_res)
Intercepts = data.frame(INLA_res2, MN_res2)

	INLA_res2	MN_res2
	0.53	0.53
general	0.24	0.25
vocation	0.22	0.22

And the coefficients too:

coef = data.frame(Multinomial_fn = summary(test2)$coefficients[, 2], INLA_result = r$summary.fixed[, 
    1])

	Multinomial_fn	INLA_result
general	-0.07	-0.07
vocation	-0.12	-0.12

As my daughter says: “Same:Same”. This may seem like a bit of a complicated exercise to achieve something almost identical to what has been done in a few lines of code in multinom(). But this simple model is now extensible to a spatio-temporal framework using the spde and random walk models described in the INLA tutorials. Stay tuned…