This is a bit of a follow-up to the thread here, where I’m now looking for feedback on the solution I’ve come up with. To recap, this is the 8schools scenario but in addition to modelling variability in the mean from school to school, I also want to model the variability in the within-school variability. That is, the standard 8schools model takes the observed standard deviation within a given school as an error-free measurement, while I want to express that there is likely some error in that measurement and partially-pool information on this quantity across schools. To do this you need to know the mean, SD and N (number of students observed) for each school (c.f. typical 8schools models, where the SD & N are collapsed to an SE before input to the model).
Below is R code to simulate some data and sample given the Stan model, which is also included below. Note that I attempt to achieve weak/data-informed priors by standardizing the means and log-SDs.
I seem to be able to get posteriors that correspond with whatever input parameters I use to simulate data, but I think there still might be something missing. Specifically, I feel like the line in the model expressing the observed SD as a sample from a population should include information about how much measurement effort was included (i.e. how many students were observed), because it makes sense to “believe” SDs from large samples more than those from small samples, no?
Edit: After some further playing, I’m sure something more needs to be done to account for measurement effort, as making K vary from school to school (replacing K = rep(100,times=N) with K = runif(N,2,100) in the R code) leads to a posterior on variability_sd_log that doesn’t match the input of sd_of_sds_log (the posterior distribution seems to overestimate the simulation parameter). Any ideas?
library(rstan)
#population parameters
mean_of_means = 100 #arbitrary value
sd_of_means = 20 #arbitrary value
mean_of_sds_log = 4 #arbitrary value
sd_of_sds_log = .1 #probably want this one to be small-ish
#measurement parameters
N = 1000 # ex. number of schools
K = rep(100,times=N) # ex. number of students observed per school
#generate true means & sds per school
true_means = rnorm(N,mean_of_means,sd_of_means)
true_sds = exp(rnorm(N,mean_of_sds_log,sd_of_sds_log))
#loop to generate raw data and consequent observed means and sds
obs_means = rep(NA,times=N)
obs_sds = rep(NA,times=N)
for(i in 1:N){
individual_obs = rnorm(K[i],true_means[i],true_sds[i])
obs_means[i] = mean(individual_obs)
obs_sds[i] = sd(individual_obs)
}
#package data for Stan
data_for_stan = list(
N = N # number of experiments
, K = K # number of observations per experiment
, obs_means = obs_means # observed mean from each experiment (priors assume that these have been standardized to a mean of 0 and sd of 1)
, obs_sds = obs_sds # observed standard deviations in each experiment
)
#compile model
meta_mod = rstan::stan_model('meta.stan')
#sample model
post = rstan::sampling(
object = meta_mod
, data = data_for_stan
, chains = 4
, cores = 4
)
check_hmc_diagnostics(post)
summary(
post
, par = c(
'centrality_intercept'
, 'centrality_sd'
, 'variability_intercept_log'
, 'variability_sd_log'
)
)$summary
Contents of meta.stan:
data {
int<lower=2> N; // number of experiments
vector<lower=2>[N] K ; // number of observations per experiment
vector[N] obs_means; // observed mean from each experiment (priors assume that these have been standardized to a mean of 0 and sd of 1)
vector<lower=0>[N] obs_sds; // observed standard deviations in each experiment
}
transformed data{
//standardize means for easy weakly-informed priors
real mean_of_means = mean(obs_means) ;
real sd_of_means = sd(obs_means) ;
vector[N] obs_means_z = (obs_means-mean_of_means)/sd_of_means ;
// log then standardize sds for easy weakly-informed priors
vector[N] obs_sds_log = log(obs_sds) ;
real mean_of_sds = mean(obs_sds_log) ;
real sd_of_sds = sd(obs_sds_log) ;
vector[N] obs_sds_z = (obs_sds_log-mean_of_sds)/sd_of_sds ;
}
parameters {
//centrality parameters
real centrality_z_intercept;
real<lower=0> centrality_z_sd;
vector[N] centrality_noncentering_helper;
//variability parameters
real variability_z_intercept;
real<lower=0> variability_z_sd;
vector[N] variability_noncentering_helper;
}
transformed parameters{
vector[N] true_experiment_mean = (centrality_z_intercept + centrality_z_sd * centrality_noncentering_helper) * sd_of_means + mean_of_means ;
vector[N] true_experiment_log_sd = (variability_z_intercept + variability_z_sd * variability_noncentering_helper) * sd_of_sds + mean_of_sds ;
}
model {
// priors for centrality
target += std_normal_lpdf(centrality_z_intercept); //peaked at zero. Implements data-driven weakly-informed prior when data is standardized to mean=0 & sd=1
target += weibull_lpdf(centrality_z_sd | 2, 1); //peaked around 1, positive and zero-avoiding. Implements data-driven weakly-informed prior when data is standardized to mean=0 & sd=1
target += std_normal_lpdf(centrality_noncentering_helper);
// priors for variability (these need tweaking!)
target += std_normal_lpdf(variability_z_intercept); //peaked at zero. Implements data-driven weakly-informed prior when data is standardized to mean=0 & sd=1
target += weibull_lpdf(variability_z_sd | 2, 1); //peaked around 1, positive and zero-avoiding. Implements data-driven weakly-informed prior when data is standardized to mean=0 & sd=1
target += std_normal_lpdf(variability_noncentering_helper);
// likelihood
target += normal_lpdf(obs_sds_z | variability_z_intercept , variability_z_sd );
target += normal_lpdf(obs_means | true_experiment_mean, exp(true_experiment_log_sd) ./ sqrt(K) );
}
generated quantities{
//transform back to original data scales
real centrality_intercept = (centrality_z_intercept * sd_of_means) + mean_of_means ;
real<lower=0> centrality_sd = centrality_z_sd * sd_of_means ;
real variability_intercept_log = (variability_z_intercept * sd_of_sds) + mean_of_sds ;
real<lower=0> variability_sd_log = variability_z_sd * sd_of_sds ;
}