Problem with dlkjcorr in varying slopes simulation

[jacksutton]

Problem with dlkjcorr in varying slopes simulation

[rmcelreath]

I guess you mean code block 13.12? Without your code, I can't help you much.

In general, the correlation matrix will read its dimension off the dimension of the multi_normal distribution it is embedded in. If it's not embedded in one, then you'll need to manually specify it's dimension.

[jacksutton]

Yes, exactly, block 13.12. I've been working through the examples in ch 12 & 13. Here is the code, including the data simulation:

simulate data

N_cafes <- 20

library(MASS)
set.seed(5) # used to replicate samples
vary_effects <- mvrnorm(N_cafes, Mu, Sigma)
a_cafe <- vary_effects[,1]
b_cafe <- vary_effects[,2]

N_visits <- 10
afternoon <- rep(0:1, N_visits*N_cafes/2)
cafe_id <- rep(1:N_cafes, each=N_visits)
mu <- a_cafe[cafe_id] + b_cafe[cafe_id] * afternoon
sigma <- 0.5 # sd within cafes
wait <- rnorm(N_visits * N_cafes, mu, sigma)
d <- data.frame(cafe=cafe_id, afternoon=afternoon, wait=wait)

fit model

m13.1 <- map2stan(
alist(
wait ~ dnorm(mu, Sigma),
mu <- a_cafe[cafe] + b_cafe[cafe] * afternoon,
c(a_cafe, b_cafe) ~ dmvnorm2(c(a,b), sigma_cafe, Rho),
a ~ dnorm(0, 10),
b ~ dnorm(0, 10),
sigma_cafe ~ dcauchy(0, 2),
sigma ~ dcauchy(0, 2),
Rho ~ dlkjcorr(2)
),
data = d, iter = 5000, warmup = 2000, chains = 2
)

[rmcelreath]

Your code is missing some variable assignments, like Mu. This is the needed code, starting from the beginning of the chapter:

library(rethinking)

## R code 13.1
a <- 3.5            # average morning wait time
b <- (-1)           # average difference afternoon wait time
sigma_a <- 1        # std dev in intercepts
sigma_b <- 0.5      # std dev in slopes
rho <- (-0.7)       # correlation between intercepts and slopes

## R code 13.2
Mu <- c( a , b )

## R code 13.3
cov_ab <- sigma_a*sigma_b*rho
Sigma <- matrix( c(sigma_a^2,cov_ab,cov_ab,sigma_b^2) , ncol=2 )

## R code 13.5
sigmas <- c(sigma_a,sigma_b) # standard deviations
Rho <- matrix( c(1,rho,rho,1) , nrow=2 ) # correlation matrix

# now matrix multiply to get covariance matrix
Sigma <- diag(sigmas) %*% Rho %*% diag(sigmas)

## R code 13.6
N_cafes <- 20

## R code 13.7
library(MASS)
set.seed(5) # used to replicate example
vary_effects <- mvrnorm( N_cafes , Mu , Sigma )

## R code 13.8
a_cafe <- vary_effects[,1]
b_cafe <- vary_effects[,2]

## R code 13.10
N_visits <- 10
afternoon <- rep(0:1,N_visits*N_cafes/2)
cafe_id <- rep( 1:N_cafes , each=N_visits )
mu <- a_cafe[cafe_id] + b_cafe[cafe_id]*afternoon
sigma <- 0.5  # std dev within cafes
wait <- rnorm( N_visits*N_cafes , mu , sigma )
d <- data.frame( cafe=cafe_id , afternoon=afternoon , wait=wait )

## R code 13.12
m13.1 <- map2stan(
    alist(
        wait ~ dnorm( mu , sigma ),
        mu <- a_cafe[cafe] + b_cafe[cafe]*afternoon,
        c(a_cafe,b_cafe)[cafe] ~ dmvnorm2(c(a,b),sigma_cafe,Rho),
        a ~ dnorm(0,10),
        b ~ dnorm(0,10),
        sigma_cafe ~ dcauchy(0,2),
        sigma ~ dcauchy(0,2),
        Rho ~ dlkjcorr(2)
    ) ,
    data=d ,
    iter=5000 , warmup=2000 , chains=2 )

If that doesn't work for you, then I suspect there is a gremlin in your global R namespace.

[jacksutton]

I forgot to say: Doh! Thank you.