Common Statistical Errors

class: center, middle, inverse, title-slide

.title[
# Common Statistical Errors
]
.author[
### Yue Jiang
]
.date[
### STA 440 / Spring 2025
]

---

### Announcements

- HW 1 (check early afternoon in class repository; due Jan. 31).
- Lab Abstract 1 (due Feb. 7).

---

### Choosing a model / choosing variables

.question[
Why do we even want to select a model? Why would we want to select variables?
Are model selection and variable selection the same thing?
]

We might want to **explain** or **describe** relationships in populations (often
by performing statistical inference on a sample from this population). In some
cases, this might give us causal or mechanistic understanding of the phenomenon
at hand; in others, we might only be able to explore associations.

We alternatively might want to **predict** values of new observations, often by
examining how well our model does using our data at hand (perhaps in terms of
performance on holdout data.)

---

### A common occurrance

```
##  height     age logfev1
##    1.52 14.1985 0.60432
##    1.68 15.9562 1.16002
##    1.42 11.1102 0.53649
##    1.42 10.8008 0.62594
##    1.47 11.1923 0.87129
##    1.63 14.6502 1.21194
##    1.30  7.3128 0.18232
##    1.66 18.0726 1.34025
##    1.63 14.3874 1.17248
##    1.24  7.1348 0.54232
##    1.22  8.2930 0.41211
##    1.35  9.9384 0.60432
##    1.27  7.4251 0.37156
##    1.63 17.0568 1.16002
##    1.63 15.4552 1.19089
##    1.71 17.7413 1.25846
##    1.43 11.0856 0.60977
##    1.57 14.5188 1.02245
##    1.49 11.5811 0.94391
##    1.63 15.7043 1.30563
##    1.60 16.1232 1.13462
##    1.26  8.7310 0.49470
##    1.68 17.2567 1.27257
##    1.51 15.6140 0.93609
##    1.33  9.2156 0.29267
##    1.45 10.6694 0.68813
##    1.63 16.9665 1.01523
##    1.69 17.6756 1.17866
##    1.62 14.3874 1.00430
##    1.57 13.2813 0.70310
##    1.55 12.6324 0.69813
```

---

### A common occurrence

``` r
m2 <- lm(logfev1 ~ height + age, data = fev)
summary(m2)
```

```
## 
## Call:
## lm(formula = logfev1 ~ height + age, data = fev)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.310946 -0.072796  0.001951  0.100668  0.239072 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.36514    0.55433  -2.463   0.0202
## height       1.09719    0.56574   1.939   0.0626
## age          0.04315    0.02472   1.746   0.0918
## 
## Residual standard error: 0.1401 on 28 degrees of freedom
## Multiple R-squared:  0.8393,	Adjusted R-squared:  0.8278 
## F-statistic: 73.11 on 2 and 28 DF,  p-value: 7.667e-12
```

.question[
Is "low collinearity" an assumption of the linear model?
]

---

### A common occurrence

---

### A common occurrance

``` r
m1 <- lm(logfev1 ~ height, data = fev)
summary(m1)
```

```
## 
## Call:
## lm(formula = logfev1 ~ height, data = fev)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27684 -0.08115  0.04251  0.09438  0.23525 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -2.2213     0.2674  -8.307 3.71e-09
## height        2.0390     0.1763  11.564 2.21e-12
## 
## Residual standard error: 0.1449 on 29 degrees of freedom
## Multiple R-squared:  0.8218,	Adjusted R-squared:  0.8157 
## F-statistic: 133.7 on 1 and 29 DF,  p-value: 2.213e-12
```

---

### A common occurrance

``` r
m2 <- lm(logfev1 ~ age, data = fev)
summary(m2)
```

```
## 
## Call:
## lm(formula = logfev1 ~ age, data = fev)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34938 -0.10021  0.01975  0.04908  0.22280 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.307985   0.105422  -2.921  0.00668
## age          0.088860   0.007791  11.405 3.08e-12
## 
## Residual standard error: 0.1466 on 29 degrees of freedom
## Multiple R-squared:  0.8177,	Adjusted R-squared:  0.8114 
## F-statistic: 130.1 on 1 and 29 DF,  p-value: 3.083e-12
```

.question[
What would you do?
]

---

### DNA methylation

---

### DNA methylation

---

### DNA methylation

---

### DNA methylation

---

### In-class activity

Load the dataset `cpg_dat`, which contains 1,000 observations of DNA methylation
M-values 20 CpG sites labeled `X1` through `X20` (because who wants to use
cg05575921 as a variable name?). A standardized measure of lung
capacity (FEV1) is presented in the outcome `y`, represented as standard
deviations above/below the age-adjusted population mean.

Using your favorite 
model/variable selection techniques for linear regression, come up with the 
"best" model (you choose what that means!).

The file is located at

`https://github.com/sta440-sp25/class-files`

---

### In-class activity

.question[
What was the form of the final model you chose? Is there sufficient evidence 
that at least one of these `$\beta$` terms is non-zero? That is, test the 
following hypotheses:

- `$H_0:\beta_1 = \beta_2 = \cdots = \beta_{20} = 0$`
- `$H_1:$` At least one `$\beta_k \neq 0$` for `$k = 1, 2, \cdots, 20$`.

What hypothesis test did you use? What might you conclude with these data?

Check on Gradescope for an additional activity regarding a *specific* CpG site
of interest. We'll discuss this together.
]

---

### Explanation vs. prediction

Is there a relationship between standardized age-adjusted lung capacity and DNA methylation at a specific site (while controlling for methylation at other CpG sites)?

Given a patient's DNA methylation at these 20 sites, can we predict their lung capacity? Can we come up with a model 
that does the best job of doing so (in some sense)?

We might use similar (or even the same) tools to both cases, but these are
**distinct** questions. With that said, we're often interested in addressing
both questions.

See also this Stack Exchange discussion [here](https://stats.stackexchange.com/questions/1194/practical-thoughts-on-explanatory-vs-predictive-modeling).

---

### Model selection

.question[
Why do we even want to select a model?
]

"Choosing the best model" is an unfortunate and perhaps misleading use of words.
In explaining our data, there might be multiple reasonable models that all do
an acceptable job explaining.

We might also consider .vocab[parsimony]. Ease of explanation and 
conceptualization might more than make up for a marginal increase in some
metric.

---

### Back to our activity...

How many people/teams did something like the following...

---

### Back to our activity...

``` r
m1 <- lm(y ~ ., data = dat); round(summary(m1)$coef[-1,], 3)
```

```
##     Estimate Std. Error t value Pr(>|t|)
## X1     0.086      0.031   2.739    0.006
## X2    -0.018      0.032  -0.555    0.579
## X3     0.000      0.032  -0.008    0.993
## X4     0.016      0.032   0.517    0.605
## X5    -0.017      0.032  -0.535    0.593
## X6    -0.031      0.032  -0.995    0.320
## X7    -0.032      0.032  -1.014    0.311
## X8    -0.006      0.030  -0.208    0.835
## X9     0.029      0.032   0.890    0.374
## X10   -0.013      0.031  -0.400    0.689
## X11   -0.064      0.032  -2.029    0.043
## X12    0.015      0.031   0.466    0.641
## X13    0.034      0.031   1.087    0.277
## X14    0.024      0.032   0.747    0.455
## X15    0.004      0.032   0.117    0.907
## X16   -0.046      0.031  -1.488    0.137
## X17   -0.025      0.031  -0.804    0.422
## X18   -0.032      0.032  -1.004    0.316
## X19    0.070      0.031   2.251    0.025
## X20    0.042      0.031   1.358    0.175
```

---

### Back to our activity...

``` r
m2 <- lm(y ~ X1 + X11 + X19, data = dat)
summary(m2)
```

```
## 
## Call:
## lm(formula = y ~ X1 + X11 + X19, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.81047 -0.65261  0.00075  0.65403  3.03093 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.01698    0.03122   0.544  0.58648
## X1           0.08442    0.03087   2.735  0.00636
## X11         -0.06220    0.03151  -1.974  0.04866
## X19          0.06947    0.03093   2.246  0.02492
## 
## Residual standard error: 0.9851 on 996 degrees of freedom
## Multiple R-squared:  0.01626,	Adjusted R-squared:  0.01329 
## F-statistic: 5.486 on 3 and 996 DF,  p-value: 0.000969
```

---

### Back to our activity...

Or this...

---

### Back to our activity...

``` r
library(MASS)
m3 <- stepAIC(m1, direction = "both", trace = F)
summary(m3)
```

```
## 
## Call:
## lm(formula = y ~ X1 + X11 + X16 + X19, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.83423 -0.65687 -0.00603  0.64656  3.06008 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.01436    0.03123   0.460  0.64570
## X1           0.08579    0.03086   2.780  0.00554
## X11         -0.06139    0.03149  -1.950  0.05149
## X16         -0.04919    0.03057  -1.609  0.10795
## X19          0.06933    0.03091   2.243  0.02510
## 
## Residual standard error: 0.9843 on 995 degrees of freedom
## Multiple R-squared:  0.01881,	Adjusted R-squared:  0.01486 
## F-statistic: 4.769 on 4 and 995 DF,  p-value: 0.0008177
```

---

### Back to our activity...

``` r
set.seed(123); n <- 1000; p <- 20; y <- rnorm(n, 0, 1)
x <- matrix(rnorm(n * p), nrow = n); dat <- data.frame(y, x)
m4 <- lm(y ~ X1 + X11 + X16 + X19, data = dat); summary(m4)
```

---

### Back to our activity...

Remember, if the null hypothesis is true, then p-values should have a uniform
distribution on (0, 1).

``` r
temp <- numeric(5000)
for(i in 1:5000){
 set.seed(i)
 y <- rnorm(n, 0, 1); x <- matrix(rnorm(n * p), nrow = n);
 dat <- data.frame(y, x)
 mod <- summary(lm(y ~ ., data = dat))
 temp[i] <- pf(mod$fstatistic[1], mod$fstatistic[2], mod$fstatistic[3], lower.tail = F)
}
```

---

### Back to our activity...

``` r
round(temp[1:25], 3)
```

```
##  [1] 0.827 0.093 0.366 0.062 0.894 0.330 0.152 0.589 0.348 0.600 0.459 0.190
## [13] 0.840 0.388 0.052 0.724 0.899 0.291 0.197 0.126 0.302 0.437 0.035 0.729
## [25] 0.405
```

``` r
hist(temp, xlab = "p-value from overall F-test", 
     ylab = "Count", main = "")
```

---

### Back to our activity...

``` r
temp <- numeric(25)
for(i in 1:25){
 set.seed(i)
 y <- rnorm(n, 0, 1); x <- matrix(rnorm(n * p), nrow = n);
 dat <- data.frame(y, x)
 mod <- summary(stepAIC(lm(y ~ ., data = dat), direction = "both", trace = F))
 temp[i] <- pf(mod$fstatistic[1], mod$fstatistic[2], mod$fstatistic[3], lower.tail = F)
}
round(temp, 3)
```

```
##  [1] 0.119 0.001 0.026 0.001 0.035 0.009 0.001 0.018 0.013 0.032 0.011 0.001
## [13] 0.075 0.003 0.001 0.055 0.120 0.002 0.001 0.003 0.006 0.002 0.000 0.012
## [25] 0.005
```

.question[
What do you notice? What's going on?
]

---

### The problem with post-selection inference

Note that we chose our final model *because it did a good job with these data*.

- In classical statistical theory, data is assumed to be generated from a known
model, with inference for parameters performed on this known model.
- In the real world, it's often the case that the generating mechanism is 
unknown, *with a model being chosen based on the data*.
- The selection process that produces a model is itself random, which affects
the sample distriutions of the post-selection parameter estimates. These may 
change dramatically after conditioning on model selection procedures.

Performing statistical inference on the same data we used to choose our model
is not appropriate - *even if* asymptotically your selection procedure correctly
chooses variables almost surely!

---

### The problem with post-selection inference

- Rhee et al. (2003) studied HIV drugs and predicted drug resistance based on
location of genetic mutations.
- Using forward selection, naive p-values suggest six significant predictors.
- However, selection-adjusted p-values only suggest two significant predictors 
(figure below from Taylor and Tibshirani, 2015, *PNAS*):

---

### The problem with post-selection inference

So what can we do?

- Split the data - perform model selection using only part of the data, and then
use that model on the remaining data.
- Methods specifically dealing with post-selection inferential methods - e.g., relying on polyhedral lemma for LASSO, etc. 
- Ignore the problem entirely.