Metric data: robustness and differences of ‘means’
Peter Ralph
5 February 2018 – Advanced Biological Statistics
Overview
Summary
So far we’ve focused on discrete data: coin flips and counts.
While doing that we came across familiar things, like Poisson regression.
This week we’ll look at continuous (i.e., “metric”) data. Different sorts of predictors will lead us to different classical statistics: linear regression, \(t\)-tests, etcetera.
But, control of the underlying model will let us easily get more sophisticated, including for instance robustness to error, and penalization for overdetermined problems.
The basic ingredients
Comparison of means
If the predictor, \(X\), is discrete then we are doing a \(t\)-test, or ANOVA, or something.
Simulate data - difference in means of 3.0:
The \(t\)-test
## user system elapsed
## 0.005 0.000 0.007##
## Welch Two Sample t-test
##
## data: y by x
## t = -44.168, df = 193.19, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.164550 -2.894004
## sample estimates:
## mean in group 0 mean in group 1
## 0.9611205 3.9903971with Stan
stt_block <- "
data {
int N;
vector[N] x; // will be a vector of 0's and 1's
vector[N] y;
}
parameters {
real b0;
real b1;
real<lower=0> sigma;
}
model {
y ~ normal(b0 + b1*x, sigma);
}"
system.time( stantt <- stan(model_code=stt_block,
data=list(N=length(x), x=x, y=y), iter=1e3) )## user system elapsed
## 40.240 0.978 42.754## $summary
## mean se_mean sd 2.5% 25% 50%
## b0 0.9584959 0.0016430080 0.05127725 0.8597816 0.9234022 0.9588988
## b1 3.0320092 0.0022753493 0.06901506 2.8978689 2.9872654 3.0331521
## sigma 0.4863959 0.0006812764 0.02493617 0.4388088 0.4691070 0.4859912
## lp__ 44.0727049 0.0419015047 1.24303922 40.8969349 43.4910938 44.4028637
## 75% 97.5% n_eff Rhat
## b0 0.9914469 1.0635215 974.0250 1.000422
## b1 3.0784012 3.1683381 920.0079 1.000391
## sigma 0.5024615 0.5372778 1339.7164 1.002857
## lp__ 44.9781269 45.4684006 880.0562 1.001635
##
## $c_summary
## , , chains = chain:1
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 0.9561642 0.05454421 0.8437029 0.9217489 0.9589110
## b1 3.0333137 0.07354142 2.8979076 2.9832068 3.0326182
## sigma 0.4887660 0.02508591 0.4407525 0.4719735 0.4878538
## lp__ 43.9842256 1.28819655 40.8466789 43.3886738 44.3590545
## stats
## parameter 75% 97.5%
## b0 0.9922897 1.0620172
## b1 3.0837845 3.1771981
## sigma 0.5044771 0.5389809
## lp__ 44.8811333 45.4592424
##
## , , chains = chain:2
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 0.9587976 0.04814115 0.8637515 0.9255311 0.9596549
## b1 3.0309562 0.06787379 2.9013071 2.9836336 3.0330067
## sigma 0.4883044 0.02411460 0.4417797 0.4722400 0.4892917
## lp__ 44.1644344 1.22483135 40.9472576 43.6595758 44.4510361
## stats
## parameter 75% 97.5%
## b0 0.9883987 1.0561342
## b1 3.0744157 3.1652904
## sigma 0.5024329 0.5374876
## lp__ 45.0591509 45.4795431
##
## , , chains = chain:3
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 0.9600256 0.05338369 0.8602954 0.9226958 0.9552302
## b1 3.0329607 0.06667857 2.9000560 2.9917661 3.0351048
## sigma 0.4838985 0.02503191 0.4359803 0.4664332 0.4831959
## lp__ 44.0170075 1.27496415 40.9000264 43.3256537 44.3764698
## stats
## parameter 75% 97.5%
## b0 0.9983666 1.068852
## b1 3.0732049 3.166966
## sigma 0.5014593 0.535272
## lp__ 45.0050491 45.486117
##
## , , chains = chain:4
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 0.9589962 0.04881229 0.8711463 0.9248723 0.9603675
## b1 3.0308062 0.06793050 2.8955776 2.9887929 3.0333363
## sigma 0.4846146 0.02519660 0.4378644 0.4666779 0.4828211
## lp__ 44.1251522 1.17581446 41.2491122 43.6387616 44.4166274
## stats
## parameter 75% 97.5%
## b0 0.9885910 1.0586022
## b1 3.0787049 3.1610643
## sigma 0.5010647 0.5335561
## lp__ 44.9715461 45.4541174Standard linear regression
Standard linear regression
## user system elapsed
## 0.003 0.000 0.003##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.39045 -0.32827 0.00434 0.33770 1.29782
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.03743 0.03511 29.55 <2e-16 ***
## x 2.99838 0.01114 269.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4964 on 198 degrees of freedom
## Multiple R-squared: 0.9973, Adjusted R-squared: 0.9973
## F-statistic: 7.24e+04 on 1 and 198 DF, p-value: < 2.2e-16with Stan
slr_block <- "
data {
int N;
vector[N] x;
vector[N] y;
}
parameters {
real b0;
real b1;
real<lower=0> sigma;
}
model {
y ~ normal(b0 + b1*x, sigma);
}"
system.time( stanlr <- stan(model_code=slr_block,
data=list(N=length(x), x=x, y=y), iter=1e3) )## user system elapsed
## 0.353 0.021 4.764## $summary
## mean se_mean sd 2.5% 25% 50%
## b0 1.037869 0.0007832024 0.03445802 0.9713526 1.014941 1.0375510
## b1 2.998214 0.0002726515 0.01091294 2.9762261 2.991127 2.9983839
## sigma 0.499058 0.0005992648 0.02528762 0.4519479 0.481801 0.4974935
## lp__ 38.915232 0.0365957577 1.15005479 35.8583462 38.361609 39.2207213
## 75% 97.5% n_eff Rhat
## b0 1.0614043 1.1055692 1935.6758 0.9986115
## b1 3.0058555 3.0190104 1602.0198 0.9990105
## sigma 0.5155453 0.5513456 1780.6489 0.9992986
## lp__ 39.7494696 40.2749670 987.5875 1.0029886
##
## $c_summary
## , , chains = chain:1
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 1.0367225 0.034279995 0.9669801 1.0117994 1.0382511
## b1 2.9979769 0.009875374 2.9770437 2.9919000 2.9980638
## sigma 0.4990248 0.024079538 0.4564641 0.4826994 0.4972312
## lp__ 39.0575615 1.022991990 36.6202467 38.5023542 39.3546764
## stats
## parameter 75% 97.5%
## b0 1.0609250 1.0998919
## b1 3.0047203 3.0161260
## sigma 0.5143986 0.5517112
## lp__ 39.8186780 40.2877270
##
## , , chains = chain:2
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 1.0373640 0.03561805 0.9728951 1.0134970 1.0373803
## b1 2.9980877 0.01138842 2.9769848 2.9901239 2.9980656
## sigma 0.4979226 0.02372709 0.4550867 0.4824173 0.4958306
## lp__ 38.9026250 1.17136079 35.9881874 38.4118545 39.2227635
## stats
## parameter 75% 97.5%
## b0 1.0610916 1.1057651
## b1 3.0066883 3.0183298
## sigma 0.5125704 0.5490354
## lp__ 39.7453365 40.2690270
##
## , , chains = chain:3
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 1.0383469 0.03298995 0.9744880 1.0163542 1.0374661
## b1 2.9985349 0.01177823 2.9757089 2.9903582 2.9984920
## sigma 0.4998314 0.02606133 0.4557211 0.4811195 0.4980518
## lp__ 38.8550151 1.18167025 35.7202653 38.3079196 39.0694001
## stats
## parameter 75% 97.5%
## b0 1.0599133 1.1053201
## b1 3.0068983 3.0215172
## sigma 0.5174916 0.5513328
## lp__ 39.7534883 40.2581306
##
## , , chains = chain:4
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## b0 1.0390416 0.03494682 0.9708443 1.0154732 1.037488 1.0639421
## b1 2.9982554 0.01053387 2.9766606 2.9917165 2.998576 3.0052957
## sigma 0.4994534 0.02715986 0.4467995 0.4794095 0.499711 0.5181878
## lp__ 38.8457266 1.20606007 35.6770810 38.3080068 39.131457 39.7021003
## stats
## parameter 97.5%
## b0 1.1060811
## b1 3.0179989
## sigma 0.5519271
## lp__ 40.2863087Multivariate linear regrssion
Simulate data:
truth <- list(b0=1.0,
b=c(3.0, -1.0, 0.0, 0.0),
sigma=0.5)
n <- 200
x <- matrix(rnorm(4*n, mean=0, sd=3), ncol=4)
y <- truth$b0 + x %*% truth$b + rnorm(n, mean=0, sd=truth$sigma)Your turn: compare standard multivariate regression to Stan.
Note: %*% in R, and * in Stan, are matrix multiplication.
Standard linear regression
## user system elapsed
## 0.002 0.000 0.001##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.41425 -0.31589 0.01571 0.33286 1.74978
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.964532 0.037556 25.683 <2e-16 ***
## x1 3.013153 0.012082 249.402 <2e-16 ***
## x2 -0.981735 0.013644 -71.956 <2e-16 ***
## x3 0.001413 0.013201 0.107 0.915
## x4 0.005444 0.011416 0.477 0.634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5269 on 195 degrees of freedom
## Multiple R-squared: 0.997, Adjusted R-squared: 0.9969
## F-statistic: 1.614e+04 on 4 and 195 DF, p-value: < 2.2e-16with Stan
mlr_block <- "
data {
int N;
matrix[N,4] x;
vector[N] y;
}
parameters {
real b0;
vector[4] b1;
real<lower=0> sigma;
}
model {
y ~ normal(to_vector(b0 + x * b1), sigma);
b0 ~ normal(0, 10);
b1 ~ normal(0, 10);
sigma ~ normal(0, 10);
}"
system.time( stanmlr <- stan(model_code=mlr_block,
data=list(N=nrow(x), x=x,
y=as.vector(y)), iter=1e3) )## user system elapsed
## 36.700 0.754 39.127## $summary
## mean se_mean sd 2.5% 25%
## b0 0.964168876 0.0008609641 0.03850349 0.88884240 0.939406253
## b1[1] 3.013312508 0.0002716845 0.01215010 2.99088002 3.004693625
## b1[2] -0.981814139 0.0003059612 0.01368300 -1.00843766 -0.990836467
## b1[3] 0.001582385 0.0002901135 0.01297427 -0.02414813 -0.007005947
## b1[4] 0.004641061 0.0002539090 0.01135515 -0.01752785 -0.003051936
## sigma 0.529757796 0.0005968461 0.02669177 0.48294737 0.510904188
## lp__ 26.980574920 0.0591634454 1.73877680 22.52377782 26.065571894
## 50% 75% 97.5% n_eff Rhat
## b0 0.963721191 0.98890566 1.04202257 2000.0000 0.9998817
## b1[1] 3.013401938 3.02162108 3.03658448 2000.0000 0.9995003
## b1[2] -0.981705463 -0.97276853 -0.95519086 2000.0000 0.9998914
## b1[3] 0.001302625 0.01037515 0.02706174 2000.0000 0.9990247
## b1[4] 0.004710580 0.01223199 0.02721156 2000.0000 0.9985843
## sigma 0.528650061 0.54710170 0.58528911 2000.0000 1.0021144
## lp__ 27.275462315 28.26884818 29.30228460 863.7355 1.0021258
##
## $c_summary
## , , chains = chain:1
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 0.965489164 0.03704322 0.89244744 0.941751021 0.965325283
## b1[1] 3.013113234 0.01218186 2.99019878 3.003942192 3.012948908
## b1[2] -0.980972706 0.01429928 -1.00831464 -0.989874077 -0.981252601
## b1[3] 0.001954729 0.01326935 -0.02370373 -0.006891810 0.001518715
## b1[4] 0.004267654 0.01248824 -0.02003368 -0.002961459 0.004026469
## sigma 0.530161630 0.02618649 0.48441958 0.511656124 0.529157268
## lp__ 26.892122057 1.86887860 22.29523081 26.018117229 27.271768512
## stats
## parameter 75% 97.5%
## b0 0.99094337 1.04145311
## b1[1] 3.02167827 3.03497667
## b1[2] -0.97236903 -0.95061024
## b1[3] 0.01051886 0.02815576
## b1[4] 0.01167012 0.02887848
## sigma 0.54668856 0.58473576
## lp__ 28.24984013 29.31057866
##
## , , chains = chain:2
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 0.960877330 0.04094817 0.88499478 0.933355962 0.9586259404
## b1[1] 3.013851692 0.01178230 2.99257299 3.005174650 3.0142191821
## b1[2] -0.982426095 0.01318563 -1.00714283 -0.991586586 -0.9817596634
## b1[3] 0.001354411 0.01353855 -0.02427455 -0.007713342 0.0008937733
## b1[4] 0.005093480 0.01076634 -0.01455499 -0.001942040 0.0049601497
## sigma 0.529517950 0.02667113 0.48416906 0.510762964 0.5290034111
## lp__ 26.972577521 1.76904808 22.47286797 26.063578522 27.2688494951
## stats
## parameter 75% 97.5%
## b0 0.98649739 1.04059419
## b1[1] 3.02175558 3.03775450
## b1[2] -0.97374190 -0.95573130
## b1[3] 0.01066009 0.02654051
## b1[4] 0.01240653 0.02722870
## sigma 0.54806488 0.58231373
## lp__ 28.28041391 29.27620668
##
## , , chains = chain:3
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 0.9648426609 0.03819984 0.88914973 0.941831465 0.9636390260
## b1[1] 3.0130886611 0.01230454 2.99055062 3.005078489 3.0129760251
## b1[2] -0.9812709921 0.01372253 -1.00814500 -0.990570389 -0.9812624414
## b1[3] 0.0009627094 0.01222275 -0.02324193 -0.007495341 0.0007426121
## b1[4] 0.0046265670 0.01131579 -0.01715998 -0.003408141 0.0046722917
## sigma 0.5322338608 0.02586258 0.48689854 0.513296973 0.5303819103
## lp__ 27.0633123377 1.73410399 22.56749594 26.152443961 27.3821452402
## stats
## parameter 75% 97.5%
## b0 0.987576156 1.04327983
## b1[1] 3.021585263 3.03674708
## b1[2] -0.971687247 -0.95550022
## b1[3] 0.009788693 0.02534619
## b1[4] 0.012226587 0.02638471
## sigma 0.550021541 0.58476783
## lp__ 28.351457029 29.35447666
##
## , , chains = chain:4
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 0.965466349 0.03763059 0.89006349 0.940383424 0.966373290
## b1[1] 3.013196445 0.01234384 2.98860738 3.004647147 3.013449878
## b1[2] -0.982586764 0.01346876 -1.00882031 -0.992128027 -0.982204177
## b1[3] 0.002057692 0.01283640 -0.02349400 -0.005717890 0.001772736
## b1[4] 0.004576544 0.01078270 -0.01722273 -0.003643438 0.004868278
## sigma 0.527117742 0.02783492 0.47745834 0.507373992 0.526057944
## lp__ 26.994287764 1.57074798 23.49262217 26.053048026 27.190228323
## stats
## parameter 75% 97.5%
## b0 0.99085690 1.04217216
## b1[1] 3.02093133 3.03652066
## b1[2] -0.97302255 -0.95851713
## b1[3] 0.01081647 0.02662559
## b1[4] 0.01237443 0.02409439
## sigma 0.54459562 0.59317977
## lp__ 28.21249435 29.28008931Testing with crossvalidation
Crossvalidation
Just because a model fits doesn’t mean that it’s any good.
Divide your data randomly into 5 pieces.
Fit your model on 4/5ths, and see how well it predicts the remaining 1/5th.
Do this for each of the 5 pieces.
A better model should have better crossvalidation accuracy.
Question: for linear regression, how do we “see how well it predicts”? Write down the math, then code it up.
five_groups <- sample(rep(1:5, each=nrow(x)/5))
do_crossval <- function (k) {
# fit using not group k
subset_lm <- lm( y ~ x, subset=(five_groups != k))
# predict on group k
yhat <- predict(subset_lm,
newdata=list(x=x[five_groups == k,]))
# compute crossvalidation accuracy
S <- sqrt( mean( ( y[five_groups == k] - yhat )^2 ) )
return(S)
}
crossvals <- sapply(1:5, do_crossval)The mean crossvalidation score for ordinary linear regression was 0.5403888, with a standard deviation of 0.070462.
What is overfitting?
Even completely independent sets of numbers will correlate a little, because of noise.
When you have a lot of variables, there may be some that correlate well with the response variable just by chance.
If you have as many variables than observations, then there is (almost) always a linear model that fits perfectly (with \(\epsilon = 0\)).
Example:
nvars <- 200
truth <- list(b0=1.0,
b=rep(0.0, nvars),
sigma=0.5)
n <- 200
x <- matrix(rnorm(nvars*n, mean=0, sd=3), ncol=nvars)
y <- truth$b0 + x %*% truth$b + rnorm(n, mean=0, sd=truth$sigma)
the_lm <- lm(y ~ x)
range(predict(the_lm) - y)## [1] -5.750955e-14 6.311618e-14Stochastic minute
If \(X \sim \Cauchy(\text{center}=\mu, \text{scale}=\sigma)\), then \(X\) has probability density \[\begin{aligned} f(x \given \mu, \sigma) = \frac{1}{\pi\left( 1 + \left( \frac{x - \mu}{\sigma} \right)^2 \right)} . \end{aligned}\]
The Cauchy is a good example of a distribution with “heavy tails”: rare, very large values.
If \(Y\) and \(Z\) are independent \(\Normal(0,1)\) then \(Y/Z \sim \Cauchy(0,1)\).
If \(X_1, X_2, \ldots, X_n\) are independent \(\Cauchy(0,1)\) then \(\max(X_1, \ldots, X_n)\) is of size \(n\).
Also, \(\frac{1}{n}(X_1 + \ldots + X_n) \sim \Cauchy(0,1)\).
- If \(X_1, X_2, \ldots, X_n\) are independent \(\Cauchy(0,1)\) then \[\begin{aligned} \frac{1}{n} \left(X_1 + \cdots + X_n\right) \sim \Cauchy(0,1) . \end{aligned}\]
Wait, what?!?
A single value has the same distribution as the mean of 1,000 of them?
Let’s look:
\(X \sim \Normal(0,1)\)
\(X \sim \Cauchy(0,1)\)
Another way to look at it: extreme values
n <- 100
plot(c(cummax(rcauchy(n))), type='l', ylab='max value so far', xlab='number of samples', col='red')
lines(c(cummax(rnorm(n))), col='black')
legend("bottomright", lty=1, col=c('black', 'red'), legend=c('normal', 'cauchy'))
Another way to look at it: extreme values
n <- 1000
plot(c(cummax(rcauchy(n))), type='l', ylab='max value so far', xlab='number of samples', col='red')
lines(c(cummax(rnorm(n))), col='black')
legend("bottomright", lty=1, col=c('black', 'red'), legend=c('normal', 'cauchy'))
Another way to look at it: extreme values
n <- 1e6
plot(c(cummax(rcauchy(n))), type='l', ylab='max value so far', xlab='number of samples', col='red')
lines(c(cummax(rnorm(n))), col='black')
legend("bottomright", lty=1, col=c('black', 'red'), legend=c('normal', 'cauchy'))
Problem #1: noise
Cauchy noise
Let’s see what happens if \[\begin{aligned} Y_i &= b_0 + b_1 X_1 + \epsilon_i \\ \epsilon_i &\sim \Cauchy(0, \sigma) . \end{aligned}\]
The \(t\)-test
## user system elapsed
## 0.005 0.000 0.005##
## Welch Two Sample t-test
##
## data: y by x
## t = -0.5875, df = 904.03, p-value = 0.557
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.779966 1.499041
## sample estimates:
## mean in group 0 mean in group 1
## 2.537384 3.177847with Stan
sct_block <- "
data {
int N;
vector[N] x; // will be a vector of 0's and 1's
vector[N] y;
}
parameters {
real b0;
real b1;
real<lower=0> sigma;
}
model {
y ~ cauchy(b0 + b1*x, sigma);
}"
system.time( stanct <- stan(model_code=sct_block,
data=list(N=length(x), x=x, y=y), iter=1e3) )## user system elapsed
## 36.315 0.833 42.321## $summary
## mean se_mean sd 2.5% 25%
## b0 1.0291922 0.0011314808 0.03548784 0.960985 1.0044436
## b1 2.9313917 0.0016139906 0.04953027 2.833088 2.8983148
## sigma 0.5246031 0.0006601783 0.02314277 0.481779 0.5086437
## lp__ -728.2095555 0.0491822653 1.29123990 -731.396073 -728.8063460
## 50% 75% 97.5% n_eff Rhat
## b0 1.0295111 1.0533330 1.0994354 983.7046 0.9991930
## b1 2.9325418 2.9645720 3.0299830 941.7582 0.9997976
## sigma 0.5245932 0.5401948 0.5702207 1228.8759 1.0033877
## lp__ -727.8744378 -727.2580625 -726.7144310 689.2818 1.0019354
##
## $c_summary
## , , chains = chain:1
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 1.0296636 0.03459815 0.9696624 1.0038920 1.0285112
## b1 2.9317775 0.04820848 2.8307235 2.9032111 2.9355758
## sigma 0.5271417 0.02222283 0.4881341 0.5122751 0.5252849
## lp__ -728.1911557 1.29963400 -731.4383352 -728.7771593 -727.8187576
## stats
## parameter 75% 97.5%
## b0 1.0526276 1.0937105
## b1 2.9618142 3.0210954
## sigma 0.5409076 0.5757164
## lp__ -727.2403446 -726.7335290
##
## , , chains = chain:2
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 1.0300621 0.03471651 0.9617799 1.0053069 1.0320335
## b1 2.9287701 0.04875061 2.8307552 2.8946702 2.9309229
## sigma 0.5252629 0.02262749 0.4834991 0.5093539 0.5248081
## lp__ -728.1085467 1.26765283 -731.3837532 -728.6215158 -727.7928074
## stats
## parameter 75% 97.5%
## b0 1.0538801 1.1024453
## b1 2.9623883 3.0256776
## sigma 0.5411029 0.5693632
## lp__ -727.2401254 -726.7140901
##
## , , chains = chain:3
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 1.0273191 0.03512203 0.9611255 1.0044036 1.0282320
## b1 2.9317931 0.04885671 2.8355257 2.8966320 2.9341094
## sigma 0.5217075 0.02383327 0.4763099 0.5053075 0.5222038
## lp__ -728.1931195 1.24596210 -731.1086065 -728.8880672 -727.8999648
## stats
## parameter 75% 97.5%
## b0 1.051544 1.0955920
## b1 2.965257 3.0294150
## sigma 0.537420 0.5680276
## lp__ -727.275904 -726.7089190
##
## , , chains = chain:4
##
## stats
## parameter mean sd 2.5% 25% 50%
## b0 1.0297237 0.03747682 0.9515017 1.0037207 1.0294222
## b1 2.9332261 0.05224475 2.8352613 2.8997922 2.9281443
## sigma 0.5243002 0.02358604 0.4819830 0.5070324 0.5253618
## lp__ -728.3453999 1.34219529 -731.6390165 -729.0326680 -728.0149964
## stats
## parameter 75% 97.5%
## b0 1.0561818 1.1026390
## b1 2.9688785 3.0367147
## sigma 0.5401849 0.5708525
## lp__ -727.3210549 -726.6999041Because we correctly model the noise, using Stan, we are not thrown off by large values.
Robust ANOVA
Metric data from groups
Suppose we have numerical observations coming from \(m\) different groups, and want to know if the means are different between groups, and by how much.
e.g., dry leaf mass after 10d of growth in 6 different conditions
dry leaf mass after 10d of growth in 6 different conditions
Each group has a different location and scale; noise about these is Cauchy; group locations are random deviations; group scale parameters are random deviations from a common scale.
\[\begin{aligned} Y_i &\sim \Cauchy(b_{g_i}, \sigma_{g_i}) \\ \sigma_g &\sim \log\Normal(\mu, 1) \\ b_g &\sim \Normal(0, \eta) \\ \eta &\sim \Normal(0, 15) \\ \mu &\sim \Normal(0, 1) \end{aligned}\]
\[\begin{aligned} Y_i &\sim \Cauchy(b_{g_i}, \sigma_{g_i}) \\ \sigma_g &\sim \log\Normal(\mu, 1) \\ b_g &\sim \Normal(z, \eta) \\ z &\sim \Normal(0, 10) \\ \eta &\sim \Normal(0, 15) \\ \mu &\sim \Exp(1) \end{aligned}\]
Changes:
Added common mean, \(z\).
Put “shrinkage” prior on \(\mu\).
anova_block <- "
data {
int N; // number of obs
int m; // number of groups
vector[N] y; // leaf mass
int g[N]; // group ID
}
parameters {
vector[m] b; // group 'means'
real z; // overall mean
// group 'SD's
vector<lower=0>[m] sigma;
// variability in group means
real<lower=0> eta;
// log-mean of sigma
real mu;
}
model {
y ~ cauchy(b[g], sigma[g]);
// lognormal, encouraging all group sigmas to be the same,
// but trying not to constrain them if they want to be different
sigma ~ lognormal(mu, 1);
// unsure of scale for b,
// so put a hyperprior on it
b ~ normal(z, eta);
z ~ normal(0, 15);
eta ~ normal(0, 15);
// encourage mu to be near zero if reasonable
mu ~ exponential(1);
}"
## note we subtract the median of y
anova_data <- list(N=length(y),
m=length(unique(x)),
y=y - median(y),
g=x)
anova_fit <- stan(model_code=anova_block,
data=anova_data,
iter=1e3,
control=list(adapt_delta=0.99))## Warning: There were 8 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: Examine the pairs() plot to diagnose sampling problems
## $summary
## mean se_mean sd 2.5% 25%
## b[1] 2.0004977 0.060226283 2.4055064 -2.548039e+00 0.5329497
## b[2] 3.8144570 0.020915938 0.9353892 2.066095e+00 3.1516557
## b[3] -5.7398166 0.030528892 1.3652936 -8.577400e+00 -6.5666481
## b[4] -6.0629232 0.113239374 3.8691468 -1.146395e+01 -9.0973767
## b[5] 2.6730416 0.073261829 2.6358666 -3.047754e+00 1.2429469
## b[6] 7.3137300 0.037505761 1.4435569 4.030945e+00 6.5016253
## b[7] -1.9352014 0.026194048 0.8962760 -3.482414e+00 -2.4739148
## b[8] -6.7625511 0.026930714 1.0678358 -8.534017e+00 -7.4462092
## b[9] -9.3529212 0.016674325 0.7456985 -1.076861e+01 -9.8430554
## b[10] -0.9128103 0.031145738 1.2559245 -3.471611e+00 -1.6279343
## b[11] -0.7444011 0.036753731 1.3239965 -3.608819e+00 -1.4710115
## b[12] -5.1873019 0.023807518 0.7926762 -6.751619e+00 -5.6460311
## b[13] 2.6825066 0.022497768 0.8931097 6.314364e-01 2.1892071
## b[14] 6.6480984 0.082148658 1.8631067 1.888808e+00 6.3080801
## b[15] 6.7591424 0.018793675 0.7710591 5.299906e+00 6.2701552
## b[16] -1.2011538 0.039757111 1.5286743 -4.259362e+00 -2.1498786
## b[17] 1.6801114 0.024140282 0.8950193 1.348572e-01 1.0918949
## b[18] 5.3025488 0.021110648 0.9440969 3.232415e+00 4.7801497
## b[19] -1.2398895 0.021077778 0.9426269 -3.156877e+00 -1.8406593
## b[20] 3.0965591 0.046225322 1.4624299 1.106025e-03 2.2716071
## z 0.1980326 0.028923300 1.2934893 -2.370657e+00 -0.6344199
## sigma[1] 6.7592524 0.065289636 2.9198413 2.633624e+00 4.6987140
## sigma[2] 1.7238744 0.019495772 0.7889384 6.991711e-01 1.1826245
## sigma[3] 3.6843377 0.027070229 1.2106175 1.827846e+00 2.8446029
## sigma[4] 5.8745465 0.097619356 3.6738025 1.371243e+00 3.4116733
## sigma[5] 6.8143446 0.120401297 4.2525121 1.915057e+00 3.7804684
## sigma[6] 2.2300302 0.026557779 1.0971283 7.607752e-01 1.4594705
## sigma[7] 1.8058863 0.030570861 1.0300583 5.323175e-01 1.1044089
## sigma[8] 2.3876930 0.025597703 1.1447641 8.645015e-01 1.5827505
## sigma[9] 1.9491344 0.014723320 0.6584469 9.432482e-01 1.4563038
## sigma[10] 2.7461600 0.045705377 1.6757890 8.587729e-01 1.6454605
## sigma[11] 2.8846814 0.031913402 1.4272107 1.040292e+00 1.9280316
## sigma[12] 1.8356603 0.026778170 0.9339195 6.702639e-01 1.1715258
## sigma[13] 1.9138968 0.024162875 0.9487308 6.743648e-01 1.2724416
## sigma[14] 1.6457907 0.082139713 2.0317448 1.867137e-01 0.5347309
## sigma[15] 2.0410611 0.016453656 0.7358298 9.389705e-01 1.5109605
## sigma[16] 4.6449526 0.044763812 2.0018985 1.966467e+00 3.2362394
## sigma[17] 2.0591803 0.024054311 0.9756655 7.397800e-01 1.3800453
## sigma[18] 2.2680606 0.022241730 0.9946804 9.099254e-01 1.5861674
## sigma[19] 2.1196747 0.019300087 0.8369106 9.327253e-01 1.5290413
## sigma[20] 2.4411677 0.033763920 1.2816655 8.159465e-01 1.5625834
## eta 5.4681789 0.027580363 1.0167004 3.846221e+00 4.7664420
## mu 0.8123159 0.006326862 0.2560255 3.210090e-01 0.6408437
## lp__ -514.5477186 0.209934288 5.2435538 -5.262728e+02 -517.9062655
## 50% 75% 97.5% n_eff Rhat
## b[1] 1.9374549 3.30142628 7.6021465 1595.2947 1.0002656
## b[2] 3.8108844 4.44782582 5.6447755 2000.0000 0.9994052
## b[3] -5.6963034 -4.87163472 -3.1231710 2000.0000 1.0008704
## b[4] -6.7569869 -3.68865706 2.4926365 1167.4425 1.0010146
## b[5] 2.7360219 4.20312675 8.1382729 1294.4681 1.0000534
## b[6] 7.4657472 8.32873044 9.6527558 1481.3982 0.9995790
## b[7] -2.0623853 -1.46335695 0.1102782 1170.7881 1.0007709
## b[8] -6.8879211 -6.15632066 -4.3667455 1572.2195 1.0008960
## b[9] -9.3910884 -8.87684560 -7.7592761 2000.0000 0.9993227
## b[10] -0.9402218 -0.17524705 1.7090401 1626.0346 1.0010109
## b[11] -0.6337882 0.07354723 1.5995252 1297.6897 1.0017385
## b[12] -5.2001326 -4.74999217 -3.6082229 1108.5706 0.9986734
## b[13] 2.7671614 3.26385210 4.2232260 1575.9077 1.0004246
## b[14] 6.9776489 7.38317924 9.3457845 514.3687 1.0077838
## b[15] 6.6982013 7.20502140 8.4121039 1683.2628 1.0003745
## b[16] -1.1906708 -0.21954563 1.8608588 1478.4285 0.9996612
## b[17] 1.5592458 2.21110085 3.7589153 1374.6120 1.0043429
## b[18] 5.3803251 5.91453708 7.0151613 2000.0000 0.9991106
## b[19] -1.1989284 -0.56860462 0.4844041 2000.0000 1.0001292
## b[20] 3.1160869 4.01238469 5.7625231 1000.8989 1.0009716
## z 0.1840209 1.05019140 2.7329651 2000.0000 0.9997611
## sigma[1] 6.3054520 8.20112088 13.9468963 2000.0000 1.0002093
## sigma[2] 1.5645541 2.12190166 3.5689688 1637.5906 0.9988775
## sigma[3] 3.5048337 4.33448789 6.6791108 2000.0000 0.9983651
## sigma[4] 5.0879729 7.46297633 14.4521364 1416.3146 1.0017588
## sigma[5] 5.6732000 8.79229985 17.4543776 1247.4662 0.9986339
## sigma[6] 1.9881258 2.73664997 4.9750285 1706.5983 0.9988829
## sigma[7] 1.5818598 2.21425079 4.4218746 1135.2937 1.0034177
## sigma[8] 2.1664996 2.92170086 5.2175857 2000.0000 0.9995263
## sigma[9] 1.8740519 2.29354331 3.4855117 2000.0000 1.0015570
## sigma[10] 2.2980578 3.35093212 7.0928144 1344.3244 1.0024416
## sigma[11] 2.5808759 3.53640301 6.4708878 2000.0000 0.9996793
## sigma[12] 1.6280221 2.25504667 4.3304266 1216.3459 0.9996215
## sigma[13] 1.7256902 2.33576459 4.2780385 1541.6606 0.9997241
## sigma[14] 0.9604672 1.99242508 6.9312344 611.8317 1.0065147
## sigma[15] 1.9361880 2.43185165 3.7498887 2000.0000 1.0020214
## sigma[16] 4.2698982 5.59210056 9.5455152 2000.0000 1.0019560
## sigma[17] 1.8664669 2.54034919 4.3180990 1645.1899 1.0015402
## sigma[18] 2.0900833 2.73847725 4.8941140 2000.0000 0.9989870
## sigma[19] 1.9762718 2.53953596 4.1536794 1880.3541 1.0004673
## sigma[20] 2.1728968 3.01403200 5.6308545 1440.9327 1.0018720
## eta 5.3259840 6.06059398 7.8726263 1358.8954 0.9994464
## mu 0.8115207 0.98114637 1.3095732 1637.5313 0.9996643
## lp__ -514.1804552 -510.82354786 -505.4932146 623.8565 1.0062171
##
## $c_summary
## , , chains = chain:1
##
## stats
## parameter mean sd 2.5% 25% 50%
## b[1] 2.0328337 2.3551052 -2.4358767 0.4699576 2.0437588
## b[2] 3.8200824 0.8889261 2.1836775 3.2624002 3.8348619
## b[3] -5.6940812 1.3165047 -8.4505179 -6.4620288 -5.6416212
## b[4] -5.9589489 3.9755278 -11.5202637 -9.0083960 -6.7448875
## b[5] 2.7089811 2.7008774 -2.8595678 1.2528529 2.7253266
## b[6] 7.3839436 1.3832035 4.0080298 6.5477630 7.5654119
## b[7] -1.9500774 0.9268734 -3.5970750 -2.4417340 -2.0462459
## b[8] -6.7861595 1.1002842 -8.5653843 -7.5484210 -6.9432811
## b[9] -9.3651755 0.7390296 -10.7725458 -9.8240939 -9.3983840
## b[10] -0.8484606 1.3256987 -3.6376498 -1.6457248 -0.8795872
## b[11] -0.7729124 1.3458731 -3.7928678 -1.5899345 -0.6642670
## b[12] -5.1921257 0.9355172 -6.8861889 -5.7810239 -5.2764717
## b[13] 2.6703621 0.8870186 0.6319409 2.1511995 2.7684584
## b[14] 6.4792372 2.1312935 1.0481047 6.3006906 6.9741965
## b[15] 6.7798336 0.7390957 5.4620110 6.2544000 6.7099665
## b[16] -1.2248053 1.4615902 -4.2469123 -2.1690054 -1.1948179
## b[17] 1.7911492 1.0119049 0.2065372 1.1153786 1.6619529
## b[18] 5.3008257 0.9199658 3.3848819 4.7596542 5.3893675
## b[19] -1.2080730 0.8570986 -2.9075219 -1.7422671 -1.1812023
## b[20] 3.1713555 1.4105407 0.2572616 2.3507805 3.1449793
## z 0.1619997 1.3182416 -2.3403580 -0.6627080 0.1632631
## sigma[1] 6.7275080 2.9573108 2.5730636 4.6292229 6.2070832
## sigma[2] 1.7058989 0.7130230 0.7225411 1.1925718 1.5627038
## sigma[3] 3.6557134 1.2390603 1.9216094 2.7481876 3.3861014
## sigma[4] 5.8321879 3.2344121 1.4099061 3.4797750 5.1465407
## sigma[5] 6.9927125 4.4869762 1.8502003 3.7131857 5.7006677
## sigma[6] 2.1908199 0.9976401 0.7764749 1.4887165 1.9973617
## sigma[7] 1.8026329 0.9920322 0.5467018 1.1298715 1.6032180
## sigma[8] 2.4724619 1.2230874 0.8081620 1.5862196 2.2052066
## sigma[9] 1.9222065 0.6063502 0.9695710 1.4537111 1.8889957
## sigma[10] 2.8951642 2.0290836 0.7127981 1.5971079 2.3920787
## sigma[11] 2.9539038 1.5283990 0.9778602 1.9192530 2.6293947
## sigma[12] 1.8399284 0.9177647 0.6648495 1.2215913 1.6172720
## sigma[13] 1.9336464 0.9712836 0.6632813 1.2373772 1.7388733
## sigma[14] 1.6325242 1.9582402 0.1766534 0.4908258 0.9413603
## sigma[15] 2.0293959 0.6877194 0.9692958 1.5370337 1.9232676
## sigma[16] 4.5288552 1.8647100 2.0137488 3.1810038 4.2040870
## sigma[17] 2.1189713 1.0442071 0.7524678 1.3981614 1.8818776
## sigma[18] 2.2715464 1.0162099 0.8672859 1.5958611 2.0707446
## sigma[19] 2.1514588 0.8804499 0.9196101 1.5315339 2.0006615
## sigma[20] 2.4740094 1.3114944 0.8220980 1.6042636 2.2072213
## eta 5.4971067 0.9576591 3.9249795 4.8027596 5.3710264
## mu 0.8094879 0.2561303 0.2777192 0.6419005 0.8070607
## lp__ -514.9778600 5.3017840 -526.4015096 -518.2497560 -514.3841990
## stats
## parameter 75% 97.5%
## b[1] 3.43747613 7.49782608
## b[2] 4.39150396 5.58275047
## b[3] -4.84186517 -3.24241094
## b[4] -3.49056460 2.66409791
## b[5] 4.31239738 8.27368598
## b[6] 8.39382867 9.49677356
## b[7] -1.46594071 -0.09586154
## b[8] -6.20187484 -4.16416834
## b[9] -8.94093792 -7.89129450
## b[10] -0.02723903 1.84007895
## b[11] 0.02251390 1.70438313
## b[12] -4.73246807 -3.41647988
## b[13] 3.28446754 4.30449479
## b[14] 7.38454769 9.00321292
## b[15] 7.26312226 8.25949094
## b[16] -0.27877346 1.79322363
## b[17] 2.34554473 4.07381866
## b[18] 5.91538670 6.92303495
## b[19] -0.59332927 0.28199010
## b[20] 4.09586295 5.86948288
## z 1.07349983 2.73444968
## sigma[1] 8.14807220 14.28830159
## sigma[2] 2.10827255 3.43783121
## sigma[3] 4.31192188 6.62891040
## sigma[4] 7.43204192 13.71099631
## sigma[5] 8.88122063 18.39688853
## sigma[6] 2.69994783 4.61834389
## sigma[7] 2.19069095 4.51634986
## sigma[8] 3.10168507 5.38929922
## sigma[9] 2.26752509 3.29376865
## sigma[10] 3.49293338 8.63586545
## sigma[11] 3.70829065 6.66842752
## sigma[12] 2.24420103 4.40330058
## sigma[13] 2.47236510 3.90571045
## sigma[14] 1.92366952 7.62755966
## sigma[15] 2.38741960 3.72876710
## sigma[16] 5.56674667 8.63968815
## sigma[17] 2.62804342 4.70028156
## sigma[18] 2.73016189 4.92785409
## sigma[19] 2.58424270 4.26622628
## sigma[20] 2.95582930 5.61860780
## eta 6.09737467 7.43442190
## mu 0.97359744 1.29630218
## lp__ -511.22419206 -506.22522632
##
## , , chains = chain:2
##
## stats
## parameter mean sd 2.5% 25% 50%
## b[1] 1.9969441 2.3075134 -2.3445752 0.6203264 1.8725834
## b[2] 3.8479973 0.9476985 2.2125113 3.1539469 3.8397003
## b[3] -5.7500117 1.3858451 -8.9407923 -6.5580871 -5.6100789
## b[4] -5.9234311 3.8632466 -11.3869300 -8.9880723 -6.5367612
## b[5] 2.6662983 2.6179892 -3.0312986 1.3516699 2.7015137
## b[6] 7.2924815 1.4061182 4.2248337 6.4725155 7.4643387
## b[7] -1.9085329 0.9325294 -3.2947648 -2.4889362 -2.1299060
## b[8] -6.7842030 1.0559300 -8.4571778 -7.4266516 -6.9109198
## b[9] -9.3099387 0.7639253 -10.7802058 -9.8067107 -9.3384813
## b[10] -0.8828697 1.3118418 -3.5259611 -1.5956960 -0.8970007
## b[11] -0.6798157 1.2939600 -3.3985567 -1.4396333 -0.6832903
## b[12] -5.1732732 0.7343941 -6.4857937 -5.6108401 -5.1829617
## b[13] 2.6629862 0.8538398 0.6625243 2.2105956 2.7637611
## b[14] 6.6261333 1.9575666 1.7506848 6.0327720 6.9256158
## b[15] 6.7962029 0.8296134 5.1934128 6.2995338 6.7464349
## b[16] -1.1779690 1.5158027 -4.2138455 -2.0847384 -1.1386707
## b[17] 1.6424071 0.8306390 0.1046742 1.0741156 1.5320774
## b[18] 5.2970062 0.9079774 3.2264158 4.7605604 5.3585823
## b[19] -1.3023842 0.9512427 -3.3062255 -1.8973540 -1.2464402
## b[20] 2.9726763 1.7422636 -1.3221462 2.1353728 3.0760897
## z 0.2716038 1.2622824 -2.1668716 -0.5219997 0.2352290
## sigma[1] 6.8097286 2.9287797 2.6571521 4.6115748 6.5931938
## sigma[2] 1.7035745 0.9143979 0.6683924 1.1287519 1.5284213
## sigma[3] 3.7122603 1.2079753 1.8670808 2.8789163 3.5032702
## sigma[4] 6.1855552 3.8792786 1.7750056 3.5654419 5.2806261
## sigma[5] 6.8635933 4.3842699 1.9412039 3.9073501 5.6859233
## sigma[6] 2.2387985 1.1126122 0.7600400 1.4472525 1.9784026
## sigma[7] 1.7459429 0.9606330 0.5126567 1.0956008 1.5653736
## sigma[8] 2.3845180 1.1387824 0.9346546 1.5859533 2.1849052
## sigma[9] 1.9605819 0.6874100 0.9505906 1.4542897 1.8506106
## sigma[10] 2.8065994 1.5244945 0.9319761 1.7434476 2.3779669
## sigma[11] 2.9364293 1.4004674 1.1139399 2.0069412 2.6408828
## sigma[12] 1.8679369 0.9823361 0.6719143 1.1499115 1.6640020
## sigma[13] 1.9212575 0.9133235 0.7806373 1.3109717 1.7406391
## sigma[14] 1.9276526 2.4529726 0.1960917 0.5856412 1.0607559
## sigma[15] 2.0757220 0.7865242 0.9542418 1.5042219 1.9891191
## sigma[16] 4.6631129 1.8260804 2.0408246 3.3432561 4.3370938
## sigma[17] 2.0483458 0.9452213 0.7977965 1.3756457 1.9147808
## sigma[18] 2.2905033 0.9969409 0.9147586 1.6157591 2.0653062
## sigma[19] 2.1128392 0.7955496 0.9010698 1.5348701 2.0012651
## sigma[20] 2.4845815 1.3163087 0.8179326 1.5559962 2.1653844
## eta 5.4218260 1.0061599 3.8377853 4.7232576 5.2463085
## mu 0.8275075 0.2710516 0.3366577 0.6434033 0.8238716
## lp__ -514.7416362 5.3620473 -527.1797270 -518.1435593 -514.3334155
## stats
## parameter 75% 97.5%
## b[1] 3.2394273 7.1540369
## b[2] 4.5046584 5.6340100
## b[3] -4.8563369 -3.3165037
## b[4] -3.7378460 2.8182033
## b[5] 4.1943422 8.0335462
## b[6] 8.2528441 9.7191028
## b[7] -1.5140689 0.4039937
## b[8] -6.1661944 -4.7222840
## b[9] -8.8126164 -7.7764848
## b[10] -0.1238670 1.7246723
## b[11] 0.1269496 1.7913320
## b[12] -4.7901716 -3.4450804
## b[13] 3.2310897 3.9807287
## b[14] 7.3630803 10.2669962
## b[15] 7.2644536 8.5677175
## b[16] -0.1390477 1.7594996
## b[17] 2.1908596 3.2933117
## b[18] 5.8792438 6.9637505
## b[19] -0.6276683 0.3752896
## b[20] 3.9855606 6.0825477
## z 1.0699799 2.7358275
## sigma[1] 8.4163391 13.8308765
## sigma[2] 2.0631840 3.7910649
## sigma[3] 4.4173413 6.5757670
## sigma[4] 7.9022639 15.3661912
## sigma[5] 8.8628592 18.4838193
## sigma[6] 2.7737641 4.9875559
## sigma[7] 2.1316923 3.9740806
## sigma[8] 2.8763866 5.2089962
## sigma[9] 2.3121339 3.6420945
## sigma[10] 3.5227675 6.9794127
## sigma[11] 3.5330862 6.6133841
## sigma[12] 2.3017262 4.4075756
## sigma[13] 2.2745081 4.3400505
## sigma[14] 2.1873732 8.4267815
## sigma[15] 2.4625319 3.8720214
## sigma[16] 5.5182722 9.0446598
## sigma[17] 2.5288280 4.2042036
## sigma[18] 2.7948363 4.7330063
## sigma[19] 2.5606354 3.7947111
## sigma[20] 3.1322195 5.6217127
## eta 5.9819495 7.8065980
## mu 1.0069409 1.3574087
## lp__ -510.8348466 -506.0501102
##
## , , chains = chain:3
##
## stats
## parameter mean sd 2.5% 25% 50%
## b[1] 2.1791042 2.5178799 -2.5321858 0.7711011 2.0170806
## b[2] 3.8000511 0.9852116 1.9434208 3.1393114 3.7810666
## b[3] -5.7087007 1.4674757 -8.6808930 -6.5628572 -5.6710763
## b[4] -6.5211822 3.7480290 -11.2743901 -9.4098615 -7.4138093
## b[5] 2.6968527 2.7816596 -4.4164012 1.3260608 2.8190104
## b[6] 7.2212197 1.6216395 3.2596299 6.4632028 7.4459603
## b[7] -1.8812000 0.8732212 -3.4681630 -2.4405759 -2.0072722
## b[8] -6.7102815 1.1714157 -8.6336479 -7.4394545 -6.8459384
## b[9] -9.3524299 0.7609646 -10.7460054 -9.8647123 -9.4305974
## b[10] -0.9161628 1.2452733 -3.4166401 -1.6422338 -0.9943315
## b[11] -0.7190106 1.2895007 -3.5935095 -1.3580885 -0.5486893
## b[12] -5.2025032 0.7742154 -6.7468765 -5.6375958 -5.1461556
## b[13] 2.6537664 0.9978123 0.4452065 2.1408854 2.7590095
## b[14] 6.6622917 1.9547254 1.6840497 6.3953089 7.0046451
## b[15] 6.7200885 0.7732623 5.3166502 6.2151897 6.6813661
## b[16] -1.2240034 1.3752872 -3.7727046 -2.1508849 -1.2498052
## b[17] 1.6346265 0.8601532 0.2142913 1.0906452 1.5173888
## b[18] 5.2547962 1.0343736 2.8323313 4.7252983 5.3303989
## b[19] -1.2049879 1.0010056 -3.2140923 -1.8802838 -1.2055224
## b[20] 3.1924410 1.2837466 0.6311140 2.4556544 3.1832640
## z 0.1740235 1.2725749 -2.2479653 -0.7007524 0.1647213
## sigma[1] 6.8042291 2.9943756 2.4632357 4.8015582 6.2396711
## sigma[2] 1.7456688 0.7790501 0.7218473 1.1963269 1.5704739
## sigma[3] 3.7048521 1.2643084 1.7612136 2.8963767 3.5161802
## sigma[4] 5.4321847 3.3415484 1.2648468 3.0022899 4.6537648
## sigma[5] 6.7651001 4.3786117 2.0293067 3.8001892 5.5523819
## sigma[6] 2.2865051 1.1882955 0.7638665 1.5124397 2.0035943
## sigma[7] 1.8726096 1.1601998 0.5064208 1.1338820 1.5940513
## sigma[8] 2.3564115 1.1148982 0.9105006 1.6084186 2.1647256
## sigma[9] 1.9081959 0.6399942 0.9050804 1.4411266 1.8477822
## sigma[10] 2.7290251 1.6033937 0.9201798 1.6830406 2.2548686
## sigma[11] 2.7701060 1.3729455 1.0864426 1.9020986 2.4573398
## sigma[12] 1.8697967 1.0230954 0.6534886 1.1589249 1.6409366
## sigma[13] 1.9615661 1.0647421 0.6430493 1.2674755 1.7542604
## sigma[14] 1.6137602 2.1347799 0.1838878 0.5313934 0.9653480
## sigma[15] 2.0393168 0.7093241 0.9629960 1.5199201 1.9697861
## sigma[16] 4.5654865 2.1168538 1.8532315 3.0735285 4.1308736
## sigma[17] 2.0535561 0.9938773 0.7706025 1.3314794 1.8434072
## sigma[18] 2.2681702 0.9920663 0.9564820 1.5386466 2.1182462
## sigma[19] 2.1096132 0.8416030 0.9518300 1.5234280 1.9786798
## sigma[20] 2.3254372 1.0832222 0.8877991 1.5718431 2.1593914
## eta 5.4720844 1.0850787 3.7463123 4.7567843 5.3024838
## mu 0.8096770 0.2486056 0.3217010 0.6434201 0.7924835
## lp__ -514.7589854 5.1158320 -525.2052548 -518.0331450 -514.3324618
## stats
## parameter 75% 97.5%
## b[1] 3.44440051 8.23011247
## b[2] 4.50766117 5.64651540
## b[3] -4.82486570 -2.89305336
## b[4] -4.34359238 2.22721349
## b[5] 4.26984141 8.01033285
## b[6] 8.28920540 9.74316013
## b[7] -1.38477632 -0.05788307
## b[8] -6.09862074 -3.90835006
## b[9] -8.83571936 -7.71231303
## b[10] -0.27702603 1.78349161
## b[11] 0.06908281 1.46071362
## b[12] -4.72625558 -3.87852520
## b[13] 3.25438844 4.45202664
## b[14] 7.38487677 9.15963622
## b[15] 7.14950307 8.33184262
## b[16] -0.32851383 1.53557399
## b[17] 2.18422437 3.56274286
## b[18] 5.88365778 7.07169789
## b[19] -0.50805517 0.60987918
## b[20] 4.04407464 5.65136678
## z 1.06178708 2.47807517
## sigma[1] 8.21080004 13.93875596
## sigma[2] 2.13086876 3.71538178
## sigma[3] 4.25657769 6.77359607
## sigma[4] 7.36103622 13.19766635
## sigma[5] 8.66206008 16.63045840
## sigma[6] 2.75550824 5.30177493
## sigma[7] 2.26863418 4.88902654
## sigma[8] 2.81383891 4.88160735
## sigma[9] 2.24739756 3.32305947
## sigma[10] 3.34843570 6.91303324
## sigma[11] 3.34878044 6.16971731
## sigma[12] 2.29972209 4.65159144
## sigma[13] 2.34412053 5.16978198
## sigma[14] 2.04868027 6.44036403
## sigma[15] 2.46929625 3.55934395
## sigma[16] 5.45771673 10.02218327
## sigma[17] 2.52342688 4.49726112
## sigma[18] 2.76571785 4.87615495
## sigma[19] 2.50136893 4.28273448
## sigma[20] 2.91978715 4.87879955
## eta 6.04982229 8.03654485
## mu 0.96125038 1.36640213
## lp__ -511.39488711 -505.19935757
##
## , , chains = chain:4
##
## stats
## parameter mean sd 2.5% 25% 50%
## b[1] 1.7931089 2.4276849 -3.0106836 0.4174966 1.6870973
## b[2] 3.7896973 0.9187639 2.0804906 3.1370758 3.7706767
## b[3] -5.8064727 1.2854783 -8.3054729 -6.6554154 -5.8192659
## b[4] -5.8481307 3.8609576 -11.6208669 -9.0077377 -6.4228204
## b[5] 2.6200342 2.4376182 -2.2651582 1.0076933 2.7105097
## b[6] 7.3572753 1.3460874 4.5004220 6.5057460 7.4399199
## b[7] -2.0009954 0.8477368 -3.5433323 -2.5262706 -2.0664679
## b[8] -6.7695604 0.9307552 -8.4440038 -7.4267244 -6.8462397
## b[9] -9.3841408 0.7182033 -10.6897880 -9.8880401 -9.3885849
## b[10] -1.0037483 1.1297839 -3.1561928 -1.6305070 -1.0032337
## b[11] -0.8058656 1.3654529 -3.8150098 -1.5119675 -0.6520619
## b[12] -5.1813056 0.7090097 -6.6067512 -5.5999983 -5.2186482
## b[13] 2.7429117 0.8239467 1.1186277 2.2491195 2.7707394
## b[14] 6.8247315 1.2843625 3.5248731 6.4000825 7.0118399
## b[15] 6.7404446 0.7386228 5.3602419 6.2966056 6.6674645
## b[16] -1.1778374 1.7418221 -4.5979159 -2.1189591 -1.1762998
## b[17] 1.6522627 0.8594472 0.1263638 1.0950926 1.5451679
## b[18] 5.3575673 0.9080715 3.4345802 4.8572487 5.4355501
## b[19] -1.2441129 0.9549160 -3.0485468 -1.9088098 -1.1733773
## b[20] 3.0497637 1.3646101 0.2462423 2.1845917 3.0159024
## z 0.1845035 1.3207826 -2.7506013 -0.5863815 0.1638500
## sigma[1] 6.6955438 2.8024576 2.8203997 4.7699337 6.3268255
## sigma[2] 1.7403555 0.7351927 0.7245925 1.1940699 1.5973815
## sigma[3] 3.6645249 1.1295230 1.8321543 2.8528704 3.5839018
## sigma[4] 6.0482582 4.1312518 1.4011802 3.6286090 5.0833703
## sigma[5] 6.6359725 3.7213600 1.8497196 3.7246595 5.7198210
## sigma[6] 2.2039973 1.0822513 0.7181505 1.4208693 1.9739216
## sigma[7] 1.8023599 0.9946695 0.5867989 1.0447189 1.5771744
## sigma[8] 2.3373805 1.0969616 0.7832446 1.5657184 2.1034985
## sigma[9] 2.0055533 0.6937603 0.9365104 1.4891288 1.9394817
## sigma[10] 2.5538515 1.4746410 0.8732389 1.6021722 2.1584343
## sigma[11] 2.8782864 1.3989269 1.0033575 1.8889067 2.6302722
## sigma[12] 1.7649792 0.7955062 0.6977436 1.1517059 1.6040382
## sigma[13] 1.8391174 0.8283834 0.6691360 1.2680019 1.6955958
## sigma[14] 1.4092258 1.4109632 0.1876624 0.5233833 0.9147096
## sigma[15] 2.0198097 0.7566523 0.9040355 1.4810239 1.8822997
## sigma[16] 4.8223559 2.1702543 2.0307211 3.4399425 4.3353261
## sigma[17] 2.0158480 0.9144391 0.7153862 1.4053645 1.8372188
## sigma[18] 2.2420223 0.9754686 0.9245454 1.5693596 2.0900833
## sigma[19] 2.1047876 0.8295296 0.9912017 1.5278918 1.9450625
## sigma[20] 2.4806425 1.3915943 0.7391739 1.5032504 2.1633497
## eta 5.4816983 1.0153221 3.8775785 4.7655737 5.3600690
## mu 0.8025913 0.2477354 0.3247655 0.6333511 0.8066669
## lp__ -513.7123930 5.1132732 -524.6108363 -516.6447569 -513.4916488
## stats
## parameter 75% 97.5%
## b[1] 3.09960928 7.4839069
## b[2] 4.42278200 5.6553634
## b[3] -4.94233838 -3.3399291
## b[4] -3.18012947 1.9908752
## b[5] 4.05158907 7.7556758
## b[6] 8.37046728 9.5510377
## b[7] -1.50033557 -0.1679667
## b[8] -6.16452792 -4.8800290
## b[9] -8.92855222 -7.7762080
## b[10] -0.25719889 1.1776447
## b[11] 0.05990573 1.4522475
## b[12] -4.74748104 -3.7595991
## b[13] 3.29937724 4.1461720
## b[14] 7.39548850 9.1221536
## b[15] 7.15154209 8.4340048
## b[16] -0.10841897 2.2791314
## b[17] 2.14754329 3.5640411
## b[18] 5.93889916 6.9382630
## b[19] -0.53785077 0.4385380
## b[20] 3.89417427 5.5803023
## z 1.03617701 2.9462455
## sigma[1] 8.04340388 13.4913715
## sigma[2] 2.17183480 3.4474533
## sigma[3] 4.31504970 6.2284109
## sigma[4] 7.38649857 15.3371979
## sigma[5] 8.80653349 15.1119497
## sigma[6] 2.71310703 4.8694557
## sigma[7] 2.22861132 4.1914147
## sigma[8] 2.91845303 5.0121214
## sigma[9] 2.38419882 3.6484986
## sigma[10] 3.08967825 6.3953078
## sigma[11] 3.57021494 6.1041079
## sigma[12] 2.23644001 3.5935179
## sigma[13] 2.19566927 3.9035105
## sigma[14] 1.74438459 5.4041139
## sigma[15] 2.42096303 3.7803624
## sigma[16] 5.79601574 10.3939511
## sigma[17] 2.48950850 4.0881436
## sigma[18] 2.64090199 4.7544966
## sigma[19] 2.50103226 3.9335365
## sigma[20] 3.04546092 6.2119047
## eta 6.07962299 7.9041645
## mu 0.98097561 1.2638856
## lp__ -509.95864231 -505.2163494What do we want to know?
Do the treatments affect mean growth?
How much does leaf mass in treatment 8 differ from group 12?
How much variation is explained by treatment?
How do we find it out?
See if any credible intervals overlap? And, see if credible interval for eta > 0.
Posterior distribution of
b[12] - b[8]: if>0thenb[12] > b[8]; if zero is in it thenb[12] = b[8].Compare
sigmatoeta. Or take the ratio of MADs before and after subtractingb[-].
Do treatments affect growth?
First let’s simply see if the credible intervals overlap. Many of these don’t overlap even close, so it sure looks like some of the treatments were different from other ones.
Now, let’s look at eta, which is the SD parameter on the prior for b, the group location parameters. If eta is zero, then all the bs are the same.
It looks like eta has a wide posterior distribution, but all the posterior support is quite far from zero, indicating that the group locations are different between the groups.
## mean se_mean sd 2.5% 25% 50% 75%
## eta 5.468179 0.02758036 1.0167 3.846221 4.766442 5.325984 6.060594
## 97.5% n_eff Rhat
## eta 7.872626 1358.895 0.9994464
How much do treatments 8 and 12 differ?
Here we look at the posterior distribution of b[12] - b[8].
bsamps <- extract(anova_fit)$b
diff_samps <- bsamps[,12] - bsamps[,8]
hist(diff_samps, main="posterior of b[12] - b[8]", xlab='b[12] - b[8]')
The posterior mean difference is 1.5752492 but a 95% credible interval goes from -1.1545359 to 3.9845409. Since zero is in the credible interval, we can’t really tell from the data if there is a real difference here, but if there is, it’s probably less than 5 or so.
How much variation is explained by treatment?
We want to compare eta (which describes how much different treatments affect leaf mass) to the values of sigma (which describe how much leaf mass varies within each treatment).
First we need to check how much sigma vary.
The posterior distributions of different groups’ sigma values overlap a lot; so we’ll look at the posterior distribution of the mean sigma value to describe how much within-group variation there is.
It looks like sigma is between 2 and 4, roughly.
We really want to know how big eta tends to be relative to sigma
eta_samples <- extract(anova_fit)$eta
ratio_samples <- eta_samples / mean_sigmas
hist(ratio_samples)
The posterior mean ratio of eta to sigma is 1.8691119, with a 95% credible interval of 1.1587616 to 2.8474252. In other words, it looks like differences in leaf mass between treatments are almost twice the typical differences between plants that had the same treatment.