## 10.4 Simulation-based inference for regression

Recall in Subsection 10.2.5 when we interpreted the third through seventh columns of a regression table, we stated that R doesn’t do simulations to compute these values. Rather R uses theory-based methods that involve mathematical formulas.

In this section, we’ll use the simulation-based methods you previously learned in Chapters 8 and 9 to recreate the values in the regression table in Table 10.1. In particular, we’ll use the `infer`

package workflow to

- Construct a 95% confidence interval for the population slope \(\beta_1\) using bootstrap resampling with replacement. We did this previously in Sections 8.4 with the
`pennies`

data and 8.6 with the`mythbusters_yawn`

data. - Conduct a hypothesis test of \(H_0: \beta_1 = 0\) versus \(H_A: \beta_1 \neq 0\) using a permutation test. We did this previously in Sections 9.3 with the
`promotions`

data and 9.5 with the`movies_sample`

IMDb data.

### 10.4.1 Confidence interval for slope

We’ll construct a 95% confidence interval for \(\beta_1\) using the `infer`

workflow outlined in Subsection 8.4.2. Specifically, we’ll first construct the bootstrap distribution for the fitted slope \(b_1\) using our single sample of 463 courses:

`specify()`

the variables of interest in`evals_ch5`

with the formula:`score ~ bty_avg`

.`generate()`

replicates by using`bootstrap`

resampling with replacement from the original sample of 463 courses. We generate`reps = 1000`

replicates using`type = "bootstrap"`

.`calculate()`

the summary statistic of interest: the fitted`slope`

\(b_1\).

Using this bootstrap distribution, we’ll construct the 95% confidence interval using the percentile method and (if appropriate) the standard error method as well. It is important to note in this case that the bootstrapping with replacement is done *row-by-row*. Thus, the original pairs of `score`

and `bty_avg`

values are always kept together, but different pairs of `score`

and `bty_avg`

values may be resampled multiple times. The resulting confidence interval will denote a range of plausible values for the unknown population slope \(\beta_1\) quantifying the relationship between teaching and “beauty” scores for *all* professors at UT Austin.

Let’s first construct the bootstrap distribution for the fitted slope \(b_1\):

```
bootstrap_distn_slope <- evals_ch5 %>%
specify(formula = score ~ bty_avg) %>%
generate(reps = 1000, type = "bootstrap") %>%
calculate(stat = "slope")
bootstrap_distn_slope
```

```
# A tibble: 1,000 x 2
replicate stat
<int> <dbl>
1 1 0.0651055
2 2 0.0382313
3 3 0.108056
4 4 0.0666601
5 5 0.0715932
6 6 0.0854565
7 7 0.0624868
8 8 0.0412859
9 9 0.0796269
10 10 0.0761299
# … with 990 more rows
```

Observe how we have 1000 values of the bootstrapped slope \(b_1\) in the `stat`

column. Let’s visualize the 1000 bootstrapped values in Figure 10.8.

Observe how the bootstrap distribution is roughly bell-shaped. Recall from Subsection 8.7.1 that the shape of the bootstrap distribution of \(b_1\) closely approximates the shape of the sampling distribution of \(b_1\).

#### Percentile-method

First, let’s compute the 95% confidence interval for \(\beta_1\) using the percentile method. We’ll do so by identifying the 2.5th and 97.5th percentiles which include the middle 95% of values. Recall that this method does not require the bootstrap distribution to be normally shaped.

```
percentile_ci <- bootstrap_distn_slope %>%
get_confidence_interval(type = "percentile", level = 0.95)
percentile_ci
```

```
# A tibble: 1 x 2
`2.5%` `97.5%`
<dbl> <dbl>
1 0.0323411 0.0990027
```

The resulting percentile-based 95% confidence interval for \(\beta_1\) of (0.032, 0.099) is similar to the confidence interval in the regression Table 10.1 of (0.035, 0.099).

#### Standard error method

Since the bootstrap distribution in Figure 10.8 appears to be roughly bell-shaped, we can also construct a 95% confidence interval for \(\beta_1\) using the standard error method.

In order to do this, we need to first compute the fitted slope \(b_1\), which will act as the center of our standard error-based confidence interval. While we saw in the regression table in Table 10.1 that this was \(b_1\) = 0.067, we can also use the `infer`

pipeline with the `generate()`

step removed to calculate it:

```
# A tibble: 1 x 1
stat
<dbl>
1 0.0666370
```

We then use the `get_ci()`

function with `level = 0.95`

to compute the 95% confidence interval for \(\beta_1\). Note that setting the `point_estimate`

argument to the `observed_slope`

of 0.067 sets the center of the confidence interval.

```
se_ci <- bootstrap_distn_slope %>%
get_ci(level = 0.95, type = "se", point_estimate = observed_slope)
se_ci
```

```
# A tibble: 1 x 2
lower upper
<dbl> <dbl>
1 0.0333767 0.0998974
```

The resulting standard error-based 95% confidence interval for \(\beta_1\) of \((0.033, 0.1)\) is slightly different than the confidence interval in the regression Table 10.1 of \((0.035, 0.099)\).

#### Comparing all three

Let’s compare all three confidence intervals in Figure 10.9, where the percentile-based confidence interval is marked with solid lines, the standard error based confidence interval is marked with dashed lines, and the theory-based confidence interval (0.035, 0.099) from the regression table in Table 10.1 is marked with dotted lines.

```
visualize(bootstrap_distn_slope) +
shade_confidence_interval(endpoints = percentile_ci, fill = NULL,
linetype = "solid", color = "grey90") +
shade_confidence_interval(endpoints = se_ci, fill = NULL,
linetype = "dashed", color = "grey60") +
shade_confidence_interval(endpoints = c(0.035, 0.099), fill = NULL,
linetype = "dotted", color = "black")
```

Observe that all three are quite similar! Furthermore, none of the three confidence intervals for \(\beta_1\) contain 0 and are entirely located above 0. This is suggesting that there is in fact a meaningful positive relationship between teaching and “beauty” scores.

### 10.4.2 Hypothesis test for slope

Let’s now conduct a hypothesis test of \(H_0: \beta_1 = 0\) vs. \(H_A: \beta_1 \neq 0\). We will use the `infer`

package, which follows the hypothesis testing paradigm in the “There is only one test” diagram in Figure 9.14.

Let’s first think about what it means for \(\beta_1\) to be zero as assumed in the null hypothesis \(H_0\). Recall we said if \(\beta_1 = 0\), then this is saying there is no relationship between the teaching and “beauty” scores. Thus assuming this particular null hypothesis \(H_0\) means that in our “hypothesized universe” there is no relationship between `score`

and `bty_avg`

. We can therefore shuffle/permute the `bty_avg`

variable to no consequence.

We construct the null distribution of the fitted slope \(b_1\) by performing the steps that follow. Recall from Section 9.2 on terminology, notation, and definitions related to hypothesis testing where we defined the *null distribution*: the sampling distribution of our test statistic \(b_1\) assuming the null hypothesis \(H_0\) is true.

`specify()`

the variables of interest in`evals_ch5`

with the formula:`score ~ bty_avg`

.`hypothesize()`

the null hypothesis of`independence`

. Recall from Section 9.3 that this is an additional step that needs to be added for hypothesis testing.`generate()`

replicates by permuting/shuffling values from the original sample of 463 courses. We generate`reps = 1000`

replicates using`type = "permute"`

here.`calculate()`

the test statistic of interest: the fitted`slope`

\(b_1\).

In this case, we `permute`

the values of `bty_avg`

across the values of `score`

1000 times. We can do this shuffling/permuting since we assumed a “hypothesized universe” of no relationship between these two variables. Then we `calculate`

the `"slope"`

coefficient for each of these 1000 `generate`

d samples.

```
null_distn_slope <- evals %>%
specify(score ~ bty_avg) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "slope")
```

Observe the resulting null distribution for the fitted slope \(b_1\) in Figure 10.10.

Notice how it is centered at \(b_1\) = 0. This is because in our hypothesized universe, there is no relationship between `score`

and `bty_avg`

and so \(\beta_1 = 0\). Thus, the most typical fitted slope \(b_1\) we observe across our simulations is 0. Observe, furthermore, how there is variation around this central value of 0.

Let’s visualize the \(p\)-value in the null distribution by comparing it to the observed test statistic of \(b_1\) = 0.067 in Figure 10.11. We’ll do this by adding a `shade_p_value()`

layer to the previous `visualize()`

code.

Since the observed fitted slope 0.067 falls far to the right of this null distribution and thus the shaded region doesn’t overlap it, we’ll have a \(p\)-value of 0. For completeness, however, let’s compute the numerical value of the \(p\)-value anyways using the `get_p_value()`

function. Recall that it takes the same inputs as the `shade_p_value()`

function:

```
# A tibble: 1 x 1
p_value
<dbl>
1 0
```

This matches the \(p\)-value of 0 in the regression table in Table 10.1. We therefore reject the null hypothesis \(H_0: \beta_1 = 0\) in favor of the alternative hypothesis \(H_A: \beta_1 \neq 0\). We thus have evidence that suggests there is a significant relationship between teaching and “beauty” scores for *all* instructors at UT Austin.

When the conditions for inference for regression are met and the null distribution has a bell shape, we are likely to see similar results between the simulation-based results we just demonstrated and the theory-based results shown in the regression table in Table 10.1.

*Learning check*

**(LC10.2)** Repeat the inference but this time for the correlation coefficient instead of the slope. Note the implementation of `stat = "correlation"`

in the `calculate()`

function of the `infer`

package.