random slope for categorical predictor

467 views
Skip to first unread message

szm...@stanford.edu

unread,
Apr 4, 2018, 8:55:31 PM4/4/18
to R-inla discussion group
Hello!

I was wondering if it is possible model random slope for categorical predictor in INLA. In lme4 syntax, it will be

lme.res <- lmer(BETA ~ 1 + POPULATION + (0 + POPULATION|VARID), data=data)

where VARID and POPULATION are both factors. I tried

inla.res <- inla(BETA ~ 1 + POPULATION + f(VARID, POPULATION, model='iid'),
                 control.compute=list(config = TRUE), data=data)

It runs, but gave a warning message:

Warning message in `[<-.factor`(`*tmp*`, is.na(www), value = 0):
“invalid factor level, NA generated”


Helpdesk

unread,
Apr 5, 2018, 2:20:12 AM4/5/18
to szm...@stanford.edu, R-inla discussion group
population has to be numeric in f(varid, population), as this term is
like

u[i] * population[i]


where u is a random effect indexed by varid
> --
> 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 post to this group, send email to r-inla-discussion-group@googlegr
> oups.com.
> Visit this group at https://groups.google.com/group/r-inla-discussion
> -group.
> For more options, visit https://groups.google.com/d/optout.

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

Thierry Onkelinx

unread,
Apr 5, 2018, 3:39:36 AM4/5/18
to he...@r-inla.org, szm...@stanford.edu, R-inla discussion group
You'll need one of the "iidnd" models to fit this if you want to
replicate the (0 + POPULATION|VARID) model from lme4. The lme4 model
also fits the correlation among the random slopes. See
https://www.math.ntnu.no/inla/r-inla.org/doc/latent/iid123d.pdf

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE
AND FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry....@inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no
more than asking him to perform a post-mortem examination: he may be
able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does
not ensure that a reasonable answer can be extracted from a given body
of data. ~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////
> To post to this group, send an email to r-inla-disc...@googlegroups.com.

Mamie Wang

unread,
Apr 5, 2018, 11:58:33 AM4/5/18
to Thierry Onkelinx, r-inla-disc...@googlegroups.com, he...@r-inla.org
Thanks, Thierry! POPULATION has four levels and I made them dummy variables using model.matrix. And here is the updated model using iid4d: 

nod2.data$VARID.int <- as.numeric(as.factor(nod2.data$VARID))
population <- model.matrix(~0+POPULATION, data=nod2.data)
nod2.data <- cbind(nod2.data, population)
n.loc <- max(nod2.data$VARID.int)
nod2.data$VARID2.int <- nod2.data$VARID.int + n.loc
nod2.data$VARID3.int <- nod2.data$VARID.int + 2 * n.loc
nod2.data$VARID4.int <- nod2.data$VARID.int + 3 * n.loc

inla.res <- inla(BETA ~ 1 + POPULATION + f(VARID.int, POPULATIONaj, model='iid4d', n=n.loc * 4) +
                                          f(VARID2.int, POPULATIONfc, copy='VARID.int') +
                                          f(VARID3.int, POPULATIONfin, copy='VARID.int') +
                                          f(VARID4.int, POPULATIONnfe, copy='VARID.int'),
                 control.compute=list(config = TRUE), 
                 data=nod2.data)

The estimates for fixed effect and the variance for the random effects are quite similar between lme4 and INLA, but the INLA estimate for the correlation seems to be more conservative compared with the lme4 as shown in the image attached below. I also compared the same model in brms which uses rstan as backend and same syntax for formula interface to fit multilevel models, which also gives lower estimate. I was wondering why is the likelihood method gives higher correlation estimate compared to the Bayesian methods?



In addition, I was wondering if it is possible to implement LKJ prior for the correlation matrix between random effects in INLA as in Stan? 

Mamie

Mamie Wang
Biomedical Informatics MS Candidate, 2018
Stanford University, School of Medicine
LinkedIn Mamie Wang | Skype szmamie






Reply all
Reply to author
Forward
0 new messages