I wasn’t able to figure out how to do this in brms, but I have made a simple sim and simple model with initial pointers on the model from @Stephen_Martin, and it seems to work as intended! Thanks @Stephen_Martin !
I checked the true_proportion sd’s from the model against the observed_counts (I have sorted them by observed_count in the .csv sheet) and sure enough the sd’s in general increase as the observed_counts decrease. Nice!
Here is my code for sim and simple model:
#define true proportions and generate observed proportions from the true, using beta distribution scaled by the observed_count
n <- 100
observed_count <- rnbinom(n, mu=50, size=1) + 10
true_proportion <- runif(n, 0.0000001, 0.9999999)
mu <- true_proportion
shape1 <- mu * (observed_count + 1)
shape2 <- (1 - mu) * (observed_count + 1)
observed_proportion <- rbeta(n, shape1, shape2)
observed_proportion <- ifelse(observed_proportion==0, (observed_proportion + 0.0001), observed_proportion)
observed_proportion <- ifelse(observed_proportion==1, (observed_proportion - 0.0001), observed_proportion)
#generate outcome variable with the true proportion as the predictor
a_b <- 0.01
b_b <- 0.7
y_p <- a_b + true_proportion*b_b
y <- rbinom(n, 1, y_p)
#data
d1 <- list(N = n,
y = y,
observed_proportion = observed_proportion,
observed_count = observed_count)
#stan model
stan_code <- "
data {
int<lower=0> N;
vector[N] observed_count;
vector[N] observed_proportion;
int <lower=0, upper=1> y [N];
}
parameters {
vector<lower=0, upper=1>[N] true_proportion;
real alpha;
real beta;
}
model {
observed_proportion ~ beta_proportion(true_proportion, observed_count + 1);
y ~ bernoulli_logit(alpha + beta * true_proportion);
alpha ~ normal(0, 5);
beta ~ normal(0, 5);
}
"
#fit model
fit1 <- stan(model_code = stan_code, data = d1,
chains = 4, iter = 2000, warmup = 1000,
thin = 1, cores = 4)
print(fit1)
And the check sheet.
count vs sd check.csv (3.8 KB)
I was thinking I should put a prior on the parameter true_proportion
, but I am not sure how, since it is a vector…unfortunately, I am not as familiar with Stan code as I am brms. The things I tried gave me various errors.