Hello all,
I’m trying to fit a three-level nested linear model. My data are confidential, so I can’t post them, but the response variable is continuous, strictly positive and heteroskedastic. There’s good reason to believe this response decreases linearly with some predictor, so I’d like to create a lognormal regression model.
I generated some data below to make a reproducible example. As an arbitrary example, say we’re measuring the time students spend on classwork each week over the course of a 15-week class, and modeling how this changes linearly with time (week). We collected data for 10 years, following 10 students each year. I’m assuming that the slope and mean intercept vary randomly across years, and that there is an individual-level deviation from the mean intercept. I’m attaching a plot of the first four years of simulated data to give you an idea of what that looks like.
I followed these two examples when writing my model. When I use a normal likelihood, the model runs fine and seems to recover the correct parameter values. However, when I use a lognormal likelihood, initialization fails with the message “Initialization between (-2, 2) failed after 100 attempts.” So, I think there is something I am not understanding about how the lognormal model is implemented. Any suggestions or resources would be appreciated!
Data generation
library(rstan)
set.seed(42)
nyears<-10 #Data from 10 years
nsubj<-10 #10 individuals sampled per year
nweeks <- 15 #15 observations per individual each year
b0_mean <- 20 #Mean main intercept across years
b0_sd <- 3 #SD main intercept across years
b1_mean <- 0.5 #Mean slope across years
b1_sd <- 0.15 #SD slope across years
subj_sd <- 2.1 #SD of individual-level intercept
pop_sd <- 1.18 #Observation level SD
#Generate slopes and intercepts for each year
b0_k <- rnorm(nyears, mean = b0_mean, sd = b0_sd)
b1_k <- rnorm(nyears, mean = b1_mean, sd = b1_sd)
#Generate individual-level intercepts each year
b0_jk <- rnorm(nyears*nsubj, mean = rep(b0_k, each = nsubj), sd = subj_sd)
b0_ijk <- rep(b0_jk, each = nweeks) #The intercept for each observation (for each year, individual, week)
#Generate the slope for each observation
b1_ijk <- rep(b1_k, each = nweeks*nsubj)
#Generate predictor - say week in class
week_ijk <- rep(rep(c(1:nweeks), nsubj), nyears)
range(b0_ijk - b1_ijk * week_ijk) #Make sure response is >0
response_ijk <- rlnorm(nweeks*nyears*nsubj, meanlog = log(b0_ijk - b1_ijk * week_ijk), sdlog = log(pop_sd))
subj_ijk <- rep(c(1:(nsubj*nyears)), each = nweeks) #The individual index
week_ijk <- rep(rep(c(1:nweeks), nsubj), nyears)
year_ijk <- rep(c(1:nyears), each = nweeks*nsubj)
data_df <- data.frame(year = year_ijk, week = week_ijk, subj = subj_ijk, response = response_ijk)
#truepars_df<-distinct(data.frame(b0=b0_ijk, b1=b1_ijk, subj=subj_ijk, year=year_ijk))
year_for_subj <- unique(data_df[c("subj","year")])[,"year"]
dat <- list(n_obs = nrow(data_df),
n_subj = length(unique(data_df$subj)),
n_years = length(unique(data_df$year)),
subj = subj_ijk,
year = year_ijk,
year_for_subj = year_for_subj,
y_ijk = response_ijk,
x_ijk = week_ijk)
Stan code
data {
// Define variables in data
// Number of level-1 observations
int<lower=0> n_obs;
// Number of level-2 clusters
int<lower=0> n_subj;
// Number of level-3 clusters
int<lower=0> n_years;
// Cluster IDs
int<lower=1, upper=n_subj> subj[n_obs];
int<lower=1, upper=n_years> year[n_obs];
// Level 3 look up vector for level 2
int<lower=1, upper=n_years> year_for_subj[n_subj];
// Continuous outcome
vector[n_obs] y_ijk;
//Predictor
vector[n_obs] x_ijk;
}
parameters {
// Population intercept and slope
real beta_0;
real beta_1;
// Level-1 errors
real<lower=0> sigma_e0;
// Level-2 random effect
vector[n_subj] u_0jk;
real<lower=0> sigma_u0jk;
// Level-3 random effects
vector[n_years] u_0k;
real<lower=0> sigma_u0k;
vector[n_years] u_1k;
real<lower=0> sigma_u1k;
}
transformed parameters {
// Individual mean
real mu[n_obs];
// Varying intercepts
vector[n_subj] beta_0jk;
vector[n_years] beta_0k;
// Varying slope
vector[n_years] beta_1k;
// Varying intercepts definition
// Level-3 with centered parametrization
beta_0k = beta_0 + u_0k * sigma_u0k;
// Level-2
beta_0jk = beta_0k[year_for_subj] + u_0jk * sigma_u0jk;
// Varying slope definition
beta_1k = beta_1 + u_1k * sigma_u1k;
// Individual mean
for (i in 1:n_obs) {
mu[i] = beta_0jk[subj[i]] - beta_1k[year[i]] * x_ijk[i];
}
}
model {
// Prior part of Bayesian inference
beta_0 ~ normal(0, 100);
beta_1 ~ normal(0,10);
sigma_e0 ~ normal(0, 10);
sigma_u0k ~ normal(0,10);
sigma_u0jk ~ normal(0,10);
sigma_u1k ~ normal(0,10);
// Random effects distribution
u_0k ~ normal(0, 1);
u_0jk ~ normal(0, 1);
u_1k ~ normal(0, 1);
// Likelihood
//y_ijk ~ normal(mu, sigma_e0);//Normal liklihood works
y_ijk ~ lognormal(log(mu), log(sigma_e0)); //lognormal doesn't work
}
}
Running the model
fileName <- "./three level lognormal model noncentered.stan"
stan_code <- readChar(fileName, file.info(fileName)$size)
cat(stan_code)
resStan_nc <- stan(model_code = stan_code, data = dat,
chains = 3, iter = 3000, warmup = 500, thin = 10)