In charlotte-ngs/asmss2022: Applied Statistical Methods in Animal Science - Spring Semester 2022
Linear Regression {#asm-regr}
met$set_this_rmd_file(ps_this_rmd_file = ifelse(rstudioapi::isAvailable(),
rstudioapi::getSourceEditorContext()$path,
whereami::thisfile()))
Introduction {#introduction-regr}
This chapter is based on the book r met$add("Searle1971") and on the course notes r met$add("Buhlmann2016"). An interesting online book on Multiple Linear Regression is given by r met$add("Lilja2016"). Regression analysis is used to assess relationships between a given variable and other measurements or observations on the same animal. The relationships between the variable and the other characteristics of the animal are estimated based on observed data.
Example {#example-regr}
A classical example of a regression analysis in animal science is the relationship between body weight and breast circumference in cattle. This example has a practical application because the results of the regression analysis of body weight on breast circumference can be used for a measuring band. With such a measuring band the breast circumference of an animal is measured. On the back side of the band, the estimated body weight can directly be determined.
At this point the question is how is it possible to determine the relationship between the values of breast circumference in centimeters and body weight in kilograms. The answer to this question can be given by a regression analysis. The most important pre-requisites for doing a regression analysis is to have a dataset. For our example, Table \@ref(tab:tab-reg-bw-bc) shows such a dataset which can be used for a regression analysis.
tbl_reg <- tibble::tibble(Animal = c(1:10),
`Breast Circumference` = c(176, 177, 178, 179, 179, 180, 181, 182,183, 184),
`Body Weight` = c(471, 463, 481, 470, 496, 491, 518, 511, 510, 541))
n_nr_obs <- nrow(tbl_reg)
knitr::kable(tbl_reg,
booktabs = TRUE,
longtable = TRUE,
caption = paste0("Breast Circumference and Body Weight for ",
nrow(tbl_reg)," Animals", collapse = ""),
escape = FALSE)
The dataset in Table \@ref(tab:tab-reg-bw-bc) contains measurements of body weight and breast circumference for $r nrow(tbl_reg)$ animals. Figure \@ref(fig:fig-reg-bw-bc) is a graphical representation of our example dataset.
library(ggplot2)
ggplot(tbl_reg, aes(x = `Breast Circumference`, y = `Body Weight`)) +
geom_point(color = "blue")
The diagram in Figure \@ref(fig:fig-reg-bw-bc) shows on the x-axis the breast circumference in cm and on the y-axis the body weight in kg. Each of the blue dots correspond to an observation of one animal in the dataset. A diagram like the one shown in Figure \@ref(fig:fig-reg-bw-bc) is also called a dot plot. From a first visual inspection of the dot plot for our dataset, we can see that there is a tendency that larger values of breast circumference of animals are related to heavier animals. The relationship is not deterministic that means there are exceptions which do not follow the rule of this relationship. One example of an exception are animals 1 and 2. Animal 2 has a larger breast circumference value compared to animal 1, but animal 2 has a lower body weight compared to animal 1. But despite such exceptions, we can still observe that on average there is a relationship between breast circumference and body weight. Furthermore, the apparent relationship between breast circumference and body weight seams to be the same for low and high values of breast circumference. Based on this last fact, we can say that the relationship between breast circumference and body weight is linear.
# see https://rpubs.com/kaz_yos/ggplot2-stat-function and
# https://ggplot2.tidyverse.org/reference/geom_function.html#examples
cn_odg_path <- file.path(here::here(), "cn", "odg")
library(ggplot2)
# quadratic function
p_quad <- ggplot() + xlim(-5,5)
fun.1 <- function(x) x^2 + x
p_quad <- p_quad + stat_function(fun = fun.1)
ggsave(filename = file.path(cn_odg_path, "quadratic_function_plot.png"))
p_log <- ggplot() + xlim(0,10)
p_log <- p_log + stat_function(fun = function(x) log(x))
ggsave(filename = file.path(cn_odg_path, "log_function_plot.png"))
#rmdhelp::use_odg_graphic(ps_path = "odg/fig-non-linear-rel.odg")
knitr::include_graphics(path = "odg/fig-non-linear-rel.png")
Figure \@ref(fig:fig-non-linear-rel) shows two examples of non-linear functions. Both examples show that the relationship between the shown variables $x$ and $y$ is not the same over the shown range of values. Hence by inspecting the diagram of the two variables breast circumference and body weight of our example in Figure \@ref(fig:fig-reg-bw-bc), we can say that the relationship between the variables of our example dataset show a linear relationship.
Regression Model {#regression-model-regr}
With the regression model we want to find a mathematical formulation to describe and to quantify the relationship between the two variables from our example, breast circumference and body weight. One possibility to find this relationship is to take an animal with $x$ cm of breast circumference. Then the question is what would be the expected value for its body weight $y$ in kg. Under the assumption of a linear relationship between the variables from our example, the expected value $E(y)$ for the body weight $y$ can be written as
\begin{equation}
E(y) = b_0 + b_1 * x
(#eq:exbwregmodelregr)
\end{equation}
The above reasoning on how the variables are related is often referred to as model building. The model shown in \@ref(eq:exbwregmodelregr) is called a linear model, because the expected body weight ($E(y)$) is a linear function of the unknown parameters $b_0$ and $b_1$. The number of possible models between two variables $x$ and $y$ is infinite. And therefore the process of model building is difficult and requires some experience. But in general, we can say that a simpler model with fewer unknown parameters is always preferable over a more complex model as long as the simpler model is able to capture the most important aspects of a relationship between two variables.
Observations {#observations-regr}
For an animal with a breast circumference of $x$ cm, the body weight ($y$) will not exactly be $b_0 + b_1*x$. It has to be noted here that the values for $b_0$ and $b_1$ will be the same for all animals. The fact of the discrepancy between the recorded body weights ($y$) and the output of the model is taken into account by writing $E(y)$ instead of $y$ in the model shown in equation \@ref(eq:exbwregmodelregr). For a given observed body weight $y_i$ of animal $i$ with a breast circumference of $x_i$, we can write
\begin{equation}
E(y_i) = b_0 + b_1 * x_i
(#eq:bwregmodelanimaliregr)
\end{equation}
where $E(y_i)$ is not the same as $y_i$. The difference $y_i - E(y_i)$ represents the difference between the observed body weight from its expected value $E(y_i)$ and is written as
\begin{equation}
e_i = y_i - E(y_i) = y_i - b_0 - b_1 * x_i
(#eq:exbwregmodelresidualregr)
\end{equation}
Hence for the body weight $y_i$ of animal $i$, we can write
\begin{equation}
y_i = b_0 + b_1 * x_i + e_i
(#eq:exbwregmodelobservedvalueregr)
\end{equation}
Equations \@ref(eq:bwregmodelanimaliregr), \@ref(eq:exbwregmodelresidualregr) and \@ref(eq:exbwregmodelobservedvalueregr) apply to all observations $y_1, y_2, \ldots, y_{r n_nr_obs}$, of our example dataset shown in Table \@ref(tab:tab-reg-bw-bc). The $e_i$ terms for all observations might take many different values. They include potential measurement errors or deficiencies of the model itself. Due to the described properties of $e_i$'s, they are considered to be random variables and are usually called random errors or random residuals.
To complete the description of our model in terms of equation \@ref(eq:exbwregmodelobservedvalueregr), further characteristics of the random errors ($e_i$) must be specified. These characteristics consist of
- the expected value $E(e_i)$ of $e_i$ and
- the variance $var(e_i)$ of $e_i$.
Useful specifications are that the expected value $E(e_i)$ is zero and all covariances between any pair of $e_i$ and $e_j$ with $i \ne j$ are also zero. Then the variance $var(e_i)$ is assumed to be a constant for all $i$ and is represented by the symbol $\sigma^2$. Summarizing all the proposed properties in a mathematical notation, we obtain
\begin{equation}
E(e_i) = 0
(#eq:expvalueerrorregr)
\end{equation}
which is obtained from the definition of $e_i$ given in \@ref(eq:exbwregmodelresidualregr). The variance $var(e_i)$
\begin{equation}
var(e_i) = E\left[ e_i - E(e_i) \right]^2 = E(e_i^2) = \sigma^2
(#eq:varerrorregr)
\end{equation}
and
\begin{equation}
cov(e_i,e_j) = E\left[ e_i - E(e_i) \right]\left[ e_j - E(e_j) \right] = E(e_ie_j) = 0
(#eq:coverrorregr)
\end{equation}
Equations \@ref(eq:bwregmodelanimaliregr) - \@ref(eq:coverrorregr) give a constitutional description of the linear model that we have designed so far for our example dataset. These properties form the basis for the procedure used to estimated the unknown parameters $b_0$ and $b_1$.
Parameter Estimation {#parameter-estimation-regr}
There are several methods to estimate the unknown parameters $b_0$ and $b_1$ of the proposed linear model. The most frequently used method which is also implemented in the R-function lm() is called least squares.
Least squares estimation is based on the idea of minimizing the sum of the squared deviations of the observations $y_i$ from their expected values. This sum can be written as
\begin{equation}
\mathbf{e}^T\mathbf{e} = \sum_{i=1}^N e_i^2 = \sum_{i=1}^N \left[ y_i - E(e_i) \right]^2 = \sum_{i=1}^N \left[ y_i - b_0 - b_1*x_i \right]^2
(#eq:sumsquareresidualregr)
\end{equation}
where $\mathbf{e}^T = \left[\begin{array}{cccc}e_1 & e_2 & \ldots & e_N \end{array}\right]$ is the transpose of the vector $\mathbf{e}$ of length $N$ containing all the $e_i$ values and $N$ stands for the number of observations.
Although $b_0$ and $b_1$ are fixed (but unknown) values, we treat them for a moment like mathematical variables. The reason for this changed view is that we want to find the values for $b_0$ and $b_1$ that minimize the expression in \@ref(eq:sumsquareresidualregr). The resulting values from the minimization of \@ref(eq:sumsquareresidualregr) will be represented by the symbols $\hat{b}_0$ and $\hat{b}_1$ and they will be called the least squares estimators of $b_0$ and $b_1$.
Minimization of \@ref(eq:sumsquareresidualregr) is done by taking partial derivatives with respect to both unknowns $b_0$ and $b_1$ and setting both derivatives to zero^[The verification of higher order derivative to confirm that the obtained extreme value is a minimum is not done here.]. This yields two equations from which solutions called $\hat{b}_0$ and $\hat{b}_1$ can be computed.
\begin{align}
\frac{\partial\mathbf{e}^T\mathbf{e}}{\partial b_0} &= -2 \sum_{i=1}^N \left[y_i - b_0 - b_1x_i\right] \notag \
&= -2\left[\sum_{i=1}^N y_i - Nb_0 - b_1\sum_{i=1}^N x_i\right]
(#eq:derivssqresidualb0regr)
\end{align}
\begin{align}
\frac{\partial\mathbf{e}^T\mathbf{e}}{\partial b_1} &= -2 \sum_{i=1}^N x_i\left[y_i - b_0 - b_1x_i\right] \notag \
&= -2 \left[\sum_{i=1}^N x_iy_i - b_0 \sum_{i=1}^N x_i - b_1 \sum_{i=1}^N x_i^2 \right]
(#eq:derivssqresidualb1regr)
\end{align}
Setting both expression in \@ref(eq:derivssqresidualb0regr) and \@ref(eq:derivssqresidualb1regr) to zero and writing them in terms of $\hat{b}_0$ and $\hat{b}_1$ gives
\begin{equation}
N\hat{b}_0 + \hat{b}_1x. = y.
(#eq:solderivssqresidualb0regr)
\end{equation}
and
\begin{equation}
\hat{b}_0x. + \hat{b}_1(x^2). = (xy).
(#eq:solderivssqresidualb1regr)
\end{equation}
using the dot notation for the following sums $x. = \sum_{i=1}^N x_i$, $y. = \sum_{i=1}^N y_i$, $(x^2). = \sum_{i=1}^N x_i^2$ and $(xy). = \sum_{i=1}^N x_iy_i$. With the bar notation for the means, we can further write
\begin{equation}
\bar{x}. = {x. \over N}
(#eq:barxdotregr)
\end{equation}
and
\begin{equation}
\bar{y}. = {y. \over N}
(#eq:barydotregr)
\end{equation}
The solutions in \@ref(eq:solderivssqresidualb0regr) and \@ref(eq:solderivssqresidualb1regr) can then be written as
\begin{equation}
\hat{b}_0 = \bar{y}. - \hat{b}_1\bar{x}.
(#eq:solderivssqresidualb0regr)
\end{equation}
and
\begin{equation}
\hat{b}_1 = \frac{(xy). - N\bar{x}.\bar{y}.}{(x^2). - N\bar{x}.^2}
(#eq:solderivssqresidualb1regr)
\end{equation}
Estimates for Example Dataset {#estimates-for-example-dataset-regr}
The estimates as shown in \@ref(eq:solderivssqresidualb0regr) and \@ref(eq:solderivssqresidualb1regr) are computed for our example dataset. We start by computing all the components of the formulas for the estimators, then we plug those components in and get the results.
# formula components
x_dot <- sum(tbl_reg$`Breast Circumference`)
x_bar <- x_dot / n_nr_obs
y_dot <- sum(tbl_reg$`Body Weight`)
y_bar <- y_dot / n_nr_obs
xy_dot <- sum(tbl_reg$`Breast Circumference`*tbl_reg$`Body Weight`)
x2_dot <- sum(tbl_reg$`Breast Circumference`^2)
# estimators
b1_hat <- (xy_dot - n_nr_obs * x_bar * y_bar) / (x2_dot - n_nr_obs * x_bar^2)
b0_hat <- y_bar - b1_hat * x_bar
$$
N = r n_nr_obs,\ \bar{x}. = r x_bar,\ \bar{y}. = r y_bar,\ (xy). = r xy_dot,\ (x^2). = r x2_dot
$$
\begin{equation}
\hat{b}_1 = \frac{r xy_dot - r n_nr_obs * r x_bar * r y_bar}{r x2_dot - r n_nr_obs * r x_bar^2} = r round(b1_hat, digits=3)
(#eq:numsolderivssqresidualb1regr)
\end{equation}
\begin{equation}
\hat{b}_0 = r y_bar - r b1_hat * r x_bar = r round(b0_hat, digits=3)
(#eq:numsolderivssqresidualb0regr)
\end{equation}
Obtain Parameter Estimates in R {#obtain-parameter-estimates-in-r-regr}
As already mentioned, the same type of computation as shown in Section \@ref(estimates-for-example-dataset-regr) is also implemented in the R-function lm(). In what follows, we show how the estimates of the linear regression model are obtained using lm(). An important pre-requisite for using the function lm() is that the dataset is assigned to a dataframe. Here we assume that we have a dataframe named tbl_reg that contains our dataset. Then the parameter estimates are obtained using the following statements in R.
lm_bw_bc <- lm(`Body Weight` ~ `Breast Circumference`, data = tbl_reg)
summary(lm_bw_bc)
The output of the summary function shows a lot of information about the data and the model. Under the section Coefficients there are two entries
(Intercept) corresponding to our $b_0$ and Breast Circumference corresponding to our $b_1$.
The values in the first column entitled Estimate correspond to the values that we have computed in the previous section.
The General Case {#the-general-case-regr}
Suppose that in the study on body weight and breast circumference, we have an additional observation for each animal consisting of the height of the animal. The new extended data set is shown in Table \@ref(tab:tab-bw-mult-reg).
tbl_bw_ext <- dplyr::bind_cols(tbl_reg,
tibble::tibble(Height = c(161, 121, 157, 165, 136, 123, 163, 149, 143, 130)))
knitr::kable(tbl_bw_ext,
booktabs = TRUE,
longtable = TRUE,
caption = paste0("Extended Dataset of Body Weight for ",
nrow(tbl_reg)," Animals", collapse = ""),
escape = FALSE)
The model developed so far is not extended to be
\begin{equation}
E(y) = b_0 + b_1x_1 + b_2x_2
(#eq:extendedbwmodelregr)
\end{equation}
where $x_1$ represents the breast circumference in cm, and $x_2$ stands for the height of the animal in cm. Thus for animal $i$ with a breast circumference of $x_{i1}$ cm and a height of $x_{i2}$ the body weight $y_i$ can be written as
\begin{equation}
y_i = b_0 + b_1x_{i1} + b_2x_{i2} + e_i
(#eq:extendedbwdecompregr)
\end{equation}
As the number of $x$ variables increase, the equations such as shown in \@ref(eq:extendedbwdecompregr) are getting longer and they are getting tedious to handle. This problem is solved by a change of notation. Let us define the matrix $\mathbf{X}$ and the vectors $\mathbf{y}$, $\mathbf{e}$ and $\mathbf{b}$ as follows.
mat_X <- matrix(c("x_{10}", "x_{11}", "x_{12}",
"x_{20}", "x_{21}", "x_{22}",
".", ".", ".",
".", ".", ".",
".", ".", ".",
"x_{N0}", "x_{N1}", "x_{N2}"), ncol = 3, byrow = TRUE)
vec_y <- c("y_1", "y_2", ".", ".", ".", "y_N")
vec_e <- c("e_1", "e_2", ".", ".", ".", "e_N")
vec_b <- c("b_0", "b_1", "b_2")
cat("$$\n")
cat(paste0(rmdhelp::bmatrix(pmat = mat_X, ps_name = "\\mathbf{X}"), collapse = ""), "\n")
cat(", \\ \n")
cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_y, ps_name = "\\mathbf{y}")))
cat(", \\ \n")
cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_e, ps_name = "\\mathbf{e}")))
cat("\\text{ and }\n")
cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_b, ps_name = "\\mathbf{b}")))
cat("$$\n")
Because all equations contain the term $b_0$, the first column of matrix $\mathbf{X}$ consists of all ones. Hence, $x_{10} = x_{20} = \ldots = x_{N0} = 1$. The complete set of equations for all animals in our extended dataset represented by \@ref(eq:extendedbwdecompregr) is
\begin{equation}
\mathbf{y} = \mathbf{X}\mathbf{b} + \mathbf{e} \text{, with } E(\mathbf{y}) = \mathbf{X}\mathbf{b}
(#eq:matvecextendedbwdecompregr)
\end{equation}
The big advantage of the matrix-vector notation is that equation \@ref(eq:matvecextendedbwdecompregr) is invariant to the number of $x$-variables. That means no matter how many $x$-variables we include in our dataset, the linear regression model can always be written as shown in equation \@ref(eq:matvecextendedbwdecompregr). The only thing that changes are the definitions of $\mathbf{X}$ and $\mathbf{b}$. In the general case with $k$ variables
mat_X_gen <- matrix(c("x_{10}", "x_{11}", ".", "x_{1k}",
"x_{20}", "x_{21}", ".", "x_{2k}",
".", ".", ".", ".",
".", ".", ".", ".",
".", ".", ".", ".",
"x_{N0}", "x_{N1}", ".", "x_{Nk}"), ncol = 4, byrow = TRUE)
vec_b_gen <- c("b_0", "b_1", ".", ".", "b_k")
cat("$$\n")
cat(paste0(rmdhelp::bmatrix(pmat = mat_X_gen, ps_name = "\\mathbf{X}"), collapse = ""), "\n")
cat(", \\ \n")
cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_b_gen, ps_name = "\\mathbf{b}")))
cat("$$\n")
The topic of parameter estimation can also be shown using matrix-vector notation. The only restriction that we are currently imposing is that the number of $x$-variables $k$ is smaller than the number of observations $N$.
The model specification is only complete, if the properties of the vector $\mathbf{e}$ of random residuals is defined. In accordance with equations \@ref(eq:expvalueerrorregr), \@ref(eq:varerrorregr) and \@ref(eq:coverrorregr), we can write
\begin{equation}
E(\mathbf{e}) = \mathbf{0}
(#eq:matvecexpvalueerrorregr)
\end{equation}
where $E(\mathbf{e})$ stands for the vector of length $N$ containing the expected values of all random residuals and $\mathbf{0}$ is a vector of $N$ zeros and
\begin{equation}
var(\mathbf{e}) = E\left[\mathbf{e} - E(\mathbf{e}) \right]\left[\mathbf{e} - E(\mathbf{e}) \right]^T
= E(\mathbf{e}\mathbf{e}^T) = \sigma^2 \mathbf{I}_N
(#eq:matvecvarerrorregr)
\end{equation}
where $var(\mathbf{e})$ is the variance-covariance matrix between all random residuals and $\mathbf{I}_N$ is the $N\times N$ identity matrix.
Derivation of the least squares estimator of $\mathbf{b}$ follows the same principles as shown in equations \@ref(eq:derivssqresidualb0regr) - \@ref(eq:solderivssqresidualb1regr). The sum of squares of the deviations of the observations from their expected values using $E(\mathbf{e}) = \mathbf{0}$ and hence $E(\mathbf{y}) = \mathbf{Xb}$, is
\begin{align}
\mathbf{e}^T\mathbf{e} &= \left[\mathbf{y} - E(\mathbf{y}) \right]^T\left[\mathbf{y} - E(\mathbf{y}) \right] \notag \
&= \left[\mathbf{y} - \mathbf{Xb} \right]^T\left[\mathbf{y} - \mathbf{Xb} \right]\notag \
&= \mathbf{y}^T\mathbf{y} - 2 \mathbf{b}^T\mathbf{X}^T\mathbf{y} + \mathbf{b}^T\mathbf{X}^T\mathbf{X}\mathbf{b}
(#eq:matvecsumsqresidualregr)
\end{align}
The least squares estimator $\hat{\mathbf{b}}$ is found by minimizing $\mathbf{e}^T\mathbf{e}$ with respect to all elements of $\mathbf{b}$. This is corresponds to $\partial{\mathbf{e}^T\mathbf{e}}/\partial \mathbf{b}$ and is also called the gradient of $\mathbf{e}^T\mathbf{e}$ with respect to $\mathbf{b}$. Setting the gradient to zero leads to the following equations
\begin{equation}
\mathbf{X}^T\mathbf{X}\hat{\mathbf{b}} = \mathbf{X}^T\mathbf{y}
(#eq:leastsquarenormalequationregr)
\end{equation}
These equations are known as the least squares normal equations. Provided $(\mathbf{X}^T\mathbf{X})$ can be inverted, the unique solution for $\hat{\mathbf{b}}$ can be written as
\begin{equation}
\hat{\mathbf{b}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}
(#eq:solleastsquareestregr)
\end{equation}
Example Dataset Continued {#example-dataset-continued-regr}
We use the original dataset shown in Table \@ref(tab:tab-reg-bw-bc) to illustrate the estimation of the least squares parameters using the matrix-vector notation. The matrix $\mathbf{X}$ and the vectors $\mathbf{y}$, $\mathbf{b}$ and $\mathbf{e}$ are defined as follows.
nr_lines_shown <- 2
nr_lines_dot <- 3
mat_X_short <- matrix(c(rep(1, nr_lines_shown),
rep('.', nr_lines_dot), 1,
tbl_reg$`Breast Circumference`[1:nr_lines_shown],
rep('.', nr_lines_dot), tbl_reg$`Breast Circumference`[n_nr_obs]), ncol = 2)
vec_y_short <- c(tbl_reg$`Body Weight`[1:nr_lines_shown], rep('.', nr_lines_dot), tbl_reg$`Body Weight`[n_nr_obs])
vec_b_short <- c("b_0", "b_1")
vec_e_short <- c(rmdhelp::vecGetVecElem(psBaseElement = "e", pnVecLen = nr_lines_shown, psResult = 'latex'),
rep('.', nr_lines_dot), paste0("e_{", n_nr_obs, "}", collapse = ""))
# show matrix X and vectors
cat("$$\n")
cat(paste0(rmdhelp::bmatrix(pmat = mat_X_short, ps_name = "\\mathbf{X}"), collapse = ""), "\n")
cat("\\text{, }\n")
cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_y_short, ps_name = "y")))
cat("\\text{, }\n")
cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_b_short, ps_name = "b")))
cat("\\text{, }\n")
cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_e_short, ps_name = "e")))
cat("$$\n")
The solution is obtained from equation \@ref(eq:solleastsquareestregr) for this we also need
mat_X_num <- matrix(c(rep(1,n_nr_obs), tbl_reg$`Breast Circumference`), ncol = 2)
mat_XTX <- crossprod(mat_X_num)
mat_XTXinv <- solve(mat_XTX)
vec_y <- tbl_reg$`Body Weight`
mat_XTy <- crossprod(mat_X_num, vec_y)
vec_sol <- solve(mat_XTX, mat_XTy)
# show the matrices
cat("$$\n")
cat(paste0(rmdhelp::bmatrix(pmat = mat_XTX, ps_name = "X^TX"), collapse = ""), "\n")
cat("\\text{, }\n")
cat(paste0(rmdhelp::bmatrix(pmat = round(mat_XTXinv, digits = 3), ps_name = "(X^TX)^{-1}"), collapse = ""), "\n")
cat("\\text{, }\n")
cat(paste0(rmdhelp::bmatrix(pmat = round(mat_XTy, digits = 3), ps_name = "X^Ty"), collapse = ""), "\n")
cat("$$\n")
The solution vector $\hat{\mathbf{b}}$ contains the two components corresponding to the estimates of the two unknown parameters.
vec_sol_sym <- c("\\hat{b}_0", "\\hat{b}_1")
# show the solutions
cat("$$\n")
cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_sol_sym, ps_name = "\\hat{\\mathbf{b}}"), collapse = ""), "\n")
cat(" = ")
cat(paste0(rmdhelp::bmatrix(pmat = round(vec_sol, digits = 3)), collapse = ""), "\n")
cat("$$\n")
Comparing the above shown solutions to the results received earlier shows that for this example, the solution using the matrix-vector notation are the same.
charlotte-ngs/asmss2022 documentation built on June 7, 2022, 1:33 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Linear Regression {#asm-regr}
met$set_this_rmd_file(ps_this_rmd_file = ifelse(rstudioapi::isAvailable(), rstudioapi::getSourceEditorContext()$path, whereami::thisfile()))
Introduction {#introduction-regr}
This chapter is based on the book r met$add("Searle1971") and on the course notes r met$add("Buhlmann2016"). An interesting online book on Multiple Linear Regression is given by r met$add("Lilja2016"). Regression analysis is used to assess relationships between a given variable and other measurements or observations on the same animal. The relationships between the variable and the other characteristics of the animal are estimated based on observed data.
Example {#example-regr}
A classical example of a regression analysis in animal science is the relationship between body weight and breast circumference in cattle. This example has a practical application because the results of the regression analysis of body weight on breast circumference can be used for a measuring band. With such a measuring band the breast circumference of an animal is measured. On the back side of the band, the estimated body weight can directly be determined.
At this point the question is how is it possible to determine the relationship between the values of breast circumference in centimeters and body weight in kilograms. The answer to this question can be given by a regression analysis. The most important pre-requisites for doing a regression analysis is to have a dataset. For our example, Table \@ref(tab:tab-reg-bw-bc) shows such a dataset which can be used for a regression analysis.
tbl_reg <- tibble::tibble(Animal = c(1:10), `Breast Circumference` = c(176, 177, 178, 179, 179, 180, 181, 182,183, 184), `Body Weight` = c(471, 463, 481, 470, 496, 491, 518, 511, 510, 541)) n_nr_obs <- nrow(tbl_reg) knitr::kable(tbl_reg, booktabs = TRUE, longtable = TRUE, caption = paste0("Breast Circumference and Body Weight for ", nrow(tbl_reg)," Animals", collapse = ""), escape = FALSE)
The dataset in Table \@ref(tab:tab-reg-bw-bc) contains measurements of body weight and breast circumference for $r nrow(tbl_reg)$ animals. Figure \@ref(fig:fig-reg-bw-bc) is a graphical representation of our example dataset.
library(ggplot2) ggplot(tbl_reg, aes(x = `Breast Circumference`, y = `Body Weight`)) + geom_point(color = "blue")
The diagram in Figure \@ref(fig:fig-reg-bw-bc) shows on the x-axis the breast circumference in cm and on the y-axis the body weight in kg. Each of the blue dots correspond to an observation of one animal in the dataset. A diagram like the one shown in Figure \@ref(fig:fig-reg-bw-bc) is also called a dot plot. From a first visual inspection of the dot plot for our dataset, we can see that there is a tendency that larger values of breast circumference of animals are related to heavier animals. The relationship is not deterministic that means there are exceptions which do not follow the rule of this relationship. One example of an exception are animals 1 and 2. Animal 2 has a larger breast circumference value compared to animal 1, but animal 2 has a lower body weight compared to animal 1. But despite such exceptions, we can still observe that on average there is a relationship between breast circumference and body weight. Furthermore, the apparent relationship between breast circumference and body weight seams to be the same for low and high values of breast circumference. Based on this last fact, we can say that the relationship between breast circumference and body weight is linear.
# see https://rpubs.com/kaz_yos/ggplot2-stat-function and # https://ggplot2.tidyverse.org/reference/geom_function.html#examples cn_odg_path <- file.path(here::here(), "cn", "odg") library(ggplot2) # quadratic function p_quad <- ggplot() + xlim(-5,5) fun.1 <- function(x) x^2 + x p_quad <- p_quad + stat_function(fun = fun.1) ggsave(filename = file.path(cn_odg_path, "quadratic_function_plot.png")) p_log <- ggplot() + xlim(0,10) p_log <- p_log + stat_function(fun = function(x) log(x)) ggsave(filename = file.path(cn_odg_path, "log_function_plot.png"))
#rmdhelp::use_odg_graphic(ps_path = "odg/fig-non-linear-rel.odg") knitr::include_graphics(path = "odg/fig-non-linear-rel.png")
Figure \@ref(fig:fig-non-linear-rel) shows two examples of non-linear functions. Both examples show that the relationship between the shown variables $x$ and $y$ is not the same over the shown range of values. Hence by inspecting the diagram of the two variables breast circumference and body weight of our example in Figure \@ref(fig:fig-reg-bw-bc), we can say that the relationship between the variables of our example dataset show a linear relationship.
Regression Model {#regression-model-regr}
With the regression model we want to find a mathematical formulation to describe and to quantify the relationship between the two variables from our example, breast circumference and body weight. One possibility to find this relationship is to take an animal with $x$ cm of breast circumference. Then the question is what would be the expected value for its body weight $y$ in kg. Under the assumption of a linear relationship between the variables from our example, the expected value $E(y)$ for the body weight $y$ can be written as
\begin{equation} E(y) = b_0 + b_1 * x (#eq:exbwregmodelregr) \end{equation}
The above reasoning on how the variables are related is often referred to as model building. The model shown in \@ref(eq:exbwregmodelregr) is called a linear model, because the expected body weight ($E(y)$) is a linear function of the unknown parameters $b_0$ and $b_1$. The number of possible models between two variables $x$ and $y$ is infinite. And therefore the process of model building is difficult and requires some experience. But in general, we can say that a simpler model with fewer unknown parameters is always preferable over a more complex model as long as the simpler model is able to capture the most important aspects of a relationship between two variables.
Observations {#observations-regr}
For an animal with a breast circumference of $x$ cm, the body weight ($y$) will not exactly be $b_0 + b_1*x$. It has to be noted here that the values for $b_0$ and $b_1$ will be the same for all animals. The fact of the discrepancy between the recorded body weights ($y$) and the output of the model is taken into account by writing $E(y)$ instead of $y$ in the model shown in equation \@ref(eq:exbwregmodelregr). For a given observed body weight $y_i$ of animal $i$ with a breast circumference of $x_i$, we can write
\begin{equation} E(y_i) = b_0 + b_1 * x_i (#eq:bwregmodelanimaliregr) \end{equation}
where $E(y_i)$ is not the same as $y_i$. The difference $y_i - E(y_i)$ represents the difference between the observed body weight from its expected value $E(y_i)$ and is written as
\begin{equation} e_i = y_i - E(y_i) = y_i - b_0 - b_1 * x_i (#eq:exbwregmodelresidualregr) \end{equation}
Hence for the body weight $y_i$ of animal $i$, we can write
\begin{equation} y_i = b_0 + b_1 * x_i + e_i (#eq:exbwregmodelobservedvalueregr) \end{equation}
Equations \@ref(eq:bwregmodelanimaliregr), \@ref(eq:exbwregmodelresidualregr) and \@ref(eq:exbwregmodelobservedvalueregr) apply to all observations $y_1, y_2, \ldots, y_{r n_nr_obs}$, of our example dataset shown in Table \@ref(tab:tab-reg-bw-bc). The $e_i$ terms for all observations might take many different values. They include potential measurement errors or deficiencies of the model itself. Due to the described properties of $e_i$'s, they are considered to be random variables and are usually called random errors or random residuals.
To complete the description of our model in terms of equation \@ref(eq:exbwregmodelobservedvalueregr), further characteristics of the random errors ($e_i$) must be specified. These characteristics consist of
- the expected value $E(e_i)$ of $e_i$ and
- the variance $var(e_i)$ of $e_i$.
Useful specifications are that the expected value $E(e_i)$ is zero and all covariances between any pair of $e_i$ and $e_j$ with $i \ne j$ are also zero. Then the variance $var(e_i)$ is assumed to be a constant for all $i$ and is represented by the symbol $\sigma^2$. Summarizing all the proposed properties in a mathematical notation, we obtain
\begin{equation} E(e_i) = 0 (#eq:expvalueerrorregr) \end{equation}
which is obtained from the definition of $e_i$ given in \@ref(eq:exbwregmodelresidualregr). The variance $var(e_i)$
\begin{equation} var(e_i) = E\left[ e_i - E(e_i) \right]^2 = E(e_i^2) = \sigma^2 (#eq:varerrorregr) \end{equation}
and
\begin{equation} cov(e_i,e_j) = E\left[ e_i - E(e_i) \right]\left[ e_j - E(e_j) \right] = E(e_ie_j) = 0 (#eq:coverrorregr) \end{equation}
Equations \@ref(eq:bwregmodelanimaliregr) - \@ref(eq:coverrorregr) give a constitutional description of the linear model that we have designed so far for our example dataset. These properties form the basis for the procedure used to estimated the unknown parameters $b_0$ and $b_1$.
Parameter Estimation {#parameter-estimation-regr}
There are several methods to estimate the unknown parameters $b_0$ and $b_1$ of the proposed linear model. The most frequently used method which is also implemented in the R-function lm() is called least squares.
Least squares estimation is based on the idea of minimizing the sum of the squared deviations of the observations $y_i$ from their expected values. This sum can be written as
\begin{equation} \mathbf{e}^T\mathbf{e} = \sum_{i=1}^N e_i^2 = \sum_{i=1}^N \left[ y_i - E(e_i) \right]^2 = \sum_{i=1}^N \left[ y_i - b_0 - b_1*x_i \right]^2 (#eq:sumsquareresidualregr) \end{equation}
where $\mathbf{e}^T = \left[\begin{array}{cccc}e_1 & e_2 & \ldots & e_N \end{array}\right]$ is the transpose of the vector $\mathbf{e}$ of length $N$ containing all the $e_i$ values and $N$ stands for the number of observations.
Although $b_0$ and $b_1$ are fixed (but unknown) values, we treat them for a moment like mathematical variables. The reason for this changed view is that we want to find the values for $b_0$ and $b_1$ that minimize the expression in \@ref(eq:sumsquareresidualregr). The resulting values from the minimization of \@ref(eq:sumsquareresidualregr) will be represented by the symbols $\hat{b}_0$ and $\hat{b}_1$ and they will be called the least squares estimators of $b_0$ and $b_1$.
Minimization of \@ref(eq:sumsquareresidualregr) is done by taking partial derivatives with respect to both unknowns $b_0$ and $b_1$ and setting both derivatives to zero^[The verification of higher order derivative to confirm that the obtained extreme value is a minimum is not done here.]. This yields two equations from which solutions called $\hat{b}_0$ and $\hat{b}_1$ can be computed.
\begin{align} \frac{\partial\mathbf{e}^T\mathbf{e}}{\partial b_0} &= -2 \sum_{i=1}^N \left[y_i - b_0 - b_1x_i\right] \notag \ &= -2\left[\sum_{i=1}^N y_i - Nb_0 - b_1\sum_{i=1}^N x_i\right] (#eq:derivssqresidualb0regr) \end{align}
\begin{align} \frac{\partial\mathbf{e}^T\mathbf{e}}{\partial b_1} &= -2 \sum_{i=1}^N x_i\left[y_i - b_0 - b_1x_i\right] \notag \ &= -2 \left[\sum_{i=1}^N x_iy_i - b_0 \sum_{i=1}^N x_i - b_1 \sum_{i=1}^N x_i^2 \right] (#eq:derivssqresidualb1regr) \end{align}
Setting both expression in \@ref(eq:derivssqresidualb0regr) and \@ref(eq:derivssqresidualb1regr) to zero and writing them in terms of $\hat{b}_0$ and $\hat{b}_1$ gives
\begin{equation} N\hat{b}_0 + \hat{b}_1x. = y. (#eq:solderivssqresidualb0regr) \end{equation}
and
\begin{equation} \hat{b}_0x. + \hat{b}_1(x^2). = (xy). (#eq:solderivssqresidualb1regr) \end{equation}
using the dot notation for the following sums $x. = \sum_{i=1}^N x_i$, $y. = \sum_{i=1}^N y_i$, $(x^2). = \sum_{i=1}^N x_i^2$ and $(xy). = \sum_{i=1}^N x_iy_i$. With the bar notation for the means, we can further write
\begin{equation} \bar{x}. = {x. \over N} (#eq:barxdotregr) \end{equation}
and
\begin{equation} \bar{y}. = {y. \over N} (#eq:barydotregr) \end{equation}
The solutions in \@ref(eq:solderivssqresidualb0regr) and \@ref(eq:solderivssqresidualb1regr) can then be written as
\begin{equation} \hat{b}_0 = \bar{y}. - \hat{b}_1\bar{x}. (#eq:solderivssqresidualb0regr) \end{equation}
and
\begin{equation} \hat{b}_1 = \frac{(xy). - N\bar{x}.\bar{y}.}{(x^2). - N\bar{x}.^2} (#eq:solderivssqresidualb1regr) \end{equation}
Estimates for Example Dataset {#estimates-for-example-dataset-regr}
The estimates as shown in \@ref(eq:solderivssqresidualb0regr) and \@ref(eq:solderivssqresidualb1regr) are computed for our example dataset. We start by computing all the components of the formulas for the estimators, then we plug those components in and get the results.
# formula components x_dot <- sum(tbl_reg$`Breast Circumference`) x_bar <- x_dot / n_nr_obs y_dot <- sum(tbl_reg$`Body Weight`) y_bar <- y_dot / n_nr_obs xy_dot <- sum(tbl_reg$`Breast Circumference`*tbl_reg$`Body Weight`) x2_dot <- sum(tbl_reg$`Breast Circumference`^2) # estimators b1_hat <- (xy_dot - n_nr_obs * x_bar * y_bar) / (x2_dot - n_nr_obs * x_bar^2) b0_hat <- y_bar - b1_hat * x_bar
$$
N = r n_nr_obs,\ \bar{x}. = r x_bar,\ \bar{y}. = r y_bar,\ (xy). = r xy_dot,\ (x^2). = r x2_dot
$$
\begin{equation}
\hat{b}_1 = \frac{r xy_dot - r n_nr_obs * r x_bar * r y_bar}{r x2_dot - r n_nr_obs * r x_bar^2} = r round(b1_hat, digits=3)
(#eq:numsolderivssqresidualb1regr)
\end{equation}
\begin{equation}
\hat{b}_0 = r y_bar - r b1_hat * r x_bar = r round(b0_hat, digits=3)
(#eq:numsolderivssqresidualb0regr)
\end{equation}
Obtain Parameter Estimates in R {#obtain-parameter-estimates-in-r-regr}
As already mentioned, the same type of computation as shown in Section \@ref(estimates-for-example-dataset-regr) is also implemented in the R-function lm(). In what follows, we show how the estimates of the linear regression model are obtained using lm(). An important pre-requisite for using the function lm() is that the dataset is assigned to a dataframe. Here we assume that we have a dataframe named tbl_reg that contains our dataset. Then the parameter estimates are obtained using the following statements in R.
lm_bw_bc <- lm(`Body Weight` ~ `Breast Circumference`, data = tbl_reg) summary(lm_bw_bc)
The output of the summary function shows a lot of information about the data and the model. Under the section Coefficients there are two entries
(Intercept)corresponding to our $b_0$ andBreast Circumferencecorresponding to our $b_1$.
The values in the first column entitled Estimate correspond to the values that we have computed in the previous section.
The General Case {#the-general-case-regr}
Suppose that in the study on body weight and breast circumference, we have an additional observation for each animal consisting of the height of the animal. The new extended data set is shown in Table \@ref(tab:tab-bw-mult-reg).
tbl_bw_ext <- dplyr::bind_cols(tbl_reg, tibble::tibble(Height = c(161, 121, 157, 165, 136, 123, 163, 149, 143, 130))) knitr::kable(tbl_bw_ext, booktabs = TRUE, longtable = TRUE, caption = paste0("Extended Dataset of Body Weight for ", nrow(tbl_reg)," Animals", collapse = ""), escape = FALSE)
The model developed so far is not extended to be
\begin{equation} E(y) = b_0 + b_1x_1 + b_2x_2 (#eq:extendedbwmodelregr) \end{equation}
where $x_1$ represents the breast circumference in cm, and $x_2$ stands for the height of the animal in cm. Thus for animal $i$ with a breast circumference of $x_{i1}$ cm and a height of $x_{i2}$ the body weight $y_i$ can be written as
\begin{equation} y_i = b_0 + b_1x_{i1} + b_2x_{i2} + e_i (#eq:extendedbwdecompregr) \end{equation}
As the number of $x$ variables increase, the equations such as shown in \@ref(eq:extendedbwdecompregr) are getting longer and they are getting tedious to handle. This problem is solved by a change of notation. Let us define the matrix $\mathbf{X}$ and the vectors $\mathbf{y}$, $\mathbf{e}$ and $\mathbf{b}$ as follows.
mat_X <- matrix(c("x_{10}", "x_{11}", "x_{12}", "x_{20}", "x_{21}", "x_{22}", ".", ".", ".", ".", ".", ".", ".", ".", ".", "x_{N0}", "x_{N1}", "x_{N2}"), ncol = 3, byrow = TRUE) vec_y <- c("y_1", "y_2", ".", ".", ".", "y_N") vec_e <- c("e_1", "e_2", ".", ".", ".", "e_N") vec_b <- c("b_0", "b_1", "b_2") cat("$$\n") cat(paste0(rmdhelp::bmatrix(pmat = mat_X, ps_name = "\\mathbf{X}"), collapse = ""), "\n") cat(", \\ \n") cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_y, ps_name = "\\mathbf{y}"))) cat(", \\ \n") cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_e, ps_name = "\\mathbf{e}"))) cat("\\text{ and }\n") cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_b, ps_name = "\\mathbf{b}"))) cat("$$\n")
Because all equations contain the term $b_0$, the first column of matrix $\mathbf{X}$ consists of all ones. Hence, $x_{10} = x_{20} = \ldots = x_{N0} = 1$. The complete set of equations for all animals in our extended dataset represented by \@ref(eq:extendedbwdecompregr) is
\begin{equation} \mathbf{y} = \mathbf{X}\mathbf{b} + \mathbf{e} \text{, with } E(\mathbf{y}) = \mathbf{X}\mathbf{b} (#eq:matvecextendedbwdecompregr) \end{equation}
The big advantage of the matrix-vector notation is that equation \@ref(eq:matvecextendedbwdecompregr) is invariant to the number of $x$-variables. That means no matter how many $x$-variables we include in our dataset, the linear regression model can always be written as shown in equation \@ref(eq:matvecextendedbwdecompregr). The only thing that changes are the definitions of $\mathbf{X}$ and $\mathbf{b}$. In the general case with $k$ variables
mat_X_gen <- matrix(c("x_{10}", "x_{11}", ".", "x_{1k}", "x_{20}", "x_{21}", ".", "x_{2k}", ".", ".", ".", ".", ".", ".", ".", ".", ".", ".", ".", ".", "x_{N0}", "x_{N1}", ".", "x_{Nk}"), ncol = 4, byrow = TRUE) vec_b_gen <- c("b_0", "b_1", ".", ".", "b_k") cat("$$\n") cat(paste0(rmdhelp::bmatrix(pmat = mat_X_gen, ps_name = "\\mathbf{X}"), collapse = ""), "\n") cat(", \\ \n") cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_b_gen, ps_name = "\\mathbf{b}"))) cat("$$\n")
The topic of parameter estimation can also be shown using matrix-vector notation. The only restriction that we are currently imposing is that the number of $x$-variables $k$ is smaller than the number of observations $N$.
The model specification is only complete, if the properties of the vector $\mathbf{e}$ of random residuals is defined. In accordance with equations \@ref(eq:expvalueerrorregr), \@ref(eq:varerrorregr) and \@ref(eq:coverrorregr), we can write
\begin{equation} E(\mathbf{e}) = \mathbf{0} (#eq:matvecexpvalueerrorregr) \end{equation}
where $E(\mathbf{e})$ stands for the vector of length $N$ containing the expected values of all random residuals and $\mathbf{0}$ is a vector of $N$ zeros and
\begin{equation} var(\mathbf{e}) = E\left[\mathbf{e} - E(\mathbf{e}) \right]\left[\mathbf{e} - E(\mathbf{e}) \right]^T = E(\mathbf{e}\mathbf{e}^T) = \sigma^2 \mathbf{I}_N (#eq:matvecvarerrorregr) \end{equation}
where $var(\mathbf{e})$ is the variance-covariance matrix between all random residuals and $\mathbf{I}_N$ is the $N\times N$ identity matrix.
Derivation of the least squares estimator of $\mathbf{b}$ follows the same principles as shown in equations \@ref(eq:derivssqresidualb0regr) - \@ref(eq:solderivssqresidualb1regr). The sum of squares of the deviations of the observations from their expected values using $E(\mathbf{e}) = \mathbf{0}$ and hence $E(\mathbf{y}) = \mathbf{Xb}$, is
\begin{align}
\mathbf{e}^T\mathbf{e} &= \left[\mathbf{y} - E(\mathbf{y}) \right]^T\left[\mathbf{y} - E(\mathbf{y}) \right] \notag \
&= \left[\mathbf{y} - \mathbf{Xb} \right]^T\left[\mathbf{y} - \mathbf{Xb} \right]\notag \
&= \mathbf{y}^T\mathbf{y} - 2 \mathbf{b}^T\mathbf{X}^T\mathbf{y} + \mathbf{b}^T\mathbf{X}^T\mathbf{X}\mathbf{b}
(#eq:matvecsumsqresidualregr)
\end{align}
The least squares estimator $\hat{\mathbf{b}}$ is found by minimizing $\mathbf{e}^T\mathbf{e}$ with respect to all elements of $\mathbf{b}$. This is corresponds to $\partial{\mathbf{e}^T\mathbf{e}}/\partial \mathbf{b}$ and is also called the gradient of $\mathbf{e}^T\mathbf{e}$ with respect to $\mathbf{b}$. Setting the gradient to zero leads to the following equations
\begin{equation} \mathbf{X}^T\mathbf{X}\hat{\mathbf{b}} = \mathbf{X}^T\mathbf{y} (#eq:leastsquarenormalequationregr) \end{equation}
These equations are known as the least squares normal equations. Provided $(\mathbf{X}^T\mathbf{X})$ can be inverted, the unique solution for $\hat{\mathbf{b}}$ can be written as
\begin{equation} \hat{\mathbf{b}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y} (#eq:solleastsquareestregr) \end{equation}
Example Dataset Continued {#example-dataset-continued-regr}
We use the original dataset shown in Table \@ref(tab:tab-reg-bw-bc) to illustrate the estimation of the least squares parameters using the matrix-vector notation. The matrix $\mathbf{X}$ and the vectors $\mathbf{y}$, $\mathbf{b}$ and $\mathbf{e}$ are defined as follows.
nr_lines_shown <- 2 nr_lines_dot <- 3 mat_X_short <- matrix(c(rep(1, nr_lines_shown), rep('.', nr_lines_dot), 1, tbl_reg$`Breast Circumference`[1:nr_lines_shown], rep('.', nr_lines_dot), tbl_reg$`Breast Circumference`[n_nr_obs]), ncol = 2) vec_y_short <- c(tbl_reg$`Body Weight`[1:nr_lines_shown], rep('.', nr_lines_dot), tbl_reg$`Body Weight`[n_nr_obs]) vec_b_short <- c("b_0", "b_1") vec_e_short <- c(rmdhelp::vecGetVecElem(psBaseElement = "e", pnVecLen = nr_lines_shown, psResult = 'latex'), rep('.', nr_lines_dot), paste0("e_{", n_nr_obs, "}", collapse = "")) # show matrix X and vectors cat("$$\n") cat(paste0(rmdhelp::bmatrix(pmat = mat_X_short, ps_name = "\\mathbf{X}"), collapse = ""), "\n") cat("\\text{, }\n") cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_y_short, ps_name = "y"))) cat("\\text{, }\n") cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_b_short, ps_name = "b"))) cat("\\text{, }\n") cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_e_short, ps_name = "e"))) cat("$$\n")
The solution is obtained from equation \@ref(eq:solleastsquareestregr) for this we also need
mat_X_num <- matrix(c(rep(1,n_nr_obs), tbl_reg$`Breast Circumference`), ncol = 2) mat_XTX <- crossprod(mat_X_num) mat_XTXinv <- solve(mat_XTX) vec_y <- tbl_reg$`Body Weight` mat_XTy <- crossprod(mat_X_num, vec_y) vec_sol <- solve(mat_XTX, mat_XTy) # show the matrices cat("$$\n") cat(paste0(rmdhelp::bmatrix(pmat = mat_XTX, ps_name = "X^TX"), collapse = ""), "\n") cat("\\text{, }\n") cat(paste0(rmdhelp::bmatrix(pmat = round(mat_XTXinv, digits = 3), ps_name = "(X^TX)^{-1}"), collapse = ""), "\n") cat("\\text{, }\n") cat(paste0(rmdhelp::bmatrix(pmat = round(mat_XTy, digits = 3), ps_name = "X^Ty"), collapse = ""), "\n") cat("$$\n")
The solution vector $\hat{\mathbf{b}}$ contains the two components corresponding to the estimates of the two unknown parameters.
vec_sol_sym <- c("\\hat{b}_0", "\\hat{b}_1") # show the solutions cat("$$\n") cat(paste0(rmdhelp::bcolumn_vector(pvec = vec_sol_sym, ps_name = "\\hat{\\mathbf{b}}"), collapse = ""), "\n") cat(" = ") cat(paste0(rmdhelp::bmatrix(pmat = round(vec_sol, digits = 3)), collapse = ""), "\n") cat("$$\n")
Comparing the above shown solutions to the results received earlier shows that for this example, the solution using the matrix-vector notation are the same.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.