Bayesian GLMs

web.stanford.edu/class/stats305b

Jonathan Taylor

Winter 2022

Bayesian methods for GLMs (7.2)

Given design matrix \(X\) and binary \(Y\), we want a Bayesian model for the law of \(Y|X\).
Model

\[ \begin{aligned} \beta|X &\sim g \\ Y|X,\beta & \sim \text{Binomial}\left(\frac{\exp(X^T\beta)}{1 + \exp(X^T\beta)}\right) \\ \end{aligned} \]

Log of posterior density is proportional to (suppressing \(X\) as is commonly done)

\[ \begin{aligned} \log ( g(\beta | Y) ) &= \log L(\beta |Y) + \log (g(\beta))\\ &= \sum_{i=1}^nY_iX_i^T\beta - \log(1 + \exp(X_i^T\beta)) + \log (g(\beta)) \end{aligned} \]

Generally is not going to be recognizable as a familiar distribution.

Metropolis-Hastings

Requires a proposal kernel \(K : \mathbb{R}^p \times \mathbb{R}^p \rightarrow \mathbb{R}\) so that \(K(\cdot | \beta)\) is a pdf for each \(\beta\).
Simple proposal: \(\beta' \overset{D}{=} \beta + N(0, \Sigma)\) (or other symmetric density) is

\[ K(\beta'|\beta) = N(\beta,\Sigma) \]

Accept proposal with probability

\[ \min\left(1, \frac{g(\beta'|Y) \cdot K(\beta|\beta')}{g(\beta|Y) \cdot K(\beta'|\beta)} \right) \]

With a symmetric proposal, the ratio is just

\[ \min\left(1, \frac{g(\beta'|Y)}{g(\beta|Y)} \right) \]

Gibbs sampling

Instead of proposing a new \(\mathbb{R}^p\) vector, moves a coordinate at a time.
Define

\[ \begin{aligned} g_j(\beta_j | \beta_{-j},Y) &\propto g(\beta|Y) \\ &= \frac{g(\beta | Y)}{\int_{-\infty}^{\infty} g(\beta|Y) \; d\beta_j} \end{aligned} \]

Each \((\beta_{-j},Y)\) (again \(X\) suppressed) determines a conditional density for \(\beta_j\).
Sometimes (but not always) this density will be “simpler” than \(g(\beta|Y)\). Even possible that this family is from a known distribution, particularly using conjugate priors.
Gibbs sampler repeatedly cycles:
1. Initialize \(\beta^{(0)}\)
2. Pick a coordinate \(1 \leq j \leq p\).
3. Draw \(U^{cur}\) from density \(u \mapsto g_j( u | \beta^{(t-1)}_{-j},Y)\). This is our update for \(j\)-th coordinate of \(\beta_j\).
4. Set \(\beta_j^{(t)} = U^{cur}, \beta_{-j}^{(t)}=\beta_{-j}^{(t-1)}\).
5. Increment \(t += 1\) and return to step 1.
6. Repeat steps 1 and 2 many times, may also randomly choose which coordinate to select.

Threshold model

Recall the latent threshold model for binary data. With (inverse) link \(F\), we model

\[ Y_i = \begin{cases} 1 & T_i - X_i^T\beta \leq 0 \\ 0 & \text{otherwise} \end{cases} \]

with \(T_i \sim F\).

In this case, posterior distribution can be realized using \(T_i\)’s instead of \(Y\)

\[ g(\beta | Y) = g(\beta | T \leq X\beta). \]

MCMC samplers draw samples \((T,\beta) | T \leq X \beta\) with \(T_i|\beta, Y_i \overset{IID}{\sim} F, \beta \sim g\).
Known as data augmentation…

Threshold probit model (7.2.6)

When \(F=\Phi\), i.e. with probit link and \(\beta \sim N(\mu, \Sigma)\), then \((T,\beta)\) are jointly Gaussian.
Posterior sampling reduces to sampling multivariate truncated Normal. A common sampling problem for which Gibbs sampling is quite tractable.

Gibbs-ish sampler for multivariate truncated Normal

Suppose we want to sample from \(Z \sim N(\mu,\Sigma) | AZ \leq b\) for fixed \((A,b)\).
For simplicity, assume \(\Sigma=I\).
What is the Gibbs move for coordinate \(j\)?
The inequalities \(AZ \leq b\), conditional on \(Z_{-j}\), reduce to \[ L(Z_{-j}, A, b) \leq Z_j \leq U(Z_{-j}, A, b) \]
Draw \(Z_j\) from \(N(\mu_j, 1)\) truncated to interval \([L, U]\).
Example of hit and run sampling.

STAN

While these are general purpose sampling schemes (with a huge literature of their own), practically speaking people use software for (at least simple) Bayesian analyses.
Sampling is not Gibbs or MH, but a variant of Hamiltonian MCMC.
Let’s do a simple example in Agresti’s favorite crabs data.

library(glmbb) # for `crabs` data
data(crabs)
M_crabs = glm(y ~ width + color, family=binomial, data=crabs)
X_crabs = model.matrix(M_crabs)
crabs_stan = list(X = X_crabs,
                  N = nrow(X_crabs),
                  P = ncol(X_crabs),
                  y = crabs$y)

Stan file `logit_uninformative.stan`

data {
  int<lower=1> N;
  int<lower=1> P;
  matrix[N,P] X;
  int<lower=0,upper=1> y[N];
}
parameters {
  vector[P] beta;
}
model {
  beta ~ normal(0, 100);
  y ~ bernoulli_logit(X * beta);
}

library(rstan)

## Loading required package: StanHeaders

## Loading required package: ggplot2

## rstan (Version 2.21.3, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)

system.time(samples <- rstan::stan("logit_uninformative.stan", 
                                   data=crabs_stan))

## 
## SAMPLING FOR MODEL 'logit_uninformative' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 4.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.42 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 1.31019 seconds (Warm-up)
## Chain 1:                0.922471 seconds (Sampling)
## Chain 1:                2.23266 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'logit_uninformative' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.3e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.27487 seconds (Warm-up)
## Chain 2:                1.03749 seconds (Sampling)
## Chain 2:                2.31236 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'logit_uninformative' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.6e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 1.30311 seconds (Warm-up)
## Chain 3:                0.997809 seconds (Sampling)
## Chain 3:                2.30091 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'logit_uninformative' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.8e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 1.37414 seconds (Warm-up)
## Chain 4:                1.07209 seconds (Sampling)
## Chain 4:                2.44624 seconds (Total)
## Chain 4:

##    user  system elapsed 
##  31.006   1.524  34.112

samples

## Inference for Stan model: logit_uninformative.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##           mean se_mean   sd    2.5%    25%    50%    75%  97.5% n_eff Rhat
## beta[1] -12.16    0.07 2.78  -17.88 -13.95 -12.08 -10.34  -6.82  1562    1
## beta[2]   0.49    0.00 0.11    0.28   0.42   0.49   0.56   0.72  1588    1
## beta[3]  -1.17    0.01 0.62   -2.41  -1.60  -1.15  -0.75   0.01  2169    1
## beta[4]   0.31    0.02 0.85   -1.27  -0.27   0.29   0.87   2.02  2195    1
## beta[5]   0.28    0.01 0.43   -0.53  -0.01   0.29   0.57   1.13  2108    1
## lp__    -96.32    0.04 1.61 -100.29 -97.14 -96.02 -95.12 -94.18  1461    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Feb  2 11:02:23 2022.
## 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).

samples

## Inference for Stan model: logit_uninformative.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##           mean se_mean   sd    2.5%    25%    50%    75%  97.5% n_eff Rhat
## beta[1] -12.16    0.07 2.78  -17.88 -13.95 -12.08 -10.34  -6.82  1562    1
## beta[2]   0.49    0.00 0.11    0.28   0.42   0.49   0.56   0.72  1588    1
## beta[3]  -1.17    0.01 0.62   -2.41  -1.60  -1.15  -0.75   0.01  2169    1
## beta[4]   0.31    0.02 0.85   -1.27  -0.27   0.29   0.87   2.02  2195    1
## beta[5]   0.28    0.01 0.43   -0.53  -0.01   0.29   0.57   1.13  2108    1
## lp__    -96.32    0.04 1.61 -100.29 -97.14 -96.02 -95.12 -94.18  1461    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Feb  2 11:02:23 2022.
## 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).

summary(M_crabs)

## 
## Call:
## glm(formula = y ~ width + color, family = binomial, data = crabs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1124  -0.9848   0.5243   0.8513   2.1413  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -11.6090     2.7070  -4.289 1.80e-05 ***
## width         0.4680     0.1055   4.434 9.26e-06 ***
## colordarker  -1.1061     0.5921  -1.868   0.0617 .  
## colorlight    0.2238     0.7771   0.288   0.7733    
## colormedium   0.2962     0.4176   0.709   0.4781    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 225.76  on 172  degrees of freedom
## Residual deviance: 187.46  on 168  degrees of freedom
## AIC: 197.46
## 
## Number of Fisher Scoring iterations: 4

confint(M_crabs)

## Waiting for profiling to be done...

##                   2.5 %      97.5 %
## (Intercept) -17.2164739 -6.55237372
## width         0.2712817  0.68704361
## colordarker  -2.3138635  0.02792233
## colorlight   -1.2396603  1.89189591
## colormedium  -0.5325502  1.11213712

post_samples = extract(samples, "beta")$beta
plot(post_samples[,2])

Informative prior

Instead of prior \(\beta \sim N(0, 100 \cdot I)\) use \(N(0,I)\)

samples_i = rstan::stan("logit_informative.stan", 
                         data=crabs_stan)

## 
## SAMPLING FOR MODEL 'logit_informative' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 3.8e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.617123 seconds (Warm-up)
## Chain 1:                0.413415 seconds (Sampling)
## Chain 1:                1.03054 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'logit_informative' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.2e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.681517 seconds (Warm-up)
## Chain 2:                0.340647 seconds (Sampling)
## Chain 2:                1.02216 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'logit_informative' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.3e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.71395 seconds (Warm-up)
## Chain 3:                0.392543 seconds (Sampling)
## Chain 3:                1.10649 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'logit_informative' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.5e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.582875 seconds (Warm-up)
## Chain 4:                0.385174 seconds (Sampling)
## Chain 4:                0.968049 seconds (Total)
## Chain 4:

post_samples_i = extract(samples_i, "beta")$beta

Posterior density for \(\beta_{width}\)

plot(density(post_samples[,2]), ylim=c(0,11))
lines(density(post_samples_i[,2]), col='blue') # with shrinkage
# Bernstein von Mises: using uninformative
# should look Gaussian centered at MLE...
lines(seq(0,1,length=101), 
      dnorm(seq(0,1,length=101), 
            mean=coef(M_crabs)[2], 
            sd=sqrt(vcov(M_crabs)[2,2])), 
      col='red')

mean(post_samples[,2]) # MCMC estimate of posterior mean

## [1] 0.4902418

plot(post_samples[,2], post_samples[,4])

cor(post_samples)

##              [,1]        [,2]        [,3]        [,4]         [,5]
## [1,]  1.000000000 -0.99227389 -0.04083923  0.08692002 -0.003829013
## [2,] -0.992273885  1.00000000 -0.03157052 -0.14152248 -0.097927918
## [3,] -0.040839227 -0.03157052  1.00000000  0.25568340  0.470530071
## [4,]  0.086920023 -0.14152248  0.25568340  1.00000000  0.352136162
## [5,] -0.003829013 -0.09792792  0.47053007  0.35213616  1.000000000

Checking convergence

Auto-correlation function of samples

acf(post_samples[,2])

Point estimates

STAN can find posterior mode / MAP estimator

rstan::optimizing(rstan::stan_model('logit_uninformative.stan'), 
                  data=crabs_stan)

## $par
##     beta[1]     beta[2]     beta[3]     beta[4]     beta[5] 
## -11.6042280   0.4677605  -1.1054456   0.2260108   0.2963959 
## 
## $value
## [1] -93.73533
## 
## $return_code
## [1] 0
## 
## $theta_tilde
##        beta[1]   beta[2]   beta[3]   beta[4]   beta[5]
## [1,] -11.60423 0.4677605 -1.105446 0.2260108 0.2963959

M_crabs$coef

## (Intercept)       width colordarker  colorlight colormedium 
## -11.6089904   0.4679560  -1.1061215   0.2237977   0.2962146