```{r mlr-setup, include=FALSE, purl=FALSE}
knitr::opts_knit$set(unnamed.chunk.label = "mlr-")
```

# Linear regression in matrix form  {#chap-mlr}
\chaptermark{Linear regression}

This chapter shows how to write linear regression models in matrix form. The purpose is to get you comfortable writing multivariate linear models in different matrix forms before we start working with time series versions of these models.  Each matrix form is an equivalent model for the data, but written in different forms.  You do not need to worry which form is better or worse at this point.  Simply get comfortable writing multivariate linear models in different matrix forms.

A script with all the R code in the chapter can be downloaded  [here](./Rcode/linear-regression-models-matrix.R).  The Rmd file of this chapter can be downloaded [here](./Rmds/linear-regression-models-matrix.Rmd).

### Data and packages {-}

This chapter uses the **stats**, **MARSS** and **datasets** packages. Install those packages, if needed, and load:
```{r mlr-loadpackages, results='hide', warning=FALSE, message=FALSE}
library(stats)
library(MARSS)
library(datasets)
```

We will work with the ```stackloss``` dataset available in the **datasets** package.  The dataset consists of 21 observations on the efficiency of a plant that produces nitric acid as a function of three explanatory variables: air flow, water temperature and acid concentration.  We are going to use just the first 4 datapoints so that it is easier to write the matrices, but the concepts extend to as many datapoints as you have.
```{r mlr-stackloss.data}
data(stackloss, package="datasets")
dat = stackloss[1:4,] #subsetted first 4 rows
dat
```


## A simple regression: one explanatory variable {#sec-mlr-lr1}

We will start by regressing stack loss against air flow.  In R using the ```lm()``` function this is
```{r mlr-stackloss1, results='hide'}
#the dat data.frame is defined on the first page of the chapter
lm(stack.loss ~ Air.Flow, data=dat)
```

This fits the following model for the $i$-th measurment:
\begin{equation}
(\#eq:stacklossi)
stack.loss_i = \alpha + \beta air_i + e_i, \text{ where } e_i \sim \text{N}(0,\sigma^2) 
\end{equation}
We will write the model for all the measurements together in two different ways, Form 1 and Form 2.


## Matrix Form 1  {#sec-mlr-form1}
In this form, we have the explanatory variables in a matrix on the left of our parameter matrix:
\begin{equation}
(\#eq:stackloss-form1)
\begin{bmatrix}stack.loss_1\\stack.loss_2\\stack.loss_3\\stack.loss_4\end{bmatrix}
= 
\begin{bmatrix}1&air_1\\1&air_2\\1&air_3\\1&air_4\end{bmatrix}
\begin{bmatrix}\alpha\\ \beta\end{bmatrix}
+
\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}
\end{equation}
You should work through the matrix algebra to make sure you understand why Equation \@ref(eq:stackloss-form1) is Equation \@ref(eq:stacklossi) for all the $i$ data points together.


We can write the first line of Equation \@ref(eq:stackloss-form1) succinctly as
\begin{equation}
(\#eq:stackloss-form1-succinct)
\mathbf{y} = \mathbf{Z}\mathbf{x} + \mathbf{e}
\end{equation}
where $\mathbf{x}$ are our parameters, $\mathbf{y}$ are our response variables, and $\mathbf{Z}$ are our explanatory variables (with a 1 column for the intercept).  The ```lm()``` function uses Form 1, and we can recover the $\mathbf{Z}$ matrix for Form 1 by using the ```model.matrix()``` function on the output from a ```lm()``` call:
```{r mlr-model.matrix}
fit=lm(stack.loss ~ Air.Flow, data=dat)
Z=model.matrix(fit)
Z[1:4,]
```


### Solving for the parameters  {#sec-mlr-solveform1}

Note: You will not need to know how to solve linear matrix equations for this course.  This section just shows you what the ```lm()``` function is doing to estimate the parameters.

Notice that $\mathbf{Z}$ is not a square matrix and its inverse does not exist but the inverse of $\mathbf{Z}^\top\mathbf{Z}$ exists---if this is a solveable problem.  We can go through the following steps to solve for $\mathbf{x}$, our parameters $\alpha$ and $\beta$.

Start with $\mathbf{y} = \mathbf{Z}\mathbf{x} + \mathbf{e}$ and multiply by $\mathbf{Z}^\top$ on the left to get
\begin{equation*}
\mathbf{Z}^\top\mathbf{y} = \mathbf{Z}^\top\mathbf{Z}\mathbf{x} + \mathbf{Z}^\top\mathbf{e}
\end{equation*}
Multiply that by $(\mathbf{Z}^\top\mathbf{Z})^{-1}$ on the left to get
\begin{equation*}
(\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{y} = (\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{Z}\mathbf{x} + (\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{e}
\end{equation*}
$(\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{Z}$  equals the identity matrix, thus
\begin{equation*}
(\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{y} = \mathbf{x} + (\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{e}\\
\end{equation*}
Move $\mathbf{x}$ to the right by itself, to get
\begin{equation*}
(\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{y} - (\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{e} = \mathbf{x}
\end{equation*}


Let's assume our errors, the $\mathbf{e}$, are i.i.d. which means that
\begin{equation*}
\mathbf{e} \sim \text{MVN}\begin{pmatrix}0,
\begin{bmatrix}
\sigma^2&0&0&0\\ 0&\sigma^2&0&0\\ 0&0&\sigma^2&0\\ 0&0&0&\sigma^2
\end{bmatrix}
\end{pmatrix}
\end{equation*}
This equation means $\mathbf{e}$ is drawn from a multivariate normal distribution with a variance-covariance matrix that is diagonal with equal variances.
Under that assumption, the expected value of $(\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{e}$ is zero.  So we can solve for $\mathbf{x}$ as
\begin{equation*}
\mathbf{x} = (\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{y}
\end{equation*}

Let's try that with R and compare to what you get with ```lm()```:
```{r mlr-stackloss2}
y=matrix(dat$stack.loss, ncol=1)
Z=cbind(1,dat$Air.Flow) #or use model.matrix() to get Z
solve(t(Z)%*%Z)%*%t(Z)%*%y
coef(lm(stack.loss ~ Air.Flow, data=dat))
```
As you see, you get the same values.


### Form 1 with multiple explanatory variables  {#sec-mlr-form1m}

We can easily extend Form 1 to multiple explanatory variables.  Let's say we wanted to fit this model:
\begin{equation}
(\#eq:stacklossi-mult)
stack.loss_i = \alpha + \beta_1 air_i + \beta_2 water_i + \beta_3 acid_i + e_i 
\end{equation}
With ```lm()```, we can fit this with
```{r mlr-form1_mult_exp_var}
fit1.mult=lm(stack.loss ~ Air.Flow + Water.Temp + Acid.Conc., data=dat)
```
Written in matrix form (Form 1), this is
\begin{equation}
(\#eq:stackloss-form1-mult)
\begin{bmatrix}stack.loss_1\\stack.loss_2\\stack.loss_3\\stack.loss_4\end{bmatrix}
= 
\begin{bmatrix}1&air_1&water_1&acid_1\\1&air_2&water_2&acid_2\\1&air_3&water_3&acid_3\\1&air_4&water_4&acid_4\end{bmatrix}
\begin{bmatrix}\alpha\\ \beta_1 \\ \beta_2 \\ \beta_3\end{bmatrix}
+
\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}
\end{equation}
Now $\mathbf{Z}$ is a matrix with 4 columns and $\mathbf{x}$ is a column vector with 4 rows.  We can show the $\mathbf{Z}$ matrix again directly from our ```lm()``` fit:
```{r mlr-model.matrix.mult.exp.var}
Z=model.matrix(fit1.mult)
Z
```


We can solve for $\mathbf{x}$ just like before and compare to what we get with ```lm()```:
```{r mlr-stackloss3}
y=matrix(dat$stack.loss, ncol=1)
Z=cbind(1,dat$Air.Flow, dat$Water.Temp, dat$Acid.Conc)
#or Z=model.matrix(fit2)
solve(t(Z)%*%Z)%*%t(Z)%*%y
coef(fit1.mult)
```
Take a look at the $\mathbf{Z}$ we made in R.  It looks exactly like what is in our model written in matrix form (Equation \@ref(eq:stackloss-form1-mult)).


### When does Form 1 arise?  {#sec-mlr-form1when}

This form of writing a regression model will come up when you work with dynamic linear models (DLMs).  With DLMs, you will be fitting models of the form $\mathbf{y}_t=\mathbf{Z}_t\mathbf{x}_t+\mathbf{e}_t$.  In these models you have multiple $\mathbf{y}$ at regular time points and you allow your regression parameters, the $\mathbf{x}$, to evolve through time as a random walk.


### Matrix Form 1b: The transpose of Form 1  {#sec-mlr-form1b}
We could also write Form 1 as follows:
\begin{equation}
(\#eq:stackloss-form1b) 
\begin{split}
\begin{bmatrix}stack.loss_1&stack.loss_2&stack.loss_3 &stack.loss_4\end{bmatrix}
= \\
\begin{bmatrix}\alpha& \beta_1 & \beta_2 & \beta_3 \end{bmatrix}
\begin{bmatrix}1&1&1&1\\air_1&air_2&air_3&air_4\\wind_1&wind_2&wind_3&wind_4\\acid_1&acid_2&acid_3&acid_4\end{bmatrix}
+
\begin{bmatrix}e_1&e_2&e_3&e_4\end{bmatrix}
\end{split}
\end{equation}
This is just the transpose of Form 1.  Work through the matrix algebra to make sure you understand why Equation \@ref(eq:stackloss-form1b) is Equation \@ref(eq:stacklossi) for all the $i$ data points together and why it is equal to the transpose of Equation \@ref(eq:stackloss-form1).  You'll need the relationship $(\mathbf{A}\mathbf{B})^\top=\mathbf{B}^\top \mathbf{A}^\top$.

Let's write Equation \@ref(eq:stackloss-form1b) as $\mathbf{y} = \mathbf{D}\mathbf{d}$, where $\mathbf{D}$ contains our parameters.  Then we can solve for $\mathbf{D}$ following the steps in Section \@ref(sec-mlr-solveform1) but multiplying from the right instead of from the left.  Work through the steps to show that 
$\mathbf{d} = \mathbf{y}\mathbf{d}^\top(\mathbf{d}\mathbf{d}^\top)^{-1}$.

```{r mlr-stack.loss.form1b.comp}
y=matrix(dat$stack.loss, nrow=1)
d=rbind(1, dat$Air.Flow, dat$Water.Temp, dat$Acid.Conc)
y%*%t(d)%*%solve(d%*%t(d))
coef(fit1.mult)
```


## Matrix Form 2  {#sec-mlr-form2}

In this form, we have the explanatory variables in a matrix on the right of our parameter matrix as in Form 1b but we arrange everything a little differently:
\begin{equation}
(\#eq:stackloss-form2)
\begin{bmatrix}stack.loss_1\\stack.loss_2\\stack.loss_3\\stack.loss_4\end{bmatrix}
= 
\begin{bmatrix}
\beta&0&0&0\\
0&\beta&0&0\\
0&0&\beta&0\\
0&0&0&\beta
\end{bmatrix}
\begin{bmatrix}air_1\\air_2\\air_3\\air_4\end{bmatrix}
+
\begin{bmatrix}
\alpha\\
\alpha\\
\alpha\\
\alpha
\end{bmatrix} + 
\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}
\end{equation}
Work through the matrix algebra to make sure you understand why Equation \@ref(eq:stackloss-form2) is the same as Equation \@ref(eq:stacklossi) for all the $i$ data points together.

We will write Form 2 succinctly as
\begin{equation}
(\#eq:form2-succinct)
\mathbf{y}=\mathbf{Z}\mathbf{x}+\mathbf{a}+\mathbf{e}
\end{equation}


### Form 2 with multiple explanatory variables  {#sec-mlr-form2m}

The $\mathbf{x}$ is a column vector of the explanatory variables. If we have more explanatory variables, we add them to the column vector at the bottom.  So if we had air flow, water temperature and acid concentration as explanatory variables, $\mathbf{x}$ looks like
\begin{equation}
(\#eq:ss2-form2)
\begin{bmatrix}air_1 \\ air_2 \\ air_3 \\ air_4 \\ water_1 \\ water_2 \\ water_3 \\ water_4 \\ acid_1 \\ acid_2 \\ acid_3 \\ acid_4 \end{bmatrix}
\end{equation}
Add columns to the  $\mathbf{Z}$ matrix for each new variable.
\begin{equation}
\begin{bmatrix}
\beta_1 & 0 & 0 & 0 & \beta_2 & 0 & 0 & 0 & \beta_3 & 0 & 0 & 0\\
0 & \beta_1 & 0 & 0 & 0 & \beta_2 & 0 & 0 & 0 & \beta_3 & 0 & 0\\
0&0&\beta_1&0&0&0&\beta_2&0&0&0&\beta_3&0\\
0&0&0&\beta_1&0&0&0&\beta_2&0&0&0&\beta_3
\end{bmatrix}
\end{equation}
The number of rows of $\mathbf{Z}$ is always $n$, the number of rows of $\mathbf{y}$, because the number of rows on the left and right of the equal sign must match.  The number of columns in $\mathbf{Z}$ is determined by the size of $\mathbf{x}$.  Each explanatory variable (like air flow and wind) appears $n$ times ($air_1$, $air_2$, $\dots$, $air_n$, etc).  So if the number of explanatory variables is $k$, the number of columns in $\mathbf{Z}$ is $k \times n$. The $\mathbf{a}$ column matrix holds the intercept terms.

### When does Form 2 arise?  {#sec-mlr-form2when}
Form 2 is similar to how multivariate time series models are typically written for reading by humans (on a whiteboard or paper).  In these models, we see equations like this:
\begin{equation}
(\#eq:ss1)
\begin{bmatrix}y_1\\y_2\\y_3\\y_4\end{bmatrix}_t
= 
\begin{bmatrix}
\beta_a&\beta_b\\
\beta_a&0.1\\
\beta_b&\beta_a\\
0&\beta_a
\end{bmatrix}
\begin{bmatrix}x_1 \\ x_2 \end{bmatrix}_t
+
\begin{bmatrix}
a\\
a\\
a\\
a
\end{bmatrix} +
\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}_t
\end{equation}
In this case, $\mathbf{y}_t$ is the set of four observations at time $t$ and $\mathbf{x}_t$ is the set of two explanatory variables at time $t$. The $\mathbf{Z}$ is showing how we are modeling the effects of $x_1$ and $x_2$ on the $y$s.  Notice that the effects are not consistent across the $x$ and $y$.  This model would not be possible to fit with ```lm()``` but will be easy to fit with ```MARSS()```.  

## Groups of intercepts  {#sec-mlr-intercepts}

Let's say that the odd numbered plants are in the north and the even numbered are in the south. We want to include this as a factor in our model that affects the intercept.  Let's go back to just having air flow be our explanatory variable.  Now if the plant is in the north our model is
\begin{equation}
(\#eq:groups-intercepts1)
stack.loss_i = \alpha_n + \beta air_i + e_i, \text{ where } e_i \sim \text{N}(0,\sigma^2) 
\end{equation}
If the plant is in the south, our model is
\begin{equation}
(\#eq:groups-intercepts2)
stack.loss_i = \alpha_s + \beta air_i + e_i, \text{ where } e_i \sim \text{N}(0,\sigma^2) 
\end{equation}
We'll add north/south as a factor called `reg' (region) to our dataframe:
```{r mlr-newstackloss}
dat = cbind(dat, reg=rep(c("n","s"),4)[1:4])
dat
```
And we can easily fit this model with ```lm()```.
```{r mlr-newstackloss.lm}
fit2 = lm(stack.loss ~ -1 + Air.Flow + reg, data=dat)
coef(fit2)
```
The -1 is added to the ```lm()``` call to get rid of $\alpha$.  We just want the $\alpha_n$ and $\alpha_s$ intercepts coming from our regions.  


### North/South intercepts in Form 1  {#sec-mlr-intercepts-form1}

Written in matrix form, Form 1 for this model is
\begin{equation}
(\#eq:stackloss-form1-ns)
\begin{bmatrix}stack.loss_1\\ stack.loss_2\\ stack.loss_3\\ stack.loss_4\end{bmatrix}
= 
\begin{bmatrix}air_1&1&0\\ air_2&0&1 \\air_3&1&0\\air_4&0&1\end{bmatrix}
\begin{bmatrix}\beta \\ \alpha_n \\ \alpha_s \end{bmatrix}
+
\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}
\end{equation}
Notice that odd plants get $\alpha_n$ and even plants get $\alpha_s$.  Use ```model.matrix()``` to see that this is the $\mathbf{Z}$ matrix that ```lm()``` formed. Notice the matrix output by ```model.matrix()``` looks exactly like $\mathbf{Z}$ in Equation \@ref(eq:stackloss-form1-ns).
```{r mlr-newstackloss-lm-model-matrix}
Z=model.matrix(fit2)
Z[1:4,]
```

We can solve for the parameters using $\mathbf{x} = (\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{y}$ as we did for Form 1 before by adding on the 1s and 0s columns we see in the $\mathbf{Z}$ matrix in Equation \@ref(eq:stackloss-form1-ns).  We could build this $\mathbf{Z}$ using the following R code:
```{r mlr-stackloss-form1-Z}
Z=cbind(dat$Air.Flow,c(1,0,1,0),c(0,1,0,1))
colnames(Z)=c("beta","regn","regs")
```
Or just use ```model.matrix()```.  This will save time when models are more complex.
```{r mlr-stackloss.form1.Z.lm}
Z=model.matrix(fit2)
Z[1:4,]
```
Now we can solve for the parameters:
```{r mlr-stackloss.form1.ns.solve}
y=matrix(dat$stack.loss, ncol=1)
solve(t(Z)%*%Z)%*%t(Z)%*%y
```
Compare to the output from ```lm()``` and you will see it is the same. 
```{r mlr-stackloss.form1.lm.coef}
coef(fit2)
```


### North/South intercepts in Form 2  {#sec-mlr-intercepts-form2}
We would write this model in Form 2 as
\begin{equation}
(\#eq:stackloss-form2-ns)
\begin{bmatrix}stack.loss_1\\ stack.loss_2\\ stack.loss_3\\ stack.loss_4\end{bmatrix}
= 
\begin{bmatrix}
\beta&0&0&0\\
0&\beta&0&0\\
0&0&\beta&0\\
0&0&0&\beta
\end{bmatrix}\begin{bmatrix}air_1\\air_2\\air_3\\air_4\end{bmatrix}
+
\begin{bmatrix}
\alpha_n\\
\alpha_s\\
\alpha_n\\
\alpha_s
\end{bmatrix}
+
\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}=\mathbf{Z}\mathbf{x}+\mathbf{a}+\mathbf{e}
\end{equation}

## Groups of $\beta$'s  {#sec-mlr-betas}

Now let's say that the plants have different owners, Sue and Aneesh, and we want to have $\beta$ for the air flow effect vary by owner.  If the plant is in the north and owned by Sue, the model is
\begin{equation}
(\#eq:groups-betas1)
stack.loss_i = \alpha_n + \beta_s air_i + e_i, \text{ where } e_i \sim \text{N}(0,\sigma^2) 
\end{equation}
If it is in the south and owned by Aneesh, the model is
\begin{equation}
(\#eq:groups-betas2)
stack.loss_i = \alpha_s + \beta_a air_i + e_i, \text{ where } e_i \sim \text{N}(0,\sigma^2) 
\end{equation}
You get the idea.

Now we need to add an operator variable as a factor in our stackloss dataframe. Plants 1,3 are run by Sue and plants 2,4 are run by Aneesh.
```{r mlr-newstackloss3}
dat = cbind(dat, owner=c("s","a"))
dat
```
Since the operator names can be replicated the length of our data set,  R fills in the operator colmun by replicating our string of operator names to the right length, conveniently (or alarmingly).

We can easily fit this model with ```lm()``` using the ":" notation.
```{r mlr-newstackloss.lm3}
coef(lm(stack.loss ~ -1 + Air.Flow:owner + reg, data=dat))
```
Notice that we have 4 datapoints and are estimating 4 parameters.  We are not going to be able to estimate any more parameters than data points.  If we want to estimate any more, we'll need to use the fuller stackflow dataset (which has 21 data points).


### Owner $\beta$'s in Form 1  {#sec-mlr-betas-form1}
Written in Form 1, this model is
\begin{equation}
(\#eq:stackloss-form1-owner)
\begin{bmatrix}stack.loss_1\\ stack.loss_2\\ stack.loss_3\\ stack.loss_4\end{bmatrix}
= 
\begin{bmatrix}1&0&0&air_1\\ 0&1&air_2&0 \\ 1&0&0&air_3\\ 0&1&air_4&0\end{bmatrix}
\begin{bmatrix}\alpha_n \\ \alpha_s \\ \beta_a \\ \beta_s \end{bmatrix}
+
\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}=\mathbf{Z}\mathbf{x}+\mathbf{e}
\end{equation}
The air data have been written to the right of the 1s and 0s for north/south intercepts because that is how ```lm()``` writes this model in Form 1 and I want to duplicate that (for teaching purposes). Also the $\beta$'s are ordered to be alphabetical because ```lm()``` writes the $\mathbf{Z}$ matrix like that.

Now our model is more complicated and using ```model.matrix()``` to get our $\mathbf{Z}$ saves us a lot tedious matrix building.
```{r mlr-stackloss.form1.Z.lm.beta}
fit3=lm(stack.loss ~ -1 + Air.Flow:owner + reg, data=dat)
Z=model.matrix(fit3)
Z[1:4,]
```
Notice the matrix output by ```model.matrix()``` looks exactly like $\mathbf{Z}$ in Equation \@ref(eq:stackloss-form1-owner) (ignore the attributes info).  Now we can solve for the parameters:
```{r mlr-stackloss.form1.owner.solve}
y=matrix(dat$stack.loss, ncol=1)
solve(t(Z)%*%Z)%*%t(Z)%*%y
```
Compare to the output from ```lm()``` and you will see it is the same. 


### Owner $\beta$'s in Form 2  {#sec-mlr-betas-form2}
To write this model in Form 2, we just add subscripts to the $\beta$'s in our Form 2 $\mathbf{Z}$ matrix:
\begin{equation}
(\#eq:stackloss-form2-owners)
\begin{bmatrix}stack.loss_1\\ stack.loss_2\\ stack.loss_3\\ stack.loss_4\end{bmatrix}
= 
\begin{bmatrix}
\beta_s&0&0&0\\
0&\beta_a&0&0\\
0&0&\beta_s&0\\
0&0&0&\beta_a
\end{bmatrix}\begin{bmatrix}air_1\\air_2\\air_3\\air_4\end{bmatrix}
+
\begin{bmatrix}
\alpha_n\\
\alpha_s\\
\alpha_n\\
\alpha_s
\end{bmatrix}
+
\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}=\mathbf{Z}\mathbf{x}+\mathbf{a}+\mathbf{e}
\end{equation}


## Seasonal effect as a factor  {#sec-mlr-season-factor}

Let's imagine that the data were taken consecutively in time by quarter.  We want to model the seasonal effect as an intercept change.  We will drop all other effects for now.
If the data were collected in quarter 1, the model is
\begin{equation}
(\#eq:season-i)
stack.loss_i = \alpha_1 + e_i, \text{ where } e_i \sim \text{N}(0,\sigma^2) 
\end{equation}
If collected in quarter 2, the model is
\begin{equation}
(\#eq:season-i2)
stack.loss_i = \alpha_2 + e_i, \text{ where } e_i \sim \text{N}(0,\sigma^2) 
\end{equation}
etc.

We add a column to our dataframe to account for season:
```{r mlr-newstackloss.qtr}
dat = cbind(dat, qtr=paste(rep("qtr",4),1:4,sep=""))
dat
```
And we can easily fit this model with ```lm()```.
```{r mlr-newstackloss.lm.qtr}
coef(lm(stack.loss ~ -1 + qtr, data=dat))
```
The -1 is added to the ```lm()``` call to get rid of $\alpha$.  We just want the $\alpha_1$, $\alpha_2$, etc. intercepts coming from our quarters.  

For comparison look at 
```{r mlr-newstackloss.lm.qtr2}
coef(lm(stack.loss ~ qtr, data=dat))
```
Why does it look like that when -1 is missing from the ```lm()``` call?  Where did the intercept for quarter 1 go and why are the other intercepts so much smaller?


### Seasonal intercepts written in Form 1  {#sec-mlr-season-form1}
Remembering that ```lm()``` puts models in Form 1, look at the $\mathbf{Z}$ matrix for Form 1:
```{r mlr-stackloss.form1.Z.lm.qtr1}
fit4=lm(stack.loss ~ -1 + qtr, data=dat)
Z=model.matrix(fit4)
Z[1:4,]
```

Written in Form 1, this model is
\begin{equation}
(\#eq:stackloss-form1-season)
\begin{bmatrix}stack.loss_1\\ stack.loss_2\\ stack.loss_3\\ stack.loss_4\end{bmatrix}
= 
\begin{bmatrix}1&0&0&0\\ 0&1&0&0 \\ 0&0&1&0\\ 0&0&0&1\end{bmatrix}
\begin{bmatrix}\alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{bmatrix}
+
\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}=\mathbf{Z}\mathbf{x}+\mathbf{e}
\end{equation}


Compare to the model that ```lm()``` is using when the intercept included.  What does this model look like written in matrix form?
```{r mlr-stackloss.form1.Z.lm.qtr2}
fit5=lm(stack.loss ~ qtr, data=dat)
Z=model.matrix(fit5)
Z[1:4,]
```


### Seasonal intercepts written in Form 2  {#sec-mlr-season-form2}
We do not need to add 1s and 0s to our $\mathbf{Z}$ matrix in Form 2; we just add subscripts to our intercepts matrix like we did when we had north-south intercepts.  In this model, we do not have any explanatory variables so $\mathbf{Z}\mathbf{x}$ does not appear.
\begin{equation}
(\#eq:stackloss-form2-season)
\begin{bmatrix}stack.loss_1\\ stack.loss_2\\ stack.loss_3\\ stack.loss_4\end{bmatrix}
= 
\begin{bmatrix}
\alpha_1\\
\alpha_2\\
\alpha_3\\
\alpha_4
\end{bmatrix}
+
\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}=\mathbf{a}+\mathbf{e}
\end{equation}


## Seasonal effect plus other explanatory variables*  {#sec-mlr-season-w-var}
With our four data points, we are limited to estimating four parameters.  Let's use the full 21 data points so we can estimate some more complex models.  We'll add an owner variable and a quarter variable to the stackloss dataset.

```{r mlr-set-up-big-dataset}
data(stackloss, package="datasets")
fulldat=stackloss
n=nrow(fulldat)
fulldat=cbind(fulldat, 
          owner=rep(c("sue","aneesh","joe"),n)[1:n], 
          qtr=paste("qtr",rep(1:4,n)[1:n],sep=""),
          reg=rep(c("n","s"),n)[1:n])
```

Let's fit a model where there is only an effect of air flow, but that effect varies by owner and by quarter.  We also want a different intercept for each quarter.  So if datapoint $i$ is from quarter $j$ on a plant owned by owner $k$, the model is
\begin{equation}
(\#eq:stackloss-mult-beta)
stack.loss_i = \alpha_j + \beta_{j,k} air_i + e_i
\end{equation}
So there there are $4 \times 3$ $\beta$'s (4 quarters and 3 owners) and 4 $\alpha$'s (4 quarters).

With ```lm()```, we fit the model as:
```{r mlr-complex.lm}
fit7 = lm(stack.loss ~ -1 + qtr + Air.Flow:qtr:owner, data=fulldat)
```

Take a look at $\mathbf{Z}$ for Form 1 using ```model.matrix(Z)```.  It's not shown since it is large:
```{r mlr-complex.lm.Z, results='hide'}
model.matrix(fit7)
```
The $\mathbf{x}$ will be
\begin{equation}
\begin{bmatrix}\alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \\ \beta_{1,a} \\ \beta_{2,a} \\ \beta_{3,a} \\ \dots \end{bmatrix}
\end{equation}

Take a look at the model matrix that ```lm()``` is using and make sure you understand how $\mathbf{Z}\mathbf{x}$ produces Equation \@ref(eq:stackloss-mult-beta).
```{r mlr-complex.form1.solve, results='hide'}
Z=model.matrix(fit7)
```

## Models with confounded parameters*  {#sec-mlr-confound}

Try adding region as another factor in your model along with quarter and fit with ```lm()```:
```{r mlr-confounded}
coef(lm(stack.loss ~ -1 + Air.Flow + reg + qtr, data=fulldat))
```
The estimate for quarter 1 is gone (actually it was set to 0) and the estimate for quarter 4 is NA.  Look at the $\mathbf{Z}$ matrix for Form 1 and see if you can figure out the problem.  Try also writing out the model for the 1st plant and you'll see what part of the problem is and why the estimate for quarter 1 is fixed at 0.  
```{r mlr-confounded2}
fit=lm(stack.loss ~ -1 + Air.Flow + reg + qtr, data=fulldat)
Z=model.matrix(fit)
```
But why is the estimate for quarter 4 equal to NA?  What if the ordering of north and south regions was different, say 1 through 4 north, 5 through 8 south, 9 through 12 north, etc?
```{r mlr-confounded23}
fulldat2=fulldat
fulldat2$reg2 = rep(c("n","n","n","n","s","s","s","s"),3)[1:21]
fit=lm(stack.loss ~ Air.Flow + reg2 + qtr, data=fulldat2)
coef(fit)
```
Now an estimate for quarter 4 appears.

The problem is two-fold.  First by having both region and quarter intercepts, we created models where 2 intercepts appear for one $i$ model and we cannot estimate both.  ```lm()``` helps us out by setting one of the factor effects to 0.  It will chose the first alphabetically. But as we saw with the model where odd numbered plants were north and even numbered were south, we can still have a situation where one of the intercepts is non-identifiable.  ```lm()``` helps us out by alerting us to the problem by setting one to NA.  

Once you start developing your own models, you will need to make sure that all your parameters are identifiable.  If they are not, your code will simply `chase its tail'.  The code will generally take forever to converge or if you did not try different starting conditions, it may look like it converged but actually the estimates for the confounded parameters are meaningless.   So you will need to think carefully about the model you are fitting and consider if there are multiple parameters measuring the same thing (for example 2 intercept parameters).

## Solving for the parameters for Form 2*  {#sec-mlr-solveform2}
Solving for the parameters when the model is written in Form 2 is not straight-forward. We could re-write the model in Form 1, or another approach is to use Kronecker products and permutation matrices.

To solve for $\alpha$ and $\beta$, we need our parameters in a column matrix like so $\left[ \begin{smallmatrix}\alpha\\\beta\end{smallmatrix} \right]$.  We start by moving the intercept matrix, $\mathbf{a}$ into $\mathbf{Z}$.
\begin{equation}
(\#eq:stackloss-form2-solve)
\begin{bmatrix}stack.loss_1\\stack.loss_2\\stack.loss_3\\stack.loss_4\end{bmatrix}
= 
\begin{bmatrix}
\alpha&\beta&0&0&0\\
\alpha&0&\beta&0&0\\
\alpha&0&0&\beta&0\\
\alpha&0&0&0&\beta
\end{bmatrix}
\begin{bmatrix}1\\air_1\\air_2\\air_3\\air_4\end{bmatrix}
+
\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}
= \mathbf{Z}\mathbf{x} + \mathbf{e}.
\end{equation}


Then we rewrite $\mathbf{Z}\mathbf{x}$ in Equation \@ref(eq:stackloss-form2-solve) in `vec' form:  if $\mathbf{Z}$ is a $n \times m$ matrix and  $\mathbf{x}$ is a matrix with 1 column and $m$ rows, then $\mathbf{Z}\mathbf{x} = (\mathbf{x}^\top \otimes \mathbf{I}_n)\,\text{vec}(\mathbf{Z})$. The symbol $\otimes$ means Kronecker product and just ignore it since you'll never see it again in our course (or google 'kronecker product' if you are curious).   The "vec" of a matrix is that matrix rearranged as a single column:
\begin{equation*}
\,\text{vec} \begin{bmatrix}
1&2\\
3&4
\end{bmatrix} = \begin{bmatrix}
1\\3\\2\\4
\end{bmatrix} 
\end{equation*}
Notice how you just take each column one by one and stack them under each other.  In R, the vec is 
```{r mlr-vec}
A=matrix(1:6,nrow=2,byrow=TRUE)
vecA = matrix(A,ncol=1)
```
$\mathbf{I}_n$ is a $n \times n$ identity matrix, a diagonal matrix with all 0s on the off-diagonals and all 1s on the diagonal.  In R, this is simply ```diag(n)```.

To show how we solve for $\alpha$ and $\beta$, let's use an example with only 3 data points so Equation \@ref(eq:stackloss-form2-solve) becomes:
\begin{equation}
(\#eq:stackloss-form2-small)
\begin{bmatrix}stack.loss_1\\stack.loss_2\\stack.loss_3\end{bmatrix}
= 
\begin{bmatrix}
\alpha&\beta&0&0\\
\alpha&0&\beta&0\\
\alpha&0&0&\beta
\end{bmatrix}
\begin{bmatrix}1\\air_1\\air_2\\air_3\end{bmatrix}
+
\begin{bmatrix}e_1\\e_2\\e_3\end{bmatrix}
\end{equation}
Using $\mathbf{Z}\mathbf{x} = (\mathbf{x}^\top \otimes \mathbf{I}_n)\,\text{vec}(\mathbf{Z})$, this means
\begin{equation}
\begin{bmatrix}
\alpha&\beta&0&0\\
\alpha&0&\beta&0\\
\alpha&0&0&\beta
\end{bmatrix}
\begin{bmatrix}1\\air_1\\air_2\\air_3\end{bmatrix}
=\big(\begin{bmatrix}1&air_1&air_2& air_3\end{bmatrix} \otimes \begin{bmatrix}1&0&0\\ 0&1&0 \\ 0&0&1 \end{bmatrix} \bigr)
\begin{bmatrix}
\alpha\\
\alpha\\
\alpha\\
\beta\\
0\\
0\\
0\\
\beta\\
0\\
0\\
0\\
\beta
\end{bmatrix}
\end{equation}
We need to rewrite the $\,\text{vec}(\mathbf{Z})$ as a `permutation' matrix times $\left[ \begin{smallmatrix}\alpha\\\beta\end{smallmatrix} \right]$:
\begin{equation}
\begin{bmatrix}
\alpha\\
\alpha\\
\alpha\\
\beta\\
0\\
0\\
0\\
\beta\\
0\\
0\\
0\\
\beta
\end{bmatrix}
=
\begin{bmatrix}
1&0\\
1&0\\
1&0\\
0&1\\
0&0\\
0&0\\
0&0\\
0&1\\
0&0\\
0&0\\
0&0\\
0&1\\
\end{bmatrix}
\begin{bmatrix}
\alpha\\
\beta
\end{bmatrix} = \mathbf{P}\mathbf{p}
\end{equation}
where $\mathbf{P}$ is the permutation matrix and $\mathbf{p}=\left[ \begin{smallmatrix}\alpha\\\beta\end{smallmatrix} \right]$.
Thus,
\begin{equation}
(\#eq:stackloss-form2-small2)
\mathbf{y}=\mathbf{Z}\mathbf{x}+\mathbf{e} = (\mathbf{x}^\top \otimes \mathbf{I}_n)\mathbf{P}\begin{bmatrix}\alpha\\ \beta\end{bmatrix} = \mathbf{M}\mathbf{p} + \mathbf{e}
\end{equation}
where $\mathbf{M}=(\mathbf{x}^\top \otimes \mathbf{I}_n)\mathbf{P}$.
We can solve for $\mathbf{p}$, the parameters, using 
$$(\mathbf{M}^\top\mathbf{M})^{-1}\mathbf{M}^\top\mathbf{y}$$
as before.  


#### Code to solve for parameters in Form 2  {#sec-mlr-solveform2code}
In the homework, you will use the R code in this section to solve for the parameters in Form 2.  

```{r mlr-stackloss.form2.solve, tidy=FALSE}
#make your y and x matrices
y=matrix(dat$stack.loss, ncol=1)
x=matrix(c(1,dat$Air.Flow),ncol=1)
#make the Z matrix
n=nrow(dat) #number of rows in our data file
k=1
#Z has n rows and 1 col for intercept, and n cols for the n air data points
#a list matrix allows us to combine "characters" and numbers
Z=matrix(list(0),n,k*n+1) 
Z[,1]="alpha"
diag(Z[1:n,1+1:n])="beta" 
#this function creates that permutation matrix for you
P=MARSS:::convert.model.mat(Z)$free[,,1]
M=kronecker(t(x),diag(n))%*%P
solve(t(M)%*%M)%*%t(M)%*%y
coef(lm(dat$stack.loss ~ dat$Air.Flow))
```
Go through this code line by line at the R command line.  Look at ```Z```. It is a list matrix that allows you to combine numbers (the 0s) with character string (names of parameters).  Notice that ```class(Z[1,3])="numeric"``` while ```class(Z[1,2])="character"```.  This is important.  ```0``` in R is a number while ```"0"``` would be a character (the name of a parameter).
Look at the permutation matrix ```P```.  Try ```MARSS:::convert.model.mat(Z)$free``` and see that it returns a 3D matrix, which is why the ```[,,1]``` appears (to get us a 2D matrix).  To use more data points, you can redefine 
```dat``` to say ```dat=stackloss``` to use all 21 data points.

Here's another example. Rewrite the model with multiple intercepts (Equation \@ref(eq:stackloss-form2-ns) ) as
\begin{equation}
(\#eq:stackloss-form2-ns-compact)
\begin{bmatrix}stack.loss_1\\ stack.loss_2\\ stack.loss_3\\ stack.loss_4\end{bmatrix}
= 
\begin{bmatrix}
\alpha_n&\beta&0&0&0\\
\alpha_s&0&\beta&0&0\\
\alpha_n&0&0&\beta&0\\
\alpha_s&0&0&0&\beta
\end{bmatrix}\begin{bmatrix}1\\air_1\\air_2\\air_3\\air_4\end{bmatrix}
+
\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}=\mathbf{Z}\mathbf{x}+\mathbf{a}+\mathbf{e}
\end{equation}
To estimate the parameters, we need to be able to write a list matrix that looks like $\mathbf{Z}$ in Equation \@ref(eq:stackloss-form2-ns-compact). We can use the same code as above with $\mathbf{Z}$ changed to look like that in Equation \@ref(eq:stackloss-form2-ns-compact).
```{r mlr-stackloss.form2.ns}
y=matrix(dat$stack.loss, ncol=1)
x=matrix(c(1,dat$Air.Flow),ncol=1)
n=nrow(dat)
k=1
#list matrix allows us to combine numbers and character strings
Z=matrix(list(0),n,k*n+1)
Z[seq(1,n,2),1]="alphanorth"
Z[seq(2,n,2),1]="alphasouth"
diag(Z[1:n,1+1:n])="beta"
P=MARSS:::convert.model.mat(Z)$free[,,1]
M=kronecker(t(x),diag(n))%*%P
solve(t(M)%*%M)%*%t(M)%*%y
```

Similarly to estimate the parameters for Equation \@ref(eq:stackloss-form2-owners), we change the $\beta$'s in our $\mathbf{Z}$ list matrix to have owner designations:
```{r mlr-stackloss.form2.owners}
Z=matrix(list(0),n,k*n+1)
Z[seq(1,n,2),1]="alphanorth"
Z[seq(2,n,2),1]="alphasouth"
diag(Z[1:n,1+1:n])=rep(c("beta.s","beta.a"),n)[1:n]
P=MARSS:::convert.model.mat(Z)$free[,,1]
M=kronecker(t(x),diag(n))%*%P
solve(t(M)%*%M)%*%t(M)%*%y
```
The parameters estimates are the same as with the model in Form 1, though $\beta$'s are given in reversed order simply due to the way ```convert.model.mat()``` is ordering the columns in Form 2's $\mathbf{Z}$.


\clearpage

## Problems  {#sec-mlr-problems}

For the homework questions, we will using part of the ```airquality``` data set in R.  Load that as
```{r mlr-homework-data, tidy=FALSE}
data(airquality, package="datasets")
#remove any rows with NAs omitted.
airquality=na.omit(airquality)
#make Month a factor (i.e., the Month number is a name rather than a number)
airquality$Month=as.factor(airquality$Month)
#add a region factor
airquality$region = rep(c("north","south"),60)[1:111]
#Only use 5 data points for the homework so you can show the matrices easily
homeworkdat = airquality[1:5,]
```


1. Using Form 1 $\mathbf{y}=\mathbf{Z}\mathbf{x}+\mathbf{e}$, write out the model, showing the $\mathbf{Z}$ and $\mathbf{x}$ matrices, being fit by this command
    ```{r mlr-hw1}
    fit=lm(Ozone ~ Wind + Temp, data=homeworkdat)
    ```

2. For the above model, write out the following R code.
    a. Create the $\mathbf{y}$ and $\mathbf{Z}$ matrices in R.
    b. Solve for $\mathbf{x}$ (the parameters).  Show that they match what you get from the first ```lm()``` call.

3. Add -1 to your ```lm()``` call in question 1: 
    ```{r mlr-hw1b}
    fit=lm(Ozone ~ -1 + Wind + Temp, data=homeworkdat)
    ```
    a. What changes in  your model? 
    b. Write out the in Form 1 as an equation. Show the new $\mathbf{Z}$ and $\mathbf{x}$ matrices. 
    c. Solve for the parameters ($\mathbf{x}$) and show they match what is returned by ```lm()```.

4.  For the model for question 1,
    a. Write in Form 2 as an equation.
    b. Adapt the code from subsection \@ref(sec-mlr-solveform2code) and construct new `Z`, `y` and `x` in R code.
    c. Solve for the parameters using the code from subsection \@ref(sec-mlr-solveform2code).

5. A model of the ozone data with only a region (north/south) effect can be written:
    ```{r mlr-hw4}
    fit=lm(Ozone ~ -1 + region, data=homeworkdat)
    ```

    a. Write this model in Form 1 as an equation.
    b. Solve for the parameter values and show that they match what you get from the ```lm()``` call.

6. Using the same model from question 5,

   a. Write the model in Form 2 as an equation. 
   b. Write out the `Z` and `x` in R code.
   c. Solve for the parameter values and show that they match what you get from the ```lm()``` call.  To do this, you adapt the code from subsection \@ref(sec-mlr-solveform2code).

7. Write the model below in Form 2 as an equation. Show the $\mathbf{Z}$, $\mathbf{y}$ and $\mathbf{x}$ matrices.
    ```{r mlr-hw5}
    fit=lm(Ozone ~ Temp:region, data=homeworkdat)
    ```

8.  Using the airquality dataset with 111 data points

    a. Write the model below in Form 2.
    ```{r mlr-hw8}
    fit=lm(Ozone ~ -1 + Temp:region + Month, data=airquality)
    ```
    b. Solve for the parameters by adapting code from subsection \@ref(sec-mlr-solveform2code).

