Multivariate Models

class: center, middle, inverse, title-slide

# Multivariate Models
### Yue Jiang
### STA 440 / Spring 2025

---

### Coffee

---

### Bean characteristics

In processing coffee beans, there are two methods prior to roasting that affect
the flavor profile of the coffee.

Dry process beans are immediately dried after picking. Since the skin of the
fruit is attached the entire time, it locks in the sugars, resulting in a
sweeter, less acidic taste.

Wet process beans have the skin and pulp of the fruit removed after picking, and
the seeds ferment in a wet mass and then soaked in water before the drying 
process. This results in a more acidic taste.

---

### Coffee

---

### Coffee

---

### Coffee

---

### Coffee

---

### Coffee

---

### Coffee

---

### Some predictors of aroma

```r
round(summary(lm(aroma ~ process + moisture + country, 
                 data = coffee))$coef, 3)
```

```
##                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)            7.537      0.037 201.172    0.000
## processWashed / Wet   -0.066      0.033  -1.999    0.046
## moisture              -0.293      0.260  -1.124    0.261
## countryColombia        0.205      0.047   4.386    0.000
## countryEthiopia        0.485      0.065   7.422    0.000
## countryGuatemala       0.088      0.047   1.861    0.063
## countryMexico          0.014      0.046   0.306    0.759
## countryTaiwan          0.083      0.055   1.494    0.136
## countryUS (Hawaii)    -0.005      0.051  -0.096    0.924
```

---

### Some predictors of flavor

```r
round(summary(lm(flavor ~ process + moisture + country, 
                 data = coffee))$coef, 3)
```

```
##                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)            7.586      0.039 192.466    0.000
## processWashed / Wet   -0.091      0.035  -2.604    0.009
## moisture              -0.199      0.274  -0.728    0.467
## countryColombia        0.116      0.049   2.368    0.018
## countryEthiopia        0.517      0.069   7.527    0.000
## countryGuatemala      -0.002      0.050  -0.038    0.970
## countryMexico         -0.110      0.049  -2.263    0.024
## countryTaiwan         -0.034      0.058  -0.585    0.559
## countryUS (Hawaii)    -0.052      0.054  -0.975    0.330
```

---

### Some predictors of aftertaste

```r
round(summary(lm(aftertaste ~ process + moisture + country, 
                 data = coffee))$coef, 3)
```

```
##                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)            7.448      0.040 187.882    0.000
## processWashed / Wet   -0.079      0.035  -2.258    0.024
## moisture              -0.176      0.276  -0.640    0.523
## countryColombia        0.187      0.049   3.794    0.000
## countryEthiopia        0.457      0.069   6.617    0.000
## countryGuatemala      -0.104      0.050  -2.077    0.038
## countryMexico         -0.132      0.049  -2.698    0.007
## countryTaiwan          0.037      0.059   0.637    0.524
## countryUS (Hawaii)     0.017      0.054   0.319    0.750
```

---

### Some predictors of body

```r
round(summary(lm(body ~ process + moisture + country, 
                 data = coffee))$coef, 3)
```

```
##                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)            7.576      0.032 238.572    0.000
## processWashed / Wet   -0.106      0.028  -3.792    0.000
## moisture              -0.299      0.221  -1.353    0.177
## countryColombia        0.186      0.040   4.713    0.000
## countryEthiopia        0.383      0.055   6.912    0.000
## countryGuatemala       0.017      0.040   0.433    0.665
## countryMexico         -0.068      0.039  -1.721    0.086
## countryTaiwan         -0.008      0.047  -0.164    0.869
## countryUS (Hawaii)     0.072      0.043   1.658    0.098
```

---

### Some predictors of balance

```r
round(summary(lm(balance ~ process + moisture + country, 
                 data = coffee))$coef, 3)
```

```
##                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)            7.557      0.038 197.894    0.000
## processWashed / Wet   -0.106      0.034  -3.146    0.002
## moisture              -0.462      0.265  -1.742    0.082
## countryColombia        0.278      0.048   5.833    0.000
## countryEthiopia        0.477      0.067   7.168    0.000
## countryGuatemala       0.053      0.048   1.104    0.270
## countryMexico         -0.089      0.047  -1.887    0.060
## countryTaiwan         -0.026      0.057  -0.461    0.645
## countryUS (Hawaii)     0.118      0.052   2.259    0.024
```

---

### Five separate models

```
##                      aroma flavor aftertaste   body balance
## (Intercept)          7.537  7.586      7.448  7.576   7.557
## processWashed / Wet -0.066 -0.091     -0.079 -0.106  -0.106
## moisture            -0.293 -0.199     -0.176 -0.299  -0.462
## countryColombia      0.205  0.116      0.187  0.186   0.278
## countryEthiopia      0.485  0.517      0.457  0.383   0.477
## countryGuatemala     0.088 -0.002     -0.104  0.017   0.053
## countryMexico        0.014 -0.110     -0.132 -0.068  -0.089
## countryTaiwan        0.083 -0.034      0.037 -0.008  -0.026
## countryUS (Hawaii)  -0.005 -0.052      0.017  0.072   0.118
```

```
##                      aroma flavor aftertaste   body balance
## (Intercept)          0.000  0.000      0.000  0.000   0.000
## processWashed / Wet  0.046  0.009      0.024  0.000   0.002
## moisture             0.261  0.467      0.523  0.177   0.082
## countryColombia      0.000  0.018      0.000  0.000   0.000
## countryEthiopia      0.000  0.000      0.000  0.000   0.000
## countryGuatemala     0.063  0.970      0.038  0.665   0.270
## countryMexico        0.759  0.024      0.007  0.086   0.060
## countryTaiwan        0.136  0.559      0.524  0.869   0.645
## countryUS (Hawaii)   0.924  0.330      0.750  0.098   0.024
```

---

### Five separate models

.question[
What do we think of these estimates?
]

---

### Modeling all at once

.question[
What if you wanted to model the outcomes *jointly* in a .vocab[multivariate] model?
]

Suppose we have the model

`\begin{align*}
\mathbf{Y}_{n \times k} = \mathbf{X}_{n \times p}\mathbf{B}_{p \times k} + \mathbf{E}_{n \times k}
\end{align*}`

where `$\mathbf{Y}$` is now a matrix of the `$k$` responses for `$n$` observations, `$\mathbf{B}$` now represents a matrix of the `$p$` parameters to be estimated for each of the `$k$` models, and `$\mathbf{E}$` is the matrix of errors for each of the `$n$` observations in the `$k$` models.

---

### Modeling all at once

`\begin{align*}
\mathbf{Y}_{n \times k} = \mathbf{X}_{n \times p}\mathbf{B}_{p \times k} + \mathbf{E}_{n \times k}
\end{align*}`

.question[
- Suppose we had some estimate `$\widehat{\mathbf{B}}$`. What are its dimensions?
- Consider the notion of a "residual" `$\mathbf{Y} - \mathbf{X}\widehat{\mathbf{B}}$`. What are its dimensions?
- What might a "least squares" criterion be defined as regarding these residuals?
- How might we find such a "least squares" estimate?
]

---

### Modeling all at once

As it turns out, the minimizer `$\widehat{\mathbf{B}} \in \mathcal{M}_{p \times k}$` of the residual sum of the squares is

`\begin{align*}
\widehat{\mathbf{B}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{Y}.
\end{align*}`

.question[
How does this compare to the individual least squares estimates `$\widehat{\mathbf{\beta}}$` from each of the `$k$` (in the case of our coffee example, 5) models?
]

---

### Modeling all at once

`\begin{align*}
\widehat{\mathbf{B}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{Y}.
\end{align*}`

.question[
What happens if `$p > n$`? What happens if `$k > n$`?
]

---

### Modeling all at once

```r
mv <- lm(cbind(aroma, flavor, aftertaste,
               body, balance) ~ 
           process + moisture + country, data = coffee)
round(mv$coefficients, 3)
```

---

### Five separate models

.question[
What do we think of these p-values?
]

---

### Assumptions for the multivariate model

1. Independence
2. Linearity: `$E(\mathbf{Y} | \mathbf{X}) = \mathbf{X}\mathbf{B}$`
3. Constant Variance: `$Var(\mathbf{E} | \mathbf{X}) = \boldsymbol{\Sigma}$`
4. Normality: `$\mathbf{E} \sim N_k(\mathbf{0}, \boldsymbol{\Sigma})$`

.question[
- What is the distribution of `$\mathbf{Y} | \mathbf{X}$`?
- What is the distribution of each `$\widehat{\beta}_k$`?
- What is the joint distribution of `$\widehat{\beta}_j$` and `$\widehat{\beta}_k$`?
]

---

### Hypothesis tests for the multivariate model

.question[
How might you perform a hypothesis test on the *joint* significance of a predictor for all `$k$` responses of interest?

How about confidence intervals for coefficients in individual models?
]

(don't actually answer this unless you *reeeeally* like linear algebra)

---

### The covariance matrix of coefficients

```r
vcov(mv), 3
```

(this thing is gigantic, so I'm not showing it here)

Remember that the estimated coefficients from the models *covary*, and so we
can't use the individual p-values for each model to see if a predictor is
*jointly* significant across the outcomes.

To formally test this, we need to take into consideration covariances between
coefficients in all `$k$` models.

---

### Joint significance testing

```r
library(car)
Manova(mv) # watch the capitalization! 
```

```
## 
## Type II MANOVA Tests: Pillai test statistic
##          Df test stat approx F num Df den Df    Pr(>F)
## process   1   0.02199   3.2834      5    730  0.006118
## moisture  1   0.00612   0.8986      5    730  0.481415
## country   6   0.39254  10.4225     30   3670 < 2.2e-16
```

---

### Confidence intervals for parameters

(in output below, we have response/predictor combos, *not* interaction terms of some sort)

```r
round(confint(mv), 3)
```

```
##                                 2.5 % 97.5 %
## aroma:(Intercept)               7.463  7.610
## aroma:processWashed / Wet      -0.131 -0.001
## aroma:moisture                 -0.804  0.218
## aroma:countryColombia           0.113  0.296
## aroma:countryEthiopia           0.357  0.613
## aroma:countryGuatemala         -0.005  0.182
## aroma:countryMexico            -0.077  0.105
## aroma:countryTaiwan            -0.026  0.192
## aroma:countryUS (Hawaii)       -0.105  0.096
## flavor:(Intercept)              7.509  7.664
## flavor:processWashed / Wet     -0.159 -0.022
## flavor:moisture                -0.737  0.338
## flavor:countryColombia          0.020  0.213
## flavor:countryEthiopia          0.383  0.652
## flavor:countryGuatemala        -0.100  0.096
## flavor:countryMexico           -0.206 -0.015
## flavor:countryTaiwan           -0.149  0.080
## flavor:countryUS (Hawaii)      -0.158  0.053
## aftertaste:(Intercept)          7.370  7.526
## aftertaste:processWashed / Wet -0.148 -0.010
## aftertaste:moisture            -0.717  0.365
## aftertaste:countryColombia      0.090  0.284
## aftertaste:countryEthiopia      0.322  0.593
## aftertaste:countryGuatemala    -0.203 -0.006
## aftertaste:countryMexico       -0.228 -0.036
## aftertaste:countryTaiwan       -0.078  0.153
## aftertaste:countryUS (Hawaii)  -0.089  0.124
## body:(Intercept)                7.513  7.638
## body:processWashed / Wet       -0.161 -0.051
## body:moisture                  -0.732  0.135
## body:countryColombia            0.109  0.264
## body:countryEthiopia            0.274  0.492
## body:countryGuatemala          -0.062  0.096
## body:countryMexico             -0.145  0.010
## body:countryTaiwan             -0.100  0.085
## body:countryUS (Hawaii)        -0.013  0.157
## balance:(Intercept)             7.482  7.632
## balance:processWashed / Wet    -0.172 -0.040
## balance:moisture               -0.983  0.059
## balance:countryColombia         0.184  0.371
## balance:countryEthiopia         0.347  0.608
## balance:countryGuatemala       -0.042  0.148
## balance:countryMexico          -0.182  0.004
## balance:countryTaiwan          -0.137  0.085
## balance:countryUS (Hawaii)      0.015  0.220
```

---

### What about multivariate GLMs?

or multivariate quantile regression models?

.question[
Seriously, what about them? Try Googling and seeing what considerations need
to be made!
]