Skip to content

Instantly share code, notes, and snippets.

@mortonjt
Created May 29, 2020 19:27
Show Gist options
  • Save mortonjt/9d336566c83541e9c80f0f78c3a97750 to your computer and use it in GitHub Desktop.
Save mortonjt/9d336566c83541e9c80f0f78c3a97750 to your computer and use it in GitHub Desktop.
Show how to fit a random effects models with negative binomial
data {
int<lower=0> N; // number of samples
int<lower=0> D; // number of dimensions
int<lower=0> J; // number of subjects
int<lower=0> p; // number of covariates
real depth[N]; // sequencing depths of microbes
matrix[N, p] x; // covariate matrix
int y[N, D]; // observed microbe abundances
int<lower=1, upper=J> subj_ids[N]; // subject ids
}
parameters {
// parameters required for linear regression on the species means
matrix[p, D-1] beta; // covariates
matrix[J, D-1] alpha; // subject differentials
real<lower=0.01> disp;
}
transformed parameters {
matrix[N, D-1] lam;
matrix[N, D] lam_clr;
matrix[N, D] prob;
vector[N] z;
z = to_vector(rep_array(0, N));
lam = x * beta;
// add batch effects
// add in subject specific effects
for (n in 1:N){
lam[n] += alpha[subj_ids[n]];
}
lam_clr = append_col(z, lam);
}
model {
// setting priors ...
disp ~ inv_gamma(1., 1.);
for (i in 1:D-1){
for (j in 1:p){
beta[j, i] ~ normal(0., 5.); // uninformed prior
}
}
// generating counts
for (n in 1:N){
for (i in 1:D){
target += neg_binomial_2_log_lpmf(y[n, i] | depth[n] + lam_clr[n, i], disp);
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment