Model Assumptions

Let:

• $y_i\in R^n$ as the measured response for the $i$th subject.
  
• $x_i=(1,x_{i1},...,x_{ip})\in R^{p+1}$, $i=1,2,...,n$, be the values of the explanatory values for the $i$th subject. The first entry is used for the intercept term.

Using ’ as transpose operator, the linear regression model assumes that:

$$\begin{array}{cc}{Y}_i& =x_i^{\prime }\beta +{ϵ}_i,i=1,...,n\\ Y& =X\beta +ϵ\end{array}$$

where:

• $\beta =({\beta }_0,{\beta }_1,...,{\beta }_p{)}^{\prime }\in R^{p+1}$ is the unknown regression coefficient vector.
• Design matrix $X\in R^{n\times p+1}$ have $x_i^{\prime }$ at the $i$th row.
  
• $Y=(Y_1,...,Y_n{)}^{\prime }$ is a random vector.
• $ϵ=({ϵ}_1,{ϵ}_2,...,{ϵ}_n{)}^{\prime }\stackrel{i.i.d.}{\sim }(0,{\sigma }^2)$, are i.i.d error terms from unknown random distribution with mean 0 and unknown fixed variance ${\sigma }^2$.
  • $cov({ϵ}_i,{ϵ}_i)=var\left({ϵ}_i\right)={\sigma }^2$
  • Due to i.i.d. assumption, $cov({ϵ}_i,{ϵ}_j)=0,i\ne j$
  • Thus, the diagonal elements of $var\left(ϵ\right)$ are equal to ${\sigma }^2$ and the off-diagonal elements are zero. $$var\left(ϵ\right)=\left[\begin{array}{cccc}{\sigma }^2& 0& \cdots & 0\\ 0& {\sigma }^2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\sigma }^2\end{array}\right]$$
  • We can write $var\left(ϵ\right)={\sigma }^2{I}_n$, where $I_n$ is a $n\times n$ identity matrix, with ones on its diagonal and zero everywhere else.

Note that no assumptions are imposed onto the design matrix X: the values can be chosen by experimenter or they could be randomly realised from any distributions.

Ordinary Least Square Estimator (OLS) of $\beta$

Assuming that $(X^{\prime }X{)}^{-1}$ exists, that is $\left(X^{\prime }X\right)$ is invertible, the ordinary least square estimator of $\beta$ is define by minimizing the error terms, the $ϵ=Y-XB$. Specifically, we are to minimizing the squared error terms.

$$\hat{\beta }=\underset{b\in R^p}{argmin}min(Y-Xb{)}^{\prime }(Y-Xb)$$

It can be shown that the solution is $$\hat{\beta }=(X^{\prime }X{)}^{-1}{X}^{\prime }Y$$ given that the $(X^{\prime }X{)}^{-1}$ exists.

Estimator of ${\sigma }^2$

An unbiased estimator of ${\sigma }^2$ is $$\begin{array}{cc}{\hat{\sigma }}_{ϵ}^2& =\frac{1}{n-p}\sum _{i=1}^n{ϵ}^2\\ & =\frac{1}{n-p}(Y-X\hat{\beta }{)}^{\prime }(Y-X\hat{\beta })\end{array}$$

Visualize data points

Let’s visualize the data.

library(plotly)
pcolor = rep("purple", length(x1))
plot_ly(df_lm,
        x = ~x1,
        y = ~x2,
        z = ~y,
        type = "scatter3d", 
        mode = "markers",
        marker = list(color = pcolor, 
                      opacity = 0.4))

Estimating $\beta$ and $ϵ$

The $\beta$s can be calculated through two ways: using formulas and using R function.

# Using formula
beta_form = qr.solve(crossprod(xmat)) %*% crossprod(xmat,y)
# USing r function
beta_r = lm.fit(x = xmat, y = y)$coefficients
compare = cbind(beta_form, beta_r, beta)
colnames(compare)= c("formula", "r_function", "actual")
compare

##      formula r_function actual
##    9.8791689  9.8791689   10.0
## x1 2.0064254  2.0064254    2.0
## x2 0.5127404  0.5127404    0.5

The calculated (estimated) value for the $\beta$s are similar to the one that was used to generate the model, that is 10, 2, and 0.5 for ${\beta }_0$, ${\beta }_1$ and ${\beta }_2$ respectively.

The realization of ${\hat{\sigma }}_{ϵ}^2$ is:

eps_hat = 1/(n - p)*sum((y - xmat %*% beta_form)^2)
eps_hat

## [1] 17.76159

The original ${\sigma }^2$ is 16 and the realized value is 17.7615874.

Visualize fitted plane

So what’s is the fitted plane on the data points? Let’s visualize it.

library(reshape2)
# predict over sensible grid of values
x1s <- seq(min(x1), max(x1), length.out = 100)
x2s <- seq(min(x2), max(x2), length.out = 100)
surface <- expand.grid(x1 = x1s,x2 = x2s,KEEP.OUT.ATTRS = F)

# Making prediction
surface$y <- apply(surface[,c("x1", "x2")], 1, function(x) cbind(1,x[1], x[2]) %*% beta_form)
surface <- acast(surface, x2 ~ x1, value.var = "y")

pcolor = rep("purple", length(x1))
plot_ly(df_lm,
        x = ~x1,
        y = ~x2,
        z = ~y,
        type = "scatter3d", 
        mode = "markers",
        marker = list(color = pcolor, 
                      opacity = 0.7)) %>%
  add_trace(z = surface,
            x = x1s,
            y = x2s,
            type = "surface",
            opacity = 0.6)

We can see that the plane fitted from $\hat{\beta }$ follows the data points pretty closely. However, note that Linear models are very rigid (i.e., a plane).