text:
# 18 Metric Predicted Variable with Multiple Metric Predictors

We will consider models in which the predicted variable is an additive combination of predictors, all of which have proportional influence on the prediction. This kind of model is called multiple linear regression. We will also consider nonadditive combinations of predictors, which are called interactions. (p. 509, emphasis in the original)

## 18.1 Multiple linear regression

If $$y \sim \operatorname{Normal} (\mu, \sigma)$$ and $$\mu = \beta_0 + \beta_1 x_1 + \beta_2 x_2$$, then it’s also the case that we can rewrite the formula for $$y$$ as

$y \sim \operatorname{Normal} (\beta_0 + \beta_1 x_1 + \beta_2 x_2, \sigma).$

As Kruschke pointed out, the basic model “assumes homogeneity of variance, which means that at all values of $$x_1$$ and $$x_2$$, the variance $$\sigma^2$$ of $$y$$ is the same” (p. 510).

If we presume the data for the two $$x$$ variables are uniformly distributed within 0 and 10, we can make the data for Figure 18.1 like this.

library(tidyverse)

n <- 300

set.seed(18)
d <-
tibble(x_1 = runif(n = n, 0, 10),
x_2 = runif(n = n, 0, 10)) %>%
mutate(y   = rnorm(n = n, mean = 10 + x_1 + 2 * x_2))

head(d)
## # A tibble: 6 x 3
##     x_1   x_2     y
##   <dbl> <dbl> <dbl>
## 1 8.23   8.62  36.9
## 2 7.10   1.33  20.3
## 3 9.66   1.08  19.5
## 4 0.786  7.09  25.1
## 5 0.536  6.67  25.6
## 6 5.75   4.72  25.7

Here are three of the scatter plots in Figure 18.1.

theme_set(theme_grey() +
theme(panel.grid = element_blank()))

d %>%
ggplot(aes(x = x_1, y = y)) +
geom_point(alpha = 2/3) +
coord_cartesian(ylim = 0:50)

d %>%
ggplot(aes(x = x_2, y = y)) +
geom_point(alpha = 2/3) +
coord_cartesian(ylim = 0:50)

d %>%
ggplot(aes(x = x_1, y = x_2)) +
geom_vline(xintercept = seq(from = 0, to = 10, by = .5), color = "grey98") +
geom_hline(yintercept = seq(from = 0, to = 10, by = .5), color = "grey98") +
geom_point(alpha = 2/3)

As in previous chapters, I’m not aware that ggplot2 allows for three-dimensional wireframe plots of the kind in the upper left panel. If you’d like to make one in base R, have at it. But with respect to the plots in the off-diagonal, I have no idea what Kruschke did with the compressed and diagonal grids. If you have that figured out, please share your code.

For Figure 18.2, the $$x$$ variables look to be multivariate normal with a correlation of about -.95. We can simulate such data with help from the MASS package.

Sven Hohenstein’s answer to this stats.stackexchange.com question provides the steps for simulating the data. First, we’ll need to specify the desired means and standard deviations for our variables. Then we’ll make a correlation matrix with 1s on the diagonal and the desired correlation coefficient, $$\rho$$ on the off-diagonal. Since the correlation matrix is symmetric, both off-diagonal positions are the same. Then we convert the correlation matrix to a covariance matrix.

mus <- c(5, 5)
sds <- c(2, 2)

cors <- matrix(c(1, -.95,
-.95, 1),
ncol = 2)
cors
##       [,1]  [,2]
## [1,]  1.00 -0.95
## [2,] -0.95  1.00
covs <- sds %*% t(sds) * cors
covs
##      [,1] [,2]
## [1,]  4.0 -3.8
## [2,] -3.8  4.0

Now we’ve defined our means, standard deviations, and covariance matrix, we’re ready to simulate the data with the MASS::mvrnorm() function.

# how many data points would you like to simulate?
n <- 300

set.seed(18.2)

d <-
MASS::mvrnorm(n         = n,
mu        = mus,
Sigma     = covs,
empirical = T) %>%
as_tibble() %>%
set_names("x_1", "x_2") %>%
mutate(y = rnorm(n = n, mean = 10 + x_1 + 2 * x_2))

Now we have our simulated data in hand, we’re ready for three of the four panels of Figure 18.2.

d %>%
ggplot(aes(x = x_1, y = y)) +
geom_point(alpha = 2/3) +
coord_cartesian(xlim = 0:10,
ylim = 0:50)

d %>%
ggplot(aes(x = x_2, y = y)) +
geom_point(alpha = 2/3) +
coord_cartesian(xlim = 0:10,
ylim = 0:50)

d %>%
ggplot(aes(x = x_1, y = x_2)) +
geom_vline(xintercept = seq(from = 0, to = 10, by = .5), color = "grey98") +
geom_hline(yintercept = seq(from = 0, to = 10, by = .5), color = "grey98") +
geom_point(alpha = 2/3) +
coord_cartesian(xlim = 0:10,
ylim = 0:10)

Pretty close.

### 18.1.1 The perils of correlated predictors.

Figures 18.1 and 18.2 show data generated from the same model. In both figures, $$\sigma = 2, \beta_0 = 10, \beta_1 = 1, \text{ and } \beta_2 = 2$$. All that differs between the two figures is the distribution of the $$\langle x_1, x_2 \rangle$$ values, which is not specified by the model. In Figure 18.1, the $$\langle x_1, x_2 \rangle$$ values are distributed independently. In Figure 18.2, the $$\langle x_1, x_2 \rangle$$ values are negatively correlated: When $$x_1$$ is small, $$x_2$$ tends to be large, and when $$x_1$$ is large, $$x_2$$ tends to be small. (p. 510)

If you look closely at our simulation code from above, you’ll see we have done so, too.

Real data often have correlated predictors. For example, consider trying to predict a state’s average high-school SAT score on the basis of the amount of money the state spends per pupil. If you plot only mean SAT against money spent, there is actually a decreasing trend… (p. 513, emphasis in the original)

Before we remake Figure 18.3 to examine that decreasing trend, we’ll need to load the data.

my_data <- read_csv("data.R/Guber1999data.csv")

glimpse(my_data)
## Observations: 50
## Variables: 8
## $State <chr> "Alabama", "Alaska", "Arizona", "Arkansas", "Californi… ##$ Spend     <dbl> 4.405, 8.963, 4.778, 4.459, 4.992, 5.443, 8.817, 7.030…
## $StuTeaRat <dbl> 17.2, 17.6, 19.3, 17.1, 24.0, 18.4, 14.4, 16.6, 19.1, … ##$ Salary    <dbl> 31.144, 47.951, 32.175, 28.934, 41.078, 34.571, 50.045…
## $PrcntTake <dbl> 8, 47, 27, 6, 45, 29, 81, 68, 48, 65, 57, 15, 13, 58, … ##$ SATV      <dbl> 491, 445, 448, 482, 417, 462, 431, 429, 420, 406, 407,…
## $SATM <dbl> 538, 489, 496, 523, 485, 518, 477, 468, 469, 448, 482,… ##$ SATT      <dbl> 1029, 934, 944, 1005, 902, 980, 908, 897, 889, 854, 88…

Now we have the Guber (1999) data, here are three of the four plots.

library(patchwork)

p1 <-
my_data %>%
ggplot(aes(x = Spend, y = SATT)) +
geom_point(alpha = 2/3) +
coord_cartesian(ylim = 800:1120)

p2 <-
my_data %>%
ggplot(aes(x = PrcntTake, y = SATT)) +
geom_point(alpha = 2/3) +
xlab("% Take") +
coord_cartesian(ylim = 800:1120)

p3 <-
my_data %>%
ggplot(aes(x = PrcntTake, y = Spend)) +
geom_vline(xintercept = seq(from = 5,   to = 80, by = 5), color = "grey98") +
geom_hline(yintercept = seq(from = 3.5, to = 10, by = .5), color = "grey98") +
geom_point(alpha = 2/3) +
xlab("% Take")

# bind them together and add a title
wrap_elements(grid::textGrob('No 3D wireframe plots for us')) +
p1 + p2 + p3 +
plot_annotation(title = "SATT ~ N(m,sd=31.5), m = 993.8 + −2.9 %Take + 12.3 Spend")

### 18.1.2 The model and implementation.

We’ll make a custom function to standardize the criterion and predictor values.

standardize <- function(x) {
(x - mean(x)) / sd(x)
}

my_data <-
my_data %>%
mutate(prcnt_take_z = standardize(PrcntTake),
spend_z      = standardize(Spend),
satt_z       = standardize(SATT))

Let’s open brms.

library(brms)

Now we’re ready to fit the model. As Kruschke pointed out, the priors on the standardized predictors are set with

an arbitrary standard deviation of 2.0. This value was chosen because standardized regression coefficients are algebraically constrained to fall between −1 and +1 in least-squares regression, [1] and therefore, the regression coefficients will not exceed those limits by much. A normal distribution with standard deviation of 2.0 is reasonably flat over the range from −1 to +1. (p. 516)

With data like this, even a prior(normal(0, 1), class = b) would be only mildly regularizing.

This is a good place to emphasize how priors in brms are given classes. If you’d like all parameters within a given class to have the prior, you can just specify one prior argument within that class. For our fit1, both parameters of class = b have a normal(0, 2) prior. So we can just include one statement to handle both. Had we wanted different priors for the coefficients for spend_z and prcnt_take_z, we’d need to include two prior() arguments with at least one including a coef argument.

fit1 <-
brm(data = my_data,
family = student,
satt_z ~ 1 + spend_z + prcnt_take_z,
prior = c(prior(normal(0, 2), class = Intercept),
prior(normal(0, 2), class = b),
prior(normal(0, 1), class = sigma),
prior(exponential(one_over_twentynine), class = nu)),
chains = 4, cores = 4,
stanvars = stanvar(1/29, name = "one_over_twentynine"),
seed = 18)

Check the model summary.

print(fit1)
##  Family: student
##   Links: mu = identity; sigma = identity; nu = identity
## Formula: satt_z ~ 1 + spend_z + prcnt_take_z
##    Data: my_data (Number of observations: 50)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
##
## Population-Level Effects:
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept        0.00      0.06    -0.12     0.13 1.00     4222     2585
## spend_z          0.23      0.08     0.08     0.38 1.00     3599     2990
## prcnt_take_z    -1.03      0.08    -1.18    -0.87 1.00     3490     3160
##
## Family Specific Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.42      0.05     0.32     0.53 1.00     3523     2413
## nu       33.07     27.90     4.27   107.31 1.00     3398     2360
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

So when we use a multivariable model, increases in spending now appear associated with increases in SAT scores.

### 18.1.3 The posterior distribution.

Based on Equation 18.1, we can convert the standardized coefficients from our multivariable model back to their original metric as follows:

\begin{align*} \beta_0 & = \operatorname{SD}_y \zeta_0 + M_y - \operatorname{SD}_y \sum_j \frac{\zeta_j M_{x_j}}{\operatorname{SD}_{x_j}} \;\;\; \text{and} \\ \beta_j & = \frac{\operatorname{SD}_y \zeta_j}{\operatorname{SD}_{x_j}}. \end{align*}

To use them, we’ll first extract the posterior samples.

post <- posterior_samples(fit1)

head(post)
##   b_Intercept b_spend_z b_prcnt_take_z     sigma        nu      lp__
## 1  0.03332938 0.2246157      -1.044109 0.4381914 58.891563 -35.40259
## 2 -0.05051637 0.2270969      -1.076683 0.3759888 13.331876 -35.92027
## 3 -0.04916057 0.2770627      -1.092591 0.3992656 18.839407 -35.58017
## 4  0.01043252 0.3638677      -1.141387 0.4167015  5.558774 -37.95336
## 5  0.02980380 0.3671874      -1.262402 0.4401791 29.336185 -39.44394
## 6  0.00317025 0.3080994      -1.259354 0.4584194 46.137031 -39.81174

Like we did in Chapter 17, let’s wrap the consequences of Equation 18.1 into two functions.

make_beta_0 <- function(zeta_0, zeta_1, zeta_2, sd_x_1, sd_x_2, sd_y, m_x_1, m_x_2, m_y) {
sd_y * zeta_0 + m_y - sd_y * ((zeta_1 * m_x_1 / sd_x_1) + (zeta_2 * m_x_2 / sd_x_2))
}

make_beta_j <- function(zeta_j, sd_j, sd_y) {
sd_y * zeta_j / sd_j
}

After saving a few values, we’re ready to use our custom functions.

sd_x_1 <- sd(my_data$Spend) sd_x_2 <- sd(my_data$PrcntTake)
sd_y   <- sd(my_data$SATT) m_x_1 <- mean(my_data$Spend)
m_x_2  <- mean(my_data$PrcntTake) m_y <- mean(my_data$SATT)

post <-
post %>%
mutate(b_0 = make_beta_0(zeta_0 = b_Intercept,
zeta_1 = b_spend_z,
zeta_2 = b_prcnt_take_z,
sd_x_1 = sd_x_1,
sd_x_2 = sd_x_2,
sd_y   = sd_y,
m_x_1  = m_x_1,
m_x_2  = m_x_2,
m_y    = m_y),
b_1 = make_beta_j(zeta_j = b_spend_z,
sd_j   = sd_x_1,
sd_y   = sd_y),
b_2 = make_beta_j(zeta_j = b_prcnt_take_z,
sd_j   = sd_x_2,
sd_y   = sd_y))

glimpse(post)
## Observations: 4,000
## Variables: 9
## $b_Intercept <dbl> 0.033329378, -0.050516372, -0.049160566, 0.010432… ##$ b_spend_z      <dbl> 0.22461570, 0.22709693, 0.27706273, 0.36386773, 0…
## $b_prcnt_take_z <dbl> -1.0441089, -1.0766835, -1.0925912, -1.1413875, -… ##$ sigma          <dbl> 0.4381914, 0.3759888, 0.3992656, 0.4167015, 0.440…
## $nu <dbl> 58.891563, 13.331876, 18.839407, 5.558774, 29.336… ##$ lp__           <dbl> -35.40259, -35.92027, -35.58017, -37.95336, -39.4…
## $b_0 <dbl> 998.4585, 994.5899, 980.0592, 961.1825, 973.4781,… ##$ b_1            <dbl> 12.331810, 12.468035, 15.211248, 19.977000, 20.15…
## $b_2 <dbl> -2.919049, -3.010119, -3.054593, -3.191014, -3.52… Here’s the top panel of Figure 18.5. library(tidybayes) # here are the primary data post %>% transmute(Intercept = b_0, Spend = b_1, Percent Take = b_2, Scale = sigma * sd_y, Normality = nu %>% log10()) %>% gather() %>% # the plot ggplot(aes(x = value)) + geom_histogram(color = "grey92", fill = "grey67", size = .2, bins = 40) + stat_pointintervalh(aes(y = 0), point_interval = mode_hdi, .width = c(.95, .5)) + scale_y_continuous(NULL, breaks = NULL) + xlab(NULL) + facet_wrap(~key, scales = "free", ncol = 3) The slope on spending has a mode of about 13, which suggests that SAT scores rise by about 13 points for every extra$1000 spent per pupil. The slope on percentage taking the exam (PrcntTake) is also credibly non-zero, with a mode around −2.8, which suggests that SAT scores fall by about 2.8 points for every additional 1% of students who take the test. (p. 517)

If you want those exact modes and, say, 50% intervals around them, you can just use tidybayes::mode_hdi().

post %>%
transmute(Spend          = b_1,
Percent Take = b_2) %>%
gather() %>%
group_by(key) %>%
mode_hdi(value, .width = .5)
## # A tibble: 2 x 7
##   key          value .lower .upper .width .point .interval
##   <chr>        <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>
## 1 Percent Take -2.86  -3.04  -2.74    0.5 mode   hdi
## 2 Spend        12.9   10.3   15.9     0.5 mode   hdi

The brms::bayes_R2() function makes it easy to compute a Bayesian $$R^2$$. Simply feed a brm() fit object into bayes_R2() and you’ll get back the posterior mean, $$SD$$, and 95% intervals.

bayes_R2(fit1)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.8139321 0.02239821 0.7572305 0.8429055

I’m not going to go into the technical details here, but you should be aware that the Bayeisan $$R^2$$ returned from the bayes_R2() function is not calculated the same as it is with OLS. If you want to dive in, check out the paper by Gelman, Goodrich, Gabry, and Vehtari. Anyway, if you’d like to view the Bayesian $$R^2$$ distribution rather than just get the summaries, specify summary = F, convert the output to a tibble, and plot as usual.

bayes_R2(fit1, summary = F) %>%
as_tibble() %>%

ggplot(aes(x = R2)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 25) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = .95) +
scale_y_continuous(NULL, breaks = NULL) +
labs(subtitle = expression(paste("Bayesian ", italic(R)^2)),
x = NULL) +
coord_cartesian(xlim = c(.6, 1))

Since the brms::bayes_R2() function is not identical with Kruschke’s method in the text, the results might differ a bit.

We can get a sense of the scatter plots with bayesplot::mcmc_pairs().

library(bayesplot)

color_scheme_set("gray")

post %>%
transmute(Intercept      = b_0,
Spend          = b_1,
Percent Take = b_2,
Scale          = sigma * sd_y,
Normality      = nu %>% log10()) %>%
mcmc_pairs(off_diag_args = list(size = 1/8, alpha = 1/8))

One way to get the Pearson’s correlation coefficients among the parameters is with psych::lowerCor().

post %>%
transmute(Intercept      = b_0,
Spend          = b_1,
Percent Take = b_2,
Scale          = sigma * sd_y,
Normality      = nu %>% log10()) %>%
psych::lowerCor(digits = 3)
##              Intrcp Spend  PrcntT Scale  Nrmlty
## Intercept     1.000
## Spend        -0.934  1.000
## Percent Take  0.349 -0.613  1.000
## Scale         0.076 -0.083  0.080  1.000
## Normality     0.099 -0.108  0.089  0.342  1.000

Kruschke finished the subsection with the observation: “Sometimes we are interested in using the linear model to predict $$y$$ values for $$x$$ values of interest. It is straight forward to generate a large sample of credible $$y$$ values for specified $$x$$ values” (p. 519).

Like we practiced with in the last chapter, the simplest way to do so in brms is with the fitted() function. For a quick example, say we wanted to know what the model would predict if we were to have a standard-score increase in spending and a simultaneous standard-score decrease in the percent taking the exam. We’d just specify those values in a tibble and feed that tibble into fitted() along with the model.

nd <-
tibble(spend_z      =  1,
prcnt_take_z = -1)

fitted(fit1,
newdata = nd)
##      Estimate Est.Error      Q2.5    Q97.5
## [1,] 1.262009 0.1547968 0.9596958 1.568147

### 18.1.4 Redundant predictors.

As a simplified example of correlated predictors, think of just two data points: Suppose $$y = 1 \text{ for } \langle x_1, x_2 \rangle = \langle 1, 1 \rangle \text{ and } y = 2 \text{ for } \langle x_1, x_2 \rangle = \langle 2, 2 \rangle$$. The linear model, $$y = \beta_1 x_1 + \beta_2 x_2$$ is supposed to satisfy both data points, and in this case both are satisfied by $$1 = \beta_1 + \beta_2$$. Therefore, many different combinations of $$\beta_1$$ and $$\beta_2$$ satisfy the data. For example, it could be that $$\beta_1 = 2$$ and $$\beta_2 = -1$$, or $$\beta_1 = 0.5$$ and $$\beta_2 = 0.5$$, or $$\beta_1 = 0$$ and $$\beta_2 = 1$$. In other words, the credible values of $$\beta_1$$ and $$\beta_2$$ are anticorrelated and trade-off to fit the data. (p. 519)

Here are what those data look like. You would not want to fit a regression model with these data.

tibble(x_1 = 1:2,
x_2 = 1:2,
y   = 1:2)
## # A tibble: 2 x 3
##     x_1   x_2     y
##   <int> <int> <int>
## 1     1     1     1
## 2     2     2     2

We can take percentages and turn them into their inverse re-expressed as a proportion.

percent_take <- 37

(100 - percent_take) / 100
## [1] 0.63

Let’s make a redundant predictor and then standardize() it.

my_data <-
my_data %>%
mutate(prop_not_take   = (100 - PrcntTake) / 100) %>%
mutate(prop_not_take_z = standardize(prop_not_take))

glimpse(my_data)
## Observations: 50
## Variables: 13
## $State <chr> "Alabama", "Alaska", "Arizona", "Arkansas", "Cal… ##$ Spend           <dbl> 4.405, 8.963, 4.778, 4.459, 4.992, 5.443, 8.817,…
## $StuTeaRat <dbl> 17.2, 17.6, 19.3, 17.1, 24.0, 18.4, 14.4, 16.6, … ##$ Salary          <dbl> 31.144, 47.951, 32.175, 28.934, 41.078, 34.571, …
## $PrcntTake <dbl> 8, 47, 27, 6, 45, 29, 81, 68, 48, 65, 57, 15, 13… ##$ SATV            <dbl> 491, 445, 448, 482, 417, 462, 431, 429, 420, 406…
## $SATM <dbl> 538, 489, 496, 523, 485, 518, 477, 468, 469, 448… ##$ SATT            <dbl> 1029, 934, 944, 1005, 902, 980, 908, 897, 889, 8…
## $prcnt_take_z <dbl> -1.0178453, 0.4394222, -0.3078945, -1.0925770, 0… ##$ spend_z         <dbl> -1.10086058, 2.24370805, -0.82716069, -1.0612364…
## $satt_z <dbl> 0.8430838, -0.4266207, -0.2929676, 0.5223163, -0… ##$ prop_not_take   <dbl> 0.92, 0.53, 0.73, 0.94, 0.55, 0.71, 0.19, 0.32, …
## $prop_not_take_z <dbl> 1.0178453, -0.4394222, 0.3078945, 1.0925770, -0.… We’re ready to fit the redundant-predictor model. fit2 <- brm(data = my_data, family = student, satt_z ~ 0 + Intercept + spend_z + prcnt_take_z + prop_not_take_z, prior = c(prior(normal(0, 2), class = b, coef = "Intercept"), prior(normal(0, 2), class = b, coef = "spend_z"), prior(normal(0, 2), class = b, coef = "prcnt_take_z"), prior(normal(0, 2), class = b, coef = "prop_not_take_z"), prior(normal(0, 1), class = sigma), prior(exponential(one_over_twentynine), class = nu)), chains = 4, cores = 4, stanvars = stanvar(1/29, name = "one_over_twentynine"), seed = 18, # this will let us use prior_samples() later on sample_prior = "yes") You might notice a few things about the brm() code. First, we have used the ~ 0 + Intercept + ... syntax instead of the default syntax for intercepts. In normal situations, we would have been in good shape using the typical ~ 1 + ... syntax for the intercept, especially given our use of standardized data. However, since brms version 2.5.0, using the sample_prior argument to draw samples from the prior distribution will no longer allow us to return samples from the typical brms intercept. Bürkner addressed the issue on the Stan forums. As he pointed out, if you want to get prior samples from an intercept, you’ll have to use the alternative syntax. The other thing to point out is that even though we used the same prior on all the predictors, including the intercept, we still explicitly spelled each out with the coef argument. If we hadn’t been explicit like this, we would only get a single b vector from the prior_samples() function. But since we want separate vectors for each of our predictors, we used the verbose code. If you’re having a difficult time understanding these two points, experiment. Fit the model in a few different ways with either the typical or the alternative intercept syntax and with either the verbose prior code or the simplified prior(normal(0, 2), class = b) code. And after each, execute prior_samples(fit2). You’ll see. Let’s move on. Kruschke mentioned high autocorrelations in the prose. Here are the autocorrelation plots for our $$\beta$$s. post <- posterior_samples(fit2, add_chain = T) mcmc_acf(post, pars = c("b_Intercept", "b_spend_z", "b_prcnt_take_z", "b_prop_not_take_z"), lags = 10) Looks like HMC made a big difference. The $$N_{eff}/N$$ ratios weren’t terrible, either. neff_ratio(fit2)[1:6] %>% mcmc_neff() + yaxis_text(hjust = 0) The brms::vcov() function returns a variance/covariance matrix–or a correlation matrix when you set correlation = T–of the population-level parameters (i.e., the fixed effects). It returns the values to a decadent level of precision, so we’ll simplify the output with round(). vcov(fit2, correlation = T) %>% round(digits = 3) ## Intercept spend_z prcnt_take_z prop_not_take_z ## Intercept 1.000 0.026 0.028 0.029 ## spend_z 0.026 1.000 0.013 0.048 ## prcnt_take_z 0.028 0.013 1.000 0.998 ## prop_not_take_z 0.029 0.048 0.998 1.000 Notice how much lower our Spend_z correlations are than those Kruschke displayed on page 520 of the text. However, it turns out the correlations among the redundant predictors were still very high. If any of the nondiagonal correlations are high (i.e., close to +1 or close to −1), be careful when interpreting the posterior distribution. Here, we can see that the correlation of PrcntTake and PropNotTake is −1.0, which is an immediate sign of redundant predictors. (p. 520) You can really get a sense of the silliness of the parameters if you plot them. We’ll use geom_halfeyeh() to get a sense of densities and summaries of the $$\beta$$s. post %>% select(b_Intercept:b_prop_not_take_z) %>% gather() %>% # this line isn't necessary, but it does allow us to arrange the parameters on the y-axis mutate(key = factor(key, levels = c("b_prop_not_take_z", "b_prcnt_take_z", "b_spend_z", "b_Intercept"))) %>% ggplot(aes(x = value, y = key)) + geom_vline(xintercept = 0, color = "white") + geom_halfeyeh(point_interval = mode_hdi, .width = .95, # the next two lines are purely aesthetic scale = "width", relative_scale = .9) + labs(x = NULL, y = NULL) + coord_cartesian(xlim = -5:5) + theme(axis.text.y = element_text(hjust = 0), axis.ticks.y = element_blank()) Yeah, on the standardized scale those are some ridiculous estimates. Let’s update our make_beta_0() function. make_beta_0 <- function(zeta_0, zeta_1, zeta_2, zeta_3, sd_x_1, sd_x_2, sd_x_3, sd_y, m_x_1, m_x_2, m_x_3, m_y) { sd_y * zeta_0 + m_y - sd_y * ((zeta_1 * m_x_1 / sd_x_1) + (zeta_2 * m_x_2 / sd_x_2) + (zeta_3 * m_x_3 / sd_x_3)) } sd_x_1 <- sd(my_data$Spend)
sd_x_2 <- sd(my_data$PrcntTake) sd_x_3 <- sd(my_data$prop_not_take)
sd_y   <- sd(my_data$SATT) m_x_1 <- mean(my_data$Spend)
m_x_2  <- mean(my_data$PrcntTake) m_x_3 <- mean(my_data$prop_not_take)
m_y    <- mean(my_data$SATT) post <- post %>% transmute(Intercept = make_beta_0(zeta_0 = b_Intercept, zeta_1 = b_spend_z, zeta_2 = b_prcnt_take_z, zeta_3 = b_prop_not_take_z, sd_x_1 = sd_x_1, sd_x_2 = sd_x_2, sd_x_3 = sd_x_3, sd_y = sd_y, m_x_1 = m_x_1, m_x_2 = m_x_2, m_x_3 = m_x_3, m_y = m_y), Spend = make_beta_j(zeta_j = b_spend_z, sd_j = sd_x_1, sd_y = sd_y), Percent Take = make_beta_j(zeta_j = b_prcnt_take_z, sd_j = sd_x_2, sd_y = sd_y), Proportion not Take = make_beta_j(zeta_j = b_prop_not_take_z, sd_j = sd_x_3, sd_y = sd_y), Scale = sigma * sd_y, Normality = nu %>% log10()) glimpse(post) ## Observations: 4,000 ## Variables: 6 ##$ Intercept             <dbl> 1061.1824, 280.0482, 873.1032, 551.0707, 6…
## $Spend <dbl> 5.454006, 17.469214, 15.816227, 12.931905,… ##$ Percent Take        <dbl> -3.1303938, 3.9409313, -1.8844443, 1.47699…
## $Proportion not Take <dbl> -26.318716, 682.212440, 105.706334, 432.37… ##$ Scale                 <dbl> 30.16235, 32.49213, 26.98000, 36.66645, 28…
## $Normality <dbl> 1.1409028, 1.3884920, 1.8217587, 1.0168820… Now we’ve done the conversions, here are the histograms of Figure 18.6. post %>% gather() %>% ggplot() + geom_histogram(aes(x = value), color = "grey92", fill = "grey67", size = .2, bins = 40) + stat_pointintervalh(aes(x = value, y = 0), point_interval = mode_hdi, .width = c(.95, .5)) + scale_y_continuous(NULL, breaks = NULL) + xlab(NULL) + facet_wrap(~key, scales = "free", ncol = 3) Their scatter plots are as follows: post %>% mcmc_pairs(off_diag_args = list(size = 1/8, alpha = 1/8)) Since brms::vcov() only computes the Pearson’s correlations among the untransformed population-level intercepts and slopes, it will be of limited use computing the correlations among all our parameters of interest. So here we’ll use psych::lowerCor() to compute the correlations among the full range of our parameters. post %>% psych::lowerCor(digits = 3) ## Intrcp Spend PrcntT PrprnT Scale Nrmlty ## Intercept 1.000 ## Spend -0.100 1.000 ## Percent Take -0.996 0.013 1.000 ## Proportion not Take -0.998 0.048 0.998 1.000 ## Scale 0.002 -0.086 0.007 0.003 1.000 ## Normality 0.033 -0.082 -0.024 -0.030 0.373 1.000 Figure 18.7 is all about the prior predictive distribution. Here we’ll extract the priors with prior_samples() and wrangle all in one step. prior <- prior_samples(fit2) %>% transmute(Intercept = make_beta_0(zeta_0 = b_Intercept, zeta_1 = b_spend_z, zeta_2 = b_prcnt_take_z, zeta_3 = b_prop_not_take_z, sd_x_1 = sd_x_1, sd_x_2 = sd_x_2, sd_x_3 = sd_x_3, sd_y = sd_y, m_x_1 = m_x_1, m_x_2 = m_x_2, m_x_3 = m_x_3, m_y = m_y), Spend = make_beta_j(zeta_j = b_spend_z, sd_j = sd_x_1, sd_y = sd_y), Percent Take = make_beta_j(zeta_j = b_prcnt_take_z, sd_j = sd_x_2, sd_y = sd_y), Proportion not Take = make_beta_j(zeta_j = b_prop_not_take_z, sd_j = sd_x_3, sd_y = sd_y), Scale = sigma * sd_y, Normality = nu %>% log10()) glimpse(prior) ## Observations: 4,000 ## Variables: 6 ##$ Intercept             <dbl> 1462.43927, 664.21947, 1289.22313, 1476.56…
## $Spend <dbl> -38.749588, 13.362327, -49.410718, -93.439… ##$ Percent Take        <dbl> 1.2899620, 6.4475454, 7.1217195, 1.7590988…
## $Proportion not Take <dbl> -310.88980, 63.95417, -280.14134, 101.8131… ##$ Scale                 <dbl> 2.822765, 6.293331, 67.317724, 11.692913, …
## Normality <dbl> 1.5678380, 1.5254245, 1.2846292, 0.6427115… Now we’ve wrangled the priors, we’re ready to make the histograms at the top of Figure 18.7. prior %>% gather() %>% ggplot(aes(x = value)) + geom_histogram(color = "grey92", fill = "grey67", size = .2, bins = 40, boundary = 0) + stat_pointintervalh(aes(y = 0), point_interval = mode_hdi, .width = c(.95, .5)) + scale_y_continuous(NULL, breaks = NULL) + xlab(NULL) + facet_wrap(~key, scales = "free", ncol = 3) Since we used the half-Gaussian prior for our $$\sigma$$, our Scale histogram looks different from Kruschke’s. Otherwise, everything’s on the up and up. Here are the scatter plots at the bottom of Figure 18.7. prior %>% mcmc_pairs(off_diag_args = list(size = 1/8, alpha = 1/8)) And finally, here are those Pearson’s correlation coefficients for the priors. prior %>% psych::lowerCor(digits = 3) ## Intrcp Spend PrcntT PrprnT Scale Nrmlty ## Intercept 1.000 ## Spend -0.826 1.000 ## Percent Take -0.256 0.003 1.000 ## Proportion not Take -0.486 0.020 0.019 1.000 ## Scale -0.007 0.005 0.024 -0.017 1.000 ## Normality 0.007 -0.006 -0.006 -0.001 -0.001 1.000 At the top of page 523, Kruschke asked us to “notice that the posterior distribution in Figure 18.6 has ranges for the redundant parameters that are only a little smaller than their priors.” With a little wrangling, we can compare the prior/posterior distributions for our redundant parameters more directly. post %>% gather(parameter, posterior) %>% bind_cols( prior %>% gather() %>% transmute(prior = value) ) %>% gather(key, value, -parameter) %>% filter(parameter %in% c("Percent Take", "Proportion not Take")) %>% ggplot(aes(x = value, fill = key)) + geom_histogram(color = "grey92", size = .2, bins = 40, boundary = 0) + stat_pointintervalh(aes(y = 0), point_interval = mode_hdi, .width = c(.95, .5)) + scale_fill_viridis_d(option = "D", begin = .35, end = .65) + scale_y_continuous(NULL, breaks = NULL) + xlab(NULL) + theme(legend.position = "none") + facet_grid(key~parameter, scales = "free") Kruschke was right. The posterior distributions are only slightly narrower than the priors for those two. With our combination of data and model, we learned virtually nothing beyond the knowledge we encoded in those priors. Kruschke mentioned SEM as a possible solution to multicollinearity. brms isn’t fully capable of SEM, at the moment, but its multivariate syntax does allow for path analysis and IRT models. However, you can currently fit a variety of Bayesian SEMs with the blavaan package. I’m not aware of any textbooks highlighting blavaan. If you’ve seen any, please share. ### 18.1.5 Informative priors, sparse data, and correlated predictors. It’s worth reproducing some of Kruschke’s prose from this subsection. The examples in this book tend to use mildly informed priors (e.g., using information about the rough magnitude and range of the data). But a benefit of Bayesian analysis is the potential for cumulative scientific progress by using priors that have been informed from previous research. Informed priors can be especially useful when the amount of data is small compared to the parameter space. A strongly informed prior essentially reduces the scope of the credible parameter space, so that a small amount of new data implies a narrow zone of credible parameter values. (p. 523) ## 18.2 Multiplicative interaction of metric predictors From page 526: Formally, interactions can have many different specific functional forms. We will consider multiplicative interaction. This means that the nonadditive interaction is expressed by multiplying the predictors. The predicted value is a weighted combination of the individual predictors and, additionally, the multiplicative product of the predictors. For two metric predictors, regression with multiplicative interaction has these algebraically equivalent expressions: \begin{align*} \mu & = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_{1 \times 2} x_1 x_2 \\ & = \beta_0 + \underbrace{(\beta_1 + \beta_{1 \times 2} x_2)}_{\text{slope of } x_1} x_1 + \beta_2 x_2 \\ & = \beta_0 + \beta_1 x_1 + \underbrace{(\beta_2 + \beta_{1 \times 2} x_1)}_{\text{slope of } x_2} x_2. \end{align*} Figure 18.8 is out of our ggplot2 repertoire. ### 18.2.1 An example. With brms, you can specify an interaction with either the x_i*x_j syntax or the x_i:x_j syntax. I typically use x_i:x_j. It’s often the case that you can just make the interaction term right in the model formula. But since we’re fitting the model with standardized predictors and then using Kruschke’s equations to convert the parameters back to the unstandardized metric, it seems easier to make the interaction term in the data, first. my_data <- my_data %>% mutate(interaction = Spend * PrcntTake) %>% mutate(interaction_z = standardize(interaction)) Now we’ll fit the model. fit3 <- brm(data = my_data, family = student, satt_z ~ 1 + spend_z + prcnt_take_z + interaction_z, prior = c(prior(normal(0, 2), class = Intercept), prior(normal(0, 2), class = b), prior(normal(0, 1), class = sigma), prior(exponential(one_over_twentynine), class = nu)), chains = 4, cores = 4, stanvars = stanvar(1/29, name = "one_over_twentynine"), seed = 18) Note that even though an interaction term might seem different kind from other regression terms, it’s just another coefficient of class = b as far as the prior() statements are concerned. Anyway, let’s inspect the summary(). summary(fit3) ## Family: student ## Links: mu = identity; sigma = identity; nu = identity ## Formula: satt_z ~ 1 + spend_z + prcnt_take_z + interaction_z ## Data: my_data (Number of observations: 50) ## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; ## total post-warmup samples = 4000 ## ## Population-Level Effects: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## Intercept -0.00 0.06 -0.12 0.11 1.00 2779 2427 ## spend_z 0.04 0.14 -0.24 0.32 1.00 1484 2038 ## prcnt_take_z -1.48 0.29 -2.03 -0.90 1.00 1336 1985 ## interaction_z 0.59 0.37 -0.12 1.30 1.00 1217 1696 ## ## Family Specific Parameters: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## sigma 0.41 0.05 0.33 0.52 1.00 2837 2503 ## nu 34.38 28.67 4.52 110.70 1.00 2685 2567 ## ## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample ## is a crude measure of effective sample size, and Rhat is the potential ## scale reduction factor on split chains (at convergence, Rhat = 1). The correlations among our parameters are about as severe as those in the text. vcov(fit3, correlation = T) %>% round(digits = 3) ## Intercept spend_z prcnt_take_z interaction_z ## Intercept 1.000 -0.009 0.003 -0.001 ## spend_z -0.009 1.000 0.713 -0.835 ## prcnt_take_z 0.003 0.713 1.000 -0.962 ## interaction_z -0.001 -0.835 -0.962 1.000 We can see that the interaction variable is strongly correlated with both predictors. Therefore, we know that there will be strong trade-offs among the regression coefficients, and the marginal distributions of single regression coefficients might be much wider than when there was no interaction included. (p. 528) Let’s convert the posterior draws to the unstandardized metric. sd_x_3 <- sd(my_datainteraction)
m_x_3  <- mean(my_data$interaction) post <- posterior_samples(fit3) %>% transmute(Intercept = make_beta_0(zeta_0 = b_Intercept, zeta_1 = b_spend_z, zeta_2 = b_prcnt_take_z, zeta_3 = b_interaction_z, sd_x_1 = sd_x_1, sd_x_2 = sd_x_2, sd_x_3 = sd_x_3, sd_y = sd_y, m_x_1 = m_x_1, m_x_2 = m_x_2, m_x_3 = m_x_3, m_y = m_y), Spend = make_beta_j(zeta_j = b_spend_z, sd_j = sd_x_1, sd_y = sd_y), Percent Take = make_beta_j(zeta_j = b_prcnt_take_z, sd_j = sd_x_2, sd_y = sd_y), Spend : Percent Take = make_beta_j(zeta_j = b_interaction_z, sd_j = sd_x_3, sd_y = sd_y), Scale = sigma * sd_y, Normality = nu %>% log10()) glimpse(post) ## Observations: 4,000 ## Variables: 6 ##$ Intercept              <dbl> 1072.273, 1052.266, 1051.125, 1068.692, 1…
## $Spend <dbl> -1.6619618, 4.1428271, 2.3548164, -1.6373… ##$ Percent Take         <dbl> -4.815437, -4.600243, -4.150411, -4.25895…
## $Spend : Percent Take <dbl> 0.31441969, 0.22608990, 0.21413414, 0.255… ##$ Scale                  <dbl> 28.15976, 32.34262, 29.09982, 29.94091, 2…
## $Normality <dbl> 1.6180165, 1.1979350, 1.5089735, 1.569651… Now we’ve done the conversions, here are the histograms of Figure 18.9. post %>% gather() %>% ggplot(aes(x = value)) + geom_histogram(color = "grey92", fill = "grey67", size = .2, bins = 40) + stat_pointintervalh(aes(y = 0), point_interval = mode_hdi, .width = c(.95, .5)) + scale_y_continuous(NULL, breaks = NULL) + xlab(NULL) + facet_wrap(~key, scales = "free", ncol = 3) “To properly understand the credible slopes on the two predictors, we must consider the credible slopes on each predictor as a function of the value of the other predictor” (p. 528). This is our motivation for the middle panel of Figure 18.9. To make it, we’ll need to expand() our post, wrangle a bit, and plot with geom_pointrange(). # this will come in handy in expand() bounds <- range(my_data$PrcntTake)

# wrangle
post %>%
expand(nesting(Spend, Spend : Percent Take),
PrcntTake = seq(from = bounds[1], to = bounds[2], length.out = 20)) %>%
mutate(slope = Spend + Spend : Percent Take * PrcntTake) %>%
group_by(PrcntTake) %>%
median_hdi(slope) %>%

# plot
ggplot(aes(x = PrcntTake, y = slope,
ymin = .lower, ymax = .upper)) +
geom_hline(yintercept = 0, color = "white") +
geom_pointrange(color = "grey50") +
labs(x     = "Value of prcnt_take",
y     = "Slope on spend",
title = expression(paste("Slope on spend is ", beta[1] + beta[3] %.% "prcnt_take")))

That worked like a charm. We’ll follow the same basic order of operations for the final panel.

# this will come in handy in expand()
bounds <- range(my_data$Spend) # wrangle post %>% expand(nesting(Percent Take, Spend : Percent Take), Spend = seq(from = bounds[1], to = bounds[2], length.out = 20)) %>% mutate(slope = Percent Take + Spend : Percent Take * Spend) %>% group_by(Spend) %>% median_hdi(slope) %>% # plot ggplot(aes(x = Spend, y = slope, ymin = .lower, ymax = .upper)) + geom_pointrange(color = "grey50") + labs(x = "Value of spend", y = "Slope on prcnt_take", title = expression(paste("Slope on prcnt_take is ", beta[2] + beta[3] %.% "spend"))) Kruschke outlined all this in the opening paragraphs of page 530. His parting words of this subsection warrant repeating: “if you include an interaction term, you cannot ignore it even if its marginal posterior distribution includes zero” (p. 530). ## 18.3 Shrinkage of regression coefficients In some research, there are many candidate predictors which we suspect could possibly be informative about the predicted variable. For example, when predicting college GPA, we might include high-school GPA, high-school SAT score, income of student, income of parents, years of education of the parents, spending per pupil at the student’s high school, student IQ, student height, weight, shoe size, hours of sleep per night, distance from home to school, amount of caffeine consumed, hours spent studying, hours spent earning a wage, blood pressure, etc. We can include all the candidate predictors in the model, with a regression coefficient for every predictor. And this is not even considering interactions, which we will ignore for now. With so many candidate predictors of noisy data, there may be some regression coefficients that are spuriously estimated to be non-zero. We would like some protection against accidentally nonzero regression coefficients. (p. 530) That’s what this section is all about. We’ll make our random noise predictors with rnorm(). set.seed(18) my_data <- my_data %>% mutate(x_rand_1 = rnorm(n = n(), 0, 1), x_rand_2 = rnorm(n = n(), 0, 1), x_rand_3 = rnorm(n = n(), 0, 1), x_rand_4 = rnorm(n = n(), 0, 1), x_rand_5 = rnorm(n = n(), 0, 1), x_rand_6 = rnorm(n = n(), 0, 1), x_rand_7 = rnorm(n = n(), 0, 1), x_rand_8 = rnorm(n = n(), 0, 1), x_rand_9 = rnorm(n = n(), 0, 1), x_rand_10 = rnorm(n = n(), 0, 1), x_rand_11 = rnorm(n = n(), 0, 1), x_rand_12 = rnorm(n = n(), 0, 1)) glimpse(my_data) ## Observations: 50 ## Variables: 27 ##$ State           <chr> "Alabama", "Alaska", "Arizona", "Arkansas", "Cal…
## $Spend <dbl> 4.405, 8.963, 4.778, 4.459, 4.992, 5.443, 8.817,… ##$ StuTeaRat       <dbl> 17.2, 17.6, 19.3, 17.1, 24.0, 18.4, 14.4, 16.6, …
## $Salary <dbl> 31.144, 47.951, 32.175, 28.934, 41.078, 34.571, … ##$ PrcntTake       <dbl> 8, 47, 27, 6, 45, 29, 81, 68, 48, 65, 57, 15, 13…
## $SATV <dbl> 491, 445, 448, 482, 417, 462, 431, 429, 420, 406… ##$ SATM            <dbl> 538, 489, 496, 523, 485, 518, 477, 468, 469, 448…
## $SATT <dbl> 1029, 934, 944, 1005, 902, 980, 908, 897, 889, 8… ##$ prcnt_take_z    <dbl> -1.0178453, 0.4394222, -0.3078945, -1.0925770, 0…
## $spend_z <dbl> -1.10086058, 2.24370805, -0.82716069, -1.0612364… ##$ satt_z          <dbl> 0.8430838, -0.4266207, -0.2929676, 0.5223163, -0…
## $prop_not_take <dbl> 0.92, 0.53, 0.73, 0.94, 0.55, 0.71, 0.19, 0.32, … ##$ prop_not_take_z <dbl> 1.0178453, -0.4394222, 0.3078945, 1.0925770, -0.…
## $interaction <dbl> 35.240, 421.261, 129.006, 26.754, 224.640, 157.8… ##$ interaction_z   <dbl> -0.94113720, 0.93111798, -0.48635915, -0.9822954…
## $x_rand_1 <dbl> 0.92645924, 1.82282117, -1.61056690, -0.28510975… ##$ x_rand_2        <dbl> -0.90258025, -1.13163679, 0.49708131, -0.5477187…
## $x_rand_3 <dbl> 0.51576102, 0.30710965, 0.66199996, 2.21990655, … ##$ x_rand_4        <dbl> 1.08730491, -1.23909473, 0.43161390, 1.06733141,…
## $x_rand_5 <dbl> -0.23846777, 0.15702031, -1.02132795, 0.75395217… ##$ x_rand_6        <dbl> 0.06014956, 1.00555800, 1.47981871, -0.82827890,…
## $x_rand_7 <dbl> 1.46961709, 0.51790320, -2.33110353, 0.11339996,… ##$ x_rand_8        <dbl> 0.03463437, -1.48737599, -0.01528284, 0.48480309…
## $x_rand_9 <dbl> -0.4556078, -0.7035475, -0.5001913, -0.6526022, … ##$ x_rand_10       <dbl> 1.2858586, -0.7474640, -0.3107255, -1.1037468, 0…
## $x_rand_11 <dbl> 0.17236599, -0.37956084, 0.31982301, 0.29678108,… ##$ x_rand_12       <dbl> -0.53048519, 0.92465424, 0.66876661, 0.30935146,…

Here’s the naïve model.

fit4 <-
update(fit1,
newdata = my_data,
formula = satt_z ~ 1 + prcnt_take_z + spend_z + x_rand_1 + x_rand_2 + x_rand_3 + x_rand_4 + x_rand_5 + x_rand_6 + x_rand_7 + x_rand_8 + x_rand_9 + x_rand_10 + x_rand_11 + x_rand_12,
seed = 18)

Here we’ll examine the posterior with posterior_summary().

posterior_summary(fit4) %>%
round(digits = 2)
##                Estimate Est.Error   Q2.5  Q97.5
## b_Intercept       -0.02      0.07  -0.16   0.10
## b_prcnt_take_z    -1.12      0.09  -1.28  -0.95
## b_spend_z          0.33      0.09   0.16   0.50
## b_x_rand_1         0.03      0.06  -0.08   0.15
## b_x_rand_2         0.02      0.09  -0.16   0.20
## b_x_rand_3         0.09      0.07  -0.06   0.22
## b_x_rand_4        -0.09      0.07  -0.22   0.05
## b_x_rand_5         0.01      0.06  -0.12   0.13
## b_x_rand_6        -0.02      0.08  -0.17   0.12
## b_x_rand_7        -0.03      0.07  -0.17   0.11
## b_x_rand_8        -0.18      0.07  -0.31  -0.04
## b_x_rand_9         0.13      0.06   0.02   0.25
## b_x_rand_10        0.00      0.05  -0.10   0.11
## b_x_rand_11        0.05      0.09  -0.13   0.22
## b_x_rand_12       -0.08      0.06  -0.21   0.04
## sigma              0.38      0.07   0.23   0.52
## nu                25.09     26.21   2.15  95.21
## lp__             -54.70      3.49 -62.56 -49.11

Before we can make Figure 18.11, we’ll need to update our make_beta_0() function to accommodate this model.

make_beta_0 <-
function(zeta_0, zeta_1, zeta_2, zeta_3, zeta_4, zeta_5, zeta_6, zeta_7, zeta_8, zeta_9, zeta_10, zeta_11, zeta_12, zeta_13, zeta_14,
sd_x_1, sd_x_2, sd_x_3, sd_x_4, sd_x_5, sd_x_6, sd_x_7, sd_x_8, sd_x_9, sd_x_10, sd_x_11, sd_x_12, sd_x_13, sd_x_14, sd_y,
m_x_1, m_x_2, m_x_3, m_x_4, m_x_5, m_x_6, m_x_7, m_x_8, m_x_9, m_x_10, m_x_11, m_x_12, m_x_13, m_x_14, m_y) {
sd_y * zeta_0 + m_y - sd_y * ((zeta_1 * m_x_1 / sd_x_1) +
(zeta_2 * m_x_2 / sd_x_2) +
(zeta_3 * m_x_3 / sd_x_3) +
(zeta_4 * m_x_4 / sd_x_4) +
(zeta_5 * m_x_5 / sd_x_5) +
(zeta_6 * m_x_6 / sd_x_6) +
(zeta_7 * m_x_7 / sd_x_7) +
(zeta_8 * m_x_8 / sd_x_8) +
(zeta_9 * m_x_9 / sd_x_9) +
(zeta_10 * m_x_10 / sd_x_10) +
(zeta_11 * m_x_11 / sd_x_11) +
(zeta_12 * m_x_12 / sd_x_12) +
(zeta_13 * m_x_13 / sd_x_13) +
(zeta_14 * m_x_14 / sd_x_14))
}

Sigh, our poor make_beta_0() and make_beta_1() code is getting obscene. I don’t have the energy to think of how to wrap this into a simpler function. But someone probably should. If that ends up as you, do share your code.

sd_x_1  <- sd(my_data$Spend) sd_x_2 <- sd(my_data$PrcntTake)
sd_x_3  <- sd(my_data$x_rand_1) sd_x_4 <- sd(my_data$x_rand_2)
sd_x_5  <- sd(my_data$x_rand_3) sd_x_6 <- sd(my_data$x_rand_4)
sd_x_7  <- sd(my_data$x_rand_5) sd_x_8 <- sd(my_data$x_rand_6)
sd_x_9  <- sd(my_data$x_rand_7) sd_x_10 <- sd(my_data$x_rand_8)
sd_x_11 <- sd(my_data$x_rand_9) sd_x_12 <- sd(my_data$x_rand_10)
sd_x_13 <- sd(my_data$x_rand_11) sd_x_14 <- sd(my_data$x_rand_12)
sd_y    <- sd(my_data$SATT) m_x_1 <- mean(my_data$Spend)
m_x_2  <- mean(my_data$PrcntTake) m_x_3 <- mean(my_data$x_rand_1)
m_x_4  <- mean(my_data$x_rand_2) m_x_5 <- mean(my_data$x_rand_3)
m_x_6  <- mean(my_data$x_rand_4) m_x_7 <- mean(my_data$x_rand_5)
m_x_8  <- mean(my_data$x_rand_6) m_x_9 <- mean(my_data$x_rand_7)
m_x_10 <- mean(my_data$x_rand_8) m_x_11 <- mean(my_data$x_rand_9)
m_x_12 <- mean(my_data$x_rand_10) m_x_13 <- mean(my_data$x_rand_11)
m_x_14 <- mean(my_data$x_rand_12) m_y <- mean(my_data$SATT)

post <-
posterior_samples(fit4) %>%
transmute(Intercept      = make_beta_0(zeta_0 = b_Intercept,
zeta_1  = b_spend_z,
zeta_2  = b_prcnt_take_z,
zeta_3  = b_x_rand_1,
zeta_4  = b_x_rand_2,
zeta_5  = b_x_rand_3,
zeta_6  = b_x_rand_4,
zeta_7  = b_x_rand_5,
zeta_8  = b_x_rand_6,
zeta_9  = b_x_rand_7,
zeta_10 = b_x_rand_8,
zeta_11 = b_x_rand_9,
zeta_12 = b_x_rand_10,
zeta_13 = b_x_rand_11,
zeta_14 = b_x_rand_12,
sd_x_1  = sd_x_1,
sd_x_2  = sd_x_2,
sd_x_3  = sd_x_3,
sd_x_4  = sd_x_4,
sd_x_5  = sd_x_5,
sd_x_6  = sd_x_6,
sd_x_7  = sd_x_7,
sd_x_8  = sd_x_8,
sd_x_9  = sd_x_9,
sd_x_10 = sd_x_10,
sd_x_11 = sd_x_11,
sd_x_12 = sd_x_12,
sd_x_13 = sd_x_13,
sd_x_14 = sd_x_14,
sd_y    = sd_y,
m_x_1   = m_x_1,
m_x_2   = m_x_2,
m_x_3   = m_x_3,
m_x_4   = m_x_4,
m_x_5   = m_x_5,
m_x_6   = m_x_6,
m_x_7   = m_x_7,
m_x_8   = m_x_8,
m_x_9   = m_x_9,
m_x_10  = m_x_10,
m_x_11  = m_x_11,
m_x_12  = m_x_12,
m_x_13  = m_x_13,
m_x_14  = m_x_14,
m_y     = m_y),
Spend          = make_beta_j(zeta_j = b_spend_z,
sd_j   = sd_x_1,
sd_y   = sd_y),
Percent Take = make_beta_j(zeta_j = b_prcnt_take_z,
sd_j   = sd_x_2,
sd_y   = sd_y),
x_rand_1       = make_beta_j(zeta_j = b_x_rand_1,
sd_j   = sd_x_3,
sd_y   = sd_y),
x_rand_2       = make_beta_j(zeta_j = b_x_rand_2,
sd_j   = sd_x_4,
sd_y   = sd_y),
x_rand_3       = make_beta_j(zeta_j = b_x_rand_3,
sd_j   = sd_x_5,
sd_y   = sd_y),
x_rand_4       = make_beta_j(zeta_j = b_x_rand_4,
sd_j   = sd_x_6,
sd_y   = sd_y),
x_rand_5       = make_beta_j(zeta_j = b_x_rand_5,
sd_j   = sd_x_7,
sd_y   = sd_y),
x_rand_6       = make_beta_j(zeta_j = b_x_rand_6,
sd_j   = sd_x_8,
sd_y   = sd_y),
x_rand_7       = make_beta_j(zeta_j = b_x_rand_7,
sd_j   = sd_x_9,
sd_y   = sd_y),
x_rand_8       = make_beta_j(zeta_j = b_x_rand_8,
sd_j   = sd_x_10,
sd_y   = sd_y),
x_rand_9       = make_beta_j(zeta_j = b_x_rand_9,
sd_j   = sd_x_11,
sd_y   = sd_y),
x_rand_10      = make_beta_j(zeta_j = b_x_rand_10,
sd_j   = sd_x_12,
sd_y   = sd_y),
x_rand_11      = make_beta_j(zeta_j = b_x_rand_11,
sd_j   = sd_x_13,
sd_y   = sd_y),
x_rand_12      = make_beta_j(zeta_j = b_x_rand_12,
sd_j   = sd_x_14,
sd_y   = sd_y),
Scale          = sigma * sd_y,
Normality      = nu %>% log10())

glimpse(post)
## Observations: 4,000
## Variables: 17
## $Intercept <dbl> 978.2629, 954.1022, 936.3345, 927.7583, 964.4377,… ##$ Spend          <dbl> 16.23816, 21.85959, 26.38032, 26.63065, 18.40972,…
## $Percent Take <dbl> -3.123617, -3.470104, -3.496485, -3.622503, -3.31… ##$ x_rand_1       <dbl> 0.37557071, 0.47578746, 2.01992361, -2.49333200, …
## $x_rand_2 <dbl> -2.747509054, -0.637063446, 3.112304494, -4.64025… ##$ x_rand_3       <dbl> 8.9258844, 1.5485414, 3.0112088, 0.4144604, 11.01…
## $x_rand_4 <dbl> -8.1747088, -9.1809834, -13.4576786, -14.0589291,… ##$ x_rand_5       <dbl> -0.5797268, 0.1277984, 0.6989087, -2.2552537, -0.…
## $x_rand_6 <dbl> -6.499016480, 2.006864887, 2.694773713, 0.7167073… ##$ x_rand_7       <dbl> -7.4236507, 2.4005544, 0.5063171, 0.9188742, 10.2…
## $x_rand_8 <dbl> -11.9493748, -12.7780330, -12.2361496, -16.807961… ##$ x_rand_9       <dbl> 4.132977, 11.727639, 12.163225, 6.259811, 13.0476…
## $x_rand_10 <dbl> 0.19896636, 4.15941187, 4.02509269, 1.58380409, -… ##$ x_rand_11      <dbl> 9.3111126, -4.5582639, -8.9587840, -9.6417867, -4…
## $x_rand_12 <dbl> -8.3675915, -2.5825711, -3.9334521, -3.8179492, -… ##$ Scale          <dbl> 26.16620, 19.98710, 20.51592, 26.75018, 32.71402,…
## $Normality <dbl> 0.7053100, 0.7184234, 0.9159173, 0.4871434, 1.752… Okay, here are the histograms of Figure 18.11. post %>% select(Intercept:x_rand_3, x_rand_10:Normality) %>% gather() %>% mutate(key = factor(key, levels = c("Intercept", "Spend", "Percent Take", "x_rand_1", "x_rand_2", "x_rand_3", "x_rand_10", "x_rand_11", "x_rand_12", "Scale", "Normality"))) %>% ggplot(aes(x = value)) + geom_histogram(color = "grey92", fill = "grey67", size = .2, bins = 40) + stat_pointintervalh(aes(y = 0), point_interval = mode_hdi, .width = c(.95, .5)) + scale_y_continuous(NULL, breaks = NULL) + xlab(NULL) + facet_wrap(~key, scales = "free", ncol = 3) And here’s the final histogram depicting the Bayesian $$R^2$$. bayes_R2(fit4, summary = F) %>% as_tibble() %>% ggplot(aes(x = R2)) + geom_histogram(color = "grey92", fill = "grey67", size = .2, bins = 25) + stat_pointintervalh(aes(y = 0), point_interval = mode_hdi, .width = .95) + scale_y_continuous(NULL, breaks = NULL) + labs(subtitle = expression(paste("Bayesian ", italic(R)^2)), x = NULL) + coord_cartesian(xlim = c(.6, 1)) Note that unlike the one Kruschke displayed in the text, our brms::bayes_R2()-based $$R^2$$ distribution did not exceed the logical right bound of 1. Sometimes when you have this many parameters you’d like to compare, it’s better to display their summaries with an ordered coefficient plot. post %>% select(Spend:x_rand_12) %>% gather() %>% ggplot(aes(x = value, y = reorder(key, value))) + geom_vline(xintercept = 0, color = "white") + stat_pointintervalh(point_interval = mode_hdi, .width = .95, color = "grey50") + labs(x = NULL, y = NULL) + theme(axis.text.y = element_text(hjust = 0), axis.ticks.y = element_blank(), panel.grid.major.y = element_line(color = "grey98", linetype = 2)) Now we can see that by chance alone, the coefficients for x_rand_8 and x_rand_9 are clearly distinct from zero. The geom_density_ridges() function from the ggridges package can be informative, too. library(ggridges) post %>% select(Spend:x_rand_12) %>% gather() %>% ggplot(aes(x = value, y = reorder(key, value), group = reorder(key, value))) + geom_vline(xintercept = 0, color = "white") + geom_density_ridges(scale = 4, fill = "grey50", color = "grey50", alpha = 3/4, size = 1/5) + labs(x = NULL, y = NULL) + theme(axis.text.y = element_text(hjust = 0), axis.ticks.y = element_blank()) With brms, we can fit something like the model Kruschke displayed in Figure 18.12 with the horseshoe() prior. From the horseshoe section of the brms (version 2.10.3) reference manual: The horseshoe prior is a special shrinkage prior initially proposed by Carvalho et al. (2009). It is symmetric around zero with fat tails and an infinitely large spike at zero. This makes it ideal for sparse models that have many regression coefficients, although only a minority of them is non-zero. The horseshoe prior can be applied on all population-level effects at once (excluding the intercept) by using set_prior("horseshoe(1)"). The 1 implies that the student-t prior of the local shrinkage parameters has 1 degrees of freedom. Based on the quote, here’s how to fit our horseshoe-prior model. fit5 <- update(fit4, newdata = my_data, formula = satt_z ~ 1 + prcnt_take_z + spend_z + x_rand_1 + x_rand_2 + x_rand_3 + x_rand_4 + x_rand_5 + x_rand_6 + x_rand_7 + x_rand_8 + x_rand_9 + x_rand_10 + x_rand_11 + x_rand_12, prior = c(prior(normal(0, 2), class = Intercept), prior(horseshoe(1), class = b), prior(normal(0, 1), class = sigma), prior(exponential(one_over_twentynine), class = nu)), seed = 18) Check the parameter summary. posterior_summary(fit5) %>% round(digits = 2) ## Estimate Est.Error Q2.5 Q97.5 ## b_Intercept -0.02 0.06 -0.14 0.11 ## b_prcnt_take_z -1.02 0.10 -1.20 -0.82 ## b_spend_z 0.20 0.11 0.00 0.39 ## b_x_rand_1 0.02 0.04 -0.05 0.10 ## b_x_rand_2 0.01 0.04 -0.08 0.11 ## b_x_rand_3 0.02 0.04 -0.06 0.12 ## b_x_rand_4 -0.03 0.05 -0.14 0.04 ## b_x_rand_5 0.00 0.04 -0.07 0.09 ## b_x_rand_6 -0.02 0.04 -0.13 0.05 ## b_x_rand_7 -0.01 0.04 -0.09 0.07 ## b_x_rand_8 -0.09 0.07 -0.24 0.02 ## b_x_rand_9 0.05 0.05 -0.02 0.17 ## b_x_rand_10 0.00 0.03 -0.06 0.07 ## b_x_rand_11 0.01 0.04 -0.07 0.12 ## b_x_rand_12 -0.02 0.04 -0.12 0.05 ## sigma 0.40 0.06 0.27 0.52 ## nu 30.80 28.89 2.86 108.14 ## lp__ -105.49 6.71 -119.76 -93.33 ## hs_c2 1.89 3.56 0.35 7.79 Our make_beta_0() and make_beta_1() code remains obscene. post <- posterior_samples(fit5) %>% transmute(Intercept = make_beta_0(zeta_0 = b_Intercept, zeta_1 = b_spend_z, zeta_2 = b_prcnt_take_z, zeta_3 = b_x_rand_1, zeta_4 = b_x_rand_2, zeta_5 = b_x_rand_3, zeta_6 = b_x_rand_4, zeta_7 = b_x_rand_5, zeta_8 = b_x_rand_6, zeta_9 = b_x_rand_7, zeta_10 = b_x_rand_8, zeta_11 = b_x_rand_9, zeta_12 = b_x_rand_10, zeta_13 = b_x_rand_11, zeta_14 = b_x_rand_12, sd_x_1 = sd_x_1, sd_x_2 = sd_x_2, sd_x_3 = sd_x_3, sd_x_4 = sd_x_4, sd_x_5 = sd_x_5, sd_x_6 = sd_x_6, sd_x_7 = sd_x_7, sd_x_8 = sd_x_8, sd_x_9 = sd_x_9, sd_x_10 = sd_x_10, sd_x_11 = sd_x_11, sd_x_12 = sd_x_12, sd_x_13 = sd_x_13, sd_x_14 = sd_x_14, sd_y = sd_y, m_x_1 = m_x_1, m_x_2 = m_x_2, m_x_3 = m_x_3, m_x_4 = m_x_4, m_x_5 = m_x_5, m_x_6 = m_x_6, m_x_7 = m_x_7, m_x_8 = m_x_8, m_x_9 = m_x_9, m_x_10 = m_x_10, m_x_11 = m_x_11, m_x_12 = m_x_12, m_x_13 = m_x_13, m_x_14 = m_x_14, m_y = m_y), Spend = make_beta_j(zeta_j = b_spend_z, sd_j = sd_x_1, sd_y = sd_y), Percent Take = make_beta_j(zeta_j = b_prcnt_take_z, sd_j = sd_x_2, sd_y = sd_y), x_rand_1 = make_beta_j(zeta_j = b_x_rand_1, sd_j = sd_x_3, sd_y = sd_y), x_rand_2 = make_beta_j(zeta_j = b_x_rand_2, sd_j = sd_x_4, sd_y = sd_y), x_rand_3 = make_beta_j(zeta_j = b_x_rand_3, sd_j = sd_x_5, sd_y = sd_y), x_rand_4 = make_beta_j(zeta_j = b_x_rand_4, sd_j = sd_x_6, sd_y = sd_y), x_rand_5 = make_beta_j(zeta_j = b_x_rand_5, sd_j = sd_x_7, sd_y = sd_y), x_rand_6 = make_beta_j(zeta_j = b_x_rand_6, sd_j = sd_x_8, sd_y = sd_y), x_rand_7 = make_beta_j(zeta_j = b_x_rand_7, sd_j = sd_x_9, sd_y = sd_y), x_rand_8 = make_beta_j(zeta_j = b_x_rand_8, sd_j = sd_x_10, sd_y = sd_y), x_rand_9 = make_beta_j(zeta_j = b_x_rand_9, sd_j = sd_x_11, sd_y = sd_y), x_rand_10 = make_beta_j(zeta_j = b_x_rand_10, sd_j = sd_x_12, sd_y = sd_y), x_rand_11 = make_beta_j(zeta_j = b_x_rand_11, sd_j = sd_x_13, sd_y = sd_y), x_rand_12 = make_beta_j(zeta_j = b_x_rand_12, sd_j = sd_x_14, sd_y = sd_y), Scale = sigma * sd_y, Normality = nu %>% log10()) glimpse(post) ## Observations: 4,000 ## Variables: 17 ##$ Intercept      <dbl> 979.2471, 1023.5298, 1005.1208, 1019.9144, 987.87…
## $Spend <dbl> 10.433664, 6.800388, 8.659539, 6.694374, 13.33200… ##$ Percent Take <dbl> -2.647076, -2.530614, -2.792637, -2.863597, -2.82…
## $x_rand_1 <dbl> -0.122140323, -1.158510122, 4.129300096, 2.557700… ##$ x_rand_2       <dbl> -0.08683436, -7.80656229, -0.29081925, -0.4379702…
## $x_rand_3 <dbl> -0.52525371, 4.52651239, -1.31939726, -1.97089650… ##$ x_rand_4       <dbl> -0.14436314, 0.16171124, 1.24412546, 1.66430716, …
## $x_rand_5 <dbl> -1.00575083, 7.17913473, -0.49749734, -0.80289945… ##$ x_rand_6       <dbl> -7.70462095, -2.23838325, -1.33121990, -2.0296176…
## $x_rand_7 <dbl> -0.52338825, -0.97967618, -1.66890548, -3.6052574… ##$ x_rand_8       <dbl> -16.30586010, -4.60620698, 0.95248757, 1.40157515…
## $x_rand_9 <dbl> -0.06430865, 14.47257381, 0.02975212, 0.24323983,… ##$ x_rand_10      <dbl> -0.693127222, 0.282968557, -0.575203802, 0.764815…
## $x_rand_11 <dbl> 0.61800434, -0.14365478, 9.95488153, 13.63182415,… ##$ x_rand_12      <dbl> 0.004654011, -9.772196168, 0.691943833, 2.6793078…
## $Scale <dbl> 32.13506, 27.57153, 36.89317, 37.52294, 37.52031,… ##$ Normality      <dbl> 1.5010222, 1.2281283, 1.4576176, 1.5595392, 1.448…

And here are the majority of the histograms of Figure 18.12.

post %>%
select(Intercept:x_rand_3, x_rand_10:Normality) %>%
gather() %>%
mutate(key = factor(key,
levels = c("Intercept", "Spend", "Percent Take",
"x_rand_1", "x_rand_2", "x_rand_3",
"x_rand_10", "x_rand_11", "x_rand_12",
"Scale", "Normality"))) %>%

ggplot(aes(x = value)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 40) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
xlab(NULL) +
facet_wrap(~key, scales = "free", ncol = 3)

Based on the distributions for the random predictors, it looks like our brms horseshoe prior regularized more aggressively than Kruschke’s hierarchical prior in the text. And interestingly, look how our marginal posterior for Spend is bimodal.

For kicks and giggles, here’s the corresponding coefficient plot for $$\beta_1$$ through $$\beta_{14}$$.

post %>%
select(Spend:x_rand_12) %>%
gather() %>%

ggplot(aes(x = value, y = reorder(key, value))) +
geom_vline(xintercept = 0, color = "white") +
stat_pointintervalh(point_interval = mode_hdi, .width = .95, color = "grey50") +
labs(x = NULL,
y = NULL) +
theme(axis.text.y  = element_text(hjust = 0),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_line(color = "grey98", linetype = 2))

But anyways, here’s that final Bayesian $$R^2$$ histogram for Figure 18.12.

bayes_R2(fit5, summary = F) %>%
as_tibble() %>%

ggplot(aes(x = R2)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 25) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = .95) +
scale_y_continuous(NULL, breaks = NULL) +
labs(subtitle = expression(paste("Bayesian ", italic(R)^2)),
x = NULL) +
coord_cartesian(xlim = c(.6, 1))

Just recall, though, that our fit5 was not exactly like Kruschke’s model. Whereas we hard coded the scale of our Student-$$t$$ horseshoe prior to be 1, Kruschke estimated it with help from the gamma distribution. I’m not aware that’s currently possible in brms. If I’m at fault and you know how to do it, please share your code.

## 18.4 Variable selection

We can rewrite the linear regression model to accommodate whether it includes a predictor as

$\mu_i = \beta_0 + \sum_j \delta_j \beta_j x_{j, i},$

where $$\delta$$ is a dummy for which 0 = not included 1 = included. I’m not aware of a way to use $$\delta$$ as an inclusion indicator in brms the way Kruschke implemented it in JAGS. And in fact, it appears this might be unfeasible within the Stan framework. But if you know of a way, please share your code. However, this issue can lead to a similar approach: information criteria. To do so, let’s follow Kruschke’s basic flow and use the first model from way back in subsection 18.1.1 as a starting point. The model formula was as follows.

sd_x_4 <- sd(my_data$Salary) ## Now we'll start extracting our posterior samples and wrangling them, by model # fit_0100 posterior_samples(fit_0100) %>% transmute(PrcntTake = make_beta_j(zeta_j = b_prcnt_take_z, sd_j = sd_x_2, sd_y = sd_y)) %>% gather() %>% mutate(fit = "fit_0100") %>% # fit_0111 bind_rows( posterior_samples(fit_0111) %>% transmute(PrcntTake = make_beta_j(zeta_j = b_prcnt_take_z, sd_j = sd_x_2, sd_y = sd_y), StuTeaRat = make_beta_j(zeta_j = b_stu_tea_rat_z, sd_j = sd_x_3, sd_y = sd_y), Salary = make_beta_j(zeta_j = b_salary_z, sd_j = sd_x_4, sd_y = sd_y)) %>% gather() %>% mutate(fit = "fit_0111") ) %>% # fit_1100 bind_rows( posterior_samples(fit_1100) %>% transmute(Spend = make_beta_j(zeta_j = b_spend_z, sd_j = sd_x_1, sd_y = sd_y), PrcntTake = make_beta_j(zeta_j = b_prcnt_take_z, sd_j = sd_x_2, sd_y = sd_y)) %>% gather() %>% mutate(fit = "fit_1100") ) %>% # fit_0001 bind_rows( posterior_samples(fit_0001) %>% transmute(Salary = make_beta_j(zeta_j = b_salary_z, sd_j = sd_x_4, sd_y = sd_y)) %>% gather() %>% mutate(fit = "fit_0001") ) %>% # fit_1000 bind_rows( posterior_samples(fit_1000) %>% transmute(Spend = make_beta_j(zeta_j = b_spend_z, sd_j = sd_x_1, sd_y = sd_y)) %>% gather() %>% mutate(fit = "fit_1000") ) %>% # fit_0101 # spend_z + prcnt_take_z + stu_tea_rat_z bind_rows( posterior_samples(fit_0101) %>% transmute(PrcntTake = make_beta_j(zeta_j = b_prcnt_take_z, sd_j = sd_x_2, sd_y = sd_y), Salary = make_beta_j(zeta_j = b_salary_z, sd_j = sd_x_4, sd_y = sd_y)) %>% gather() %>% mutate(fit = "fit_0101") ) %>% # the next two lines just help order the grid the plots appear in mutate(key = factor(key, levels = c("Spend", "PrcntTake", "StuTeaRat", "Salary")), fit = factor(fit, levels = c("fit_0100", "fit_0111", "fit_1100", "fit_0001", "fit_1000", "fit_0101"))) %>% # finally, the plot! ggplot(aes(x = value, fill = fit)) + geom_histogram(color = "grey92", size = .2, bins = 40) + stat_pointintervalh(aes(y = 0), point_interval = mode_hdi, .width = c(.95, .5)) + scale_fill_viridis_d(option = "D", begin = .2, end = .8) + scale_y_continuous(NULL, breaks = NULL) + xlab(NULL) + theme(legend.position = "none") + facet_grid(fit~key, scales = "free") Like Kruschke’s results in the text, PrcntTake was the most prominant predictor. ### 18.4.3 What to report and what to conclude. Kruschke made a great point in the opening paragraph of this subsection. It might make sense to use the single most credible model, especially if it is notably more credible than the runner up, and if the goal is to have a parsimonious explanatory description of the data. But it is important to recognize that using the single best model, when it excludes some predictors, is concluding that the regression coefficients on the excluded predictors are exactly zero. (p. 546) Later he added “A forthright report should state the posterior probabilities of the several top models. Additionally, it can be useful to report, for each model, the ratio of its posterior probability relative to that of the best model” (p. 546). With our information criteria and model weights approach, we don’t have posterior probabilities for the models themselves. But we can report on their information criteria comparisons and weights. In the final paragraph of the subsection, Kruschke wrote: When the goal is prediction of $$y$$ for interesting values of the predictors, as opposed to parsimonious explanation, then it is usually not appropriate to use only the single most probable model. Instead, predictions should be based on as much information as possible, using all models to the extent that they are credible. This approach is called Bayesian model averaging (BMA). (p. 547) It’s worth it to walk this out a bit. With brms, one can use brms::pp_average() to get the weighted posterior distributions for the model parameters. This is a natural extension of our model weights comparisons. # how many points on the x-axis? n_points <- 30 # what vales of the predictors would we like to evaluate the weighted posterior over? nd <- tibble(spend_z = seq(from = -3, to = 3, length.out = n_points), prcnt_take_z = 0, stu_tea_rat_z = 0, salary_z = 0) pp <- # the first things we feed into pp_average() are the brm() fits we'd like to average over pp_average(fit_1111, fit_0111, fit_1011, fit_1101, fit_1110, fit_0011, fit_0101, fit_0110, fit_1001, fit_1010, fit_1100, fit_0001, fit_0010, fit_0100, fit_1000, fit_0000, # here we tell it to evaluate the posterior over these predictor values newdata = nd, # we can get the mean trends using the "fitted" method method = "fitted", # by robust, we mean we'd like the Estimate in terms of posterior medians, rather than means robust = T) str(pp) ## num [1:30, 1:4] -0.0164 -0.0164 -0.0164 -0.0163 -0.0164 ... ## - attr(*, "dimnames")=List of 2 ## ..$ : NULL
##   ..$: chr [1:4] "Estimate" "Est.Error" "Q2.5" "Q97.5" ## - attr(*, "weights")= Named num [1:16] 6.86e-06 2.49e-01 1.81e-05 3.95e-06 2.70e-03 ... ## ..- attr(*, "names")= chr [1:16] "fit_1111" "fit_0111" "fit_1011" "fit_1101" ... ## - attr(*, "nsamples")= Named num [1:16] 0 995 0 0 11 0 36 0 0 0 ... ## ..- attr(*, "names")= chr [1:16] "fit_1111" "fit_0111" "fit_1011" "fit_1101" ... The pp object will require a little wrangling before it’s of use for ggplot2. pp %>% as_tibble() %>% bind_cols(nd) %>% ggplot(aes(x = spend_z, y = Estimate, ymin = Q2.5, ymax = Q97.5)) + geom_pointrange(color = "grey50") + labs(x = "Value of Spend_z", y = "Standardized SATT") We can build on this to make a plot considering each of the four predictors. But first that requires we make a new nd tibble to feed into pp_average(). # how many points on the x-axis? n_points <- 30 # what vales of the predictors would we like to evaluate the weighted posterior over? nd <- tibble(spend_z = c(seq(from = -3, to = 3, length.out = n_points), rep(0, times = n_points * 3)), prcnt_take_z = c(rep(0, times = n_points), seq(from = -3, to = 3, length.out = n_points), rep(0, times = n_points * 2)), stu_tea_rat_z = c(rep(0, times = n_points * 2), seq(from = -3, to = 3, length.out = n_points), rep(0, times = n_points)), salary_z = c(rep(0, times = n_points * 3), seq(from = -3, to = 3, length.out = n_points))) pp <- pp_average(fit_1111, fit_0111, fit_1011, fit_1101, fit_1110, fit_0011, fit_0101, fit_0110, fit_1001, fit_1010, fit_1100, fit_0001, fit_0010, fit_0100, fit_1000, fit_0000, newdata = nd, method = "fitted", robust = T, # Note the probs argument probs = c(.025, .975, .1, .9, .25, .75)) str(pp) ## num [1:120, 1:8] -0.0163 -0.0163 -0.0163 -0.0163 -0.0163 ... ## - attr(*, "dimnames")=List of 2 ## ..$ : NULL
##   ..\$ : chr [1:8] "Estimate" "Est.Error" "Q2.5" "Q97.5" ...
##  - attr(*, "weights")= Named num [1:16] 6.86e-06 2.49e-01 1.81e-05 3.95e-06 2.70e-03 ...
##   ..- attr(*, "names")= chr [1:16] "fit_1111" "fit_0111" "fit_1011" "fit_1101" ...
##  - attr(*, "nsamples")= Named num [1:16] 0 995 0 0 11 0 36 0 0 0 ...
##   ..- attr(*, "names")= chr [1:16] "fit_1111" "fit_0111" "fit_1011" "fit_1101" ...

In each panel of the plot, below, we focus on one predictor. For that predictor, we hold all other three at their mean, which, since they are all standardized, is zero. We consider the posterior predictions for standardized SAT scores across a range of values each focal predictor. The posterior predictions are depicted in terms of 95%, 80%, and 50% percentile-based interval bands and a line at the median.

pp %>%
as_tibble() %>%
mutate(x = seq(from = -3, to = 3, length.out = n_points) %>% rep(., times = 4),
predictor = rep(c("Spend_z", "PrcntTake_z", "StuTeaRat_z", "Salary_z"), each = n_points)) %>%

ggplot(aes(x = x)) +
geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5, fill = predictor),
alpha = 1/5) +
geom_ribbon(aes(ymin = Q10, ymax = Q90, fill = predictor),
alpha = 1/4) +
geom_ribbon(aes(ymin = Q25, ymax = Q75, fill = predictor),
alpha = 1/3) +
geom_line(aes(y = Estimate, color = predictor),
size = 1) +
scale_fill_viridis_d(option = "D", begin = .1, end = .6) +
scale_color_viridis_d(option = "D", begin = .1, end = .6) +
labs(x = "Standardized value of the focal predictor",
y = "Standardized SATT") +
theme(legend.position = "none") +
facet_grid(predictor ~ ., scales = "free")

Based on the weighted average across the models, the PrcntTake_z predictor was the pmost potent.

### 18.4.4 Caution: Computational methods.

To conclude this section regarding variable selection, it is appropriate to recapitulate the considerations at the beginning of the section. Variable selection is a reasonable approach only if it is genuinely plausible and meaningful that candidate predictors have zero relation to the predicted variable. The results can be surprisingly sensitive to the seemingly innocuous choice of prior for the regression coefficients, and, of course, the prior for the inclusion probability. Because of these limitations, hierarchical shrinkage priors may be a more meaningful approach. (p. 548)

### 18.4.5 Caution: Interaction variables.

When interaction terms are included in a model that also has hierarchical shrinkage on regression coefficients, the interaction coefficients should not be put under the same higher-level prior distribution as the individual component coefficients, because interaction coefficients are conceptually from a different class of variables than individual components. (pp. 548–549)

## Session info

sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base
##
## other attached packages:
##  [1] ggridges_0.5.1  bayesplot_1.7.0 tidybayes_1.1.0 brms_2.10.3
##  [5] Rcpp_1.0.2      patchwork_1.0.0 forcats_0.4.0   stringr_1.4.0
##  [9] dplyr_0.8.3     purrr_0.3.3     readr_1.3.1     tidyr_1.0.0
## [13] tibble_2.1.3    ggplot2_3.2.1   tidyverse_1.2.1
##
## loaded via a namespace (and not attached):
##   [1] colorspace_1.4-1          ellipsis_0.3.0
##   [3] rsconnect_0.8.15          ggstance_0.3.2
##   [5] markdown_1.1              base64enc_0.1-3
##   [7] rstudioapi_0.10           rstan_2.19.2
##   [9] svUnit_0.7-12             DT_0.9
##  [11] fansi_0.4.0               lubridate_1.7.4
##  [13] xml2_1.2.0                bridgesampling_0.7-2
##  [15] mnormt_1.5-5              knitr_1.23
##  [17] shinythemes_1.1.2         zeallot_0.1.0
##  [19] jsonlite_1.6              broom_0.5.2
##  [21] shiny_1.3.2               compiler_3.6.0
##  [23] httr_1.4.0                backports_1.1.5
##  [25] assertthat_0.2.1          Matrix_1.2-17
##  [27] lazyeval_0.2.2            cli_1.1.0
##  [29] later_1.0.0               htmltools_0.4.0
##  [31] prettyunits_1.0.2         tools_3.6.0
##  [33] igraph_1.2.4.1            coda_0.19-3
##  [35] gtable_0.3.0              glue_1.3.1.9000
##  [37] reshape2_1.4.3            cellranger_1.1.0
##  [39] vctrs_0.2.0               nlme_3.1-139
##  [41] crosstalk_1.0.0           psych_1.8.12
##  [43] xfun_0.10                 ps_1.3.0
##  [45] rvest_0.3.4               mime_0.7
##  [47] miniUI_0.1.1.1            lifecycle_0.1.0
##  [49] gtools_3.8.1              MASS_7.3-51.4
##  [51] zoo_1.8-6                 scales_1.0.0
##  [53] colourpicker_1.0          hms_0.4.2
##  [55] promises_1.1.0            Brobdingnag_1.2-6
##  [57] parallel_3.6.0            inline_0.3.15
##  [59] shinystan_2.5.0           yaml_2.2.0
##  [61] gridExtra_2.3             loo_2.1.0
##  [65] dygraphs_1.1.1.6          pkgbuild_1.0.5
##  [67] rlang_0.4.1               pkgconfig_2.0.3
##  [69] matrixStats_0.55.0        HDInterval_0.2.0
##  [71] evaluate_0.14             lattice_0.20-38
##  [73] rstantools_2.0.0          htmlwidgets_1.5
##  [75] labeling_0.3              tidyselect_0.2.5
##  [77] processx_3.4.1            plyr_1.8.4
##  [79] magrittr_1.5              R6_2.4.0
##  [81] generics_0.0.2            foreign_0.8-71
##  [83] pillar_1.4.2              haven_2.1.0
##  [85] withr_2.1.2               xts_0.11-2
##  [87] abind_1.4-5               modelr_0.1.4
##  [89] crayon_1.3.4              arrayhelpers_1.0-20160527
##  [91] utf8_1.1.4                rmarkdown_1.13
## [103] shinyjs_1.0


text:
# Is corrosion a physical change?

Oct 17, 2016

Corrosion is undoubtedly an example of $\text{chemical change}$.

#### Explanation:

$\text{Chemical change}$ is characterized by the formation of new substances and the making and breaking of chemical bonds. Both processess occur upon corrosion.

See this older answer regarding the importance of rust in the building of monumental structures.

Oct 17, 2016

Corrosion is a chemical change.

#### Explanation:

The most sure way to determine whether a physical or chemical change has occurred is to take a before/after "snapshot" around the change.

If the chemical composition of a material changes, then a chemical change has occurred.

Corrosion could be defined as, "... damage or destruction of living tissue or material due to a chemical or electrochemical reaction" (Chemistry.about.com).

When we boil water, the snapshot would reveal that the structure has not changed, just the distance between the atoms. Therefore, only a physical change has occurred.

When we look at an iron plate left outside, we can tell by the snapshot that oxygen has somehow worked itself into bonds with the iron. The oxygen wasn't there before! Hello rust! iron bonding with oxygen is different from iron bonding with iron, so a chemical change has occurred.

format: plain


text:
, 28.10.2019 16:29 rhaineandreirefuerzo

# In a subtraction problem, which is smaller, subtrahend or minuend?

### Another question on Math

Math, 28.10.2019 16:29
Write a mathematical model in the form of an algebraic expression reprensenting the word statements.use x as the variable 1.the sum of a certain number and 23. 2.46 less thanvdome number. 3.the product of a certain number and 45. 4.101 fewer tthan a certain number. 5.5 less than 3 times a number. 6.17 more than 7 times a number 7.15 less that twice a certain number 8.the product of 17 and the square of the difference when 7 is subtracted from a number 9.the quotient when the fourth power of the sum of a number and 16 divided by 5. 10.the sum of a number and 11,mutiplied by the difference when 14 is subtracted from the number. 11.the total number and twice 5,divided by the total of the number and 5. 12.the perimeter (p) of a square is for times its side (s)
Math, 28.10.2019 17:28
$$\sqrt[2]{8} {x}^{2} \times \sqrt[2]{12x {y}^{2} }$$


text:
# The field and landscape of affordances: Koffka’s two environments revisited

Synthese, Mar 2019

The smooth integration of the natural sciences with everyday lived experience is an important ambition of radical embodied cognitive science. In this paper we start from Koffka’s recommendation in his Principles of Gestalt Psychology that to realize this ambition psychology should be a “science of molar behaviour”. Molar behavior refers to the purposeful behaviour of the whole organism directed at an environment that is meaningfully structured for the animal. Koffka made a sharp distinction between the “behavioural environment” and the “geographical environment”. We show how this distinction picks out the difference between the environment as perceived by an individual organism, and the shared publicly available environment. The ecological psychologist James Gibson was later critical of Koffka for inserting a private phenomenal reality in between animals and the shared environment. Gibson tried to make do with just the concept of affordances in his explanation of molar behaviour. We argue however that psychology as a science of molar behaviour will need to make appeal both to the concepts of shared publicly available affordances, and of the multiplicity of relevant affordances that invite an individual to act. A version of Koffka’s distinction between the two environments remains alive today in a distinction we have made between the field and landscape of affordances. Having distinguished the two environments, we go on to provide an account of how the two environments are related. Koffka suggested that the behavioural environment forms out of the causal interaction of the individual with a pre-existing, ready-made geographical environment. We argue that such an account of the relation between the two environments fails to do justice to the complex entanglement of the social with the material aspects of the geographical environment. To better account for this sociomaterial reality of the geographical environment, we propose a process-perspective on our distinction between the landscape and field of affordances. While the two environments can be conceptually distinguished, we argue they should also be viewed as standing in a relation of reciprocal and mutual dependence.

This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs11229-019-02123-x.pdf

Julian Kiverstein, Ludger van Dijk, Erik Rietveld. The field and landscape of affordances: Koffka’s two environments revisited, Synthese, 2019, 1-18, DOI: 10.1007/s11229-019-02123-x


text:
Is global warming even real?

Discussion in 'Earth Science' started by Ilikeponies579, Dec 16, 2014.

1. billvonValued Senior Member

Messages:
18,112
You give the denier political movement a lot of ammunition to claim that alarmists make absurd, unsupported assertions and should be ignored.

3. iceauraValued Senior Member

Messages:
30,244
The denier claims have nothing to do with "ammunition".

Trying to edit or curb discussions of climate change risks and possibilities so that deniers cannot somehow find "ammunition" in them is wasted effort.

Messages:
10,890

7. billvonValued Senior Member

Messages:
18,112
I disagree. For example, the sensationalism with which the "global cooling" studies were reported in the 1970's has damaged, and continues to damage, the cause of public education on the topic of climate change. Newsweek, to their credit, apologized for their sensationalism later - and admitted it led many to draw the wrong conclusions about climate change.

8. iceauraValued Senior Member

Messages:
30,244
Which goes to show that nothing you do can prevent the denial operation from finding "ammunition" in whatever you say.

Are supposed to ignore the incoming glaciation prospects, or the methane burst hazard, or the increased odds of a short-term limited runaway hydropowered warming, or the chances of a Gulf Stream shutdown leading to arctic weather in Ireland, and the like, because such concerns can be misrepresented by the unscrupulous and dingbat?

Are we supposed to edit our statistical discussions to avoid using the word "trick", so that people who steal our emails cannot lie about what we meant?

Let's rather leave the propaganda considerations to the propagandists.

9. billvonValued Senior Member

Messages:
18,112
Agreed. And in fact I don't think there's anything you can say to a true denier that will convince him. (Indeed, at least one study has shown that they become more entrenched the more information is presented to them.)

However, the vast majority of Americans are not deniers - nor are they the sort of alarmists that think that 90%+ of humanity will be killed by climate change. These are the people who hear the argument "well, scientists once predicted an ice age - and now they are predicting a heat wave! They'll probably change their minds again in a few years" - and think that's a reasonable argument. And Newsweek, Time et al gave deniers the ammunition to make that argument.
Not at all. We should, rather, present the science that we can support, rather than the more spectacular/sensationalist guesses about our future.
Agreed. And predicting events where "half of humanity will die in less than a month" or "the temperature rises to the extinction level" is propaganda rather than science.

10. Billy TUse Sugar Cane Alcohol car FuelValued Senior Member

Messages:
23,198
NE Brail has 56 cities, most relatively small with out water, except that delivered by potable water trucks.

Here is Sao Paulo things are looking better. The month of March dumped more water into our main reservoir than any year since 2008 and it holds 19% of capacity now. - way up from the ~5% low of a few months ago. There have been many large floods in the low areas of Sao Paulo - worst I have seen in 22 years. - Cars rolled down the street by the run off waters. Yesterday a man died in one.

I have started to use some of the 90 liters of bottled water I bought.

Last edited by a moderator: Apr 2, 2015
11. Billy TUse Sugar Cane Alcohol car FuelValued Senior Member

Messages:
23,198
No I don't. I give information extremely well supported with EVIDENCE. I admit evidence will not open well closed mind, but did not realize yours was.
See the evidenc in this less than 10 minutes, then comment on the evidence or correct my comments on it below:

"Boiled" refers to video's violent rise of CH4 bubbles, not 100C.
This is the most complete SCIENTIFIC investigation yet done of the cause of the "Great Extinction"

See highly qualified author extracting from sediment cores tiny pieces from the worse ever known mass extinction period and how he processed them to come to his conclusion. Did you see the huge storage rooms with tens of thousands of racks of long core tubes? He is an expert - otherwise would never have been allowed to extract almost irreplaceable data in these ancient sea bottom cores.

The conclusions ARE "substantiated" - It is your claim they are "absurd" that is without any basis than OPINION. Fact that they are not commonly held by other experts is probably due to fact these "other experts" have not had access to the ancient mud cores, they are based on.

This conclusion is his, not mine. The most complete EVIDENCE yet analysized, does indicate we are starting the same surge up in CH4 that made extinct 95% of all creatures, - main difference is with at least ten fold more rapid CO2 release that has already driven the CH4 concentration levels three times higher that at any time in the last 800,000 years. (That comes from ice core data, not the much more distant past recorded in the sea sediment mud.)

Not only does his chemical and fossil changes (found in sediment cores) analysis lead to the same conclusion many others have reached: that the "Siberian traps" were bad for life on earth, but not the primary cause of the 95% extinction, which as most experts now agree via a CH4 surge lasting a few thousands of year with rapid and about a 10C temperature increase.

He, with his better time resolution, was able to show it was really three separate surges. The first mass of CH4 released died way, making CO2 at least one molecule of CO2 for each of CH4 that was destroyed. That surge of longer lasting CO2, was probably, like the present being initially absorbed in the Oceans (my guess, he does not say). But the Global Warming of that CO2 surge made a second CH4 surge and process repeated to make a third. If history repeats, our "2nd CH4 surge" should begin about 2050 as more of the oceans stored CO2 comes back into the air.

Last edited by a moderator: Apr 2, 2015
12. billvonValued Senior Member

Messages:
18,112
So did Time Magazine and Newsweek back in the 1970's. It took them decades to apologize, but they finally did.

13. Billy TUse Sugar Cane Alcohol car FuelValued Senior Member

Messages:
23,198
Thanks. I agree. Here is small part of the "causes" section:

"Intergovernmental Panel on Climate Change, a group of 1,300 independent scientific experts from countries all over the world under the auspices of the United Nations, concluded there's a more than 90 percent probability that human activities over the past 250 years have warmed our planet. ... The panel also concluded there's a better than 90 percent probability that human-produced greenhouse gases such as carbon dioxide, methane and nitrous oxide have caused much of the observed increase in Earth's temperatures over the past 50 years.

They said the rate of increase in global warming due to these gases is very likely to be unprecedented within the past 10,000 years or more. The panel's full Summary for Policymakers report is online at http://www.ipcc.ch/pdf/assessment-report/ar4/syr/ar4_syr_spm.pdf."

It is my opinion that the ICCP has strong tendency to understate the AGW problem, in part because governments get to review their report, even edit, before they are released. Slowly AGW will force the government to spend more than they like. For example, Governor Brown of CA, yesterday announce a billion dollar program to address CA's water problem - today noting that the snow levels in mountains from which half* of CA's H2O consumption comes from has never been so low in the 65 years that records have been keep.

When ICCP first's report estimated when Arctic would be ice free at end summer / early fall it was a few decades into the 2100s. I think that will have happen before 2020, as do some others. Their "consensus" track recorded is terrible!

Also IMO, they are grossly underestimating the longer term danger of CH4 - seem even to be ignorant of fact the half-life of this very powerful GHG is now increasing at the rate of 0.3years each year; in large part because now its release rate is greater than rate harsh solar UV can can make the OH- radical which is mainly what removes CH4 from the air (mutual destruction reactions) Now the CH4 is rising - is 3 times greater in the last few decade than even the highest concentrationin the last 800,000 year as concentration of OH- is falling. For more than 800,000 year solar UV could make OH- faster than CH4 was release, so the "surplus" OH- held CH4 concentration low - now the reverse is true - greater release of CH4 is driving the OH- concentration down. Now each CH4 molecule remains in the air longer before it finds an OH- to react with. Soon, 1 Kg of CH4 will cause more GW than 100 Kg of CO2 during the decade following "puffs" of them released simultaneously.

* The other half of CA's water consumption is being "mined" from the ground (paired satellite data - one trailing the other in LEO, so relative motion give measure of mass concentrations in near surface ground - They have been "seeing" the ground water decrease quantatively for more than 5 years.) This falling "water table" is no news to many now whose wells have gone dry.

Last edited by a moderator: Apr 2, 2015
14. Billy TUse Sugar Cane Alcohol car FuelValued Senior Member

Messages:
23,198
That is not a comment on the video or my comments on it. Yes 45 years go man did not have either the data or the understanding of global warming he does now. So this is just a "duck and weave" to avoid commenting on the evidence or substance of the conclusions, and not even a good one.

15. billvonValued Senior Member

Messages:
18,112
Right. But they did understand global cooling from high altitude aerosols, and accurately described the phenomenon. Journalists then extrapolated from those facts and predicted calamity.

"But the scientists see few signs that government leaders anywhere are even prepared to take the simple measures of stockpiling food or of introducing the variables of climatic uncertainty into economic projections of future food supplies...The longer the planners (politicians) delay, the more difficult will they find it to cope with climatic change once the results become grim reality. . . . the resulting famines could be catastrophic" -"drought and desolation" - "the most devastating outbreak of tornadoes ever recorded" - "droughts, floods, extended dry spells, long freezes, delayed monsoons" - "impossible for starving peoples to migrate" - "the present decline has taken the planet about a sixth of the way toward the Ice Age."

You are doing something very similar by posting facts about the forcing properties of methane, then extrapolating to state that "the temperature rises to the extinction level" and the possibility that "half of humanity will die in less than a month."

16. zgmcRegistered Senior Member

Messages:
831

Messages:
2,527
I'll ask some of my buds. Hell, Trippy might be able to get a copy, but he wouldn't be able to share it here.

18. milkweedValued Senior Member

Messages:
1,654

Or you could email the author himself and see if a copy is available.

One of my fav comments so far:

19. Billy TUse Sugar Cane Alcohol car FuelValued Senior Member

Messages:
23,198
Yes, I am saying: "History* can repeat." And the less than 10 minute video of post 638, the latest and most complete analysis of the "Great Extinction" (~95% of all creatures dead) tells why that is to be expected, without drastic changes in man's energy systems.

Switching liquid fuel requiring cars to sugar cane based alcohol, being an excellent way to start as:

(1) It is renewable, not finite.
(2) All needed can be grown only using abandoned pasture land.
Much of that is in poor countries and cutting cane is a low skill job. Many now with no salary would get jobs and become customers for first world products.
A win/win change for all but big oil.
(3) Switching to sugar cane based fuel, not only greatly reduces the amount of CO2 annually released, it actually REMOVES some from the air. (All the carbon in cane's roots, ocean tankers filled with alcohol, large storage tanks at ports, and millions of smaller ones in cars, was earlier removed from the air by the growing cane.) I. e. sugar cane alcohol is a slightly "net negative CO2 release" fuel.
(4) Cane based alcohol is cheaper per mile driven that gasoline when petroleum cost >$70/ barrel without ANY government subsidy. (or >$65/barrel if big oil's "depletion allowance" tax subsidy were ended. Big oil has other subsidies too. Their lobbyist are well paid and good.) Sugar cane alcohol has more than 100 independent producers, competing with each other, with no government subsidy**
(5) Sugar cane alcohol gives slightly more power than gasoline, in same IC motor and significantly more in one designed with higher compression as alcohol does not "knock" as easily.
(6) A minor advantage is it is clearer burning - and reduces need to change spark plugs, get other repairs.

Only disadvantage vs. gasoline is that a full tank will only have 70% of the driving range that gasoline does. So you will need to fill tank more often.

* I never said "half of humanity will die in a month." - That is non-scientific "alarmist talk" - I don't do that; however, there will be locally heat waves where thousands die in a month - Wet bulb temperature of 35C (95F) is lethal to humans in less than an hour. For every 1 degree C of ocean surface temperature rise there is 7% more water vapor in the air in equilibrium with salt water's increased vapor pressure.

At time of the great extinction, the temperature was already about 5C higher than now and in less than 1000 years went up about 8 degrees more at the peak, but most were already extinct. There are several ways to estimate the temperature back then and not all agree on these values, given in the video, but all do agree that the rise, if three surges of CH4 concentration do again occur - there will be a lage extinction again. 99+ percent of humans now are part of some social structure, and may not be able to survive with out it.

I like to note that the Bible is correct: "The meek will inherit the Earth."
I.e. tiny mice with large surface to volume ratios for cooling and habit of staying in their burros during the heat of the day, are the most probable "meek" - The most advanced survivors inheriting the earth, but humans living north of the arctic circle, probably will live decades longer than most. Heat & humidity will not kill them - Acidified oceans and lack of food locally grown probably will be the cause of their deaths, if they too go extinct. However, there will be thousand of un-buried bodies (not just human) and plagues may get them instead of lack of food.

** Actually in Brazil, until a few days ago alcohol competed and won, with gasoline, even with price of gasoline suppressed by the populist government, which owns more than half of PetroBrass. They could not suppress the price of alcohol too as it is produced by dozens of independent producers. All the government could do was remove some alcohol supply by requiring gasoline to have 25% alcohol content.

Brazil's socialistic government has done so much damage to the economy that collapse is a real concern of many and finally price of electricity and gasoline has been increased, both by more than 10%. This of course has dropped the president's approval rating to ~20% despite having won re-election just a couple of months ago. ~60% think she is doing a "terrible job" now. I was one of the more than one million who marched down Sao Paulo's most important street (Av. Paulista) three Sundays ago in massive protest against the socialistic government (~3 million total marching in all of Brazil that Sunday.) The presidential "spin" on the marches was: "This shows the world how strong Brazil's democracy is."

Last edited by a moderator: Apr 2, 2015
20. billvonValued Senior Member

Messages:
18,112
The Great Extinction (P-Tr) took over a million years. If your claim is that, over the next million years, 95% of all currently-existing creatures will become extinct (with new ones evolving to fill the gap) then I would not only agree with that, but I'd say we are hastening it. Indeed, mankind as we know it will almost certainly be extinct, replaced by whatever humanity evolves into. (Keep in mind that a million years ago, Homo Sapiens did not exist.)

However generally you are not saying that. You are saying that it could happen in a month or so. You say that forests are drying out, and that might lead to large forest fires. And from the resulting soot "about half of humanity will die in less than a month." (And since you don't believe you said that, here's the link to the post: http://www.sciforums.com/threads/runaway-global-warming.142102/page-26#post-3250191)

THAT is the kind of non-scientific "alarmist talk" that will make you (and by extension, other people promoting both climate change education and altenative fuels) a laughingstock. In a few years, if your posts were to get any sort of real-world traction, a denier publication would state "remember when Billy T claimed that forest fires from the drought would have half of humanity dead in a month? Check out the weather in Brazil - floods! And now climate scientists are saying that the storms might be due to climate change. What a bunch of clueless alarmists!" And while that will do nothing to convince people on either side, those statements (backed up by quotes from both you and climate scientists) may convince some moderates that climate change scientists are essentially always wrong.

21. Billy TUse Sugar Cane Alcohol car FuelValued Senior Member

Messages:
23,198
I don't know the time scale for different creatures to go extinct, but think for most the time between there population be reduced by 50% until they were extinct was on the order of a few thousand years, at most, but that is just my guess.
Yes I have noted that a much more rapid and wide spread death than normally is considered possible MIGHT BE POSSIBLE, via Hadley cell circulation (between + & - 30 degree latitude) pumping huge amount of soot into high clouds, cutting their reflectivity in half (effectively same a "step function" in increased of solar absorption heating of the atmosphere). I have speculated that IF for example Amazon's routine fires did not burn out as they normally do but spread in the increasing drought conditions* and merge to burn 50% of the Amazon rain forest in a week or so, many at least in the tropics would die in a month or so. That is a SPECULATIVE POSSIBILITY, I have invented, mainly to show, how uncertain and serious AGW may be in the near term - not 30 to 90 years form now. This is very different from the near certainty (without drastic change in man's energy system,) of wide spread deaths, if not eventual extinction (a repeat of history with same mechanism) as the video makes clear.

*56 Brazilian cities, quite a few in the Amazon rain forest, I think, now have no water to distribute in the distribution system pipes. Potable water trucks are preventing wide spread deaths.

Here quote of the full paragraph from your link where I discussed the "Hadley cell / large forest burning" theat. with two parts you seem to have missed now made bold and larger:

"The forests are drying, and trees are dying - more CO2 released than taken up, but the worst thing is that large part of the Rain Forests may burn - sending soot high up via Hadley cell, that cover half the earth (From -30 S to +30 N latitude). Currently at least 2/3 of the solar energy on them is reflected back into space but with only about 1 soot particle in ~1000 water droplets their absorption would double - 100% increase in solar energy being absorbed by them! If that happens, about half of humanity will die in less than a month. "

That is an indication of a POSSIBILITY not any prediction. I other posts I have suggested the AGW problem is so complex, the mankind may not have even recognized the main threats. Hell, I was able to think of this one, which as far as I know, no one else has mentioned the potential of the Hadley cells to greatly amplify the danger of wide spread forest fires in the tropics. Most discussion of soot being a positive feed back mechanism refers to it speeding the melting of ice on Greenland etc.

Last edited by a moderator: Apr 2, 2015
22. TrippyALEA IACTA ESTStaff Member

Messages:
10,890
No.

Benton has suggested that there were extinction events at the beginning and the end of the Guadaloupian epoch leading to a decline in biodiversity, or that biodiversity failed to recover in the leadup to the end permian event - in the same vein as Erwin and Teichert who seperately suggested that the extinction was prolonged throughout the Lopingian epoch. Rampino et al suggest, using high resolution cyclostratigraphy, time scales on the order of 10^3 - 10^4 for width. Wang and Everson suggest, based on Chinese Ostracod data that there were two pulses seperated by 720,000 - 1,200,000 years. The width proposed by Rampino seems to be supported by Twitchet who says in addition to that that it took a further 10^5 years full the full faunal effects to play out, and that the collapse of terrestrial and marine ecosystems began at the same time and precedded the delta[sup]13[/sup]C excursion. There is an older theory that there was a pulse at the end of the guadloupian followed by the larger pulse at the end of the permian that recently seems to have been regaining some traction.

So, to me at any rate, it looks like the way the evidence is playing out is that the end permian extinction was the result of a sequence of events that played out over a period of at least ten million years starting with at least one extinction pulse at the end of the Guadaloupian (which may or may not have been preceeded by one ten million years earlier at the start of the Guadaloupian), followed by an enhanced extinction rate in the last ten million years of the Permian (during the Lopingian) which terminated in two more extinction pulses of duration 10^3-10^4 years seperated by a period of 10^5-10^6 years.

23. iceauraValued Senior Member

Messages:
30,244
And those people will hear such arguments from deniers regardless of how people genuinely discussing issues edit themselves.

So editing actual discussion, in an attempt to prevent misrepresentation by deniers and the like, is vain.

So you are sure that such events are so very, very unlikely the possibility of them need not be even considered? You are completely sure of that?