coding random slope model (with or without an uncorrelated random intercept)

[Julien]

Hi,

 

Is there a simple way of coding random slope model with INLA?

 

For exemple, if I want to include a random slope in my model for a contious variable "x" that varies in different "region".

 

With the package LME4, a code like this would do the job:

 

y ~ x + (0+x|region) 

 

If I wanted to include an uncorrelated random intercept, it will become y ~ x + (1|region) + (0+x|region)

 

Could I reproduce these models with INLA? if yes, how?

 

Best regards,

 

Julien

[Håvard Rue]

On Wed, 2011-04-06 at 12:13 -0700, Julien wrote:
> Hi,
>
> Is there a simple way of coding random slope model with INLA?

yes

> For exemple, if I want to include a random slope in my model for a
> contious variable "x" that varies in different "region".
>
> With the package LME4, a code like this would do the job:
>
> y ~ x + (0+x|region)
>
> If I wanted to include an uncorrelated random intercept, it will
> become y ~ x + (1|region) + (0+x|region)
>
> Could I reproduce these models with INLA? if yes, how?

one example is here

http://www.r-inla.org/models/tools#TOC-Copying-a-model

you need to copy a model, and then its like normal.

just make sure you get the ordering of the indices correct in the iid2d
model; see

http://www.r-inla.org/models/latent-models

and the iid2d model.

hope this helps.
H

--
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
Fax : +47-7359-3524 Email: havar...@math.ntnu.no

This message was created in a Microsoft-free computing environment.

[Håvard Rue]

On Wed, 2011-04-06 at 12:13 -0700, Julien wrote:

I presumed you ment that you wanted to correlated the intercept and the
slope, but if you do not want that, then

f(region, covariate, model ="iid")

is a random slope model, where

region[i] is the region for data y[i]
covariate[i] is the covariate for data y[i]

so that this gives with z = covariate,

+ \beta_region * z_i ....

the second argument in f() is an optional weight for the f-term.

[Julien]

I am going to try that.

Thanks

Julien

Håvard Rue a écrit :

[Julien]

Thank you Havard for your previous response. It is now clear to me how to code a random slope model. However, I had trouble to obtain convergence with a simple random intercept model (see below) whereas it worked fine using lme4 package (see below). Also, I have tried to code a combined random intercept-slope model (see below) in INLA using the copy function that you mentioned in your previous message but I did not have any success. Could you indicate how to replicate the following models (fitted with lme4 packages) using the simulated data set below?

Thank you for your time,

Julien   

########################

## Simulated data set ##

########################

 

x <- sort(unique(runif(100, 0, 20)))

y1 <- 60 - 2*x + rnorm(length(x), 0, 4)

y2 <- 70 - 0.6*x + rnorm(length(x), 0, 4)

y3 <- 45 + 1*x + rnorm(length(x), 0, 4)

y4 <- 35 + 1.3*x + rnorm(length(x), 0, 4)

plot(x,y1, ylim=c(20,80), xlim=c(0,30))

points(x, y2, col="red")

points(x, y3, col="green")

points(x, y4, col="orange")

abline(a=coef(lm(y1~x))[1], b=coef(lm(y1~x))[2])

abline(a=coef(lm(y2~x))[1], b=coef(lm(y2~x))[2])

abline(a=coef(lm(y3~x))[1], b=coef(lm(y3~x))[2])

abline(a=coef(lm(y4~x))[1], b=coef(lm(y4~x))[2])

mydata <- data.frame(x1 = rep(x, times = 4), y = c(y1,y2,y3,y4), bloc = rep(c(1,2,3,4), each = length(x)), fuzzy = rep(1, 400))

 

##########

## LME4 ##

##########

require(lme4)

 

# no random effect

model_no_random <- lmer( y ~ x1 + (1|fuzzy), data = mydata)

summary(model_no_random)

 

# random intercept only

model_intercept <- lmer( y ~ x1 + (1|bloc), data = mydata)

summary(model_intercept)

 

# random slope only

model_slope <- lmer( y ~ x1 + (0 + x1|bloc), data = mydata)

summary(model_slope)

 

# random slope of x1 within bloc with correlated intercept

model_slope_intercept_correlated <- lmer( y ~ x1 + (x1|bloc), data = mydata)

summary(model_slope_intercept_correlated)

 

# uncorrelated random intercept and random slope within bloc

model_slope_intercept_uncorrelated <- lmer(y ~ x1 + (1|bloc) + (0 + x1|bloc), data = mydata)

summary(model_slope_intercept_uncorrelated)

 

##########

## INLA ##

##########

require(INLA)

 

# no random effect

formule_no_random = y ~ x1

model_no_random_INLA = inla.cpo(formula=formule_no_random, data=mydata, family="gaussian", control.compute=list(dic=TRUE))

summary(model_no_random_INLA)

 

# random intercept only

formule_random_intercept_only = y ~ x1 + f(bloc, model = "iid")

model_random_intecept_only_INLA = inla.cpo(formula=formule_random_intercept_only, data=mydata, family="gaussian", control.compute=list(dic=TRUE))

summary(model_random_intecept_only_INLA)

 

# random slope only

formule_random_slope_only = y ~ x1 + f(bloc, x1, model = "iid")

model_random_slope_only_INLA = inla.cpo(formula=formule_random_slope_only, data=mydata, family="gaussian", control.compute=list(dic=TRUE))

summary(model_random_slope_only_INLA)

 

# random slope of x1 within bloc with correlated intercept

formule_random_correlated_intercept_slope = y ~ x1 + ?

model_random_correlated_intercept_slope_INLA = inla.cpo(formula=formule_random_correlated_intercept_slope, data=mydata, family="gaussian", control.compute=list(dic=TRUE))

summary(model_random_correlated_intercept_slope_INLA)

 

# uncorrelated random intercept and random slope within bloc

formule_random_uncorrelated_intercept_slope = y ~ x1 + ?

model_random_uncorrelated_intercept_slope_INLA = inla.cpo(formula=formule_random_uncorrelated_intercept_slope, data=mydata, family="gaussian", control.compute=list(dic=TRUE))

summary(model_random_uncorrelated_intercept_slope_INLA)

 

[Håvard Rue]

On Tue, 2011-04-19 at 12:56 -0700, Julien wrote:

I think this is what you want. I hope I have understoof the models
correct; I'm not familiar to the (x1|bloc) notation...

results will surely will prior dependent.

also, I needed to make the overall intercept have a proper prior using
prec.intercept=0.0001 (default 0).

Best
H

********************************

x <- sort(unique(runif(100, 0, 20)))
y1 <- 60 - 2*x + rnorm(length(x), 0, 4)
y2 <- 70 - 0.6*x + rnorm(length(x), 0, 4)
y3 <- 45 + 1*x + rnorm(length(x), 0, 4)
y4 <- 35 + 1.3*x + rnorm(length(x), 0, 4)

plot(x,y1, ylim=c(20,80), xlim=c(0,30))
points(x, y2, col="red")
points(x, y3, col="green")
points(x, y4, col="orange")
abline(a=coef(lm(y1~x))[1], b=coef(lm(y1~x))[2])
abline(a=coef(lm(y2~x))[1], b=coef(lm(y2~x))[2])
abline(a=coef(lm(y3~x))[1], b=coef(lm(y3~x))[2])
abline(a=coef(lm(y4~x))[1], b=coef(lm(y4~x))[2])

mydata <- data.frame(x1 = rep(x, times = 4),
y = c(y1,y2,y3,y4),
bloc = rep(c(1,2,3,4),
each = length(x)),
fuzzy = rep(1, 400))

require(lme4)
require(INLA)

## no random effect
##

model_no_random <- lmer( y ~ x1 + (1|fuzzy), data = mydata)
summary(model_no_random)

formule_no_random = y ~ x1

model_no_random_INLA = inla.cpo(formula=formule_no_random,
data=mydata, family="gaussian", control.compute=list(dic=TRUE))
summary(model_no_random_INLA)

## random intercept only
##

model_intercept <- lmer( y ~ x1 + (1|bloc), data = mydata)
summary(model_intercept)

formule_random_intercept_only = y ~ x1 + f(bloc, model = "iid")
model_random_intecept_only_INLA =
inla.cpo(formula=formule_random_intercept_only,
data=mydata, family="gaussian", control.compute=list(dic=TRUE),
control.fixed = list(prec.intercept = 0.0001))
summary(model_random_intecept_only_INLA)

## random slope only
##

model_slope <- lmer( y ~ x1 + (0 + x1|bloc), data = mydata)
summary(model_slope)

formule_random_slope_only = y ~ x1 + f(bloc, x1, model = "iid")

model_random_slope_only_INLA =
inla.cpo(formula=formule_random_slope_only,
data=mydata, family="gaussian", control.compute=list(dic=TRUE),
control.fixed = list(prec.intercept = 0.0001))
summary(model_random_slope_only_INLA)

## random slope of x1 within bloc with correlated intercept
##

model_slope_intercept_correlated <- lmer( y ~ x1 + (x1|bloc), data =
mydata)
summary(model_slope_intercept_correlated)

n.block = max(mydata$block) ## = 4
mydata$i.intercept = mydata$bloc
mydata$j.intercept = mydata$bloc + n.block ## see doc for iid2d

formule_random_correlated_intercept_slope = y ~ x1 +

f(i.intercept, model="iid2d", n=2*n.block) +
f(j.intercept, x1, copy="i.intercept")
model_random_correlated_intercept_slope_INLA = inla.cpo(
formula=formule_random_correlated_intercept_slope,

data=mydata, family="gaussian", control.compute=list(dic=TRUE))
summary(model_random_correlated_intercept_slope_INLA)

## uncorrelated random intercept and random slope within bloc
##

model_slope_intercept_uncorrelated <- lmer(y ~ x1 + (1|bloc) + (0 + x1|
bloc), data = mydata)
summary(model_slope_intercept_uncorrelated)

formule_random_uncorrelated_intercept_slope = y ~ x1 +
f(i.intercept, model="iid") + f(j.intercept, x1, model="iid")
model_random_uncorrelated_intercept_slope_INLA =
inla.cpo(formula=formule_random_uncorrelated_intercept_slope,
data=mydata, family="gaussian", control.compute=list(dic=TRUE),
control.fixed = list(prec.intercept = 0.0001))
summary(model_random_uncorrelated_intercept_slope_INLA)

[inti pedroso]

Hi,

I am writing in within this post because my question is related to the previous one, I hope that is OK.

The material from the previous answer is quite useful but I still have a doubt. 

In my case I am interested on a model like:

Y ~ N(mu_y, sd_y)

mu_y = b0*1 + b1*A

b1_{i} = b2*1 + B + C + D

In my data Y has 6100 values and A ~ 1300 levels. B,C and D are variables with few levels (~5) that in the design are nested within each A value. My primary interest is to arrive to a "average" Y value for each A level after controlling for the effect of B,C and D. I could do it as a regression for each A level but the design is not balanced and I want to see if I can exploit a multi-level model to improve the estimates. 

I cannot get my head around writing the model in INLA. In lme4 language it would be something like> lmer( Y ~ 1 + (1 + B + C + D|A) ). I am not interested on the correlation between the  B, C and D effects. 

I think my problem is that I cannot workout how to write the index variables (the i and j.intercept in the previous response). 

Your help will be very much appreciated. 

I am sorry if this is a simple answers, I am new to both multilevel models and INLA.

Cheers,

Inti

[Håvard Rue]

On Wed, 2011-07-20 at 08:35 -0700, inti pedroso wrote:
> Hi,
> I am writing in within this post because my question is related to the
> previous one, I hope that is OK.
> The material from the previous answer is quite useful but I still have
> a doubt.
> In my case I am interested on a model like:
>
>
> Y ~ N(mu_y, sd_y)
> mu_y = b0*1 + b1*A
> b1_{i} = b2*1 + B + C + D

A B C and D are factors with one coof for one level; sorry.

[Håvard Rue]

On Wed, 2011-07-20 at 08:35 -0700, inti pedroso wrote:

> Hi,
> I am writing in within this post because my question is related to the
> previous one, I hope that is OK.
> The material from the previous answer is quite useful but I still have
> a doubt.
> In my case I am interested on a model like:
>
>
> Y ~ N(mu_y, sd_y)
> mu_y = b0*1 + b1*A
> b1_{i} = b2*1 + B + C + D
>

I wrote down the model as I understood it; its in the attachment.
Here is how to implement it:

data$A = as.numeric(data$A)
data$B = as.factor(data$B)
data$C = as.factor(data$C)
data$D = as.factor(data$D)

formula = y ~ 1 + A + B:A + C:A + D:A

## since A is numeric and B is a factor, it will expand into
## \beta_{B_i}^B * A_i

## small test-example

y = 1:100
B = as.factor(rep(1:5, each = 20))
A = as.numeric(y)
formula = y ~ 1 + B:A
r = inla(formula, data = data.frame(y,A,B))

## so that for each of the 5 \beta_{B}'s, one for each factor, should
## be 1.

> r$summary.fixed[,"mean"]
(Intercept) B1:A B2:A B3:A
B4:A
5.293578909e-10 1.000000000e+00 1.000000000e+00 1.000000000e+00
1.000000000e+00
B5:A
1.000000000e+00

best
H

[Julien Beguin]

It would be easier if you could send us a simulated data set representing your own data structure? 

 

Tell me if I am wrong, you would like to treat variables B, C, D as cross-random slopes in your model, right? 

 

Julien

---

De : r-inla-disc...@googlegroups.com [r-inla-disc...@googlegroups.com] de la part de inti pedroso [intip...@gmail.com]
Date d'envoi : 20 juillet 2011 11:35
À : r-inla-disc...@googlegroups.com
Cc : Julien Beguin
Objet : Re: Re : Re: [r-inla] coding random slope model (with or without an uncorrelated random intercept)

[inti pedroso]

Hi, Thanks for your quick response.

I have attached an R binary dataset with some real data. I have just renamed the variables.

I am sorry if I do not get this quickly. I am just starting with this type of analyses.

this is the code

##################### START R CODE

library(INLA)

load("example_data.RData") # this is the attached file, all variables are on data frame "z"

data=z

#### what variables do we have?

colnames(data)

# [1] "Y" "B" "C" "D" "A" 

# Y is the response variable, 

# A are treatments 

# B, C and D are variables that indicate variation on the treatments, e.g., doses. These variations are common for different treatments, I think you usually refer to this as crossed with treatments?

# I am interested on the mean effect of each A after accounting for these variations B, C and D. Or said in another way. 

####  What do I want?

# I want the intercept of the model to vary for each treatment A. Each treatments is different so it should be coded as a factor in the model.

data$A = as.factor(data$A)

data$B = as.factor(data$B)
data$C = as.factor(data$C)
data$D = as.factor(data$D)

# As I understand the fitting treatments A should be treated as coming from a common distribution. So the simplest model, not accounting for the effect of BCD would be to pool all observations for each treatment A independent of BCD variables but also will shrink the estimates towards the grand mean. 

# I am more familiar with the coding on the lme4 package. The simplest model using lmer would be: 

library(lme4)

(r_lmer_A_random<-lmer(Y ~ 1 + (1|A), data=data))

# where my value of interest would be the grand mean + the mean under each treatment

fixef(r_lmer_A_random) + head(ranef(r_lmer_A_random)[["A"]]) # fixef extract the fixed-effects estimates and ranef the random effects

#     (Intercept)

# 1  0.005945035

# 2 -0.035491279

# 3 -0.034360372

# 4  0.003965238

# 5 -0.002544253

# 6  0.015596639

# coding this on INLA-R would be

formula_A_random = Y ~ 1 + f(A,model="iid")

r_A_random = inla(formula_A_random, data = data)

# inspecting the results

head(r_A_random$summary.random$A)

#  ID          mean         sd  0.025quant      0.5quant 0.975quant          kld

# 1  1  0.0018349540 0.04284262 -0.08240764  0.0018187206 0.08606898 1.153688e-19

# 2  2 -0.0388059098 0.04313830 -0.12451523 -0.0384904150 0.04502659 0.000000e+00

# 3  3 -0.0377066768 0.04143557 -0.11990624 -0.0374445043 0.04290222 0.000000e+00

# 4  4 -0.0001072996 0.04472216 -0.08812133 -0.0001076699 0.08780380 2.465190e-31

# 5  5 -0.0064911496 0.04577419 -0.09677893 -0.0064236470 0.08329513 1.617034e-18

# 6  6  0.0112986571 0.04579466 -0.07850862  0.0111775134 0.10171264 6.462351e-18

# a quick look suggest the results are similar.

# Now, what I do not understand is how to code the effect of BCD.

# in the lme4 coding my model of interest would be 

(r_lmer_full<-lmer(Y ~ 1 + (1 + B + C + D|A), data=data)) # it does fit but with a warning.

# inspect the results, again I am interested on the intercept.

fixef(r_lmer_full) + head(ranef(r_lmer_full)[["A"]])[,"(Intercept)"]

# 1   1.694557e-03

# 2   -1.974743e-04

# 3   -7.863843e-06

# 4   2.102188e-03

# 5   2.032537e-03

# 6   2.327871e-03

# My doubt is on how to code the above model on inla

### END R CODE

Again, thanks a lot for your time.

I hope the problem is question now.

Best Regards,

Inti

[Julien]

Hi,

 

I tried your exemple but could not make it working in lme4 because of insufficient memory... (while it converged more easely in removing the factor B). I also have to admit that the analysis of several nested-crossed random effects involving only factors is at the edge of my knowledge so the advice of people more familiar with such data structure will be more relevant. Sorry for not being more helpful. If this model is workable, it would be very informative to know how to code it in INLA...

 

Julien    

[Helpdesk]

in that case you would do

f(j.intercept,Vist_10, model="iid",
hyper = list(prec = list(prior = "loggamma",
param = c(0.001, 0.001))))


the Gamma(eps,eps) is not a good choice, and there are better options
available; see


https://projecteuclid.org/euclid.ss/1491465621


Best
H

On Wed, 2018-11-14 at 11:48 -0800, Belay Birlie wrote:
> Dear Prof,
>

> Following your suggestion some time ago, If I want to change the
> default prior to lets say Gamma(0.001,0.001) in the uncorrelated
> random effect model below, what should I do? I tried this
> formula2 = y ~ ...+Vist_10+f(i.intercept, model="iid",param
> =c(0.001,0.001))
> +
> f(j.intercept,Vist_10, model="iid",param =c(0.001,0.001))
> and got the following error message
> Error in +f(j.intercept, Vist_10, model = "iid", param = c(0.001,
> 0.001)) :
> invalid argument to unary operator
>
> Kind regards,
> Belay

Helpdesk
he...@r-inla.org