Poisson regression, sparsification, and model selection
Peter Ralph
29 January 2018 – Advanced Biological Statistics
Overview
Summary
So far we looked at fitting probabilities of binary events predicted by various sorts of explanatory variables - for instance, “logistic regression”.
Grouping means together using hyperpriors induced shrinkage, allowing sharing of information between groups.
Simulation is also useful for debugging and, vaguely, checking model fit.
Next up
Next we’ll hit many of the same topics in a slightly different context: when the data are counts rather than proportions. (“How many birds are there?” instead of “How many of the birds are red?”)
We’ll also see how to use simulation to do formal model comparisons.
Count data
A hypothetical situation:
We have counts of transcript numbers,
from some individuals of different ages and past exposures to solar irradiation,
of two genotypes.
Model:
Counts are Poisson,
with mean that depends on age and exposure,
but effect of exposure depends on genotype.
- But actually, counts are overdispersed, so make the means random, and lognormally distributed.
Counts are Poisson,
with mean that depends on age and exposure,
but effect of exposure depends on genotype.
But actually, counts are overdispersed, so make the means random, and lognormally distributed.
\[\begin{aligned} Z_i &\sim \Poisson(\mu_i) \\ \end{aligned}\]
Counts are Poisson,
with mean that depends on age and exposure,
but effect of exposure depends on genotype.
But actually, counts are overdispersed, so make the means random, and lognormally distributed.
\[\begin{aligned} Z_i &\sim \Poisson(\mu_i) \\ \mu_i &= a + b \times \text{age}_i + c \times \text{exposure}_i \end{aligned}\]
Counts are Poisson,
with mean that depends on age and exposure,
but effect of exposure depends on genotype.
But actually, counts are overdispersed, so make the means random, and lognormally distributed.
\[\begin{aligned} Z_i &\sim \Poisson(\mu_i) \\ \mu_i &= \exp\left( a_{g_i} + b \times \text{age}_i + c_{g_i} \times \text{exposure}_i \right) \end{aligned}\]
Counts are Poisson,
with mean that depends on age and exposure,
but effect of exposure depends on genotype.
But actually, counts are overdispersed, so make the means random, and lognormally distributed.
\[\begin{aligned} Z_i &\sim \Poisson(\mu_i) \\ \mu_i &= \exp(W_i) \\ W_i &\sim \Normal(y_i, \sigma) \\ y_i &= a_{g_i} + b \times \text{age}_i + c_{g_i} \times \text{exposure}_i \end{aligned}\]
Counts are Poisson,
with mean that depends on age and exposure,
but effect of exposure depends on genotype.
But actually, counts are overdispersed, so make the means random, and lognormally distributed.
\[\begin{aligned} Z_i &\sim \Poisson(\mu_i) \\ \mu_i &\sim \log\Normal(y_i, \sigma) \\ y_i &= a_{g_i} + b \times \text{age}_i + c_{g_i} \times \text{exposure}_i \end{aligned}\]
Simulate these and compare.
\[\begin{aligned} Z_i &\sim \Poisson(\mu_i) \\ \mu_i &= \exp\left(a_{g_i} + b \times \text{age}_i \right.\\ &\qquad \left. {} + c_{g_i} \times \text{exposure}_i \right) \end{aligned}\]
\[\begin{aligned} Z_i &\sim \Poisson(\mu_i) \\ \mu_i &= \exp(X_i) \\ X_i &\sim \Normal(y_i, \sigma) \\ y_i &= a_{g_i} + b \times \text{age}_i \\ &\qquad {} + c_{g_i} \times \text{exposure}_i \end{aligned}\]
Simulate this.
true_params <- list(a=c(0, 0.2),
b=1/20,
c=c(1/30, -1/15),
sigma=1.0)
nsamples <- 500
data <- data.frame(genotype=sample(c(1,2), nsamples,
replace=TRUE),
age = rgamma(nsamples, 3, 0.1),
exposure = rexp(nsamples, 0.2))
data$y <- with(data, true_params$a[genotype] +
true_params$b * age +
true_params$c[genotype] * exposure)
data$mu <- exp(rnorm(nrow(data), mean=data$y,
sd=true_params$sigma))
data$counts <- rpois(nsamples, data$mu)\[\begin{aligned} Z_i &\sim \Poisson(\mu_i) \\ \mu_i &\sim \log\Normal(y_i, \sigma) \\ y_i &= \exp\left(a_{g_i} + b \times \text{age}_i \right.\\ &\qquad \left. {} + c_{g_i} \times \text{exposure}_i \right) \end{aligned}\]
Write a Stan block
Counts are Poisson,
with mean that depends on age and exposure,
but effect of exposure depends on genotype;
means are random, and lognormally distributed.
\[\begin{aligned} Z_i &\sim \Poisson(\mu_i) \\ \mu_i &\sim \log\Normal(y_i, \sigma) \\ y_i &= a_{g_i} + b \times \text{age}_i \\ &\qquad {} + c_{g_i} \times \text{exposure}_i \end{aligned}\]
data {
// what we know
}
parameters {
// what we want to find out
}
model {
// how they relate
}The result
Counts are Poisson,
with mean that depends on age and exposure,
but effect of exposure depends on genotype;
means are random, and lognormally distributed.
\[\begin{aligned} Z_i &\sim \Poisson(\mu_i) \\ \mu_i &\sim \log\Normal(y_i, \sigma) \\ y_i &= a_{g_i} + b \times \text{age}_i \\ &\qquad {} + c_{g_i} \times \text{exposure}_i \end{aligned}\]
new_block <- "data {
int N; // number of obs
vector[N] age;
int geno[N];
vector[N] expo;
int counts[N];
int ngeno;
}
parameters {
real<lower=0> mu[N]; // per-indiv 'mean'
vector[ngeno] a;
real b;
vector[ngeno] c;
real<lower=0> sigma;
}
model {
vector[N] y;
y = a[geno] + b * age + c[geno] .* expo;
// for (i in 1:N) {
// y[i] = a[geno[i]] + b * age[i] + c[geno[i]] * expo[i];
// }
counts ~ poisson(mu);
mu ~ lognormal(y, sigma);
a ~ normal(0, 10);
b ~ normal(0, 10);
c ~ normal(0, 10);
sigma ~ normal(0, 10);
}"
data_list <- list(
N = nrow(data),
age = data$age,
geno = data$genotype,
expo = data$exposure,
counts = data$counts,
ngeno = length(unique(data$genotype)))
new_fit <- stan(model_code=new_block,
data=data_list,
control=list(max_treedepth=12),
iter=100, chains=3)Stochastic minute
The Poisson distribution
If \(N \sim \Poisson(\mu)\) then \(N \ge 0\) and \[\begin{aligned} \P\{N = k\} = \frac{\mu^k}{k!} e^{-\mu} \end{aligned}\]
\(N\) is a nonnegative integer (i.e., a count)
\(\E[N] = \var[N] = \mu\)
If a machine makes widgets very fast, producing on average one broken widget per minute (and many good ones), each breaking independent of the others, then the number of broken widgets in \(\mu\) minutes is \(\Poisson(\mu)\).
If busses arrive randomly every \(\Exp(1)\) minutes, then the number of busses to arrive in \(\mu\) minutes is \(\Poisson(\mu)\).
Two models
Step back: a simple model
Forget that we know how the data were generated.
Let’s fit a standard Poisson model.
simple_block <- "
data {
int N;
int counts[N];
vector[N] age;
vector[N] exposure;
int genotype[N]; // between 0 and 1
}
parameters {
vector[2] a;
real b;
vector[2] c;
}
model {
vector[N] mu;
mu = exp(a[genotype] + b * age
+ c[genotype] .* exposure);
counts ~ poisson(mu);
a ~ normal(0, 5);
b ~ normal(0, 5);
c ~ normal(0, 5);
}"Note: scaling the data helps Stan find the right scale to move on.
## Inference for Stan model: e04ed23635a1fc5cf9233f0866bdf5cf.
## 3 chains, each with iter=1000; warmup=500; thin=1;
## post-warmup draws per chain=500, total post-warmup draws=1500.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5%
## a[1] 2.24 0.00 0.02 2.20 2.22 2.24 2.25 2.28
## a[2] 1.58 0.00 0.03 1.53 1.56 1.58 1.60 1.64
## b 0.91 0.00 0.01 0.89 0.90 0.91 0.92 0.93
## c[1] 0.13 0.00 0.01 0.10 0.12 0.13 0.14 0.16
## c[2] -0.38 0.00 0.03 -0.44 -0.40 -0.38 -0.36 -0.32
## lp__ 13445.48 0.08 1.67 13441.29 13444.65 13445.85 13446.71 13447.62
## n_eff Rhat
## a[1] 1079 1
## a[2] 1260 1
## b 1298 1
## c[1] 1500 1
## c[2] 1500 1
## lp__ 489 1
##
## Samples were drawn using NUTS(diag_e) at Mon Jan 29 22:36:19 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).Aside: what happens during “warmup”?
The usual plot (without warmup)
How’d we do?
Here are posterior distributions of the parameters, with the true values in red. 
What happened?
Goodness of fit
Let’s simulate up data under this model to check for goodness of fit.
We expect to not see a good fit. (Why?)
100 datasets from the posterior distribution
model {
vector[N] mu;
mu = exp(a[genotype]
+ b * age
+ c[genotype]
.* exposure);
counts ~ poisson(mu);The true data are overdispersed relative to our simulations
plot(data$counts[order(mu1)], ylab="counts", ylim=range(sim1), type='n')
segments(x0=seq_len(nrow(data)),
y0=rowMins(sim1)[order(mu1)],
y1=rowMaxs(sim1)[order(mu1)])
points(data$counts[order(mu1)], pch=20, col='red')
legend("topleft", pch=c(20,NA), lty=c(NA,1), legend=c("observed", "simulated range"), col=c('red', 'black'))
Including overdispersion
A random mean adds a level of randomness
full_block <- "
data {
int N;
int counts[N];
vector[N] age;
vector[N] exposure;
int genotype[N];
}
parameters {
vector[2] a;
real b;
vector[2] c;
vector<lower=0>[N] mu;
real<lower=0> sigma;
}
model {
vector[N] y;
y = a[genotype] + b * age + c[genotype] .* exposure;
counts ~ poisson(mu);
mu ~ lognormal(y, sigma);
a ~ normal(0, 5);
b ~ normal(0, 5);
c ~ normal(0, 5);
sigma ~ normal(0, 5);
}"## Inference for Stan model: cc2fed94df9d2f003676a260a7adbf8e.
## 3 chains, each with iter=1000; warmup=500; thin=1;
## post-warmup draws per chain=500, total post-warmup draws=1500.
##
## mean se_mean sd 2.5% 25% 50% 75%
## a[1] 1.45 0.00 0.08 1.29 1.39 1.44 1.50
## a[2] 1.27 0.00 0.08 1.12 1.22 1.27 1.33
## b 0.91 0.00 0.05 0.81 0.87 0.91 0.94
## c[1] 0.13 0.00 0.07 -0.01 0.08 0.12 0.17
## c[2] -0.40 0.00 0.09 -0.58 -0.46 -0.40 -0.33
## mu[1] 3.81 0.04 1.68 1.23 2.60 3.58 4.70
## mu[2] 3.00 0.04 1.62 0.82 1.84 2.68 3.80
## mu[3] 5.43 0.06 2.20 1.97 3.89 5.14 6.67
## mu[4] 12.72 0.09 3.43 6.82 10.32 12.43 14.77
## mu[5] 1.90 0.03 1.13 0.43 1.09 1.64 2.44
## mu[6] 45.07 0.17 6.69 33.06 40.68 44.60 49.57
## mu[7] 3.90 0.05 1.92 1.30 2.51 3.55 4.90
## mu[8] 1.12 0.02 0.78 0.20 0.55 0.93 1.46
## mu[9] 2.55 0.03 1.32 0.72 1.60 2.29 3.25
## mu[10] 20.51 0.11 4.39 12.82 17.31 20.30 23.26
## mu[11] 0.66 0.01 0.56 0.08 0.29 0.51 0.83
## mu[12] 0.64 0.01 0.54 0.08 0.26 0.49 0.84
## mu[13] 1.44 0.02 0.95 0.29 0.76 1.22 1.88
## mu[14] 2.03 0.03 1.18 0.48 1.18 1.76 2.62
## mu[15] 16.23 0.10 3.94 9.33 13.46 15.82 18.82
## mu[16] 642.68 0.66 25.72 595.43 625.25 641.75 659.83
## mu[17] 1.30 0.02 0.86 0.23 0.69 1.09 1.68
## mu[18] 1.48 0.02 0.94 0.31 0.80 1.27 1.96
## mu[19] 6.30 0.06 2.34 2.48 4.66 6.04 7.59
## mu[20] 1.86 0.03 1.15 0.38 1.01 1.62 2.45
## mu[21] 5.58 0.06 2.23 2.20 3.91 5.24 7.00
## mu[22] 3.11 0.04 1.52 0.96 2.03 2.84 3.95
## mu[23] 3.15 0.04 1.51 0.93 2.02 2.91 4.06
## mu[24] 38.38 0.16 6.22 27.45 33.92 37.99 42.45
## mu[25] 2.46 0.04 1.41 0.61 1.46 2.17 3.10
## mu[26] 1.32 0.02 0.90 0.25 0.68 1.10 1.76
## mu[27] 5.05 0.06 2.14 1.87 3.43 4.71 6.31
## mu[28] 4.29 0.05 1.95 1.37 2.90 4.01 5.27
## mu[29] 3.14 0.04 1.59 0.94 1.91 2.85 4.07
## mu[30] 3.40 0.04 1.65 1.09 2.20 3.09 4.20
## mu[31] 2.57 0.04 1.42 0.65 1.55 2.31 3.24
## mu[32] 1.04 0.02 0.74 0.17 0.52 0.85 1.37
## mu[33] 3.76 0.05 1.79 1.13 2.49 3.46 4.76
## mu[34] 1.17 0.02 0.79 0.22 0.58 0.99 1.51
## mu[35] 2.38 0.03 1.29 0.64 1.47 2.10 3.02
## mu[36] 2.06 0.03 1.17 0.50 1.21 1.84 2.65
## mu[37] 3.94 0.05 1.88 1.25 2.55 3.69 4.94
## mu[38] 1.75 0.03 1.07 0.38 0.97 1.52 2.29
## mu[39] 1.58 0.03 1.04 0.32 0.84 1.33 2.00
## mu[40] 1.48 0.02 0.97 0.28 0.81 1.27 1.92
## mu[41] 4.54 0.05 1.85 1.69 3.23 4.31 5.57
## mu[42] 5.58 0.05 2.12 2.33 4.03 5.25 6.84
## mu[43] 10.14 0.08 3.13 4.80 7.99 9.81 11.90
## mu[44] 21.70 0.12 4.52 13.73 18.39 21.38 24.66
## mu[45] 1.18 0.02 0.86 0.20 0.59 0.98 1.50
## mu[46] 1.08 0.02 0.74 0.17 0.55 0.90 1.41
## mu[47] 5.56 0.06 2.19 2.18 3.98 5.26 6.86
## mu[48] 4.21 0.05 1.78 1.53 2.95 3.97 5.20
## mu[49] 1.87 0.03 1.11 0.44 1.09 1.60 2.39
## mu[50] 9.36 0.08 2.95 4.74 7.23 8.99 11.09
## mu[51] 47.43 0.18 6.82 35.90 42.50 46.94 51.65
## mu[52] 1.43 0.02 0.95 0.30 0.72 1.19 1.88
## mu[53] 5.35 0.05 2.13 2.07 3.80 5.10 6.51
## mu[54] 33.75 0.14 5.39 24.30 29.81 33.58 37.25
## mu[55] 23.93 0.12 4.79 15.15 20.63 23.69 27.00
## mu[56] 84.01 0.22 8.49 68.44 78.25 83.75 89.33
## mu[57] 12.09 0.09 3.47 6.16 9.61 11.75 14.17
## mu[58] 1.49 0.03 0.97 0.28 0.81 1.27 1.94
## mu[59] 4.10 0.04 1.74 1.55 2.81 3.80 5.10
## mu[60] 0.91 0.02 0.71 0.15 0.42 0.70 1.20
## mu[61] 11.78 0.09 3.38 6.10 9.27 11.59 13.70
## mu[62] 3.41 0.04 1.64 1.05 2.22 3.15 4.27
## mu[63] 0.47 0.01 0.42 0.05 0.18 0.34 0.63
## mu[64] 1.18 0.02 0.84 0.21 0.61 0.97 1.52
## mu[65] 12.65 0.09 3.35 6.96 10.22 12.39 14.67
## mu[66] 1.21 0.02 0.81 0.22 0.64 1.02 1.63
## mu[67] 17.26 0.10 3.98 10.15 14.28 17.06 19.84
## mu[68] 2.10 0.03 1.29 0.48 1.18 1.83 2.71
## mu[69] 1.78 0.03 1.13 0.35 0.98 1.50 2.30
## mu[70] 1.48 0.03 1.00 0.29 0.77 1.24 1.95
## mu[71] 1.85 0.03 1.10 0.43 1.04 1.62 2.37
## mu[72] 46.88 0.17 6.72 35.04 42.11 46.37 51.23
## mu[73] 1.53 0.03 1.05 0.28 0.80 1.29 1.98
## mu[74] 2.84 0.04 1.53 0.66 1.72 2.58 3.66
## mu[75] 14.41 0.10 3.71 8.07 11.76 13.99 16.66
## mu[76] 3.74 0.05 1.74 1.24 2.48 3.39 4.72
## mu[77] 32.16 0.14 5.61 22.13 28.34 31.92 35.72
## mu[78] 11.76 0.09 3.31 6.24 9.32 11.47 13.64
## mu[79] 9.13 0.07 2.81 4.42 7.10 8.90 10.86
## mu[80] 1.74 0.03 1.10 0.39 0.97 1.52 2.24
## mu[81] 6.71 0.07 2.54 2.80 4.91 6.34 8.10
## mu[82] 1.85 0.03 1.08 0.42 1.07 1.63 2.34
## mu[83] 2.14 0.03 1.24 0.51 1.27 1.89 2.73
## mu[84] 6.78 0.06 2.50 2.86 4.94 6.51 8.32
## mu[85] 3.28 0.04 1.58 1.05 2.18 3.05 4.07
## mu[86] 18.33 0.11 4.29 11.11 15.21 17.93 21.02
## mu[87] 11.32 0.08 3.23 5.86 9.04 11.01 13.42
## mu[88] 3.31 0.04 1.72 0.97 2.05 2.96 4.24
## mu[89] 21.93 0.11 4.44 14.29 18.79 21.65 24.84
## mu[90] 4.09 0.05 1.90 1.38 2.68 3.75 5.10
## mu[91] 4.49 0.05 2.06 1.53 2.98 4.11 5.57
## mu[92] 5.46 0.06 2.23 2.02 3.87 5.11 6.77
## mu[93] 2.70 0.04 1.41 0.76 1.66 2.45 3.47
## mu[94] 21.54 0.12 4.48 13.63 18.29 21.26 24.33
## mu[95] 1.35 0.02 0.90 0.24 0.69 1.13 1.76
## mu[96] 195.76 0.37 14.19 169.01 186.34 195.35 204.67
## mu[97] 3.88 0.05 1.79 1.31 2.59 3.57 4.93
## mu[98] 3.16 0.04 1.60 0.88 1.99 2.89 4.00
## mu[99] 4.37 0.05 1.93 1.52 2.92 4.06 5.56
## mu[100] 0.74 0.02 0.63 0.11 0.32 0.57 0.99
## mu[101] 6.21 0.06 2.17 2.78 4.60 6.01 7.63
## mu[102] 13.58 0.10 3.81 7.13 10.90 13.22 15.79
## mu[103] 2.85 0.04 1.47 0.88 1.74 2.57 3.61
## mu[104] 29.54 0.14 5.42 20.18 25.86 29.30 32.91
## mu[105] 3.18 0.04 1.65 0.86 2.03 2.88 4.02
## mu[106] 2.01 0.03 1.23 0.43 1.12 1.75 2.59
## mu[107] 15.69 0.10 3.72 9.26 13.09 15.38 18.04
## mu[108] 1.75 0.03 1.13 0.34 0.93 1.48 2.30
## mu[109] 39.67 0.16 6.16 28.33 35.44 39.34 43.55
## mu[110] 1.70 0.03 1.11 0.34 0.91 1.44 2.24
## mu[111] 4.90 0.05 2.12 1.76 3.33 4.59 6.05
## mu[112] 4.97 0.05 1.99 1.92 3.56 4.69 6.07
## mu[113] 0.67 0.01 0.56 0.08 0.28 0.51 0.86
## mu[114] 4.24 0.05 1.96 1.43 2.79 3.92 5.35
## mu[115] 54.60 0.18 7.06 42.21 49.38 54.24 59.23
## mu[116] 2.75 0.04 1.40 0.74 1.71 2.51 3.54
## mu[117] 1.06 0.02 0.77 0.20 0.52 0.88 1.36
## mu[118] 13.82 0.10 3.73 7.52 11.06 13.53 16.16
## mu[119] 1.96 0.03 1.09 0.45 1.16 1.75 2.49
## mu[120] 2.00 0.03 1.25 0.44 1.15 1.71 2.54
## mu[121] 2.98 0.04 1.50 0.87 1.85 2.71 3.79
## mu[122] 1.27 0.02 0.95 0.21 0.60 1.02 1.68
## mu[123] 4.52 0.05 1.93 1.63 3.13 4.22 5.55
## mu[124] 2.43 0.04 1.38 0.56 1.47 2.15 3.06
## mu[125] 34.55 0.15 5.76 24.17 30.52 34.20 38.27
## mu[126] 1.35 0.02 0.93 0.26 0.69 1.15 1.76
## mu[127] 3.60 0.04 1.72 1.12 2.34 3.26 4.50
## mu[128] 7.41 0.07 2.53 3.34 5.68 7.07 8.82
## mu[129] 7.08 0.06 2.48 3.19 5.27 6.81 8.57
## mu[130] 3.23 0.04 1.60 0.94 2.02 2.94 4.11
## mu[131] 86.62 0.23 8.98 70.42 80.06 86.25 92.60
## mu[132] 11.36 0.09 3.35 5.74 8.95 11.04 13.51
## mu[133] 1.44 0.03 0.99 0.25 0.75 1.18 1.87
## mu[134] 18.26 0.11 4.18 11.14 15.32 17.95 20.85
## mu[135] 10.76 0.08 3.16 5.54 8.54 10.41 12.70
## mu[136] 2.60 0.04 1.41 0.67 1.56 2.35 3.35
## mu[137] 7.95 0.07 2.67 3.80 5.93 7.54 9.74
## mu[138] 0.87 0.02 0.70 0.11 0.39 0.68 1.15
## mu[139] 4.85 0.05 2.03 1.73 3.36 4.53 5.98
## mu[140] 3.94 0.05 1.76 1.43 2.63 3.66 4.93
## mu[141] 0.72 0.02 0.58 0.10 0.33 0.57 0.96
## mu[142] 3.95 0.05 1.85 1.27 2.59 3.66 4.97
## mu[143] 2.99 0.04 1.44 0.89 1.94 2.74 3.73
## mu[144] 0.70 0.01 0.55 0.09 0.31 0.55 0.93
## mu[145] 1.56 0.03 1.06 0.26 0.85 1.32 2.00
## mu[146] 0.98 0.02 0.76 0.15 0.45 0.78 1.28
## mu[147] 8.87 0.07 2.79 4.36 6.79 8.55 10.60
## mu[148] 17.23 0.11 4.14 10.06 14.37 16.94 19.73
## mu[149] 23.08 0.12 4.65 15.19 19.57 22.77 26.09
## mu[150] 31.48 0.14 5.42 21.83 27.54 31.15 34.84
## mu[151] 1.22 0.02 0.88 0.19 0.56 0.99 1.58
## mu[152] 0.87 0.02 0.67 0.13 0.40 0.68 1.12
## mu[153] 1.75 0.03 1.10 0.35 0.96 1.51 2.29
## mu[154] 59.18 0.19 7.54 44.86 53.97 59.08 63.71
## mu[155] 1.59 0.03 1.01 0.33 0.83 1.36 2.08
## mu[156] 1.23 0.02 0.82 0.24 0.64 1.02 1.60
## mu[157] 10.65 0.08 3.10 5.33 8.59 10.37 12.44
## mu[158] 0.89 0.02 0.75 0.13 0.40 0.67 1.13
## mu[159] 11.85 0.09 3.39 6.28 9.33 11.58 13.97
## mu[160] 3.19 0.04 1.54 0.96 2.01 2.92 4.02
## mu[161] 19.71 0.11 4.17 12.45 16.80 19.39 22.23
## mu[162] 1.97 0.03 1.15 0.41 1.13 1.73 2.56
## mu[163] 4.92 0.05 2.03 1.85 3.50 4.61 6.04
## mu[164] 1.86 0.03 1.14 0.40 1.04 1.61 2.40
## mu[165] 1.29 0.02 0.91 0.22 0.66 1.07 1.66
## mu[166] 1.59 0.03 1.02 0.34 0.86 1.36 2.10
## mu[167] 16.15 0.10 3.89 9.46 13.46 15.96 18.55
## mu[168] 3.65 0.04 1.67 1.18 2.48 3.38 4.50
## mu[169] 2.16 0.03 1.18 0.60 1.28 1.93 2.75
## mu[170] 49.22 0.18 7.09 35.90 44.43 48.95 53.40
## mu[171] 2.95 0.04 1.48 0.87 1.89 2.68 3.69
## mu[172] 1.31 0.02 0.92 0.23 0.65 1.09 1.72
## mu[173] 8.92 0.07 2.86 4.22 6.85 8.60 10.72
## mu[174] 29.94 0.14 5.43 20.39 26.04 29.57 33.60
## mu[175] 12.17 0.09 3.32 6.67 9.80 11.80 14.15
## mu[176] 28.93 0.14 5.34 19.85 25.13 28.56 32.35
## mu[177] 0.95 0.02 0.73 0.15 0.44 0.75 1.29
## mu[178] 1.06 0.02 0.79 0.16 0.50 0.86 1.42
## mu[179] 1.50 0.02 0.94 0.30 0.79 1.28 1.99
## mu[180] 0.48 0.01 0.43 0.05 0.19 0.35 0.64
## mu[181] 8.90 0.07 2.86 4.18 6.94 8.50 10.60
## mu[182] 2.07 0.03 1.25 0.45 1.17 1.81 2.67
## mu[183] 1.56 0.03 1.02 0.31 0.82 1.33 1.98
## mu[184] 14.28 0.09 3.64 8.37 11.76 13.92 16.53
## mu[185] 3.16 0.04 1.56 0.86 2.01 2.93 4.02
## mu[186] 4.83 0.05 2.08 1.84 3.28 4.46 6.03
## mu[187] 1.46 0.03 1.04 0.23 0.76 1.21 1.87
## mu[188] 121.93 0.29 11.12 101.37 114.00 121.31 129.74
## mu[189] 6.03 0.06 2.31 2.50 4.34 5.68 7.42
## mu[190] 18.66 0.11 4.33 11.29 15.71 18.23 21.27
## mu[191] 1.91 0.03 1.17 0.42 1.06 1.66 2.46
## mu[192] 18.60 0.11 4.12 11.41 15.70 18.28 21.20
## mu[193] 2.52 0.04 1.37 0.69 1.52 2.25 3.26
## mu[194] 1.63 0.03 1.01 0.35 0.90 1.39 2.15
## mu[195] 0.69 0.01 0.57 0.10 0.30 0.53 0.90
## mu[196] 3.21 0.04 1.65 1.02 1.96 2.89 4.11
## mu[197] 66.11 0.21 8.33 50.09 60.41 65.75 71.30
## mu[198] 3.24 0.04 1.59 0.99 2.07 2.98 4.09
## mu[199] 4.09 0.05 1.83 1.35 2.74 3.83 5.12
## mu[200] 1.77 0.03 1.08 0.38 0.99 1.55 2.28
## mu[201] 6.98 0.07 2.54 2.95 5.11 6.61 8.70
## mu[202] 3.61 0.04 1.68 1.15 2.39 3.33 4.50
## mu[203] 1.56 0.03 0.99 0.31 0.81 1.35 2.05
## mu[204] 3.65 0.04 1.66 1.22 2.46 3.40 4.52
## mu[205] 7.43 0.07 2.61 3.22 5.51 7.22 8.93
## mu[206] 1.06 0.02 0.74 0.18 0.54 0.87 1.35
## mu[207] 38.37 0.16 6.13 27.09 34.20 38.04 42.18
## mu[208] 1.64 0.03 1.04 0.35 0.89 1.41 2.13
## mu[209] 2.66 0.04 1.40 0.72 1.67 2.40 3.36
## mu[210] 1.29 0.02 0.93 0.22 0.64 1.06 1.66
## mu[211] 9.20 0.07 2.87 4.54 7.07 8.88 10.96
## mu[212] 5.00 0.05 2.08 1.92 3.53 4.66 6.20
## mu[213] 16.16 0.10 3.89 9.40 13.44 15.93 18.45
## mu[214] 2.85 0.04 1.43 0.86 1.79 2.58 3.65
## mu[215] 0.82 0.02 0.67 0.12 0.36 0.64 1.10
## mu[216] 1.92 0.03 1.19 0.39 1.07 1.62 2.47
## mu[217] 88.94 0.24 9.31 72.02 82.49 88.66 94.98
## mu[218] 3.49 0.04 1.70 1.02 2.27 3.23 4.37
## mu[219] 2.16 0.03 1.20 0.54 1.28 1.92 2.79
## mu[220] 1.05 0.02 0.75 0.18 0.51 0.86 1.37
## mu[221] 2.52 0.04 1.44 0.65 1.49 2.24 3.19
## mu[222] 10.49 0.08 3.20 5.35 8.14 10.08 12.60
## mu[223] 1.56 0.03 0.99 0.29 0.86 1.32 1.99
## mu[224] 5.54 0.06 2.15 2.15 4.01 5.24 6.79
## mu[225] 5.88 0.06 2.30 2.25 4.23 5.57 7.28
## mu[226] 1.18 0.02 0.83 0.20 0.57 0.98 1.54
## mu[227] 26.85 0.12 4.80 18.05 23.52 26.55 29.91
## mu[228] 16.14 0.10 3.84 9.44 13.36 15.92 18.38
## mu[229] 73.98 0.24 9.13 56.73 67.83 73.54 79.87
## mu[230] 1.16 0.02 0.78 0.21 0.62 0.97 1.50
## mu[231] 4.32 0.05 1.85 1.53 2.95 4.02 5.45
## mu[232] 19.25 0.11 4.36 11.82 16.15 18.98 22.15
## mu[233] 37.84 0.16 6.36 26.71 33.35 37.31 41.89
## mu[234] 13.41 0.09 3.60 7.34 10.84 12.97 15.52
## mu[235] 4.06 0.04 1.72 1.55 2.79 3.80 5.04
## mu[236] 2.37 0.03 1.32 0.62 1.42 2.09 3.04
## mu[237] 1.44 0.02 0.95 0.30 0.76 1.22 1.88
## mu[238] 1.56 0.03 1.04 0.27 0.82 1.31 2.03
## mu[239] 4.95 0.06 2.18 1.70 3.36 4.62 6.10
## mu[240] 9.91 0.08 3.19 4.74 7.64 9.61 11.78
## mu[241] 5.15 0.05 2.12 1.97 3.59 4.88 6.41
## mu[242] 11.74 0.09 3.31 6.35 9.38 11.37 13.69
## mu[243] 4.50 0.05 2.03 1.39 3.10 4.19 5.58
## mu[244] 4.12 0.05 1.84 1.43 2.81 3.82 5.06
## mu[245] 16.08 0.10 3.88 9.48 13.31 15.67 18.57
## mu[246] 4.43 0.05 1.90 1.59 3.03 4.16 5.57
## mu[247] 2.05 0.03 1.25 0.46 1.16 1.77 2.61
## mu[248] 0.59 0.01 0.52 0.06 0.23 0.43 0.77
## mu[249] 4.43 0.05 2.01 1.51 2.96 4.15 5.55
## mu[250] 2.18 0.03 1.22 0.58 1.31 1.94 2.83
## mu[251] 3.56 0.04 1.64 1.21 2.39 3.24 4.49
## mu[252] 9.32 0.07 2.82 4.76 7.28 9.09 10.89
## mu[253] 8.27 0.07 2.75 3.84 6.25 7.98 9.94
## mu[254] 7.76 0.07 2.69 3.49 5.74 7.48 9.49
## mu[255] 3.33 0.04 1.60 1.06 2.14 3.03 4.26
## mu[256] 1.37 0.02 0.91 0.25 0.71 1.15 1.77
## mu[257] 8.44 0.07 2.62 4.21 6.62 8.16 9.83
## mu[258] 0.63 0.01 0.54 0.08 0.26 0.47 0.81
## mu[259] 1.99 0.03 1.17 0.50 1.11 1.71 2.60
## mu[260] 1.40 0.03 0.98 0.25 0.73 1.17 1.78
## mu[261] 0.77 0.02 0.60 0.11 0.35 0.60 1.04
## mu[262] 1.96 0.03 1.19 0.42 1.05 1.70 2.61
## mu[263] 3.47 0.04 1.58 1.12 2.30 3.23 4.43
## mu[264] 1.55 0.03 1.00 0.29 0.84 1.33 2.05
## mu[265] 6.40 0.06 2.33 2.60 4.75 6.12 7.70
## mu[266] 6.50 0.06 2.48 2.76 4.61 6.21 8.09
## mu[267] 79.08 0.24 9.22 61.68 72.76 78.78 85.14
## mu[268] 5.45 0.06 2.16 2.27 3.87 5.15 6.76
## mu[269] 0.88 0.02 0.67 0.14 0.40 0.69 1.16
## mu[270] 1.92 0.03 1.17 0.44 1.09 1.67 2.45
## mu[271] 3.93 0.05 1.78 1.39 2.61 3.68 4.88
## mu[272] 16.80 0.10 4.05 9.91 13.75 16.55 19.43
## mu[273] 1.08 0.02 0.83 0.16 0.48 0.84 1.39
## mu[274] 4.34 0.05 1.83 1.53 2.99 4.10 5.42
## mu[275] 0.82 0.02 0.65 0.11 0.35 0.64 1.08
## mu[276] 3.40 0.04 1.63 1.04 2.19 3.14 4.33
## mu[277] 1.14 0.02 0.83 0.18 0.56 0.93 1.49
## mu[278] 2.76 0.04 1.41 0.83 1.72 2.47 3.55
## mu[279] 3.50 0.04 1.71 1.13 2.33 3.19 4.40
## mu[280] 5.25 0.05 2.07 2.11 3.81 4.97 6.36
## mu[281] 0.81 0.02 0.70 0.11 0.33 0.61 1.06
## mu[282] 4.94 0.05 1.99 1.88 3.47 4.64 6.16
## mu[283] 4.76 0.05 2.02 1.71 3.26 4.47 6.02
## mu[284] 1.26 0.02 0.92 0.19 0.63 1.02 1.64
## mu[285] 6.34 0.06 2.41 2.59 4.59 6.03 7.73
## mu[286] 3.66 0.04 1.74 1.18 2.36 3.36 4.70
## mu[287] 3.43 0.04 1.67 1.06 2.21 3.08 4.34
## mu[288] 15.30 0.10 3.93 8.74 12.53 14.92 17.62
## mu[289] 15.54 0.10 3.90 8.88 12.86 15.19 17.92
## mu[290] 6.15 0.06 2.29 2.52 4.54 5.89 7.49
## mu[291] 3.00 0.04 1.52 0.89 1.87 2.73 3.81
## mu[292] 3.60 0.04 1.55 1.30 2.48 3.39 4.45
## mu[293] 2.67 0.04 1.44 0.74 1.61 2.37 3.49
## mu[294] 86.55 0.24 9.31 70.55 79.83 86.02 92.46
## mu[295] 16.12 0.10 3.98 9.24 13.18 15.76 18.67
## mu[296] 2.02 0.03 1.15 0.52 1.17 1.78 2.63
## mu[297] 4.82 0.05 2.09 1.74 3.31 4.56 5.91
## mu[298] 1.82 0.03 1.16 0.35 1.01 1.56 2.34
## mu[299] 65.25 0.20 7.89 50.31 59.89 64.99 70.29
## mu[300] 1.62 0.03 1.03 0.32 0.90 1.40 2.09
## mu[301] 1.98 0.03 1.23 0.44 1.09 1.71 2.59
## mu[302] 5.31 0.06 2.35 1.75 3.64 4.89 6.74
## mu[303] 1.78 0.03 1.11 0.37 0.97 1.53 2.31
## mu[304] 1.92 0.03 1.16 0.40 1.06 1.68 2.53
## mu[305] 0.96 0.02 0.74 0.14 0.45 0.76 1.31
## mu[306] 15.69 0.10 3.82 9.06 13.12 15.43 17.94
## mu[307] 1.13 0.02 0.74 0.20 0.60 0.96 1.49
## mu[308] 5.95 0.06 2.19 2.45 4.37 5.71 7.14
## mu[309] 66.51 0.21 8.11 51.34 60.90 65.86 71.96
## mu[310] 2.35 0.03 1.32 0.59 1.38 2.12 3.00
## mu[311] 11.16 0.09 3.30 5.70 8.78 10.83 13.16
## mu[312] 1.88 0.03 1.14 0.39 1.04 1.63 2.44
## mu[313] 36.87 0.16 6.21 26.09 32.47 36.34 40.63
## mu[314] 11.64 0.09 3.33 6.15 9.21 11.28 13.60
## mu[315] 5.14 0.06 2.14 1.90 3.60 4.86 6.35
## mu[316] 2.82 0.04 1.40 0.85 1.85 2.57 3.53
## mu[317] 259.16 0.42 16.16 228.05 248.29 259.13 269.73
## mu[318] 4.96 0.05 2.08 1.78 3.36 4.67 6.20
## mu[319] 1.64 0.03 1.00 0.35 0.92 1.41 2.13
## mu[320] 1.81 0.03 1.08 0.38 1.03 1.60 2.33
## mu[321] 1.83 0.03 1.07 0.45 1.06 1.61 2.35
## mu[322] 1.13 0.02 0.89 0.16 0.49 0.90 1.46
## mu[323] 6.54 0.06 2.16 2.95 4.93 6.32 7.93
## mu[324] 4.39 0.05 2.04 1.49 2.92 4.05 5.49
## mu[325] 3.76 0.04 1.71 1.27 2.45 3.45 4.83
## mu[326] 26.79 0.13 5.11 17.43 23.26 26.52 29.98
## mu[327] 3.74 0.05 1.80 1.11 2.47 3.44 4.72
## mu[328] 1.05 0.02 0.81 0.14 0.50 0.82 1.36
## mu[329] 1.43 0.03 1.00 0.25 0.72 1.17 1.91
## mu[330] 3.29 0.04 1.60 0.93 2.17 3.00 4.15
## mu[331] 2.97 0.04 1.57 0.80 1.87 2.66 3.79
## mu[332] 13.30 0.09 3.47 7.39 10.74 13.02 15.47
## mu[333] 55.20 0.19 7.29 41.79 50.36 54.94 60.14
## mu[334] 1.95 0.03 1.14 0.44 1.14 1.70 2.50
## mu[335] 79.89 0.22 8.53 64.32 73.77 79.65 85.67
## mu[336] 1.96 0.03 1.16 0.45 1.14 1.71 2.51
## mu[337] 4.12 0.05 1.97 1.32 2.67 3.79 5.20
## mu[338] 1.52 0.03 0.99 0.32 0.85 1.29 1.92
## mu[339] 2.73 0.04 1.51 0.69 1.62 2.43 3.46
## mu[340] 1.86 0.03 1.10 0.44 1.05 1.60 2.43
## mu[341] 7.97 0.07 2.68 3.71 6.10 7.65 9.53
## mu[342] 5.04 0.05 2.12 1.77 3.44 4.80 6.29
## mu[343] 2.98 0.04 1.60 0.79 1.82 2.63 3.78
## mu[344] 1.87 0.03 1.14 0.41 1.05 1.63 2.41
## mu[345] 3.43 0.04 1.67 0.98 2.27 3.10 4.28
## mu[346] 14.10 0.09 3.63 8.01 11.49 13.82 16.24
## mu[347] 1.62 0.03 1.04 0.35 0.90 1.35 2.09
## mu[348] 4.27 0.05 1.93 1.44 2.87 3.98 5.31
## mu[349] 7.90 0.07 2.67 3.48 6.02 7.61 9.42
## mu[350] 2.61 0.04 1.41 0.74 1.62 2.34 3.31
## mu[351] 1.32 0.02 0.93 0.23 0.64 1.07 1.72
## mu[352] 6.96 0.06 2.47 2.94 5.20 6.71 8.32
## mu[353] 5.57 0.06 2.16 2.23 3.97 5.29 6.93
## mu[354] 1.40 0.02 0.93 0.25 0.73 1.18 1.83
## mu[355] 20.94 0.12 4.67 13.03 17.64 20.58 23.84
## mu[356] 5.03 0.05 2.06 1.99 3.54 4.70 6.25
## mu[357] 3.46 0.04 1.66 1.01 2.27 3.16 4.35
## mu[358] 1.13 0.02 0.78 0.22 0.57 0.93 1.49
## mu[359] 3.58 0.04 1.68 1.14 2.33 3.30 4.49
## mu[360] 0.76 0.02 0.58 0.13 0.35 0.61 1.00
## mu[361] 2.78 0.04 1.44 0.81 1.71 2.49 3.48
## mu[362] 1.35 0.02 0.92 0.26 0.69 1.13 1.72
## mu[363] 4.02 0.05 1.81 1.38 2.69 3.76 5.01
## mu[364] 4.95 0.05 1.98 1.85 3.53 4.66 6.06
## mu[365] 3.83 0.05 1.78 1.21 2.54 3.52 4.82
## mu[366] 3.12 0.04 1.57 0.93 1.95 2.81 3.99
## mu[367] 46.06 0.17 6.39 34.50 41.68 45.82 49.85
## mu[368] 1.36 0.02 0.87 0.31 0.72 1.15 1.77
## mu[369] 0.88 0.02 0.73 0.13 0.40 0.67 1.13
## mu[370] 0.95 0.02 0.73 0.14 0.45 0.75 1.23
## mu[371] 8.31 0.07 2.83 3.91 6.19 7.99 10.01
## mu[372] 1.13 0.02 0.83 0.15 0.54 0.90 1.46
## mu[373] 14.00 0.10 3.78 7.60 11.29 13.58 16.43
## mu[374] 60.14 0.20 7.77 45.93 55.05 59.89 65.13
## mu[375] 1.62 0.03 1.02 0.33 0.90 1.40 2.12
## mu[376] 3.90 0.05 1.79 1.31 2.58 3.61 4.89
## mu[377] 3.40 0.04 1.60 1.12 2.20 3.07 4.32
## mu[378] 19.92 0.12 4.53 12.28 16.75 19.57 22.73
## mu[379] 19.25 0.11 4.27 12.19 16.12 19.01 21.90
## mu[380] 0.75 0.02 0.64 0.09 0.31 0.58 0.99
## mu[381] 3.12 0.04 1.59 0.87 1.93 2.85 3.94
## mu[382] 2.54 0.04 1.43 0.63 1.53 2.23 3.31
## mu[383] 13.18 0.09 3.61 7.20 10.59 12.82 15.53
## mu[384] 3.30 0.04 1.66 0.95 2.10 2.99 4.22
## mu[385] 21.39 0.12 4.61 13.48 18.16 21.03 24.21
## mu[386] 11.42 0.08 3.11 6.28 9.16 11.08 13.26
## mu[387] 198.27 0.35 13.72 172.32 188.61 197.41 207.50
## mu[388] 13.20 0.09 3.39 7.31 10.85 13.02 15.26
## mu[389] 5.84 0.06 2.24 2.29 4.19 5.55 7.14
## mu[390] 40.71 0.16 6.26 29.33 36.54 40.29 44.79
## mu[391] 15.02 0.10 3.70 8.61 12.25 14.68 17.44
## mu[392] 6.09 0.06 2.22 2.69 4.43 5.81 7.43
## mu[393] 7.10 0.06 2.47 3.30 5.30 6.77 8.53
## mu[394] 1.24 0.02 0.81 0.25 0.64 1.04 1.63
## mu[395] 17.08 0.11 4.19 10.19 14.30 16.64 19.57
## mu[396] 8.08 0.06 2.51 4.06 6.22 7.80 9.60
## mu[397] 9.59 0.08 3.03 4.60 7.35 9.28 11.33
## mu[398] 3.22 0.04 1.65 0.88 1.99 2.96 4.15
## mu[399] 2.20 0.03 1.36 0.49 1.19 1.86 2.84
## mu[400] 5.28 0.05 2.11 2.01 3.73 4.92 6.52
## mu[401] 8.79 0.08 3.09 3.95 6.48 8.42 10.57
## mu[402] 1.00 0.02 0.74 0.16 0.48 0.80 1.31
## mu[403] 9.45 0.08 2.98 4.58 7.33 9.07 11.21
## mu[404] 1.58 0.03 0.99 0.34 0.86 1.36 2.07
## mu[405] 62.10 0.21 7.95 47.35 56.69 61.65 67.09
## mu[406] 7.70 0.07 2.66 3.42 5.83 7.35 9.28
## mu[407] 1.70 0.03 1.03 0.40 0.95 1.48 2.20
## mu[408] 11.65 0.09 3.35 6.07 9.32 11.23 13.69
## mu[409] 1.08 0.02 0.77 0.18 0.53 0.89 1.41
## mu[410] 1.38 0.02 0.94 0.27 0.70 1.15 1.78
## mu[411] 33.23 0.15 5.67 23.47 28.94 33.03 37.23
## mu[412] 36.58 0.15 5.97 25.66 32.40 36.13 40.52
## mu[413] 6.52 0.06 2.43 2.64 4.80 6.19 7.99
## mu[414] 1.16 0.02 0.82 0.22 0.57 0.95 1.55
## mu[415] 0.75 0.02 0.58 0.10 0.34 0.59 1.00
## mu[416] 36.94 0.14 5.55 26.51 33.33 36.70 40.29
## mu[417] 9.53 0.08 3.07 4.58 7.34 9.17 11.40
## mu[418] 15.97 0.10 3.97 9.01 13.23 15.77 18.32
## mu[419] 6.63 0.06 2.40 2.83 4.80 6.32 8.08
## mu[420] 1.58 0.02 0.96 0.35 0.87 1.39 2.02
## mu[421] 11.48 0.09 3.31 5.91 9.07 11.11 13.54
## mu[422] 1.95 0.03 1.12 0.46 1.17 1.72 2.49
## mu[423] 3.19 0.04 1.57 0.94 2.07 2.96 3.98
## mu[424] 19.68 0.12 4.64 11.78 16.30 19.26 22.65
## mu[425] 7.01 0.07 2.55 3.04 5.15 6.63 8.48
## mu[426] 3.08 0.04 1.60 0.77 1.90 2.82 4.01
## mu[427] 3.68 0.05 1.75 1.21 2.43 3.36 4.59
## mu[428] 1.12 0.02 0.82 0.20 0.53 0.88 1.47
## mu[429] 1.84 0.03 1.11 0.42 1.02 1.59 2.37
## mu[430] 5.95 0.06 2.40 2.34 4.17 5.55 7.34
## mu[431] 1.48 0.02 0.96 0.30 0.78 1.23 1.92
## mu[432] 5.96 0.06 2.31 2.47 4.35 5.62 7.20
## mu[433] 6.73 0.06 2.40 2.90 4.96 6.47 8.21
## mu[434] 3.24 0.04 1.65 0.92 2.00 2.95 4.14
## mu[435] 1.44 0.02 0.92 0.30 0.78 1.26 1.87
## mu[436] 1.12 0.02 0.81 0.16 0.54 0.91 1.51
## mu[437] 2.15 0.03 1.18 0.53 1.32 1.92 2.77
## mu[438] 1.36 0.02 0.94 0.25 0.71 1.15 1.74
## mu[439] 0.59 0.01 0.46 0.09 0.26 0.46 0.77
## mu[440] 3.62 0.05 1.81 1.08 2.31 3.28 4.61
## mu[441] 1.33 0.02 0.86 0.24 0.70 1.13 1.76
## mu[442] 2.38 0.03 1.32 0.61 1.44 2.14 3.00
## mu[443] 8.31 0.07 2.67 4.07 6.43 7.98 9.91
## mu[444] 0.52 0.01 0.48 0.05 0.19 0.38 0.69
## mu[445] 1.79 0.03 1.09 0.41 0.99 1.54 2.38
## mu[446] 2.16 0.03 1.33 0.44 1.17 1.84 2.84
## mu[447] 4.08 0.05 1.79 1.39 2.77 3.85 5.06
## mu[448] 1.61 0.03 0.99 0.32 0.92 1.39 2.06
## mu[449] 1.57 0.03 0.98 0.35 0.86 1.33 2.04
## mu[450] 18.79 0.11 4.10 11.82 15.78 18.42 21.39
## mu[451] 1.76 0.03 1.11 0.38 0.99 1.53 2.26
## mu[452] 3.13 0.04 1.62 0.89 1.95 2.80 3.99
## mu[453] 49.06 0.18 6.83 36.87 44.34 48.82 53.59
## mu[454] 147.76 0.31 11.83 125.41 139.77 147.44 155.28
## mu[455] 15.13 0.10 3.82 8.46 12.49 14.87 17.33
## mu[456] 2.97 0.04 1.52 0.82 1.86 2.69 3.76
## mu[457] 1.47 0.02 0.92 0.28 0.81 1.26 1.92
## mu[458] 3.44 0.04 1.74 1.03 2.16 3.14 4.37
## mu[459] 7.89 0.07 2.56 3.73 6.01 7.59 9.44
## mu[460] 1.48 0.02 0.96 0.30 0.78 1.24 1.93
## mu[461] 1.34 0.02 0.95 0.23 0.64 1.12 1.77
## mu[462] 44.02 0.17 6.56 31.76 39.77 43.60 48.12
## mu[463] 7.27 0.06 2.50 3.28 5.41 6.93 8.88
## mu[464] 0.81 0.02 0.68 0.10 0.35 0.62 1.05
## mu[465] 16.68 0.10 4.00 9.86 13.79 16.34 19.15
## mu[466] 12.92 0.09 3.50 6.99 10.41 12.69 15.04
## mu[467] 2.77 0.04 1.38 0.80 1.76 2.54 3.53
## mu[468] 0.97 0.02 0.74 0.15 0.46 0.76 1.26
## mu[469] 3.15 0.04 1.58 0.84 2.09 2.90 3.90
## mu[470] 1.43 0.02 0.95 0.27 0.73 1.20 1.86
## mu[471] 4.14 0.05 1.84 1.34 2.88 3.90 5.05
## mu[472] 12.92 0.09 3.36 6.92 10.57 12.74 14.95
## mu[473] 10.00 0.08 3.04 4.93 7.75 9.73 12.02
## mu[474] 4.48 0.05 1.91 1.63 3.10 4.20 5.52
## mu[475] 129.46 0.30 11.48 107.62 121.70 129.25 137.07
## mu[476] 10.39 0.08 3.19 5.34 8.15 9.99 12.21
## mu[477] 20.93 0.11 4.44 12.99 17.78 20.72 23.56
## mu[478] 30.70 0.15 5.69 20.59 26.76 30.33 34.10
## mu[479] 5.24 0.05 2.06 2.07 3.73 4.95 6.46
## mu[480] 1.18 0.02 0.81 0.21 0.61 0.98 1.54
## mu[481] 1.52 0.02 0.91 0.30 0.87 1.33 2.01
## mu[482] 1.94 0.03 1.16 0.43 1.11 1.69 2.48
## mu[483] 7.07 0.06 2.46 3.07 5.28 6.81 8.50
## mu[484] 6.64 0.07 2.57 2.79 4.65 6.31 8.18
## mu[485] 1.45 0.03 0.99 0.28 0.72 1.22 1.91
## mu[486] 31.65 0.14 5.38 21.78 28.02 31.37 35.08
## mu[487] 5.42 0.06 2.14 2.03 3.93 5.11 6.61
## mu[488] 1.32 0.02 0.91 0.23 0.68 1.10 1.65
## mu[489] 7.66 0.07 2.65 3.36 5.77 7.30 9.34
## mu[490] 0.91 0.02 0.73 0.13 0.42 0.70 1.15
## mu[491] 2.32 0.03 1.30 0.58 1.41 2.05 2.95
## mu[492] 2.36 0.03 1.24 0.68 1.48 2.14 2.99
## mu[493] 2.81 0.04 1.54 0.73 1.66 2.52 3.60
## mu[494] 0.94 0.02 0.69 0.15 0.45 0.76 1.26
## mu[495] 11.43 0.09 3.46 5.97 8.88 11.06 13.36
## mu[496] 0.68 0.01 0.56 0.09 0.29 0.52 0.90
## mu[497] 80.41 0.23 9.05 63.64 73.46 80.34 86.90
## mu[498] 3.90 0.05 1.77 1.33 2.57 3.66 4.85
## mu[499] 1.18 0.02 0.83 0.21 0.57 0.95 1.57
## mu[500] 27.20 0.14 5.40 17.78 23.27 26.90 30.69
## sigma 1.05 0.00 0.05 0.96 1.02 1.05 1.08
## lp__ 16677.31 0.94 17.86 16642.11 16665.10 16677.25 16689.85
## 97.5% n_eff Rhat
## a[1] 1.61 1500 1.00
## a[2] 1.43 1500 1.00
## b 1.01 1500 1.00
## c[1] 0.26 1500 1.00
## c[2] -0.21 1500 1.00
## mu[1] 7.79 1500 1.00
## mu[2] 6.76 1500 1.00
## mu[3] 10.71 1500 1.00
## mu[4] 20.24 1500 1.00
## mu[5] 4.82 1500 1.00
## mu[6] 58.81 1500 1.00
## mu[7] 8.74 1500 1.00
## mu[8] 3.14 1500 1.00
## mu[9] 5.93 1500 1.00
## mu[10] 29.71 1500 1.00
## mu[11] 2.06 1500 1.00
## mu[12] 2.09 1500 1.00
## mu[13] 3.84 1500 1.00
## mu[14] 5.05 1500 1.00
## mu[15] 24.69 1500 1.00
## mu[16] 692.29 1500 1.00
## mu[17] 3.58 1500 1.00
## mu[18] 3.72 1500 1.00
## mu[19] 11.72 1500 1.00
## mu[20] 4.71 1500 1.00
## mu[21] 10.48 1500 1.00
## mu[22] 6.69 1500 1.00
## mu[23] 6.69 1500 1.00
## mu[24] 51.11 1500 1.00
## mu[25] 6.28 1500 1.00
## mu[26] 3.75 1500 1.00
## mu[27] 9.92 1500 1.00
## mu[28] 9.20 1500 1.00
## mu[29] 7.02 1500 1.00
## mu[30] 7.50 1500 1.00
## mu[31] 5.99 1500 1.00
## mu[32] 2.94 1500 1.00
## mu[33] 7.99 1500 1.00
## mu[34] 3.25 1500 1.00
## mu[35] 5.46 1500 1.00
## mu[36] 4.92 1500 1.00
## mu[37] 8.57 1500 1.00
## mu[38] 4.46 1500 1.00
## mu[39] 4.21 1500 1.00
## mu[40] 3.84 1500 1.00
## mu[41] 9.07 1500 1.00
## mu[42] 10.52 1500 1.00
## mu[43] 17.25 1500 1.00
## mu[44] 30.99 1500 1.00
## mu[45] 3.58 1500 1.00
## mu[46] 3.03 1500 1.00
## mu[47] 10.93 1500 1.00
## mu[48] 8.37 1500 1.00
## mu[49] 4.62 1500 1.00
## mu[50] 16.05 1500 1.00
## mu[51] 61.59 1500 1.00
## mu[52] 3.65 1500 1.00
## mu[53] 10.45 1500 1.00
## mu[54] 44.65 1500 1.00
## mu[55] 34.12 1500 1.00
## mu[56] 101.70 1500 1.00
## mu[57] 20.24 1500 1.00
## mu[58] 4.00 1500 1.00
## mu[59] 8.27 1500 1.00
## mu[60] 2.68 1500 1.00
## mu[61] 19.57 1500 1.00
## mu[62] 7.42 1500 1.00
## mu[63] 1.62 1500 1.00
## mu[64] 3.35 1500 1.00
## mu[65] 20.00 1500 1.00
## mu[66] 3.16 1500 1.00
## mu[67] 25.62 1500 1.00
## mu[68] 5.40 1500 1.00
## mu[69] 4.74 1500 1.00
## mu[70] 4.12 1500 1.00
## mu[71] 4.55 1500 1.00
## mu[72] 61.76 1500 1.00
## mu[73] 4.22 1500 1.00
## mu[74] 6.77 1500 1.00
## mu[75] 22.63 1500 1.00
## mu[76] 7.74 1500 1.00
## mu[77] 44.07 1500 1.00
## mu[78] 19.51 1500 1.00
## mu[79] 15.41 1500 1.00
## mu[80] 4.40 1500 1.00
## mu[81] 12.79 1500 1.00
## mu[82] 4.59 1500 1.00
## mu[83] 5.32 1500 1.00
## mu[84] 12.32 1500 1.00
## mu[85] 7.17 1500 1.00
## mu[86] 27.99 1500 1.00
## mu[87] 18.58 1500 1.00
## mu[88] 7.59 1500 1.00
## mu[89] 31.32 1500 1.00
## mu[90] 8.88 1500 1.00
## mu[91] 9.60 1500 1.00
## mu[92] 10.80 1500 1.00
## mu[93] 6.17 1500 1.00
## mu[94] 30.81 1500 1.00
## mu[95] 3.65 1500 1.00
## mu[96] 225.92 1500 1.00
## mu[97] 7.92 1500 1.00
## mu[98] 6.90 1500 1.00
## mu[99] 8.85 1500 1.00
## mu[100] 2.45 1500 1.00
## mu[101] 10.83 1500 1.00
## mu[102] 22.22 1500 1.00
## mu[103] 6.30 1500 1.00
## mu[104] 40.91 1500 1.00
## mu[105] 7.18 1500 1.00
## mu[106] 5.12 1500 1.00
## mu[107] 23.56 1500 1.00
## mu[108] 4.65 1500 1.00
## mu[109] 52.34 1500 1.00
## mu[110] 4.68 1500 1.00
## mu[111] 9.71 1500 1.00
## mu[112] 9.59 1500 1.00
## mu[113] 2.10 1500 1.00
## mu[114] 8.73 1500 1.00
## mu[115] 69.03 1500 1.00
## mu[116] 5.92 1500 1.00
## mu[117] 3.05 1500 1.00
## mu[118] 21.46 1500 1.00
## mu[119] 4.58 1500 1.00
## mu[120] 5.04 1500 1.00
## mu[121] 6.53 1500 1.00
## mu[122] 3.56 1500 1.00
## mu[123] 8.83 1500 1.00
## mu[124] 5.98 1500 1.00
## mu[125] 46.83 1500 1.00
## mu[126] 3.63 1500 1.00
## mu[127] 7.83 1500 1.00
## mu[128] 12.99 1500 1.00
## mu[129] 12.73 1500 1.00
## mu[130] 7.20 1500 1.00
## mu[131] 104.70 1500 1.00
## mu[132] 18.90 1500 1.00
## mu[133] 4.02 1500 1.00
## mu[134] 27.30 1500 1.00
## mu[135] 17.67 1500 1.00
## mu[136] 6.02 1500 1.00
## mu[137] 13.99 1500 1.00
## mu[138] 2.68 1500 1.00
## mu[139] 9.46 1500 1.00
## mu[140] 8.01 1500 1.00
## mu[141] 2.10 1500 1.00
## mu[142] 8.41 1500 1.00
## mu[143] 6.64 1500 1.00
## mu[144] 2.19 1500 1.00
## mu[145] 4.36 1500 1.00
## mu[146] 3.17 1500 1.00
## mu[147] 15.01 1500 1.00
## mu[148] 26.11 1500 1.00
## mu[149] 33.00 1500 1.00
## mu[150] 42.95 1500 1.00
## mu[151] 3.57 1500 1.00
## mu[152] 2.66 1500 1.00
## mu[153] 4.48 1500 1.00
## mu[154] 75.35 1500 1.00
## mu[155] 4.10 1500 1.00
## mu[156] 3.40 1500 1.00
## mu[157] 17.90 1500 1.00
## mu[158] 2.84 1500 1.00
## mu[159] 19.34 1500 1.00
## mu[160] 6.80 1500 1.00
## mu[161] 28.92 1500 1.00
## mu[162] 4.79 1500 1.00
## mu[163] 9.72 1500 1.00
## mu[164] 4.76 1500 1.00
## mu[165] 3.73 1500 1.00
## mu[166] 4.12 1500 1.00
## mu[167] 24.22 1500 1.00
## mu[168] 7.85 1500 1.00
## mu[169] 5.19 1500 1.00
## mu[170] 64.16 1500 1.00
## mu[171] 6.82 1500 1.00
## mu[172] 3.73 1500 1.00
## mu[173] 15.07 1500 1.00
## mu[174] 41.16 1500 1.00
## mu[175] 19.25 1500 1.00
## mu[176] 39.88 1500 1.00
## mu[177] 2.85 1500 1.00
## mu[178] 3.10 1500 1.00
## mu[179] 3.86 1500 1.00
## mu[180] 1.68 1500 1.00
## mu[181] 15.52 1500 1.00
## mu[182] 5.25 1500 1.00
## mu[183] 4.01 1500 1.00
## mu[184] 22.40 1500 1.00
## mu[185] 7.07 1500 1.00
## mu[186] 9.54 1500 1.00
## mu[187] 4.03 1500 1.00
## mu[188] 143.73 1500 1.00
## mu[189] 11.41 1500 1.00
## mu[190] 27.73 1500 1.00
## mu[191] 4.90 1500 1.00
## mu[192] 27.76 1500 1.00
## mu[193] 6.07 1500 1.00
## mu[194] 4.22 1500 1.00
## mu[195] 2.07 1500 1.00
## mu[196] 7.28 1500 1.00
## mu[197] 82.71 1500 1.00
## mu[198] 6.95 1500 1.00
## mu[199] 8.44 1500 1.00
## mu[200] 4.58 1500 1.00
## mu[201] 12.64 1500 1.00
## mu[202] 7.66 1500 1.00
## mu[203] 4.02 1500 1.00
## mu[204] 7.59 1500 1.00
## mu[205] 13.31 1500 1.00
## mu[206] 3.14 1500 1.00
## mu[207] 50.73 1500 1.00
## mu[208] 4.24 1500 1.00
## mu[209] 6.12 1500 1.00
## mu[210] 3.77 1500 1.00
## mu[211] 15.46 1500 1.00
## mu[212] 9.96 1500 1.00
## mu[213] 24.48 1500 1.00
## mu[214] 6.11 1500 1.00
## mu[215] 2.61 1500 1.00
## mu[216] 4.94 1500 1.00
## mu[217] 108.01 1500 1.00
## mu[218] 7.43 1500 1.00
## mu[219] 5.16 1500 1.00
## mu[220] 3.07 1500 1.00
## mu[221] 6.28 1500 1.00
## mu[222] 17.46 1500 1.00
## mu[223] 4.16 1500 1.00
## mu[224] 10.52 1500 1.00
## mu[225] 11.12 1500 1.00
## mu[226] 3.26 1500 1.00
## mu[227] 37.18 1500 1.00
## mu[228] 24.56 1500 1.00
## mu[229] 92.69 1500 1.00
## mu[230] 3.18 1500 1.00
## mu[231] 8.76 1500 1.00
## mu[232] 28.46 1500 1.00
## mu[233] 51.52 1500 1.00
## mu[234] 21.02 1500 1.00
## mu[235] 7.98 1500 1.00
## mu[236] 5.75 1500 1.00
## mu[237] 3.80 1500 1.00
## mu[238] 4.36 1500 1.00
## mu[239] 10.14 1500 1.00
## mu[240] 17.38 1500 1.00
## mu[241] 9.97 1500 1.00
## mu[242] 19.22 1500 1.00
## mu[243] 9.50 1500 1.00
## mu[244] 8.24 1500 1.00
## mu[245] 24.24 1500 1.00
## mu[246] 8.74 1500 1.00
## mu[247] 5.16 1500 1.00
## mu[248] 2.01 1500 1.00
## mu[249] 9.09 1500 1.00
## mu[250] 5.07 1500 1.00
## mu[251] 7.44 1500 1.00
## mu[252] 15.80 1500 1.00
## mu[253] 14.31 1500 1.00
## mu[254] 13.74 1500 1.00
## mu[255] 6.84 1500 1.00
## mu[256] 3.72 1500 1.00
## mu[257] 14.44 1500 1.00
## mu[258] 2.03 1500 1.00
## mu[259] 5.06 1500 1.00
## mu[260] 3.81 1500 1.00
## mu[261] 2.27 1500 1.00
## mu[262] 4.84 1500 1.00
## mu[263] 7.12 1500 1.00
## mu[264] 4.02 1500 1.00
## mu[265] 11.91 1500 1.00
## mu[266] 12.12 1500 1.00
## mu[267] 98.00 1500 1.00
## mu[268] 10.52 1500 1.00
## mu[269] 2.55 1500 1.00
## mu[270] 4.91 1500 1.00
## mu[271] 8.19 1500 1.00
## mu[272] 25.20 1500 1.00
## mu[273] 3.43 1500 1.00
## mu[274] 8.49 1500 1.00
## mu[275] 2.61 1500 1.00
## mu[276] 7.31 1500 1.00
## mu[277] 3.28 1500 1.00
## mu[278] 6.09 1500 1.00
## mu[279] 7.70 1500 1.00
## mu[280] 10.13 1500 1.00
## mu[281] 2.84 1500 1.00
## mu[282] 9.27 1500 1.00
## mu[283] 9.32 1500 1.00
## mu[284] 3.67 1500 1.00
## mu[285] 11.75 1500 1.00
## mu[286] 7.81 1500 1.00
## mu[287] 7.35 1500 1.00
## mu[288] 23.91 1500 1.00
## mu[289] 23.68 1500 1.00
## mu[290] 11.52 1500 1.00
## mu[291] 6.83 1500 1.00
## mu[292] 7.09 1500 1.00
## mu[293] 6.15 1500 1.00
## mu[294] 106.50 1500 1.00
## mu[295] 24.53 1500 1.00
## mu[296] 4.90 1500 1.00
## mu[297] 9.91 1500 1.00
## mu[298] 4.98 1500 1.00
## mu[299] 81.69 1500 1.00
## mu[300] 4.15 1500 1.00
## mu[301] 5.17 1500 1.00
## mu[302] 10.82 1500 1.00
## mu[303] 4.68 1500 1.00
## mu[304] 4.72 1500 1.00
## mu[305] 2.78 1500 1.00
## mu[306] 24.13 1500 1.00
## mu[307] 2.99 1500 1.00
## mu[308] 11.06 1500 1.00
## mu[309] 82.94 1500 1.00
## mu[310] 5.51 1500 1.00
## mu[311] 18.54 1500 1.00
## mu[312] 4.62 1500 1.00
## mu[313] 50.46 1500 1.00
## mu[314] 19.08 1500 1.00
## mu[315] 10.10 1500 1.00
## mu[316] 6.21 1500 1.00
## mu[317] 291.11 1500 1.00
## mu[318] 9.69 1500 1.00
## mu[319] 4.15 1500 1.00
## mu[320] 4.50 1500 1.00
## mu[321] 4.54 1500 1.00
## mu[322] 3.54 1500 1.00
## mu[323] 11.17 1500 1.00
## mu[324] 9.36 1500 1.00
## mu[325] 7.67 1500 1.00
## mu[326] 37.63 1500 1.00
## mu[327] 7.96 1500 1.00
## mu[328] 3.24 1500 1.00
## mu[329] 4.10 1500 1.00
## mu[330] 7.16 1500 1.00
## mu[331] 6.90 1500 1.00
## mu[332] 20.81 1500 1.00
## mu[333] 70.08 1500 1.00
## mu[334] 4.81 1500 1.00
## mu[335] 97.02 1500 1.00
## mu[336] 4.90 1500 1.00
## mu[337] 8.73 1500 1.00
## mu[338] 4.17 1500 1.00
## mu[339] 6.49 1500 1.00
## mu[340] 4.49 1500 1.00
## mu[341] 14.19 1500 1.00
## mu[342] 9.78 1500 1.00
## mu[343] 6.90 1500 1.00
## mu[344] 4.68 1500 1.00
## mu[345] 7.43 1500 1.00
## mu[346] 22.05 1500 1.00
## mu[347] 4.44 1500 1.00
## mu[348] 8.92 1500 1.00
## mu[349] 13.85 1500 1.00
## mu[350] 6.00 1500 1.00
## mu[351] 3.65 1500 1.00
## mu[352] 12.70 1500 1.00
## mu[353] 10.39 1500 1.00
## mu[354] 3.77 1500 1.00
## mu[355] 31.00 1500 1.00
## mu[356] 9.88 1500 1.00
## mu[357] 7.45 1500 1.00
## mu[358] 3.14 1500 1.00
## mu[359] 7.66 1500 1.00
## mu[360] 2.32 1500 1.00
## mu[361] 6.26 1500 1.00
## mu[362] 3.66 1500 1.00
## mu[363] 8.22 1500 1.00
## mu[364] 9.91 1500 1.00
## mu[365] 7.94 1500 1.00
## mu[366] 6.92 1500 1.00
## mu[367] 59.24 1500 1.00
## mu[368] 3.52 1500 1.00
## mu[369] 2.71 1500 1.00
## mu[370] 2.88 1500 1.00
## mu[371] 14.70 1500 1.00
## mu[372] 3.36 1500 1.00
## mu[373] 22.12 1500 1.00
## mu[374] 76.40 1500 1.00
## mu[375] 4.31 1500 1.00
## mu[376] 7.95 1500 1.00
## mu[377] 7.23 1500 1.00
## mu[378] 29.91 1500 1.00
## mu[379] 28.38 1500 1.00
## mu[380] 2.50 1500 1.00
## mu[381] 7.14 1500 1.00
## mu[382] 6.29 1500 1.00
## mu[383] 20.80 1500 1.00
## mu[384] 7.20 1500 1.00
## mu[385] 31.52 1500 1.00
## mu[386] 18.53 1500 1.00
## mu[387] 226.63 1500 1.00
## mu[388] 20.71 1500 1.00
## mu[389] 10.93 1500 1.00
## mu[390] 54.21 1500 1.00
## mu[391] 22.97 1500 1.00
## mu[392] 11.33 1500 1.00
## mu[393] 12.63 1500 1.00
## mu[394] 3.27 1500 1.00
## mu[395] 26.29 1500 1.00
## mu[396] 13.68 1500 1.00
## mu[397] 16.74 1500 1.00
## mu[398] 7.10 1500 1.00
## mu[399] 5.58 1500 1.00
## mu[400] 10.15 1500 1.00
## mu[401] 15.69 1500 1.00
## mu[402] 2.86 1500 1.00
## mu[403] 16.19 1500 1.00
## mu[404] 3.98 1500 1.00
## mu[405] 79.09 1500 1.00
## mu[406] 13.89 1500 1.00
## mu[407] 4.02 1500 1.00
## mu[408] 18.73 1500 1.00
## mu[409] 3.11 1500 1.00
## mu[410] 3.85 1500 1.00
## mu[411] 45.00 1500 1.00
## mu[412] 49.19 1500 1.00
## mu[413] 12.12 1500 1.00
## mu[414] 3.20 1500 1.00
## mu[415] 2.33 1500 1.00
## mu[416] 48.35 1500 1.00
## mu[417] 16.62 1500 1.00
## mu[418] 24.59 1500 1.00
## mu[419] 12.10 1500 1.00
## mu[420] 4.16 1500 1.00
## mu[421] 18.74 1500 1.00
## mu[422] 4.87 1500 1.00
## mu[423] 6.95 1500 1.00
## mu[424] 30.30 1500 1.00
## mu[425] 12.83 1500 1.00
## mu[426] 7.07 1500 1.00
## mu[427] 7.73 1500 1.00
## mu[428] 3.24 1500 1.00
## mu[429] 4.68 1500 1.00
## mu[430] 11.66 1500 1.00
## mu[431] 3.78 1500 1.00
## mu[432] 11.30 1500 1.00
## mu[433] 12.34 1500 1.00
## mu[434] 7.29 1500 1.00
## mu[435] 3.83 1500 1.00
## mu[436] 3.19 1500 1.00
## mu[437] 5.05 1500 1.00
## mu[438] 3.95 1500 1.00
## mu[439] 1.81 1500 1.00
## mu[440] 7.95 1500 1.00
## mu[441] 3.43 1500 1.00
## mu[442] 5.70 1500 1.00
## mu[443] 14.16 1500 1.00
## mu[444] 1.84 1500 1.00
## mu[445] 4.52 1500 1.00
## mu[446] 5.48 1500 1.00
## mu[447] 8.06 1500 1.00
## mu[448] 4.15 1500 1.00
## mu[449] 4.05 1500 1.00
## mu[450] 27.60 1500 1.00
## mu[451] 4.65 1500 1.00
## mu[452] 7.21 1500 1.00
## mu[453] 63.03 1500 1.00
## mu[454] 171.70 1500 1.00
## mu[455] 23.62 1500 1.00
## mu[456] 6.73 1500 1.00
## mu[457] 3.71 1500 1.00
## mu[458] 7.52 1500 1.00
## mu[459] 13.42 1500 1.00
## mu[460] 3.94 1500 1.00
## mu[461] 3.63 1500 1.00
## mu[462] 57.39 1500 1.00
## mu[463] 12.88 1500 1.00
## mu[464] 2.60 1500 1.00
## mu[465] 25.26 1500 1.00
## mu[466] 20.70 1500 1.00
## mu[467] 6.13 1500 1.00
## mu[468] 2.99 1500 1.00
## mu[469] 6.98 1500 1.00
## mu[470] 3.94 1500 1.00
## mu[471] 8.74 1500 1.00
## mu[472] 20.33 1500 1.00
## mu[473] 16.77 1500 1.00
## mu[474] 9.03 1500 1.00
## mu[475] 152.79 1500 1.00
## mu[476] 17.47 1500 1.00
## mu[477] 31.12 1500 1.00
## mu[478] 43.32 1500 1.00
## mu[479] 9.95 1500 1.00
## mu[480] 3.27 1500 1.00
## mu[481] 3.69 1500 1.00
## mu[482] 5.05 1500 1.00
## mu[483] 12.77 1500 1.00
## mu[484] 12.63 1500 1.00
## mu[485] 4.07 1500 1.00
## mu[486] 43.12 1500 1.00
## mu[487] 10.51 1500 1.00
## mu[488] 3.92 1500 1.00
## mu[489] 13.56 1500 1.00
## mu[490] 2.76 1500 1.00
## mu[491] 5.81 1500 1.00
## mu[492] 5.38 1500 1.00
## mu[493] 6.36 1500 1.00
## mu[494] 2.73 1500 1.00
## mu[495] 19.05 1500 1.00
## mu[496] 2.13 1500 1.00
## mu[497] 98.00 1500 1.00
## mu[498] 7.92 1500 1.00
## mu[499] 3.33 1500 1.00
## mu[500] 38.58 1500 1.00
## sigma 1.14 606 1.01
## lp__ 16711.80 361 1.01
##
## Samples were drawn using NUTS(diag_e) at Mon Jan 29 22:17:03 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
How’d we do now?
We still haven’t looked to see how close the true parameters were to the inferred values. Here are posterior distributions from the full model, with the true values in red. 
Posterior predictive simulations, again
# 100 datasets:
kk <- sample.int(nrow(post2$a), 100)
sims <- lapply(kk, function (k) {
a <- post2$a[k,]
b <- post2$b[k]
c <- post2$c[k,]
sigma <- post2$sigma[k]
y <- with(list2env(scaled_data),
a[genotype] +
b * age +
c[genotype] * exposure)
mu <- exp(rnorm(length(y),
mean=y, sd=sigma))
rpois(length(mu), mu)
})
sim2 <- do.call(cbind, sims)model {
vector[N] y;
y = a[genotype] + b * age + c[genotype] .* exposure;
counts ~ poisson(mu);
mu ~ lognormal(y, sigma);
}Now we cover the true data
plot(data$counts[order(mu1)], ylab="counts", ylim=range(sim2), type='n')
segments(x0=seq_len(nrow(data)),
y0=rowMins(sim2)[order(mu1)],
y1=rowMaxs(sim2)[order(mu1)])
points(data$counts[order(mu1)], pch=20, col='red')
legend("topleft", pch=c(20,NA), lty=c(NA,1), legend=c("observed", "simulated range"), col=c('red', 'black'))
Model comparison
How to compare the two models?
Two models:
counts ~ poisson(exp(a + b * age + c * exposure))counts ~ poisson(logNormal(a + b * age + c * exposure))
We just saw some plots that showed that the true data lay outside the range of the simulated data from (1) but not (2).
That was not a formal test.
We need a statistic!
Brainstorm: how can we quantify what we just saw?
Goal: come up with a single number that quantifies how much the observed data “looks like” the posterior predictive samples.
Then, the model with a better score fits better.
Same plot, zoomed in:
Implementation
gof <- function (z, x) {
# z is a vector of observed counts
# x is a (length(z) x N) matrix of simulated data
sqrt(mean( ((z - rowMeans(x)) / rowSds(x))^2 ))
}
# TEST THIS
z0 <- rnorm(5); x0 <- matrix(rnorm(50), nrow=5)
ans0 <- gof(z0, x0)
ans1 <- 0
for (i in 1:5) {
ans1 <- ans1 + (z0[i] - mean(x0[i,]))^2 / var(x0[i,])
}
ans1 <- sqrt(ans1/length(z0))
stopifnot(abs(ans0 - ans1) < 5*.Machine$double.eps)
gof_simple <- gof(data$counts, sim1)
gof_full <- gof(data$counts, sim2)The “simple” model, that did not model overdispersion has a goodness of fit score of 4.6615927, substantially larger than the lognormal model, which had a score of 1.0362712. The score communicates how far the observed data are from the mean of the simulated data, in units of standard deviations.
But does the model fit?
To see if model 2 actually does describe the data well, we’ll find the distribution of goodness-of-fit scores based on simulation: simulate some more datasets, compute the GOF score for each, and see where the observed value lies within that distribution.
On the next slide is a histogram of 100 GOF scores from simulated data; the red line shows the observed value.
more_sim2 <- do.call(cbind, sims)
gof_scores <- sapply(1:ncol(more_sim2), function (k) {
z <- more_sim2[,k]
x <- sim2
gof(z, x)
} )
hist(gof_scores)
abline(v=gof_full, col='red', lwd=2)
legend("topright", lty=2, col='red', legend="observed value")
Stan interlude
The important program blocks
data {
// what we know: the input
// declarations only
}
parameters {
// what we want to know about:
// defines the space Stan random walks in
// declarations only
}
model {
// stuff to calculate the priors and the likelihoods
// happens every step
}The program blocks
functions {
// user-defined functions
}
data {
// what we know: the input
// declarations only
}
transformed data {
// calculations to do once, at the start
}
parameters {
// what we want to know about:
// defines the space Stan random walks in
// declarations only
}
transformed parameters {
// other things that we want the posterior distribution of
// happens every step
}
model {
// stuff to calculate the priors and the likelihoods
// happens every step
}
generated quantities {
// more things we want the posterior distribution of
// but that don't affect the random walk
}On priors
Under the hood,
z ~ poisson(mu);is equivalent to
target += poisson_lpdf(z | mu);(lpdf = log posterior density function).
So, if you don’t put a prior on something, it implicitly has a uniform prior (i.e., a flat prior).
Error messages
These are important. Pay attention to them, and fix the problems.