You can download the Rmd file for this key from here.
library(datasets)
data(airquality)
#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,]Using Form 1 \(\mathbf{y}=\mathbf{Z}\mathbf{x}+\mathbf{e}\), write the model being fit by this command
fit=lm(Ozone ~ Wind + Temp, data=homeworkdat)The model is \[Ozone_i = \alpha + \beta_w Wind_i + \beta_t Temp_i + e_i.\] Form 1 for this model is: \[ \begin{bmatrix} Ozone_1 \\ Ozone_2 \\ Ozone_3 \\ Ozone_4 \\ Ozone_5 \end{bmatrix}= \begin{bmatrix} 1 & Wind_1 & Temp_1\\ 1 & Wind_2 & Temp_2 \\ 1 & Wind_3 & Temp_3 \\ 1 & Wind_4 & Temp_4 \\ 1 & Wind_5 & Temp_5 \end{bmatrix} \begin{bmatrix} \alpha \\ \beta_w \\ \beta_t \end{bmatrix} + \begin{bmatrix} e_1 \\ e_2 \\ e_3 \\ e_4 \\ e_5 \end{bmatrix} \]
Build the \(\mathbf{y}\) and \(\mathbf{Z}\) matrices for the above model in R and solve for \(\mathbf{x}\) (the parameters). Show that they match what you get from the lm call.
Here are the \(\mathbf{y}\) and \(\mathbf{Z}\). You can see they match the \(\mathbf{y}\) and \(\mathbf{Z}\) in the equation above:
y=matrix(homeworkdat$Ozone, ncol=1)
Z=cbind(1, homeworkdat$Wind, homeworkdat$Temp)
y## [,1]
## [1,] 41
## [2,] 36
## [3,] 12
## [4,] 18
## [5,] 23Z## [,1] [,2] [,3]
## [1,] 1 7.4 67
## [2,] 1 8.0 72
## [3,] 1 12.6 74
## [4,] 1 11.5 62
## [5,] 1 8.6 65Next we solve for \(\mathbf{x}\) and show it matches what we get from ‘lm’. This uses the code in Chapter 2.
solve(t(Z)%*%Z)%*%t(Z)%*%y## [,1]
## [1,] 56.6219201
## [2,] -4.9776707
## [3,] 0.2538717coef(lm(Ozone ~ Wind + Temp, data=homeworkdat))## (Intercept) Wind Temp
## 56.6219201 -4.9776707 0.2538717a. If you added -1 to your ‘lm’ call in question 1, what changes in your model?
First run the ‘lm’ call and see what changed. The intercept (\(\alpha\)) is dropped.
fit=lm(Ozone ~ -1 + Wind + Temp, data=homeworkdat)
fit##
## Call:
## lm(formula = Ozone ~ -1 + Wind + Temp, data = homeworkdat)
##
## Coefficients:
## Wind Temp
## -4.650 1.037b. Write out the new \(\mathbf{Z}\) and \(\mathbf{x}\) matrices.
The model is now \[Ozone_i = \beta_w Wind_i + \beta_t Temp_i + e_i.\] To get rid of the \(\alpha\) we drop the 1s column: \[
\begin{bmatrix} Ozone_1 \\ Ozone_2 \\ Ozone_3 \\ Ozone_4 \\ Ozone_5 \end{bmatrix}=
\begin{bmatrix} Wind_1 & Temp_1\\ Wind_2 & Temp_2 \\ Wind_3 & Temp_3 \\ Wind_4 & Temp_4 \\ Wind_5 & Temp_5 \end{bmatrix}
\begin{bmatrix} \beta_w \\ \beta_t \end{bmatrix} +
\begin{bmatrix} e_1 \\ e_2 \\ e_3 \\ e_4 \\ e_5 \end{bmatrix}
\] The \(\mathbf{Z}\) matrix is \[
\begin{bmatrix} Wind_1 & Temp_1\\ Wind_2 & Temp_2 \\ Wind_3 & Temp_3 \\ Wind_4 & Temp_4 \\ Wind_5 & Temp_5 \end{bmatrix}.
\] The \(\mathbf{x}\) matrix is \[
\begin{bmatrix} \beta_w \\ \beta_t \end{bmatrix}.
\] c. Solve and show that the result matches that returned by lm().
y=matrix(homeworkdat$Ozone, ncol=1)
Z=cbind(homeworkdat$Wind, homeworkdat$Temp)
solve(t(Z)%*%Z)%*%t(Z)%*%y## [,1]
## [1,] -4.650309
## [2,] 1.036840coef(lm(Ozone ~ -1 + Wind + Temp, data=homeworkdat))## Wind Temp
## -4.650309 1.036840a. Write the model for question 1 in Form 2.
In Form 2, \(\mathbf{y}=\mathbf{Z}\mathbf{x}+\mathbf{e}\) and the explanatory variables appear in the \(\mathbf{x}\) as a column vector (a matrix with one column). So \(\mathbf{x}\) looks like this \[\begin{bmatrix} 1 \\ Wind_1 \\ \dots \\ Wind_5 \\ Temp_1 \\ \dots \\ Temp_5 \end{bmatrix}\]
Once you get that, then you know that \(\mathbf{Z}\) is a \(5 \times (1+5+5)\) matrix. The first column of \(\mathbf{Z}\) is the \(\alpha\). The next 5 columns of \(\mathbf{Z}\) is a \(5 \times 5\) diagonal matrix with \(\beta_w\) on the diagonal, just like in Equation 1.12. For the next explanatory variable, Temperature, we tack on another \(5 \times 5\) diagonal matrix; this time with \(\beta_t\) on the diagonal. That’s all you needed to say for homework; just to show that you figured out what the form of the \(\mathbf{y}\), \(\mathbf{Z}\) and \(\mathbf{x}\) look like.
If you wanted to write it out in math form, you could show \(\mathrm{Z}\) as: \[ \begin{bmatrix} \text{column of} & 5\times5\text{ diagonal matrix}& 5\times5\text{ diagonal matrix} \\ \alpha & \text{with }\beta_w\text{ on the diagonal} & \text{with }\beta_t\text{ on the diagonal} \end{bmatrix} \] or \[ \begin{bmatrix} \alpha & \beta_w & \dots & 0 & \beta_t & \dots & 0 \\ \dots & \dots & \ddots & \dots & \dots & \ddots & \dots \\ \alpha & 0 & \dots & \beta_w & 0 & \dots & \beta_t \end{bmatrix} \]
b. Construct the \(\mathbf{Z}\), \(\mathbf{y}\) and \(\mathbf{x}\) in R code.
To do the 2nd part, you adapt the code from subsection to solve for the parameters. You will need to contruct new \(\mathbf{Z}\), \(\mathbf{y}\) and \(\mathbf{x}\) in the code.
The \(\mathbf{y}\) and \(\mathbf{x}\) are easy:
y=matrix(homeworkdat$Ozone, ncol=1)
x=matrix(c(1, homeworkdat$Wind, homeworkdat$Temp), ncol=1)We adapt the code for making \(\mathbf{Z}\) from Section 2.3.4:
#we know that Z is a 5 x (1+5+5) matrix
#n is the number of data points
n=5
#nrows = n; what about ncol?, ncol is 1 (alpha) + n (Wind) + n (Temp)
Z=matrix(list(0),n,1+n+n)
#the first column is alpha
Z[,1]="alpha"
#columns 2:112 are a diagonal matrix with betaw on the diagonal
diag(Z[,2:(n+1)])="betaw"
#columns 113:223 are a diagonal matrix with betat on the diagonal
diag(Z[,(n+2):(2*n+1)])="betat"
Z## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] "alpha" "betaw" 0 0 0 0 "betat" 0
## [2,] "alpha" 0 "betaw" 0 0 0 0 "betat"
## [3,] "alpha" 0 0 "betaw" 0 0 0 0
## [4,] "alpha" 0 0 0 "betaw" 0 0 0
## [5,] "alpha" 0 0 0 0 "betaw" 0 0
## [,9] [,10] [,11]
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] "betat" 0 0
## [4,] 0 "betat" 0
## [5,] 0 0 "betat"c. Solve for the parameters.
Now we can solve for \(\mathrm{Z}\):
require(MARSS)## Loading required package: MARSSP=MARSS:::convert.model.mat(Z)$free[,,1]
M=kronecker(t(x),diag(n))%*%P
solve(t(M)%*%M)%*%t(M)%*%y## [,1]
## alpha 56.6219201
## betaw -4.9776707
## betat 0.2538717coef(lm(Ozone ~ Wind + Temp, data=homeworkdat))## (Intercept) Wind Temp
## 56.6219201 -4.9776707 0.2538717a. Model the ozone data with only a region effect. Write this in Form 1.
First make sure you understand what model is being fit by the ‘lm’ call:
fit=lm(Ozone ~ -1 + region, data=homeworkdat)
fit##
## Call:
## lm(formula = Ozone ~ -1 + region, data = homeworkdat)
##
## Coefficients:
## regionnorth regionsouth
## 25.33 27.00The model is \[Ozone_i = \alpha_j + e_i\] where \(j\) is the region the measurement was taken in. This is an intercept (or level) only model where the intercept is determined by the region (north or south).
We want to write that model in matrix form using Form 1, \(\mathbf{y}=\mathbf{Z}\mathbf{x}+\mathbf{e}\). Eqn 2.2 shows you how to do this for Form 1, except we do not have explanatory variables besides region so we do not have a column with something like ‘air’ in it.
Matrix \(\mathbf{y}\) is as before: \[\mathbf{y} = \begin{bmatrix} Ozone_1 \\ Ozone_2 \\ \dots \\ Ozone_5\end{bmatrix}\]
Matrix \(\mathbf{x}\) are our parameters: \[\mathbf{x} = \begin{bmatrix} \alpha_n \\ \alpha_s \end{bmatrix}\]
For \(\mathbf{Z}\), we need to know that each row of \(\mathbf{x}\) is for a different data point \(i\) and tells us what region that data point is from. So \(\mathbf{Z}\) has 2 column, one for each region. If there is a 1 in column 1, it means that data point came from the north. If there is a 1 in column 2, it means that data point came from the south.
Let’s look at the regions
homeworkdat$region## [1] "north" "south" "north" "south" "north"So the first measurement is from the north, next from south, then north, … \(\mathbf{Z}\) looks like this \[\begin{bmatrix} 1&0\\ 0&1\\ 1&0\\ 0&1\\ 1&0 \end{bmatrix}\]
You can also let R show you what \(\mathbf{Z}\) is:
model.matrix(fit)## regionnorth regionsouth
## 1 1 0
## 2 0 1
## 3 1 0
## 4 0 1
## 7 1 0
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$region
## [1] "contr.treatment"The model in matrix form is \[\begin{bmatrix} Ozone_1 \\ Ozone_2 \\ \dots \\ Ozone_5\end{bmatrix}=\begin{bmatrix} 1&0\\ 0&1\\ 1&0\\ 0&1\\ 1&0 \end{bmatrix}\begin{bmatrix} \alpha_n \\ \alpha_s \end{bmatrix}+\begin{bmatrix} e_1 \\ e_2 \\ \dots \\ e_5\end{bmatrix}\]
b. Show that you can solve for the parameters and show they match what you get from lm().
Section 2.2.1 shows you how to solve for the parameters. We need to make \(\mathbf{y}\) and \(\mathbf{Z}\) in R.
y=matrix(homeworkdat$Ozone, ncol=1)We could form \(\mathbf{Z}\) like so
ndatapoints=5
Z=matrix(0,ndatapoints,2)
for(i in 1:ndatapoints) Z[i,ifelse(airquality$region[i]=="north",1,2)]=1or this
Z=cbind(
as.numeric(airquality$region=="north"),
as.numeric(airquality$region=="south")
)or just use the output from our lm call since R’s ‘lm’ is forming the \(\mathbf{Z}\) matrix too:
Z=model.matrix(fit)Now we solve for the parameters:
solve(t(Z)%*%Z)%*%t(Z)%*%y## [,1]
## regionnorth 25.33333
## regionsouth 27.00000coef(lm(Ozone ~ -1 + region, data=homeworkdat))## regionnorth regionsouth
## 25.33333 27.00000a. Write the model in question 5 in Form 2
The reason we had to use a matrix with 1s and 0s was to tell our matrix math what \(\alpha\) to use with what data point. Form 1 has the parameters appearing once in a column vector. So we need a matrix with 1s and 0s to say ‘where’ to put the \(\alpha\)s.
In Form 2, we can have the parameters just repeat in our matrix. This is how we’d write the model on the whiteboard. We just have the \(\alpha\) (intercept) column have different \(\alpha\)s in it.
The model in Form 2 is:
\[ \begin{bmatrix} Ozone_1 \\ Ozone_2 \\ Ozone_3 \\ Ozone_4 \\ Ozone_5 \end{bmatrix} = \begin{bmatrix} \alpha_n\\ \alpha_s\\ \alpha_n \\ \alpha_s \\ \alpha_n \end{bmatrix}\begin{bmatrix}1\end{bmatrix} + \begin{bmatrix} e_1\\ e_2\\ e_3 \\ e_4 \\ e_5\end{bmatrix}=\mathbf{Z}\mathbf{x}+\mathbf{e} \]
b. Solve for the parameters and show they match lm() output.
To solve for the parameters, we write \(\mathbf{y}\), \(\mathbf{Z}\), and \(\mathbf{x}\) in R and use the code in section 2.3.3.
#the number of data points
n=5
y=matrix(homeworkdat$Ozone, ncol=1)
x=matrix(1)
Z=matrix(paste("alpha",homeworkdat$region, sep="."),ncol=1)Then we solve for the parameters:
require(MARSS)
P=MARSS:::convert.model.mat(Z)$free[,,1]
M=kronecker(t(x),diag(n))%*%P
solve(t(M)%*%M)%*%t(M)%*%y## [,1]
## alpha.north 25.33333
## alpha.south 27.00000coef(lm(Ozone ~ -1 + region, data=homeworkdat))## regionnorth regionsouth
## 25.33333 27.00000Write the model below in Form 2 and show the \(\mathbf{Z}\), \(\mathbf{y}\) and \(\mathbf{x}\) parameters.
fit=lm(Ozone ~ Temp:region, data=homeworkdat)The first step is to write out what model is being fit. The model is \[ Ozone_i = \alpha + \beta_n Temp_i + e_i\] if \(i\) is from the north, and is \[ Ozone_i = \alpha + \beta_s Temp_i + e_i\] if \(i\) is from the south. So each region has the same intercept \(\alpha\) but we are including a temperature effect that is different for each region.
The model written in Form 2 is
\[\begin{bmatrix} Ozone_1 \\ Ozone_2 \\ Ozone_3 \\ Ozone_4 \\ Ozone_5 \end{bmatrix} = \begin{bmatrix} \alpha & \beta_n & 0 & 0 & 0 & 0\\ \alpha & 0 & \beta_s & 0 & 0 & 0\\ \alpha & 0 & 0 & \beta_n & 0 & 0\\ \alpha & 0 & 0 & 0 & \beta_s & 0\\ \alpha & 0 & 0 & 0 & 0 & \beta_n \end{bmatrix}\begin{bmatrix}1 \\ Temp_1 \\ Temp_2 \\ Temp_3 \\ Temp_4 \\ Temp_5 \end{bmatrix} + \begin{bmatrix} e_1\\ e_2\\ e_3 \\ e_4 \\ e_5\end{bmatrix}=\mathbf{Z}\mathbf{x}+\mathbf{e}\]
The \(\mathbf{Z}\) matrix has the parameter names and the \(\mathbf{x}\) matrix has the covariates with the row a 1 for the intercept.
Optional. Solve for the parameters.
Solving for this is a matter of writing \(\mathbf{y}\), \(\mathbf{Z}\), and \(\mathbf{x}\) in R.
n=5 #the number of data points
y=matrix(homeworkdat$Ozone, ncol=1)
x=matrix(c(1, homeworkdat$Temp),ncol=1)
#Set up Z; it is 5 x 6
Z=matrix(list(0),n,n+1)
#make the alpha column
Z[,1]="alpha"
#make the diagonal of columns 2:6 equal to the betas
diag(Z[,2:(1+n)])=paste("beta",homeworkdat$region, sep=".")Then we solve for the parameters:
require(MARSS)
P=MARSS:::convert.model.mat(Z)$free[,,1]
M=kronecker(t(x),diag(n))%*%P
solve(t(M)%*%M)%*%t(M)%*%y## [,1]
## alpha 23.86230424
## beta.north 0.01510679
## beta.south 0.05654090coef(lm(Ozone ~ Temp:region, data=homeworkdat))## (Intercept) Temp:regionnorth Temp:regionsouth
## 23.86230424 0.01510679 0.05654090a. Using the airquality dataset with 111 data points, write the model below in Form 2.
fit=lm(Ozone ~ -1 + Temp:region + Month, data=airquality)The first step is to write out what model is being fit.
\[ Ozone_i = \alpha_j + \beta_k Temp_i + e_i\] If \(i\) is from the \(j\)-th month, the intercept is \(\alpha_j\). If \(i\) is from the north, \(beta_k\) is \(\beta_n\). If it is from the south, \(beta_k\) is \(\beta_s\).
Next, let’s look at region and Month.
airquality$region[1:10]## [1] "north" "south" "north" "south" "north" "south" "north" "south"
## [9] "north" "south"airquality$Month## [1] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 7 7
## [36] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8
## [71] 8 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
## [106] 9 9 9 9 9 9
## Levels: 5 6 7 8 9Region is alternating north/south and month is grouped.
So in matrix form our model looks like \[\begin{bmatrix}Ozone_1\\ Ozone_2\\ \dots \\ Ozone_{111}\end{bmatrix} = \begin{bmatrix} \alpha_5 & \beta_n & 0 & \dots & 0 & 0\\ \alpha_5 & 0 &\beta_s & \dots & 0 & 0\\ \dots & \dots & \dots & \ddots & \dots & 0\\ \alpha_9& 0 & 0 & \dots & \beta_s & 0\\ \alpha_9& 0 & 0 & \dots & 0 & \beta_n \end{bmatrix}\begin{bmatrix}1 \\ Temp_1 \\ Temp_2 \\ \dots \\ Temp_{111} \end{bmatrix} + \begin{bmatrix} e_1\\ e_2\\ \dots \\ e_{111}\end{bmatrix}=\mathbf{Z}\mathbf{x}+\mathbf{e}\]
b. Solve for the parameters and show they match what lm() returns.*
Solving for this is a matter of writing \(\mathbf{y}\), \(\mathbf{Z}\), and \(\mathbf{x}\) in R. Notice how we don’t have to tedious create character strings for our parameter names. We just use the text in our data.frame and use the paste() function.
n=111 #the number of data points
y=matrix(airquality$Ozone, ncol=1)
x=matrix(c(1, airquality$Temp),ncol=1)
#Set up Z; it is 111 x 112
Z=matrix(list(0),n,n+1)
#make the alpha column
Z[,1]=paste("alpha",airquality$Month, sep="")
#make the diagonal of columns 2:112 equal to the betas
diag(Z[,2:(1+n)])=paste("beta",airquality$region, sep=".")Then we solve for the parameters:
require(MARSS)
P=MARSS:::convert.model.mat(Z)$free[,,1]
M=kronecker(t(x),diag(n))%*%P
solve(t(M)%*%M)%*%t(M)%*%y## [,1]
## alpha5 -160.254743
## alpha6 -187.708278
## alpha7 -173.611725
## alpha8 -172.152400
## alpha9 -181.929911
## beta.north 2.790616
## beta.south 2.758379coef(lm(Ozone ~ -1 + Temp:region + Month, data=airquality))## Month5 Month6 Month7 Month8
## -160.254743 -187.708278 -173.611725 -172.152400
## Month9 Temp:regionnorth Temp:regionsouth
## -181.929911 2.790616 2.758379