GLMs, overdispersion, and goodness-of-fit
Peter Ralph
14 January 2018 – Advanced Biological Statistics
Generalized Linear Models
Ingredients of a GLM:
- A probability distribution, \(Y\).
(the “family”; describes the output)
- A linear predictor, \(X \beta\).
(connects input to output)
- A link function, \(\E[Y] = h(X\beta)\).
(connects linear predictor to probability distribution)
Common choices:
Linear regression: \[ Y \sim \Normal(X \beta, \sigma) .\]
Logistic regression: \[ Y \sim \Binom(n, \logistic(X\beta)) .\]
Gamma regression: \[ Y \sim \Gam(\text{scale}=\exp(X\beta)/k, \text{shape}=k) .\]
Logistic regression, in R
glm(y ~ x, family=binomial)family {stats}
Family objects provide a convenient way to specify the details of the models used by functions such as glm.
binomial(link = "logit")
gaussian(link = "identity")
Gamma(link = "inverse")
inverse.gaussian(link = "1/mu^2")
poisson(link = "log")
Arguments
link: a specification for the model link function. This can be a name/expression, a literal character string, a length-one character vector, or an object of class "link-glm" (such as generated by make.link) provided it is not specified via one of the standard names given next.
The gaussian family accepts the links (as names) identity, log and inverse; the binomial family the links logit, probit, cauchit, (corresponding to logistic, normal and Cauchy CDFs respectively) log and cloglog (complementary log-log); the Gamma family the links inverse, identity and log; the poisson family the links log, identity, and sqrt; and the inverse.gaussian family the links 1/mu^2, inverse, identity and log.Unpacking logistic regression
Simulate data
100 trials, where probability of success depends on \(x\):
alpha <- -7; beta <- 1.2
x <- runif(100, 0, 10)
y <- alpha + beta * x
p <- 1 / (1 + exp(-y))
z <- rbinom(100, size=1, prob=p)
plot(z ~ x)
curve(1/(1+exp(-(-7 + 1.2 *x))), 0, 10, col='red', add=TRUE)
glm()
##
## Call:
## glm(formula = zz ~ x, family = "binomial")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.79584 -0.14355 0.01418 0.07181 2.58115
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.7512 2.4701 -3.948 7.89e-05 ***
## x 1.9531 0.4849 4.028 5.62e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 138.269 on 99 degrees of freedom
## Residual deviance: 28.807 on 98 degrees of freedom
## AIC: 32.807
##
## Number of Fisher Scoring iterations: 7Stan
\[\begin{aligned} Z &\sim \Binom(1, P) \\ P &= \logistic(\alpha + \beta X) \end{aligned}\]
Stan: fit the model
## $summary
## mean se_mean sd 2.5% 25% 50%
## alpha -11.24967 0.10744784 2.6860590 -17.262303 -12.906310 -10.916415
## beta 2.24945 0.02105587 0.5264088 1.383928 1.870049 2.190549
## lp__ -15.43677 0.03147576 0.9731003 -18.089596 -15.870265 -15.147178
## 75% 97.5% n_eff Rhat
## alpha -9.299297 -6.835742 624.9362 1.003205
## beta 2.571353 3.415786 625.0289 1.003971
## lp__ -14.714485 -14.426820 955.7903 1.003256
##
## $c_summary
## , , chains = chain:1
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## alpha -11.157536 2.5408466 -16.865943 -12.797324 -10.970090 -9.227163
## beta 2.228253 0.5085301 1.386013 1.853273 2.190844 2.543871
## lp__ -15.452112 0.9889446 -18.098891 -15.900315 -15.171180 -14.699103
## stats
## parameter 97.5%
## alpha -6.814512
## beta 3.356737
## lp__ -14.429367
##
## , , chains = chain:2
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## alpha -11.463242 2.7812442 -17.865908 -13.021270 -11.042255 -9.489456
## beta 2.297488 0.5393182 1.420122 1.902149 2.231015 2.622173
## lp__ -15.415057 0.9809461 -18.181455 -15.834569 -15.146324 -14.674684
## stats
## parameter 97.5%
## alpha -6.920368
## beta 3.513544
## lp__ -14.417750
##
## , , chains = chain:3
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## alpha -11.304196 2.7955471 -17.575864 -13.184845 -10.89553 -9.256249
## beta 2.259235 0.5482089 1.298729 1.896051 2.20113 2.616746
## lp__ -15.503305 1.0125950 -18.171424 -15.971176 -15.19399 -14.739261
## stats
## parameter 97.5%
## alpha -6.510975
## beta 3.466496
## lp__ -14.432854
##
## , , chains = chain:4
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## alpha -11.073698 2.6051648 -16.940830 -12.582058 -10.759517 -9.194292
## beta 2.212824 0.5050335 1.422149 1.843107 2.151902 2.508836
## lp__ -15.376602 0.9034170 -17.754912 -15.735517 -15.096816 -14.728251
## stats
## parameter 97.5%
## alpha -7.147461
## beta 3.297399
## lp__ -14.432523Stan: posterior distributions
samples <- extract(fit)
layout(t(1:2))
hist(samples$alpha, main=expression(alpha))
abline(v=alpha, col='red'); abline(v=coef(glm_fit)[1], col='green')
hist(samples$beta, main=expression(beta))
abline(v=beta, col='red'); abline(v=coef(glm_fit)[2], col='green')
legend("topright", lty=1, col=c('red', 'green'), legend=c('truth', 'glm'))
Identify the GLM
data {
int N;
vector[N] X;
vector[N] Y;
vector<lower=0> Z[N];
}
parameters {
real beta[2];
}
model {
Z ~ gamma(1, exp(- beta[1] * X - beta[2] * Y);
}What is
the probability distribution? (describes the output)
the linear predictor? (connects input to output)
the link function? (connects linear predictor to probability distribution)
Then, simulate from it.
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)\).
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.
Counts are Poisson,
with mean that depends on age and exposure,
but effect of exposure depends on genotype.
\[\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.
\[\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.
\[\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}\]
Poisson regression, in practice
The data
Write a Stan block
Counts are Poisson,
with mean that depends on age and exposure,
but effect of exposure depends on genotype.
\[\begin{aligned} Z_i &\sim \Poisson(y_i) \\ y_i &= \exp(a_{g_i} + b \times \text{age}_i \\ &\qquad {} + c_{g_i} \times \text{exposure}_i ) \end{aligned}\]
poisson_model <- stan_model(model_code="
data {
int N; // number of data points
vector[N] age;
vector[N] exposure;
int counts[N];
int genotype[N];
int ngenotypes;
}
parameters {
vector[ngenotypes] a; // intercepts
real b; // slope for age
vector[ngenotypes] c; // slopes for exposure
}
model {
vector[N] mu; // mean of the poissons
mu = exp(a[genotype] + b * age + c[genotype] .* exposure);
counts ~ poisson(mu);
a ~ normal(0, 100);
b ~ normal(0, 10);
c ~ normal(0, 20);
}")The result
Note: scaling the data helps Stan find the right scale to move on.
## Inference for Stan model: simple_poisson.
## 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.06 0.00 0.02 2.02 2.05 2.06 2.07 2.10
## a[2] 1.90 0.00 0.02 1.86 1.89 1.90 1.92 1.95
## b 0.84 0.00 0.01 0.82 0.84 0.84 0.85 0.87
## c[1] 0.26 0.00 0.02 0.23 0.25 0.26 0.27 0.29
## c[2] -0.14 0.00 0.02 -0.18 -0.16 -0.14 -0.13 -0.10
## lp__ 10937.60 0.06 1.64 10933.70 10936.84 10937.90 10938.81 10939.70
## n_eff Rhat
## a[1] 1218 1
## a[2] 1260 1
## b 1202 1
## c[1] 1487 1
## c[2] 1668 1
## lp__ 727 1
##
## Samples were drawn using NUTS(diag_e) at Sat Jan 5 12:59:06 2019.
## 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: a look at “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
Posterior predictive simulations
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[geno]
+ b * age
+ c[geno]
.* expo);
counts ~ poisson(mu);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'))
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.
Run code in the console. Rstudio often hides useful information.
Parameterization matters
Including overdispersion
How can we include overdispersion?
Counts are Poisson,
with mean that depends on age and exposure,
but effect of exposure depends on genotype.
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.
Actually, counts are overdispersed, so the means are 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.
Actually, counts are overdispersed, so the 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 + c_{g_i} \times \text{exposure}_i \end{aligned}\]
Your turn: add overdispersion
\[\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}\]
You can download the data here.
overdispersed_model <- stan_model(model_code="
data {
int N; // number of data points
vector[N] age;
vector[N] exposure;
int counts[N];
int genotype[N];
int ngenotypes;
}
parameters {
vector[ngenotypes] a; // intercepts
real b; // slope for age
vector[ngenotypes] c; // slopes for exposure
real<lower=0> sigma; // SD on lognormal
vector[N] mu; // mean of the poissons
}
model {
vector[N] y; // mean of the lognormals
y = a[genotype] + b * age + c[genotype] .* exposure;
mu ~ lognormal(y, sigma);
counts ~ poisson(mu);
a ~ normal(0, 100);
b ~ normal(0, 10);
c ~ normal(0, 20);
sigma ~ normal(0, 10);
}")## Inference for Stan model: lognormal_poisson.
## 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%
## mu[1] 0.96 0.02 0.72 0.16 0.47 0.75 1.23
## mu[2] 9.10 0.06 2.85 4.51 7.03 8.75 10.88
## mu[3] 26.22 0.10 5.18 17.09 22.44 25.86 29.44
## mu[4] 2.11 0.02 1.16 0.56 1.26 1.89 2.68
## mu[5] 12.39 0.06 3.25 6.83 10.06 12.18 14.32
## mu[6] 3.32 0.03 1.54 1.06 2.25 3.04 4.16
## mu[7] 2.53 0.02 1.31 0.74 1.60 2.26 3.18
## mu[8] 3.76 0.03 1.77 1.20 2.49 3.46 4.68
## mu[9] 0.93 0.01 0.68 0.17 0.46 0.76 1.20
## mu[10] 1.35 0.02 0.84 0.27 0.74 1.15 1.75
## mu[11] 2.56 0.03 1.38 0.73 1.54 2.31 3.31
## mu[12] 10.87 0.06 3.05 5.67 8.74 10.56 12.70
## mu[13] 3.25 0.03 1.57 0.98 2.18 2.96 4.00
## mu[14] 11.57 0.06 3.29 6.21 9.36 11.21 13.48
## mu[15] 0.51 0.01 0.47 0.07 0.20 0.37 0.66
## mu[16] 2.96 0.03 1.56 0.84 1.78 2.64 3.79
## mu[17] 0.86 0.01 0.64 0.16 0.40 0.70 1.11
## mu[18] 67.14 0.14 8.04 52.32 61.60 66.87 72.58
## mu[19] 1.49 0.02 0.94 0.34 0.83 1.27 1.90
## mu[20] 1.74 0.02 1.07 0.43 0.97 1.49 2.21
## mu[21] 2.30 0.03 1.33 0.53 1.33 2.03 2.99
## mu[22] 3.81 0.04 1.81 1.33 2.55 3.46 4.68
## mu[23] 7.40 0.05 2.64 3.13 5.45 7.10 8.97
## mu[24] 28.42 0.11 5.24 19.44 24.81 28.00 31.81
## mu[25] 2.57 0.02 1.32 0.74 1.59 2.32 3.30
## mu[26] 3.23 0.03 1.58 0.98 2.09 2.95 4.11
## mu[27] 0.91 0.01 0.67 0.17 0.46 0.74 1.17
## mu[28] 1.41 0.02 0.94 0.30 0.74 1.18 1.86
## mu[29] 1.69 0.02 1.01 0.37 0.95 1.49 2.19
## mu[30] 2.29 0.02 1.26 0.58 1.34 2.07 3.01
## mu[31] 8.93 0.05 2.91 4.30 6.81 8.55 10.56
## mu[32] 65.07 0.14 8.14 50.57 59.23 64.70 70.35
## mu[33] 105.66 0.19 10.31 86.37 98.51 105.62 112.76
## mu[34] 3.74 0.04 1.72 1.22 2.52 3.45 4.65
## mu[35] 23.71 0.09 4.88 15.12 20.40 23.24 26.69
## mu[36] 3.42 0.03 1.63 1.14 2.22 3.12 4.30
## mu[37] 7.12 0.04 2.40 3.12 5.37 6.84 8.76
## mu[38] 7.39 0.05 2.71 3.17 5.39 7.07 8.91
## mu[39] 2.93 0.03 1.44 0.88 1.88 2.67 3.72
## mu[40] 16.40 0.07 4.22 9.38 13.34 16.02 19.16
## mu[41] 2.11 0.02 1.14 0.56 1.28 1.89 2.66
## mu[42] 9.98 0.06 3.03 4.89 7.84 9.64 11.75
## mu[43] 4.67 0.03 1.92 1.76 3.32 4.37 5.78
## mu[44] 75.26 0.13 8.56 59.67 69.00 75.07 81.19
## mu[45] 1.02 0.01 0.71 0.18 0.52 0.84 1.30
## mu[46] 0.88 0.01 0.62 0.17 0.44 0.72 1.17
## mu[47] 5.96 0.04 2.26 2.44 4.38 5.70 7.13
## mu[48] 1.23 0.02 0.80 0.24 0.67 1.06 1.58
## mu[49] 1.61 0.02 1.04 0.34 0.86 1.38 2.07
## mu[50] 3.82 0.03 1.70 1.30 2.55 3.51 4.86
## mu[51] 3.28 0.03 1.62 1.01 2.14 2.97 4.16
## mu[52] 1.48 0.02 0.93 0.34 0.82 1.27 1.90
## mu[53] 131.28 0.22 11.13 110.32 123.66 131.10 138.17
## mu[54] 17.98 0.07 4.27 10.75 14.95 17.73 20.56
## mu[55] 68.78 0.13 8.06 53.75 63.03 68.55 74.22
## mu[56] 1.34 0.02 0.84 0.30 0.75 1.14 1.71
## mu[57] 2.90 0.03 1.54 0.83 1.78 2.64 3.70
## mu[58] 1.28 0.02 0.80 0.27 0.69 1.10 1.68
## mu[59] 0.62 0.01 0.48 0.10 0.28 0.50 0.82
## mu[60] 6.15 0.05 2.39 2.46 4.43 5.78 7.48
## mu[61] 59.15 0.15 7.69 44.69 53.77 58.71 64.24
## mu[62] 2.21 0.02 1.25 0.56 1.31 1.93 2.81
## mu[63] 1.68 0.02 1.01 0.39 0.92 1.46 2.21
## mu[64] 30.04 0.09 5.22 20.40 26.37 29.73 33.54
## mu[65] 1.66 0.02 1.01 0.40 0.92 1.42 2.19
## mu[66] 1.63 0.02 1.02 0.34 0.91 1.38 2.12
## mu[67] 17.82 0.07 4.03 10.93 14.87 17.52 20.30
## mu[68] 2.04 0.02 1.21 0.53 1.18 1.82 2.61
## mu[69] 19.25 0.07 4.13 11.99 16.42 18.96 21.72
## mu[70] 2.58 0.03 1.45 0.66 1.55 2.26 3.32
## mu[71] 7.78 0.05 2.82 3.50 5.67 7.35 9.58
## mu[72] 0.95 0.01 0.68 0.17 0.45 0.77 1.25
## mu[73] 4.08 0.03 1.77 1.49 2.80 3.77 5.04
## mu[74] 9.08 0.06 2.85 4.48 7.03 8.72 10.90
## mu[75] 9.69 0.05 2.84 4.98 7.72 9.39 11.50
## mu[76] 5.46 0.04 2.09 2.12 3.91 5.17 6.75
## mu[77] 22.56 0.08 4.47 14.71 19.38 22.36 25.53
## mu[78] 2.36 0.03 1.30 0.65 1.36 2.09 3.12
## mu[79] 2.91 0.03 1.50 0.83 1.81 2.63 3.71
## mu[80] 36.38 0.10 6.02 25.65 32.12 36.19 40.32
## mu[81] 2.33 0.02 1.26 0.65 1.41 2.10 2.96
## mu[82] 5.97 0.04 2.25 2.52 4.33 5.66 7.37
## mu[83] 7.37 0.04 2.55 3.39 5.49 7.01 8.93
## mu[84] 24.43 0.09 4.66 16.20 21.14 24.17 27.27
## mu[85] 10.94 0.05 3.32 5.67 8.46 10.61 13.14
## mu[86] 2.10 0.02 1.15 0.55 1.27 1.85 2.66
## mu[87] 0.96 0.02 0.69 0.18 0.47 0.78 1.25
## mu[88] 4.07 0.03 1.74 1.48 2.81 3.80 4.98
## mu[89] 3.12 0.03 1.56 0.97 1.92 2.85 3.93
## mu[90] 2.20 0.02 1.19 0.62 1.33 1.97 2.84
## mu[91] 2.37 0.03 1.23 0.70 1.46 2.15 3.02
## mu[92] 2.21 0.02 1.26 0.58 1.32 1.94 2.84
## mu[93] 5.34 0.04 2.05 2.20 3.84 5.07 6.50
## mu[94] 20.01 0.07 4.42 12.49 16.79 19.66 22.64
## mu[95] 3.24 0.03 1.55 1.05 2.11 2.96 4.03
## mu[96] 9.35 0.05 2.80 4.71 7.34 9.00 10.90
## mu[97] 3.67 0.03 1.70 1.27 2.47 3.37 4.52
## mu[98] 2.21 0.02 1.23 0.59 1.32 1.96 2.80
## mu[99] 1.75 0.02 1.05 0.41 1.00 1.51 2.23
## mu[100] 39.75 0.11 5.87 29.32 35.83 39.37 43.27
## mu[101] 7.32 0.05 2.55 3.14 5.53 7.03 8.76
## mu[102] 26.36 0.10 5.14 17.20 22.90 26.01 29.54
## mu[103] 2.82 0.03 1.48 0.77 1.77 2.51 3.58
## mu[104] 1.97 0.02 1.11 0.52 1.20 1.75 2.50
## mu[105] 1.30 0.02 0.88 0.24 0.68 1.08 1.68
## mu[106] 1.93 0.02 1.12 0.45 1.14 1.72 2.44
## mu[107] 8.81 0.04 2.80 4.15 6.75 8.50 10.72
## mu[108] 16.86 0.07 4.04 9.98 13.95 16.51 19.39
## mu[109] 8.99 0.05 2.79 4.47 7.05 8.66 10.63
## mu[110] 1.96 0.03 1.14 0.50 1.16 1.73 2.46
## mu[111] 2.57 0.03 1.42 0.63 1.58 2.33 3.20
## mu[112] 2.90 0.03 1.48 0.89 1.79 2.64 3.65
## mu[113] 14.08 0.07 3.71 7.99 11.46 13.70 16.28
## mu[114] 7.50 0.04 2.66 3.26 5.58 7.19 9.12
## mu[115] 5.31 0.04 2.04 2.10 3.83 5.06 6.44
## mu[116] 1.79 0.02 1.07 0.39 1.02 1.55 2.29
## mu[117] 14.41 0.07 3.71 8.11 11.70 14.07 16.94
## mu[118] 1.44 0.02 0.98 0.28 0.77 1.21 1.81
## mu[119] 8.21 0.05 2.72 3.78 6.28 8.00 9.72
## mu[120] 20.06 0.08 4.50 12.14 16.95 19.69 22.88
## mu[121] 2.81 0.03 1.36 0.90 1.79 2.56 3.59
## mu[122] 15.18 0.06 3.66 8.91 12.72 14.81 17.30
## mu[123] 2.16 0.02 1.23 0.57 1.25 1.91 2.73
## mu[124] 1.15 0.02 0.80 0.23 0.59 0.94 1.47
## mu[125] 2.17 0.02 1.15 0.59 1.30 1.94 2.82
## mu[126] 3.39 0.03 1.68 1.12 2.17 3.09 4.29
## mu[127] 22.16 0.09 4.67 14.18 18.84 21.64 25.10
## mu[128] 2.06 0.02 1.13 0.53 1.22 1.80 2.63
## mu[129] 5.42 0.04 2.16 2.15 3.80 5.12 6.70
## mu[130] 8.93 0.05 2.66 4.50 7.00 8.69 10.50
## mu[131] 1.68 0.02 1.04 0.38 0.94 1.44 2.18
## mu[132] 6.66 0.05 2.49 2.79 4.85 6.32 8.11
## mu[133] 1.14 0.02 0.80 0.21 0.56 0.94 1.49
## mu[134] 27.17 0.09 5.40 17.96 23.22 26.75 30.74
## mu[135] 13.39 0.06 3.44 7.55 10.90 13.12 15.37
## mu[136] 5.08 0.04 1.93 2.13 3.68 4.80 6.15
## mu[137] 8.30 0.05 2.75 4.04 6.29 7.97 10.05
## mu[138] 9.29 0.05 2.81 4.78 7.23 8.93 11.09
## mu[139] 4.80 0.03 1.94 1.84 3.39 4.54 5.89
## mu[140] 7.76 0.05 2.67 3.47 5.84 7.54 9.36
## mu[141] 2.17 0.02 1.20 0.57 1.28 1.93 2.77
## mu[142] 0.92 0.01 0.67 0.15 0.46 0.76 1.19
## mu[143] 27.86 0.09 5.03 19.07 24.18 27.46 31.03
## mu[144] 2.65 0.03 1.42 0.71 1.61 2.43 3.37
## mu[145] 18.58 0.08 4.08 11.28 15.73 18.27 21.15
## mu[146] 2.27 0.03 1.28 0.58 1.34 2.01 2.91
## mu[147] 2.87 0.03 1.44 0.81 1.82 2.64 3.61
## mu[148] 1.30 0.02 0.82 0.28 0.70 1.13 1.71
## mu[149] 5.89 0.04 2.22 2.36 4.24 5.66 7.21
## mu[150] 5.09 0.04 1.94 2.01 3.71 4.80 6.21
## mu[151] 4.75 0.04 2.05 1.76 3.18 4.46 5.94
## mu[152] 3.81 0.03 1.68 1.36 2.59 3.56 4.70
## mu[153] 2.23 0.03 1.27 0.54 1.31 1.97 2.84
## mu[154] 2.89 0.03 1.45 0.84 1.83 2.62 3.63
## mu[155] 6.73 0.04 2.34 3.02 5.02 6.47 8.09
## mu[156] 0.77 0.01 0.54 0.15 0.38 0.62 1.00
## mu[157] 23.07 0.09 4.71 14.42 19.74 22.77 26.01
## mu[158] 45.53 0.12 7.08 32.53 40.80 45.03 49.83
## mu[159] 0.96 0.02 0.71 0.15 0.49 0.80 1.22
## mu[160] 3.63 0.03 1.72 1.17 2.40 3.32 4.50
## mu[161] 15.99 0.07 3.94 9.19 13.24 15.65 18.47
## mu[162] 2.03 0.02 1.16 0.47 1.17 1.79 2.65
## mu[163] 107.59 0.16 10.63 87.13 100.34 107.25 114.64
## mu[164] 2.53 0.03 1.40 0.63 1.51 2.20 3.27
## mu[165] 11.25 0.05 3.22 5.94 9.01 10.98 13.33
## mu[166] 0.79 0.01 0.55 0.16 0.41 0.66 1.01
## mu[167] 11.81 0.06 3.29 6.35 9.41 11.45 13.69
## mu[168] 9.43 0.05 2.84 4.86 7.26 9.17 11.20
## mu[169] 3.98 0.03 1.82 1.33 2.70 3.67 4.95
## mu[170] 6.89 0.05 2.53 2.95 5.10 6.56 8.36
## mu[171] 2.13 0.02 1.19 0.55 1.28 1.90 2.76
## mu[172] 2.23 0.02 1.24 0.59 1.33 1.99 2.83
## mu[173] 11.31 0.07 3.26 5.87 8.97 10.98 13.38
## mu[174] 3.06 0.03 1.53 0.90 1.95 2.76 3.91
## mu[175] 100.98 0.20 10.15 81.60 94.21 100.45 107.14
## mu[176] 7.50 0.04 2.50 3.48 5.65 7.17 9.09
## mu[177] 8.39 0.05 2.77 3.86 6.41 8.02 10.01
## mu[178] 3.09 0.03 1.54 0.93 1.94 2.81 3.89
## mu[179] 7.83 0.04 2.48 3.72 6.00 7.59 9.41
## mu[180] 2.29 0.02 1.36 0.51 1.31 1.98 2.98
## mu[181] 13.01 0.06 3.44 7.40 10.56 12.73 15.09
## mu[182] 1.33 0.02 0.85 0.28 0.74 1.12 1.71
## mu[183] 1.77 0.02 1.13 0.35 0.98 1.55 2.29
## mu[184] 12.90 0.06 3.56 7.03 10.38 12.51 15.21
## mu[185] 2.27 0.02 1.21 0.63 1.43 2.06 2.86
## mu[186] 5.24 0.04 2.18 1.98 3.63 4.93 6.53
## mu[187] 1.81 0.02 1.05 0.40 1.06 1.58 2.38
## mu[188] 9.24 0.05 2.87 4.36 7.25 8.97 10.91
## mu[189] 27.48 0.09 5.09 18.37 23.84 27.21 30.87
## mu[190] 1.60 0.02 0.98 0.36 0.89 1.37 2.08
## mu[191] 8.02 0.04 2.67 3.78 6.15 7.68 9.56
## mu[192] 6.11 0.04 2.21 2.68 4.47 5.80 7.43
## mu[193] 0.86 0.01 0.65 0.13 0.40 0.69 1.13
## mu[194] 3.82 0.03 1.68 1.25 2.58 3.56 4.81
## mu[195] 10.83 0.05 3.05 5.81 8.54 10.49 12.89
## mu[196] 3.78 0.04 1.77 1.22 2.48 3.48 4.72
## mu[197] 2.19 0.02 1.26 0.53 1.30 1.94 2.79
## mu[198] 6.51 0.04 2.28 2.89 4.85 6.23 7.86
## mu[199] 1.64 0.02 1.00 0.38 0.93 1.41 2.09
## mu[200] 1.05 0.02 0.76 0.18 0.50 0.85 1.43
## mu[201] 17.49 0.07 3.93 10.65 14.72 17.22 19.90
## mu[202] 2.96 0.03 1.51 0.87 1.87 2.69 3.75
## mu[203] 5.02 0.04 2.03 1.96 3.57 4.74 6.21
## mu[204] 1.13 0.02 0.80 0.19 0.55 0.92 1.50
## mu[205] 4.78 0.04 1.96 1.74 3.38 4.45 5.94
## mu[206] 0.70 0.01 0.54 0.10 0.33 0.56 0.93
## mu[207] 2.64 0.03 1.34 0.77 1.65 2.42 3.37
## mu[208] 2.75 0.03 1.41 0.79 1.75 2.49 3.47
## mu[209] 1.71 0.02 1.00 0.41 0.97 1.52 2.22
## mu[210] 4.79 0.04 1.91 2.02 3.36 4.53 5.88
## mu[211] 18.32 0.07 4.03 11.48 15.31 17.96 21.03
## mu[212] 5.38 0.04 2.18 2.12 3.68 5.03 6.72
## mu[213] 27.56 0.10 5.14 18.51 23.91 27.17 30.91
## mu[214] 4.93 0.03 2.04 1.84 3.45 4.61 6.08
## mu[215] 1.72 0.02 1.02 0.38 0.97 1.49 2.24
## mu[216] 0.90 0.02 0.71 0.13 0.41 0.71 1.15
## mu[217] 17.26 0.07 4.18 10.53 14.22 16.84 19.90
## mu[218] 28.34 0.09 5.30 18.63 24.60 27.94 31.77
## mu[219] 1.74 0.02 1.01 0.43 1.01 1.53 2.26
## mu[220] 3.73 0.03 1.73 1.25 2.47 3.44 4.58
## mu[221] 1.92 0.02 1.11 0.47 1.09 1.68 2.47
## mu[222] 46.64 0.11 6.95 34.11 41.93 46.40 51.23
## mu[223] 5.63 0.04 2.20 2.18 4.00 5.36 6.86
## mu[224] 5.15 0.04 2.08 2.07 3.65 4.84 6.36
## mu[225] 5.76 0.03 2.16 2.41 4.13 5.45 7.08
## mu[226] 39.90 0.11 5.92 29.01 35.68 39.72 43.76
## mu[227] 5.04 0.04 2.02 2.02 3.53 4.76 6.21
## mu[228] 3.41 0.03 1.59 1.16 2.20 3.11 4.33
## mu[229] 53.09 0.13 6.99 40.27 48.12 52.71 57.98
## mu[230] 2.04 0.02 1.18 0.51 1.18 1.79 2.61
## mu[231] 103.12 0.17 9.89 84.88 96.30 102.53 109.58
## mu[232] 7.20 0.04 2.56 3.20 5.34 6.88 8.75
## mu[233] 16.05 0.07 3.82 9.60 13.25 15.81 18.55
## mu[234] 79.09 0.17 8.99 62.10 73.08 78.63 84.78
## mu[235] 6.17 0.04 2.30 2.47 4.53 5.92 7.41
## mu[236] 5.53 0.04 2.12 2.28 3.98 5.26 6.68
## mu[237] 12.88 0.06 3.48 7.07 10.29 12.52 15.12
## mu[238] 18.61 0.07 4.22 11.14 15.60 18.39 21.19
## mu[239] 2.66 0.03 1.37 0.83 1.63 2.38 3.36
## mu[240] 74.36 0.14 8.39 58.75 68.69 74.16 79.42
## mu[241] 4.31 0.04 1.88 1.50 2.89 4.01 5.42
## mu[242] 4.38 0.03 1.85 1.45 3.09 4.09 5.41
## mu[243] 2.63 0.03 1.42 0.69 1.61 2.33 3.36
## mu[244] 1.41 0.02 0.94 0.28 0.76 1.20 1.84
## mu[245] 11.23 0.06 3.19 5.76 9.01 10.83 13.14
## mu[246] 6.61 0.05 2.45 2.85 4.85 6.34 8.09
## mu[247] 5.05 0.04 2.05 1.99 3.52 4.66 6.30
## mu[248] 19.47 0.08 4.33 11.96 16.45 19.05 22.05
## mu[249] 1.07 0.02 0.74 0.19 0.54 0.91 1.42
## mu[250] 4.39 0.03 1.91 1.57 3.00 4.07 5.51
## mu[251] 1.59 0.02 1.00 0.34 0.88 1.38 2.05
## mu[252] 8.44 0.05 2.76 4.14 6.43 8.10 10.15
## mu[253] 0.82 0.01 0.60 0.14 0.41 0.67 1.06
## mu[254] 34.13 0.10 6.08 23.13 29.88 33.81 37.82
## mu[255] 16.87 0.06 3.91 10.08 13.97 16.60 19.35
## mu[256] 1.22 0.02 0.84 0.22 0.63 1.00 1.62
## mu[257] 4.46 0.04 1.79 1.78 3.09 4.20 5.58
## mu[258] 8.72 0.05 2.77 4.16 6.81 8.37 10.46
## mu[259] 23.19 0.08 4.82 14.89 19.73 22.85 26.33
## mu[260] 22.33 0.09 4.65 14.27 18.96 22.09 25.33
## mu[261] 4.22 0.03 1.81 1.58 2.85 3.98 5.26
## mu[262] 3.25 0.03 1.58 0.97 2.13 2.97 4.10
## mu[263] 4.73 0.03 2.00 1.67 3.26 4.45 5.86
## mu[264] 27.48 0.09 4.97 18.77 24.15 27.21 30.54
## mu[265] 1.86 0.02 1.11 0.41 1.04 1.64 2.39
## mu[266] 5.43 0.04 2.08 2.15 3.94 5.08 6.55
## mu[267] 13.63 0.07 3.70 7.11 11.07 13.31 15.69
## mu[268] 7.98 0.05 2.67 3.68 6.09 7.76 9.54
## mu[269] 1.32 0.02 0.87 0.27 0.69 1.11 1.72
## mu[270] 1.08 0.02 0.75 0.21 0.56 0.89 1.39
## mu[271] 3.74 0.03 1.64 1.33 2.55 3.47 4.64
## mu[272] 10.71 0.06 3.02 5.75 8.49 10.51 12.44
## mu[273] 3.14 0.03 1.60 0.93 2.01 2.85 3.91
## mu[274] 0.83 0.01 0.59 0.16 0.42 0.66 1.10
## mu[275] 2.39 0.03 1.29 0.63 1.45 2.11 3.05
## mu[276] 0.88 0.01 0.62 0.16 0.44 0.72 1.16
## mu[277] 4.37 0.03 1.91 1.58 2.94 4.11 5.41
## mu[278] 6.96 0.05 2.42 3.16 5.17 6.63 8.39
## mu[279] 15.05 0.06 3.67 8.53 12.46 14.82 17.25
## mu[280] 1.07 0.02 0.79 0.17 0.47 0.86 1.46
## mu[281] 76.29 0.14 8.24 61.72 70.34 76.00 81.67
## mu[282] 2.19 0.02 1.23 0.54 1.27 1.97 2.85
## mu[283] 8.57 0.05 2.84 4.16 6.53 8.12 10.35
## mu[284] 9.68 0.05 2.94 4.86 7.51 9.37 11.55
## mu[285] 1.48 0.02 0.93 0.35 0.83 1.27 1.92
## mu[286] 8.66 0.05 2.84 4.15 6.75 8.31 10.32
## mu[287] 2.21 0.02 1.17 0.60 1.34 2.00 2.85
## mu[288] 4.36 0.03 1.80 1.67 3.08 4.13 5.34
## mu[289] 0.92 0.02 0.70 0.15 0.42 0.73 1.21
## mu[290] 2.85 0.02 1.44 0.78 1.79 2.59 3.66
## mu[291] 131.27 0.18 11.43 109.36 123.85 130.86 138.78
## mu[292] 15.98 0.06 3.94 9.02 13.14 15.64 18.56
## mu[293] 2.28 0.02 1.26 0.57 1.34 2.05 2.90
## mu[294] 4.30 0.03 1.80 1.57 2.98 4.06 5.28
## mu[295] 1.30 0.02 0.86 0.26 0.69 1.10 1.66
## mu[296] 22.09 0.09 4.68 14.09 18.90 21.81 25.04
## mu[297] 2.57 0.02 1.30 0.81 1.65 2.34 3.22
## mu[298] 6.48 0.04 2.29 2.88 4.83 6.23 7.75
## mu[299] 10.75 0.06 3.28 5.26 8.41 10.48 12.61
## mu[300] 6.22 0.05 2.28 2.61 4.53 5.96 7.59
## mu[301] 2.08 0.02 1.14 0.54 1.29 1.86 2.60
## mu[302] 6.01 0.04 2.31 2.43 4.34 5.67 7.33
## mu[303] 3.10 0.03 1.55 0.95 1.99 2.78 4.03
## mu[304] 9.05 0.05 2.96 4.27 6.93 8.70 10.78
## mu[305] 2.21 0.03 1.24 0.57 1.29 1.96 2.86
## mu[306] 2.47 0.02 1.31 0.65 1.50 2.23 3.20
## mu[307] 17.79 0.07 3.98 10.94 15.01 17.47 20.12
## mu[308] 1.46 0.02 0.96 0.29 0.76 1.21 1.92
## mu[309] 4.28 0.04 1.86 1.58 2.89 3.98 5.32
## mu[310] 3.20 0.03 1.44 1.06 2.19 2.99 4.01
## mu[311] 4.11 0.03 1.76 1.49 2.83 3.83 5.11
## mu[312] 4.67 0.03 1.89 1.81 3.27 4.42 5.77
## mu[313] 13.21 0.06 3.52 7.36 10.69 12.89 15.35
## mu[314] 1.53 0.02 1.00 0.31 0.82 1.28 1.99
## mu[315] 1.64 0.02 0.98 0.36 0.91 1.41 2.15
## mu[316] 0.74 0.01 0.54 0.13 0.35 0.59 0.98
## mu[317] 2.76 0.03 1.48 0.74 1.69 2.47 3.47
## mu[318] 3.28 0.03 1.56 1.02 2.18 3.00 4.14
## mu[319] 39.54 0.11 5.90 28.61 35.34 39.11 43.63
## mu[320] 23.65 0.08 4.33 15.51 20.76 23.47 26.36
## mu[321] 4.54 0.04 1.94 1.67 3.14 4.26 5.71
## mu[322] 1.26 0.02 0.85 0.24 0.65 1.03 1.67
## mu[323] 34.18 0.09 5.83 23.42 30.27 33.81 37.55
## mu[324] 4.64 0.03 1.91 1.73 3.26 4.34 5.75
## mu[325] 10.42 0.05 3.15 5.11 8.12 10.20 12.37
## mu[326] 2.52 0.03 1.37 0.70 1.52 2.23 3.24
## mu[327] 5.64 0.04 2.12 2.32 4.13 5.39 6.90
## mu[328] 1.95 0.02 1.15 0.44 1.14 1.71 2.50
## mu[329] 26.77 0.08 5.04 17.76 23.17 26.49 30.01
## mu[330] 3.28 0.03 1.58 1.04 2.14 3.02 4.09
## mu[331] 2.52 0.02 1.32 0.71 1.55 2.25 3.24
## mu[332] 4.29 0.03 1.80 1.58 3.02 4.03 5.29
## mu[333] 21.07 0.07 4.34 13.24 18.05 20.76 23.81
## mu[334] 14.28 0.07 3.55 8.16 11.73 13.99 16.58
## mu[335] 2.45 0.02 1.28 0.68 1.54 2.22 3.08
## mu[336] 26.28 0.10 4.77 17.71 22.85 26.04 29.38
## mu[337] 30.09 0.10 5.44 20.59 26.27 29.79 33.72
## mu[338] 27.23 0.08 4.87 18.62 23.68 26.96 30.33
## mu[339] 12.39 0.05 3.19 6.83 10.22 12.12 14.34
## mu[340] 3.88 0.03 1.74 1.28 2.60 3.60 4.88
## mu[341] 1.44 0.02 0.92 0.31 0.79 1.21 1.89
## mu[342] 2.92 0.03 1.51 0.91 1.79 2.64 3.78
## mu[343] 2.23 0.02 1.23 0.60 1.34 1.94 2.82
## mu[344] 1.36 0.02 0.85 0.30 0.77 1.16 1.74
## mu[345] 0.83 0.01 0.65 0.12 0.38 0.65 1.08
## mu[346] 1.07 0.01 0.76 0.18 0.51 0.85 1.42
## mu[347] 6.81 0.04 2.41 2.98 5.15 6.50 8.18
## mu[348] 17.86 0.07 3.88 11.11 15.21 17.66 20.26
## mu[349] 10.98 0.05 3.20 5.79 8.58 10.57 12.87
## mu[350] 3.78 0.03 1.72 1.20 2.56 3.54 4.63
## mu[351] 18.33 0.08 4.31 10.62 15.32 18.05 21.01
## mu[352] 9.50 0.05 2.94 4.82 7.33 9.16 11.30
## mu[353] 3.08 0.03 1.52 0.95 2.01 2.77 3.85
## mu[354] 172.93 0.23 13.84 147.04 163.51 172.50 181.90
## mu[355] 3.15 0.03 1.66 0.89 1.93 2.79 4.08
## mu[356] 2.04 0.03 1.19 0.52 1.16 1.76 2.64
## mu[357] 5.71 0.04 2.20 2.24 4.12 5.44 6.96
## mu[358] 75.75 0.14 8.83 59.21 69.61 75.47 81.32
## mu[359] 4.15 0.03 1.88 1.42 2.78 3.83 5.16
## mu[360] 4.02 0.03 1.75 1.37 2.75 3.75 5.01
## mu[361] 0.57 0.01 0.45 0.07 0.26 0.44 0.74
## mu[362] 29.14 0.09 5.28 19.37 25.57 28.78 32.24
## mu[363] 8.26 0.05 2.54 4.05 6.45 7.98 9.72
## mu[364] 1.40 0.02 0.91 0.28 0.78 1.19 1.79
## mu[365] 9.13 0.05 2.88 4.60 6.98 8.79 10.91
## mu[366] 8.17 0.05 2.76 3.79 6.04 7.91 9.89
## mu[367] 16.01 0.08 3.93 9.26 13.29 15.65 18.33
## mu[368] 1.55 0.02 0.97 0.31 0.85 1.35 1.99
## mu[369] 8.96 0.05 2.90 4.37 6.85 8.58 10.72
## mu[370] 21.45 0.09 4.44 13.68 18.35 21.12 24.11
## mu[371] 3.37 0.03 1.58 1.06 2.22 3.10 4.24
## mu[372] 12.78 0.06 3.31 7.15 10.45 12.54 14.71
## mu[373] 9.36 0.05 2.88 4.68 7.19 9.06 11.23
## mu[374] 1.58 0.02 0.95 0.40 0.89 1.37 2.07
## mu[375] 5.01 0.04 2.02 1.88 3.60 4.71 6.17
## mu[376] 11.92 0.06 3.37 6.30 9.57 11.63 13.90
## mu[377] 6.95 0.04 2.56 2.93 5.15 6.63 8.38
## mu[378] 9.35 0.05 2.91 4.60 7.27 9.04 11.06
## mu[379] 21.38 0.07 4.45 13.30 18.30 21.05 24.20
## mu[380] 4.53 0.03 1.89 1.68 3.19 4.22 5.57
## mu[381] 1.74 0.02 1.05 0.39 0.96 1.49 2.27
## mu[382] 1.92 0.03 1.16 0.45 1.08 1.65 2.45
## mu[383] 1.17 0.02 0.82 0.23 0.58 0.98 1.50
## mu[384] 2.41 0.02 1.23 0.71 1.54 2.21 2.97
## mu[385] 9.87 0.05 3.02 4.97 7.81 9.54 11.64
## mu[386] 1.00 0.02 0.73 0.17 0.49 0.80 1.32
## mu[387] 2.04 0.02 1.14 0.48 1.24 1.81 2.61
## mu[388] 6.49 0.04 2.34 2.88 4.86 6.17 7.79
## mu[389] 9.47 0.05 2.72 4.84 7.54 9.19 11.20
## mu[390] 3.28 0.03 1.64 0.96 2.11 2.97 4.10
## mu[391] 2.06 0.02 1.13 0.56 1.28 1.87 2.62
## mu[392] 4.63 0.03 1.83 1.74 3.32 4.38 5.72
## mu[393] 21.01 0.08 4.50 13.16 17.90 20.73 23.83
## mu[394] 1.60 0.02 0.94 0.38 0.90 1.39 2.09
## mu[395] 11.39 0.05 3.15 6.15 9.01 11.15 13.45
## mu[396] 5.40 0.04 1.99 2.27 3.98 5.12 6.55
## mu[397] 1.25 0.02 0.81 0.26 0.69 1.04 1.59
## mu[398] 86.39 0.15 9.01 69.43 80.02 86.10 92.52
## mu[399] 2.33 0.02 1.26 0.59 1.41 2.09 3.00
## mu[400] 13.64 0.06 3.60 7.52 11.08 13.35 15.84
## mu[401] 4.89 0.04 1.98 1.83 3.42 4.64 6.03
## mu[402] 13.22 0.06 3.44 7.42 10.90 12.86 15.31
## mu[403] 1.21 0.02 0.84 0.23 0.62 1.00 1.58
## mu[404] 15.48 0.07 3.83 9.09 12.65 15.12 17.97
## mu[405] 40.40 0.09 5.92 29.74 36.52 40.08 44.11
## mu[406] 5.30 0.04 2.09 2.07 3.75 5.02 6.54
## mu[407] 1.47 0.02 0.91 0.31 0.83 1.25 1.90
## mu[408] 56.33 0.13 7.70 42.76 51.07 56.04 61.37
## mu[409] 7.13 0.04 2.48 3.42 5.31 6.85 8.60
## mu[410] 4.78 0.04 1.94 1.90 3.30 4.50 5.93
## mu[411] 6.53 0.05 2.40 2.80 4.82 6.24 7.77
## mu[412] 6.35 0.04 2.20 2.71 4.77 6.13 7.68
## mu[413] 1.85 0.02 1.14 0.41 1.01 1.60 2.42
## mu[414] 1.07 0.02 0.77 0.18 0.52 0.88 1.38
## mu[415] 4.49 0.04 1.97 1.57 3.07 4.20 5.62
## mu[416] 5.18 0.03 2.04 1.95 3.73 4.92 6.34
## mu[417] 3.40 0.03 1.67 1.05 2.13 3.06 4.40
## mu[418] 9.91 0.05 3.01 4.80 7.77 9.65 11.70
## mu[419] 3.20 0.03 1.62 0.91 2.08 2.89 4.02
## mu[420] 27.13 0.09 5.11 18.09 23.53 26.84 30.44
## mu[421] 7.03 0.05 2.47 3.13 5.19 6.71 8.48
## mu[422] 11.08 0.06 3.14 6.01 8.78 10.69 13.12
## mu[423] 2.08 0.02 1.17 0.53 1.23 1.85 2.63
## mu[424] 5.66 0.04 2.16 2.21 4.05 5.38 6.90
## mu[425] 10.64 0.05 3.09 5.70 8.35 10.29 12.72
## mu[426] 15.48 0.08 3.63 9.21 12.97 15.22 17.70
## mu[427] 1.08 0.02 0.76 0.18 0.54 0.89 1.42
## mu[428] 5.97 0.04 2.40 2.29 4.21 5.60 7.42
## mu[429] 1.71 0.02 1.04 0.38 0.93 1.48 2.25
## mu[430] 1.96 0.02 1.11 0.46 1.16 1.74 2.54
## mu[431] 2.94 0.03 1.50 0.91 1.83 2.70 3.66
## mu[432] 8.94 0.05 2.81 4.41 6.84 8.67 10.72
## mu[433] 16.93 0.08 4.39 9.43 13.92 16.57 19.43
## mu[434] 3.35 0.03 1.63 1.02 2.13 3.04 4.29
## mu[435] 6.64 0.04 2.35 2.89 4.95 6.34 8.00
## mu[436] 6.35 0.05 2.36 2.51 4.69 6.15 7.66
## mu[437] 6.21 0.05 2.24 2.75 4.57 5.96 7.47
## mu[438] 1.78 0.02 1.08 0.39 1.01 1.53 2.33
## mu[439] 6.09 0.04 2.26 2.59 4.52 5.77 7.38
## mu[440] 18.78 0.08 4.22 11.75 15.75 18.48 21.47
## mu[441] 4.26 0.03 1.83 1.49 2.90 4.01 5.28
## mu[442] 2.32 0.03 1.24 0.64 1.42 2.09 2.98
## mu[443] 81.23 0.15 8.67 65.43 75.08 81.07 87.03
## mu[444] 4.24 0.03 1.86 1.46 2.90 3.99 5.21
## mu[445] 3.35 0.03 1.66 1.04 2.13 3.08 4.17
## mu[446] 3.02 0.03 1.48 0.84 1.98 2.79 3.72
## mu[447] 3.36 0.03 1.68 1.04 2.11 3.03 4.27
## mu[448] 3.45 0.03 1.59 1.20 2.31 3.14 4.39
## mu[449] 13.02 0.07 3.38 7.55 10.49 12.76 15.15
## mu[450] 1.48 0.02 0.96 0.30 0.77 1.25 1.97
## mu[451] 3.65 0.03 1.69 1.19 2.40 3.39 4.64
## mu[452] 2.54 0.03 1.36 0.73 1.53 2.28 3.24
## mu[453] 1.27 0.02 0.85 0.24 0.66 1.04 1.72
## mu[454] 1.26 0.02 0.85 0.25 0.63 1.04 1.66
## mu[455] 0.83 0.02 0.67 0.12 0.38 0.63 1.09
## mu[456] 5.00 0.04 2.05 1.90 3.48 4.73 6.21
## mu[457] 1.48 0.02 0.87 0.38 0.84 1.29 1.90
## mu[458] 4.98 0.04 2.04 1.88 3.49 4.70 6.07
## mu[459] 31.96 0.12 5.64 22.05 27.92 31.61 35.71
## mu[460] 2.76 0.02 1.38 0.75 1.80 2.53 3.50
## mu[461] 5.77 0.04 2.28 2.24 4.16 5.44 6.98
## mu[462] 5.51 0.04 2.09 2.24 4.00 5.23 6.68
## mu[463] 1.09 0.01 0.73 0.20 0.56 0.93 1.39
## mu[464] 1.40 0.02 0.84 0.30 0.79 1.23 1.80
## mu[465] 4.42 0.03 1.89 1.68 3.06 4.13 5.45
## mu[466] 15.87 0.07 3.93 9.01 13.04 15.54 18.44
## mu[467] 2.37 0.02 1.23 0.67 1.44 2.15 3.06
## mu[468] 1.25 0.02 0.84 0.24 0.66 1.04 1.66
## mu[469] 3.39 0.03 1.53 1.17 2.27 3.14 4.17
## mu[470] 29.17 0.08 4.91 20.43 25.81 28.77 32.27
## mu[471] 1.55 0.02 0.92 0.36 0.87 1.35 2.01
## mu[472] 19.06 0.08 4.32 11.86 15.98 18.78 21.76
## mu[473] 3.37 0.03 1.53 1.12 2.26 3.14 4.23
## mu[474] 8.25 0.05 2.74 3.91 6.33 7.87 9.80
## mu[475] 17.28 0.07 3.92 10.53 14.50 16.97 19.78
## mu[476] 16.99 0.06 3.87 10.13 14.27 16.76 19.49
## mu[477] 8.81 0.05 2.79 4.27 6.84 8.41 10.47
## mu[478] 20.02 0.09 4.48 12.28 16.86 19.78 22.73
## mu[479] 16.28 0.08 4.01 9.47 13.43 15.93 18.87
## mu[480] 12.78 0.06 3.39 7.07 10.34 12.47 14.84
## mu[481] 2.58 0.02 1.37 0.73 1.53 2.32 3.37
## mu[482] 7.45 0.05 2.61 3.28 5.54 7.19 9.07
## mu[483] 4.68 0.04 1.92 1.74 3.36 4.36 5.77
## mu[484] 7.69 0.05 2.69 3.45 5.71 7.32 9.32
## mu[485] 2.41 0.02 1.25 0.67 1.46 2.18 3.10
## mu[486] 6.30 0.04 2.27 2.76 4.61 5.96 7.79
## mu[487] 92.48 0.18 9.41 75.12 85.99 92.34 98.24
## mu[488] 14.85 0.07 3.79 8.47 12.16 14.49 17.12
## mu[489] 1.03 0.02 0.74 0.19 0.51 0.84 1.32
## mu[490] 3.76 0.03 1.65 1.24 2.53 3.52 4.80
## mu[491] 64.02 0.16 7.74 49.64 58.93 63.60 68.81
## mu[492] 7.19 0.05 2.57 3.14 5.37 6.90 8.71
## mu[493] 1.44 0.02 0.91 0.29 0.77 1.21 1.88
## mu[494] 5.98 0.04 2.24 2.59 4.30 5.67 7.34
## mu[495] 2.35 0.03 1.27 0.58 1.48 2.12 2.98
## mu[496] 2.60 0.03 1.47 0.69 1.53 2.33 3.26
## mu[497] 20.04 0.08 4.27 12.64 17.04 19.72 22.67
## mu[498] 183.63 0.23 13.38 158.40 174.53 183.46 192.32
## mu[499] 4.97 0.04 1.99 1.89 3.55 4.70 6.01
## mu[500] 3.45 0.03 1.68 1.08 2.22 3.14 4.37
## a[1] 1.67 0.00 0.07 1.53 1.63 1.67 1.72
## a[2] 1.43 0.00 0.07 1.28 1.38 1.43 1.47
## b 0.87 0.00 0.05 0.77 0.84 0.87 0.90
## c[1] 0.27 0.00 0.06 0.14 0.22 0.26 0.31
## c[2] -0.19 0.00 0.08 -0.34 -0.24 -0.18 -0.13
## sigma 0.95 0.00 0.04 0.87 0.92 0.95 0.98
## lp__ 12747.44 0.87 18.07 12711.84 12734.64 12748.34 12760.16
## 97.5% n_eff Rhat
## mu[1] 2.96 2009 1.00
## mu[2] 15.39 2380 1.00
## mu[3] 36.91 2842 1.00
## mu[4] 4.76 2703 1.00
## mu[5] 19.59 2661 1.00
## mu[6] 6.89 3424 1.00
## mu[7] 5.72 3525 1.00
## mu[8] 7.91 3754 1.00
## mu[9] 2.70 2234 1.00
## mu[10] 3.38 2530 1.00
## mu[11] 6.01 2947 1.00
## mu[12] 17.74 2367 1.00
## mu[13] 7.16 2422 1.00
## mu[14] 18.93 2921 1.00
## mu[15] 1.69 1803 1.00
## mu[16] 6.71 2167 1.00
## mu[17] 2.49 1943 1.00
## mu[18] 83.27 3413 1.00
## mu[19] 3.90 2225 1.00
## mu[20] 4.37 3731 1.00
## mu[21] 5.57 2315 1.00
## mu[22] 8.11 2213 1.00
## mu[23] 13.73 3153 1.00
## mu[24] 39.00 2168 1.00
## mu[25] 5.75 3549 1.00
## mu[26] 7.05 3036 1.00
## mu[27] 2.72 2406 1.00
## mu[28] 3.73 2689 1.00
## mu[29] 4.14 2188 1.00
## mu[30] 5.30 3609 1.00
## mu[31] 15.90 2991 1.00
## mu[32] 81.49 3433 1.00
## mu[33] 126.87 2990 1.00
## mu[34] 7.93 2361 1.00
## mu[35] 34.72 2798 1.00
## mu[36] 7.38 3863 1.00
## mu[37] 12.49 3497 1.00
## mu[38] 13.38 2784 1.00
## mu[39] 6.54 2795 1.00
## mu[40] 25.75 3469 1.00
## mu[41] 4.88 3000 1.00
## mu[42] 17.11 2434 1.00
## mu[43] 9.15 3714 1.00
## mu[44] 92.39 4166 1.00
## mu[45] 2.89 2345 1.00
## mu[46] 2.43 2650 1.00
## mu[47] 11.28 2857 1.00
## mu[48] 3.20 2260 1.00
## mu[49] 4.25 2852 1.00
## mu[50] 7.82 2906 1.00
## mu[51] 7.12 3287 1.00
## mu[52] 3.82 2448 1.00
## mu[53] 153.18 2526 1.00
## mu[54] 27.37 3519 1.00
## mu[55] 85.19 3790 1.00
## mu[56] 3.37 2347 1.00
## mu[57] 6.58 2483 1.00
## mu[58] 3.38 2792 1.00
## mu[59] 1.86 2779 1.00
## mu[60] 11.60 2625 1.00
## mu[61] 75.18 2730 1.00
## mu[62] 5.31 2602 1.00
## mu[63] 4.20 2882 1.00
## mu[64] 40.68 3728 1.00
## mu[65] 4.11 3145 1.00
## mu[66] 4.36 2656 1.00
## mu[67] 26.40 3755 1.00
## mu[68] 4.99 2525 1.00
## mu[69] 28.49 3191 1.00
## mu[70] 6.38 2747 1.00
## mu[71] 14.23 3445 1.00
## mu[72] 2.76 2280 1.00
## mu[73] 8.32 2590 1.00
## mu[74] 15.70 2325 1.00
## mu[75] 16.08 3485 1.00
## mu[76] 10.06 2709 1.00
## mu[77] 31.87 3375 1.00
## mu[78] 5.44 2670 1.00
## mu[79] 6.42 2632 1.00
## mu[80] 48.39 3910 1.00
## mu[81] 5.36 2825 1.00
## mu[82] 11.03 3463 1.00
## mu[83] 13.40 3216 1.00
## mu[84] 34.41 2682 1.00
## mu[85] 18.26 4602 1.00
## mu[86] 5.23 2505 1.00
## mu[87] 2.85 1814 1.00
## mu[88] 8.16 4026 1.00
## mu[89] 6.87 2809 1.00
## mu[90] 5.00 2575 1.00
## mu[91] 5.48 2383 1.00
## mu[92] 5.25 2747 1.00
## mu[93] 9.97 3371 1.00
## mu[94] 29.70 3500 1.00
## mu[95] 6.89 2950 1.00
## mu[96] 15.91 2815 1.00
## mu[97] 7.87 2840 1.00
## mu[98] 5.33 2530 1.00
## mu[99] 4.36 3296 1.00
## mu[100] 51.71 2960 1.00
## mu[101] 13.22 2719 1.00
## mu[102] 37.46 2842 1.00
## mu[103] 6.33 3176 1.00
## mu[104] 4.76 2437 1.00
## mu[105] 3.51 2396 1.00
## mu[106] 4.86 2659 1.00
## mu[107] 15.10 4274 1.00
## mu[108] 25.67 3294 1.00
## mu[109] 15.42 3222 1.00
## mu[110] 4.70 1917 1.00
## mu[111] 6.01 2423 1.00
## mu[112] 6.47 2445 1.00
## mu[113] 22.56 2857 1.00
## mu[114] 13.15 3775 1.00
## mu[115] 10.02 2888 1.00
## mu[116] 4.43 3035 1.00
## mu[117] 22.33 2604 1.00
## mu[118] 4.05 2351 1.00
## mu[119] 14.44 3200 1.00
## mu[120] 29.96 3414 1.00
## mu[121] 5.85 2886 1.00
## mu[122] 23.21 3924 1.00
## mu[123] 5.16 3008 1.00
## mu[124] 3.27 2183 1.00
## mu[125] 4.80 3381 1.00
## mu[126] 7.35 2453 1.00
## mu[127] 32.33 2933 1.00
## mu[128] 4.74 2622 1.00
## mu[129] 10.31 3401 1.00
## mu[130] 14.79 3251 1.00
## mu[131] 4.34 2704 1.00
## mu[132] 12.74 3055 1.00
## mu[133] 3.18 2558 1.00
## mu[134] 38.62 3874 1.00
## mu[135] 21.42 3236 1.00
## mu[136] 9.49 2964 1.00
## mu[137] 14.34 3662 1.00
## mu[138] 15.36 3565 1.00
## mu[139] 9.03 3641 1.00
## mu[140] 13.88 2601 1.00
## mu[141] 5.06 3218 1.00
## mu[142] 2.70 2427 1.00
## mu[143] 38.91 3194 1.00
## mu[144] 6.33 2327 1.00
## mu[145] 27.31 2552 1.00
## mu[146] 5.63 2518 1.00
## mu[147] 6.08 3174 1.00
## mu[148] 3.35 2696 1.00
## mu[149] 10.91 2510 1.00
## mu[150] 9.75 2324 1.00
## mu[151] 9.48 2226 1.00
## mu[152] 7.70 3328 1.00
## mu[153] 5.46 2523 1.00
## mu[154] 6.68 2900 1.00
## mu[155] 11.99 3196 1.00
## mu[156] 2.30 2264 1.00
## mu[157] 32.99 2919 1.00
## mu[158] 60.40 3669 1.00
## mu[159] 2.84 2238 1.00
## mu[160] 7.65 3119 1.00
## mu[161] 24.83 3283 1.00
## mu[162] 4.70 2869 1.00
## mu[163] 128.49 4274 1.00
## mu[164] 5.84 2677 1.00
## mu[165] 18.20 4096 1.00
## mu[166] 2.21 2013 1.00
## mu[167] 19.36 3353 1.00
## mu[168] 15.73 2886 1.00
## mu[169] 8.37 3911 1.00
## mu[170] 12.64 2853 1.00
## mu[171] 5.02 2677 1.00
## mu[172] 5.27 2541 1.00
## mu[173] 18.39 2112 1.00
## mu[174] 6.71 3648 1.00
## mu[175] 123.17 2526 1.00
## mu[176] 13.29 3563 1.00
## mu[177] 14.73 3158 1.00
## mu[178] 6.83 2587 1.00
## mu[179] 13.29 3160 1.00
## mu[180] 5.54 3242 1.00
## mu[181] 20.89 3569 1.00
## mu[182] 3.63 1735 1.00
## mu[183] 4.52 2061 1.00
## mu[184] 20.52 3057 1.00
## mu[185] 5.39 2958 1.00
## mu[186] 10.19 2434 1.00
## mu[187] 4.34 2610 1.00
## mu[188] 15.72 3140 1.00
## mu[189] 38.07 3318 1.00
## mu[190] 4.18 3393 1.00
## mu[191] 14.11 3551 1.00
## mu[192] 11.30 3796 1.00
## mu[193] 2.52 1952 1.00
## mu[194] 7.79 2453 1.00
## mu[195] 17.37 3271 1.00
## mu[196] 8.10 2211 1.00
## mu[197] 5.31 2606 1.00
## mu[198] 11.61 2896 1.00
## mu[199] 4.28 2329 1.00
## mu[200] 2.96 2094 1.00
## mu[201] 26.00 2900 1.00
## mu[202] 6.75 2913 1.00
## mu[203] 9.73 3278 1.00
## mu[204] 3.13 2077 1.00
## mu[205] 9.30 2392 1.00
## mu[206] 2.16 2669 1.00
## mu[207] 5.72 2754 1.00
## mu[208] 6.12 2929 1.00
## mu[209] 4.15 2719 1.00
## mu[210] 9.36 2334 1.00
## mu[211] 26.59 3763 1.00
## mu[212] 10.51 2602 1.00
## mu[213] 38.63 2666 1.00
## mu[214] 9.65 3489 1.00
## mu[215] 4.31 2178 1.00
## mu[216] 2.84 2233 1.00
## mu[217] 26.37 3686 1.00
## mu[218] 39.61 3154 1.00
## mu[219] 4.02 2785 1.00
## mu[220] 7.98 2609 1.00
## mu[221] 4.69 3221 1.00
## mu[222] 60.93 4082 1.00
## mu[223] 10.83 3000 1.00
## mu[224] 9.74 3014 1.00
## mu[225] 10.52 4205 1.00
## mu[226] 52.17 2912 1.00
## mu[227] 9.74 2993 1.00
## mu[228] 7.19 3532 1.00
## mu[229] 67.11 3056 1.00
## mu[230] 5.01 2985 1.00
## mu[231] 123.48 3368 1.00
## mu[232] 13.02 3344 1.00
## mu[233] 24.15 2987 1.00
## mu[234] 97.24 2965 1.00
## mu[235] 11.71 2902 1.00
## mu[236] 10.50 2497 1.00
## mu[237] 20.40 3630 1.00
## mu[238] 27.63 3658 1.00
## mu[239] 6.01 2841 1.00
## mu[240] 91.54 3353 1.00
## mu[241] 8.73 2728 1.00
## mu[242] 8.76 3165 1.00
## mu[243] 6.22 2056 1.00
## mu[244] 3.74 2387 1.00
## mu[245] 18.50 3281 1.00
## mu[246] 11.99 2404 1.00
## mu[247] 9.87 3254 1.00
## mu[248] 28.97 3091 1.00
## mu[249] 2.95 1933 1.00
## mu[250] 8.56 3319 1.00
## mu[251] 4.09 2478 1.00
## mu[252] 14.80 3183 1.00
## mu[253] 2.37 1935 1.00
## mu[254] 46.58 3671 1.00
## mu[255] 25.34 3793 1.00
## mu[256] 3.39 2508 1.00
## mu[257] 8.44 2473 1.00
## mu[258] 15.00 3154 1.00
## mu[259] 33.21 3265 1.00
## mu[260] 32.21 2983 1.00
## mu[261] 8.42 3947 1.00
## mu[262] 7.18 2239 1.00
## mu[263] 9.39 3383 1.00
## mu[264] 38.07 3346 1.00
## mu[265] 4.68 2066 1.00
## mu[266] 10.47 2512 1.00
## mu[267] 21.98 2982 1.00
## mu[268] 14.09 3423 1.00
## mu[269] 3.68 2521 1.00
## mu[270] 3.09 2408 1.00
## mu[271] 7.70 3096 1.00
## mu[272] 17.40 2375 1.00
## mu[273] 7.19 3301 1.00
## mu[274] 2.24 2508 1.00
## mu[275] 5.54 2495 1.00
## mu[276] 2.49 2879 1.00
## mu[277] 8.78 4507 1.00
## mu[278] 12.47 2385 1.00
## mu[279] 22.66 3201 1.00
## mu[280] 3.07 2197 1.00
## mu[281] 93.33 3522 1.00
## mu[282] 5.13 2752 1.00
## mu[283] 14.87 2704 1.00
## mu[284] 16.03 3706 1.00
## mu[285] 3.65 3054 1.00
## mu[286] 15.08 2792 1.00
## mu[287] 5.07 2536 1.00
## mu[288] 8.64 3175 1.00
## mu[289] 2.63 2054 1.00
## mu[290] 6.35 3467 1.00
## mu[291] 153.91 4009 1.00
## mu[292] 24.30 4533 1.00
## mu[293] 5.56 2702 1.00
## mu[294] 8.34 3178 1.00
## mu[295] 3.68 2115 1.00
## mu[296] 32.22 2691 1.00
## mu[297] 5.75 2710 1.00
## mu[298] 11.88 3573 1.00
## mu[299] 18.31 3207 1.00
## mu[300] 11.21 2515 1.00
## mu[301] 4.89 2672 1.00
## mu[302] 11.27 2762 1.00
## mu[303] 6.76 3327 1.00
## mu[304] 15.67 3008 1.00
## mu[305] 5.27 2253 1.00
## mu[306] 5.49 3324 1.00
## mu[307] 27.00 2964 1.00
## mu[308] 3.80 1952 1.00
## mu[309] 8.63 2791 1.00
## mu[310] 6.74 3062 1.00
## mu[311] 8.22 2675 1.00
## mu[312] 9.13 3798 1.00
## mu[313] 21.04 3656 1.00
## mu[314] 4.01 2231 1.00
## mu[315] 3.88 2319 1.00
## mu[316] 2.19 2244 1.00
## mu[317] 6.39 2823 1.00
## mu[318] 6.99 2753 1.00
## mu[319] 51.60 3023 1.00
## mu[320] 33.09 3252 1.00
## mu[321] 9.01 2716 1.00
## mu[322] 3.42 1796 1.00
## mu[323] 46.82 4007 1.00
## mu[324] 9.25 3125 1.00
## mu[325] 17.83 3479 1.00
## mu[326] 5.95 2369 1.00
## mu[327] 10.23 3467 1.00
## mu[328] 4.96 3266 1.00
## mu[329] 37.93 3927 1.00
## mu[330] 7.14 2424 1.00
## mu[331] 5.78 2965 1.00
## mu[332] 8.66 3241 1.00
## mu[333] 30.26 3735 1.00
## mu[334] 21.55 2543 1.00
## mu[335] 5.68 2819 1.00
## mu[336] 36.44 2339 1.00
## mu[337] 41.35 2968 1.00
## mu[338] 37.34 4016 1.00
## mu[339] 19.51 4107 1.00
## mu[340] 7.65 3370 1.00
## mu[341] 3.77 2781 1.00
## mu[342] 6.60 3368 1.00
## mu[343] 5.23 2839 1.00
## mu[344] 3.57 2433 1.00
## mu[345] 2.62 2130 1.00
## mu[346] 2.96 2760 1.00
## mu[347] 12.03 3705 1.00
## mu[348] 25.95 2700 1.00
## mu[349] 18.06 3784 1.00
## mu[350] 7.83 3249 1.00
## mu[351] 27.38 3060 1.00
## mu[352] 15.74 3482 1.00
## mu[353] 6.64 2828 1.00
## mu[354] 202.46 3639 1.00
## mu[355] 7.28 2553 1.00
## mu[356] 4.96 2226 1.00
## mu[357] 10.75 2614 1.00
## mu[358] 94.69 3859 1.00
## mu[359] 8.86 3484 1.00
## mu[360] 8.22 2711 1.00
## mu[361] 1.73 2221 1.00
## mu[362] 41.11 3615 1.00
## mu[363] 14.00 3006 1.00
## mu[364] 3.61 1985 1.00
## mu[365] 15.37 3495 1.00
## mu[366] 14.18 2534 1.00
## mu[367] 25.07 2675 1.00
## mu[368] 4.03 2580 1.00
## mu[369] 15.42 3114 1.00
## mu[370] 30.91 2444 1.00
## mu[371] 7.11 2494 1.00
## mu[372] 20.20 2644 1.00
## mu[373] 15.67 3239 1.00
## mu[374] 3.90 2229 1.00
## mu[375] 9.39 2207 1.00
## mu[376] 19.37 3107 1.00
## mu[377] 12.66 3347 1.00
## mu[378] 15.74 2805 1.00
## mu[379] 30.95 4512 1.00
## mu[380] 9.07 3036 1.00
## mu[381] 4.46 3358 1.00
## mu[382] 4.80 2016 1.00
## mu[383] 3.35 2234 1.00
## mu[384] 5.56 2534 1.00
## mu[385] 16.94 3084 1.00
## mu[386] 3.01 2094 1.00
## mu[387] 4.76 2637 1.00
## mu[388] 11.75 3059 1.00
## mu[389] 15.27 3323 1.00
## mu[390] 7.55 2226 1.00
## mu[391] 5.04 3457 1.00
## mu[392] 8.69 2993 1.00
## mu[393] 30.82 3189 1.00
## mu[394] 3.94 1718 1.00
## mu[395] 18.15 3384 1.00
## mu[396] 9.96 2723 1.00
## mu[397] 3.40 2471 1.00
## mu[398] 104.45 3539 1.00
## mu[399] 5.54 2975 1.00
## mu[400] 21.97 3358 1.00
## mu[401] 9.35 2643 1.00
## mu[402] 20.85 3475 1.00
## mu[403] 3.32 3076 1.00
## mu[404] 23.92 3152 1.00
## mu[405] 52.31 4241 1.00
## mu[406] 10.04 3393 1.00
## mu[407] 3.76 2618 1.00
## mu[408] 72.06 3377 1.00
## mu[409] 12.55 4181 1.00
## mu[410] 9.33 2509 1.00
## mu[411] 12.15 2666 1.00
## mu[412] 11.29 2876 1.00
## mu[413] 4.65 3061 1.00
## mu[414] 3.13 2196 1.00
## mu[415] 9.11 2885 1.00
## mu[416] 9.87 4072 1.00
## mu[417] 7.28 2557 1.00
## mu[418] 16.62 3477 1.00
## mu[419] 7.03 2708 1.00
## mu[420] 37.89 3023 1.00
## mu[421] 12.80 2957 1.00
## mu[422] 17.90 2571 1.00
## mu[423] 5.08 2698 1.00
## mu[424] 10.56 3479 1.00
## mu[425] 17.34 3202 1.00
## mu[426] 23.47 2210 1.00
## mu[427] 3.06 1948 1.00
## mu[428] 11.75 3632 1.00
## mu[429] 4.38 2423 1.00
## mu[430] 4.58 2246 1.00
## mu[431] 6.58 2756 1.00
## mu[432] 15.29 3026 1.00
## mu[433] 26.55 3110 1.00
## mu[434] 7.20 2292 1.00
## mu[435] 11.83 3065 1.00
## mu[436] 11.39 2670 1.00
## mu[437] 11.38 1992 1.00
## mu[438] 4.55 2068 1.00
## mu[439] 11.35 2691 1.00
## mu[440] 28.01 2583 1.00
## mu[441] 8.51 3029 1.00
## mu[442] 5.25 2372 1.00
## mu[443] 98.66 3529 1.00
## mu[444] 8.63 3355 1.00
## mu[445] 7.15 2521 1.00
## mu[446] 6.62 2388 1.00
## mu[447] 7.63 3108 1.00
## mu[448] 7.15 2158 1.00
## mu[449] 20.45 2643 1.00
## mu[450] 3.96 2657 1.00
## mu[451] 7.89 3593 1.00
## mu[452] 5.80 2780 1.00
## mu[453] 3.41 1569 1.00
## mu[454] 3.36 2492 1.00
## mu[455] 2.55 1644 1.00
## mu[456] 9.87 3100 1.00
## mu[457] 3.70 2086 1.00
## mu[458] 9.85 2759 1.00
## mu[459] 43.93 2213 1.00
## mu[460] 6.18 3545 1.00
## mu[461] 11.06 2675 1.00
## mu[462] 10.32 2867 1.00
## mu[463] 2.98 2884 1.00
## mu[464] 3.55 3039 1.00
## mu[465] 8.57 2940 1.00
## mu[466] 24.48 3596 1.00
## mu[467] 5.23 2596 1.00
## mu[468] 3.38 2607 1.00
## mu[469] 7.14 3395 1.00
## mu[470] 39.82 3443 1.00
## mu[471] 3.88 2552 1.00
## mu[472] 28.81 2955 1.00
## mu[473] 7.05 3186 1.00
## mu[474] 14.45 3113 1.00
## mu[475] 25.71 3010 1.00
## mu[476] 25.13 3945 1.00
## mu[477] 15.30 2776 1.00
## mu[478] 29.70 2340 1.00
## mu[479] 25.10 2772 1.00
## mu[480] 20.15 3059 1.00
## mu[481] 5.83 3416 1.00
## mu[482] 13.19 2842 1.00
## mu[483] 9.20 2639 1.00
## mu[484] 13.89 2887 1.00
## mu[485] 5.34 2769 1.00
## mu[486] 11.31 3057 1.00
## mu[487] 111.87 2655 1.00
## mu[488] 23.25 2907 1.00
## mu[489] 2.95 2181 1.00
## mu[490] 7.44 2940 1.00
## mu[491] 80.47 2224 1.00
## mu[492] 13.28 2813 1.00
## mu[493] 3.77 2611 1.00
## mu[494] 11.21 2784 1.00
## mu[495] 5.42 2190 1.00
## mu[496] 6.08 2310 1.00
## mu[497] 29.13 2752 1.00
## mu[498] 211.38 3301 1.00
## mu[499] 9.80 2040 1.00
## mu[500] 7.58 3233 1.00
## a[1] 1.81 1562 1.00
## a[2] 1.57 1632 1.00
## b 0.96 2223 1.00
## c[1] 0.40 2681 1.00
## c[2] -0.03 1518 1.00
## sigma 1.04 891 1.00
## lp__ 12781.31 429 1.01
##
## Samples were drawn using NUTS(diag_e) at Sat Jan 5 13:00:41 2019.
## 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?
Here are posterior distributions from the full model, with the true values in red. 
Posterior predictive simulations, again
# 100 datasets:
params2 <- list(a=colMeans(post2$a),
b=mean(post2$b),
c=colMeans(post2$c),
sigma=mean(post2$sigma))
y2 <- with(list2env(scaled_data),
params2$a[geno]
+ params2$b * age
+ params2$c[geno] * expo)
sim2 <- replicate(100, {
mu = exp(rnorm(length(y2),
mean=y2,
sd=params2$sigma))
rpois(length(mu), mu)
})model {
vector[N] y;
y = a[geno] + b * age + c[geno] .* expo;
mu ~ lognormal(y, sigma);
counts ~ poisson(mu);
}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'))
Now we cover the true data
plot(data$counts[order(mu1)], ylab="counts", ylim=c(0,500), 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:
A toy problem
Is it plausible that \[ X \sim \Normal(A, A)?\] What about \(Y\)?
data.frame(
A = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10),
X = c(-0.45, -0.16, 0.15, 0.43, 0.59, 0.87, 1.1, 1.51, 2.3, 2.42,
-7.98, -2.49, 5.15, 5.86, 7.36, 8.52, 10.53, 13, 15.51, 18.8),
Y = c(-2.18, -1.47, 0.18, 0.64, 0.81, 1.18, 2.52, 2.87, 3.08, 3.1,
-49.03, -38.01, -28.81, 0.16, 9.89, 25.92, 39.1, 49.3, 50.87, 92.11))Goodness-of-fit: design
We’ll use the fact that if the simulated data “look just like” the real data, then the rank of the real data in the simulated data should be uniform.
Since we have 100 simulated datasets, that means that the probability that the real data is smaller than all the simulated data is \(1/101\), and the probability that the real data is bigger than only one of the simulated data points is \(1/101\), etcetera.
This is a distribution-free measure of goodness of fit.
Warm-up
To check that we understand what’s going on here, we’re doing to simulate data where the “sample” fits the model, so that the “sample” and the “simulations” are from the same distribution.
# 1. Simulate 99 draws from a distribution, then one more "sample".
x <- runif(99 + 1) # the "sample" will be the first one
# 2. Find the rank of the "sample" in all 100.
r <- sum(x <= x[1])
# 3. Do that a lot of times: we do it in a matrix to be efficient.
nreps <- 1000
X <- matrix(runif(100 * nreps), ncol=100)
R <- rowSums(X <= X[,1])
Chi-squraed for goodness of fit.
If \(N_i \sim \Poisson(P_i)\), then the mean and variance of \(N_i\) is \(P_i\). This motivates the chi-squared statistic: \[ \chi^2 = \sum_i \frac{(N_i - P_i)^2}{P_i} . \]
# 5. Compute the chi-squared statistic.
gof <- function (x, sim) {
R <- rowSums(x <= sim) # this is the rank minus 1
predicted <- length(x) / 10
chisq <- sum((table(R %/% 10) - predicted)^2 / predicted)
return(chisq)
}
# 6. Find the null distribution of this statistic
null_distrn <- replicate(1000, {
nreps <- 1000
X <- matrix(runif(100 * nreps), ncol=100)
gof(X[,1], X[,2:100]) })Goodness-of-fit: implementation
The goodness-of-fitscores obtained by our models were
- Original model: 978
- Second model, with more randomness: 33
Here is the “null distribution” of the goodness-of-fit scores. The first model is clearly outside this range, and the second one is on the edge.
