**Model:** I am working with a model that can be described, at the bird-eye level, as a random slope logistic model. I have 0/1 outcomes, some predictors, random intercepts and random slopes for my individuals. I am fitting the models I am satisfied with, goodness of fit and convergence and all.

**Problem**: suppose I am given a new individual with their predictors and their trajectory of the 0/1 outcomes. Can I run the sampled chains on these to obtain the random effects? Section 27.7 provides pointers in the â€śflatâ€ť case of i.i.d. data on how to pick the parameters that have been simulated, and feed them to the `generated quantities`

without the model to sample from. But in my case, I have to deal with the random effects that I am declaring as parameters â€“ and the random effects for the new data, obviously, have not been sampled. I cannot simply add that to the model and the sampler as the stuff runs for about 50 core-hours, so I cannot afford to re-run it for another 2000 new individuals that I need to produce the model predictions for.

**Reprex** as an Rmd fragment, so that the Stan model is in the middle. The goal is to get the vector `u_i[2]`

for the new data.

##
title: â€śRandom effects reprexâ€ť

output: html_document

```
knitr::opts_chunk$set(echo = TRUE,
warning = FALSE,
message = FALSE,
dev = "svg") ## export images as SVG
library(ggplot2)
library(magrittr)
library(here)
library(dplyr)
library(tidyr)
library(ggplot2)
library(tibble)
library(rstan)
library(bayesplot)
library(knitr)
```

# Simulate data

```
N <- 10
set.seed(1)
my_df <- data.frame( id=1:N, interc=rnorm(n=N,m=0.5,s=0.3), slope=rnorm(n=N, m=-0.5, s=0.3), n=sample(5:8, N, replace=T) )
my_df <- my_df[ rep(1:N, my_df$n), ]
my_df %>% group_by(id) %>% mutate(t=row_number()) %>% ungroup() %>%
mutate(y=as.integer(interc+slope*t/3+rnorm(nrow(my_df),m=0,s=0.3)>0)) -> my_df
```

## Stan code

```
data {
int<lower=1> N; // total number of observations
int<lower=0,upper=1> outcome[N]; // outcome
int<lower=1> N_people; // number of unique individuals
int<lower=1,upper=N_people> seq_person_id[N]; // SEQUENTIAL person id -- important for indexing
vector[N] predictor; // predictor that goes with random slope
}
parameters {
// variance components
vector<lower = 0>[2] tau_u;
// random effects for people
vector[2] u_i[N_people];
// correlation matrix
corr_matrix[2] rho_matrix;
}
transformed parameters {
cov_matrix[2] Sigma;
vector[N] theta;
// covariance matrix
Sigma = quad_form_diag(rho_matrix,tau_u);
// logit of outcome
for (t in 1:N)
theta[t] = u_i[seq_person_id[t],1]
+ u_i[seq_person_id[t],2]*predictor[t];
}
model {
// outcome
outcome ~ bernoulli(inv_logit(theta));
// variances
tau_u ~ cauchy(0,1);
// correlation
rho_matrix ~ lkj_corr(1);
// random effects
u_i ~ multi_normal( [0, 0]', Sigma );
}
```

# Bayesian computation

## Prepare the data for Stan

```
stan_data <- list(
N = nrow(my_df), # int<lower=1> N;
outcome = unlist(my_df[,"y"]), # int<lower=0,upper=1> outcome[N];
N_people = (my_df %>% distinct(id) %>% nrow()), # int<lower=1> N_people;
seq_person_id = unlist(my_df[,"id"]), # int<lower=1,upper=N_people> seq_person_id[N];
predictor = unlist(my_df[,"t"]) # vector[N] predictor;
)
```

## Sample from the model

```
stan_mcmc_rslope_logit <- rstan::sampling(
random_slope_logit, data=stan_data,
chains = 2, seed = 1234,
warmup = 1000, iter = 2000,
control=list(adapt_delta=0.9, stepsize_jitter = 0.1),
refresh = 100,
init = list( list(tau_u=c(0.2,0.2)), list(tau_u=c(0.3,0.3)) )
)
```

## New observations

```
new_df <- data.frame( y_new = c(0,1,1,0,1), t_new=1:5)
```

What now??