Краткое введение в мир матричной алгебры

Линейные модели…

Вадим Хайтов, Анастасия Лянгузова, Наталья Арбузова, Валентина Домашкина, Юта Тамберг, Марина Варфоломеева

Кафедра Зоологии беспозвоночных, СПбГУ

Вы сможете

  • Объяснить что такое матрицы и какие бывают их основные разновидности
  • Выполнить базовые операции с матрицами с использованием функций R
  • Применить в среде R методы матричной алгебры для решения простейших задач

Зачем понадобилась эта лекция?

Решите задачу

В популяции некоторого организма в момент времени \(t_1\) представлено четыре возрастные группы:

Возрастная группа Количество особей
0+ 100
1+ 10
2+ 5
3+ 1

Какова будет численность возрастных групп в момент времени \(t_2\) если известно, что

  1. каждая особь возрастом 2+ оставляет 10 потомков, а каждая особь возрастом 3+ - 20 потомков.
  2. вероятность выживания особей при переходе от \(t_1\) к \(t_2\) описывается следующими величинами
Переход из младшей группы в старшую Вероятность выживания
0+ - 1+ 0.1
1+ - 2+ 0.5
2+ - 3+ 0.3

Простое, но нерациональное решение

Возрастная группа Количество особей в момент времени \(t_2\)
0+ \(5 \cdot 10 + 1 \cdot 20 = 70\)
1+ \(100 \cdot 0.1 = 10\)
2+ \(10 \cdot 0.5 = 5\)
3+ \(5 \cdot 0.3 = 1.5\)

Рациональное решение

Демографический вектор в момент времени \(t_1\)

t1 <- c(100, 10, 5, 1)
t1
[1] 100  10   5   1

Рациональное решение

Демографическая матрица Лесли

\[ \begin{pmatrix} F_1 & F_2 & F_3 & F_4 \\ P_{0+ -1+}& 0 & 0 & 0 \\ 0 & P_{1+-2+} & 0 & 0 \\ 0 & 0 & P_{2+-3+} & 0 \\ \end{pmatrix} \]

lesley <- matrix(c(c(0, 0, 10, 20), c(0.1, 0, 0, 0), c(0, 0.5, 0, 0), c(0, 0, 0.3, 0)), nrow = 4, byrow = T)

lesley
     [,1] [,2] [,3] [,4]
[1,]  0.0  0.0 10.0   20
[2,]  0.1  0.0  0.0    0
[3,]  0.0  0.5  0.0    0
[4,]  0.0  0.0  0.3    0

Рациональное решение - решение средствами матричной алгебры

lesley %*% t1
     [,1]
[1,] 70.0
[2,] 10.0
[3,]  5.0
[4,]  1.5

Зачем нужны матрицы?

Матричные объекты

  • Есть много типов объектов, для которых такое выражение оказывается наиболее естественным (изображения, описания многомерных объектов и т.д.)
  • В матрицах, как и в обычных числах, скрыта информация, которую можно извлекать и преобразовывать по определенным правилам

Структура матриц

\[\begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1c} \\ a_{21} & a_{22} & \cdots & a_{2c} \\ \vdots & \vdots & \ddots & \vdots \\ a_{r1} & a_{r2} & \cdots & a_{rc} \end{pmatrix} \]

Размер (порядок) матрицы \(r \times c\)

Разновидности матриц

Вектор-строка (row matrix)

\[ \textbf {A} = \begin{pmatrix} 1 & 2 & 3 \end{pmatrix} \]

Вектор-столбец (column matrix)

\[ \textbf {B} = \begin{pmatrix} 1 \\ 4 \\ 7 \\ 10 \end{pmatrix} \]

Разновидности матриц

Прямоугольные матрицы (rectangular matrices)

\[ \textbf {C} = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ 10 & 11 & 12 \end{pmatrix} \]

\[ \textbf {D} = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{pmatrix} \]

В таком виде обычно представляются исходные данные

Квадратные матрицы (square matrices)

Это наиболее “операбельные” матрицы

\[ \textbf {E} = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix} \]

Диагональные матрицы (diagonal matrix)

\[ \textbf {F} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 5 & 0 & 0 \\ 0 & 0 & 9 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \]

Квадратные матрицы (square matrices)

Единичная матрица (identity matrix)

\[ \textbf {I} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \]

Единичная матрица (обозначение \(\textbf{I}\)) занимают особое место в матричной алгебре.
Она выполняет ту же роль, которую выполняет единица в обычной алгебре.

Особенность квадратных матриц

Для квадратных матриц могут быть найдены (но не обязательно существуют) некоторые важные для матричной алгебры показатели: определитель, инверсия, собственные значения и собственные векторы

Задание

Создайте с помощью R следующую матрицу

\[ \begin{pmatrix} 1 & 2 & 3 & 4 \\ 5 & 6 & 7 & 8 \\ 9 & 10 & 11 & 12 \\ \end{pmatrix} \]

Решение

matrix(1:12, ncol = 3, byrow = TRUE)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9
[4,]   10   11   12

Задание

Создайте с помощью R единичную матрицу

\[ \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \]

Решение

diag(rep(1,5))
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    0    1    0    0
[4,]    0    0    0    1    0
[5,]    0    0    0    0    1

Операции с матрицами

Транспонирование матриц

     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12

Транспонированная матрица \(\textbf{A}\)

     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12

Сложение матриц

     [,1] [,2] [,3]
[1,]    5    9   13
[2,]    6   10   14
[3,]    7   11   15
[4,]    8   12   16
     [,1] [,2] [,3]
[1,]    2   10   18
[2,]    4   12   20
[3,]    6   14   22
[4,]    8   16   24

Но! Нельзя складывать матрицы разных размеров

Error in `A + B`:
! неподобные многомерные матрицы

Биологическое приложение

Предположим, что мы подсчитывали двумя разными методами крупных и мелких животных трех видов в одних и тех же пробах

Крупные животные

        Sp1 Sp2 Sp3
Sample1  11   7  12
Sample2   6  11  10
Sample3  11  14  10
Sample4  11   6  12
Sample5  12  14   9

Мелкие животные

        Sp1 Sp2 Sp3
Sample1  55  49  63
Sample2  51  41  58
Sample3  43  54  51
Sample4  49  47  49
Sample5  48  48  48

Биологическое приложение

Общее обилие

        Sp1 Sp2 Sp3
Sample1  66  56  75
Sample2  57  52  68
Sample3  54  68  61
Sample4  60  53  61
Sample5  60  62  57

Простое умножение

Простое умножение на число

     [,1] [,2] [,3]
[1,]    4   20   36
[2,]    8   24   40
[3,]   12   28   44
[4,]   16   32   48

Умножение на матрицы на вектор

Простое умножение матрицы на вектор возможно только если число элементов в векторе равно числу строк в матрице

A * c(10, 11, 12, 13)
     [,1] [,2] [,3]
[1,]   10   50   90
[2,]   22   66  110
[3,]   36   84  132
[4,]   52  104  156

Все элементы первой строки матрицы умножаются на первый элемент вектора, все элементы второй строки на второй элемент вектора и т.д.

Биологическое применение

Допустим, учет организмов в части описаний проходил не на всей выборке, а лишь в ее части.

Rpocessed_portion <- c(1, 1, 1/2, 1/3, 1/4)
Processed_Factor <- 1/Rpocessed_portion
Small * Processed_Factor
        Sp1 Sp2 Sp3
Sample1  55  49  63
Sample2  51  41  58
Sample3  86 108 102
Sample4 147 141 147
Sample5 192 192 192

Скалярное произведение векторов

Это уже настоящая матричная операция

Допустимо только для векторов одинаковой размерности

\[ \textbf{a} \bullet \textbf{b} = \begin{pmatrix} a_1 \\ a_2 \\ a_3 \\ a_4 \\ a_5 \\ a_6 \\ a_7 \end{pmatrix} \times \begin{pmatrix} b_1 & b_2 & b_3 & b_4 & b_5 & b_6 & b_7 \end{pmatrix} = a_1 \cdot b_1 + a_2 \cdot b_2 + \ldots + a_7 \cdot b_7 = x \]

Результат этой операции - число (скаляр)

Биологическое применение

Сколько особей родится в популяции, если мы знаем репродуктивные характеристики всех возрастных групп?

\[ \begin{pmatrix} N_1 \\ N_2 \\ N_3 \\ N_4 \\ N_5 \\ N_6 \\ N_7 \end{pmatrix} \times \begin{pmatrix} F_1 & F_2 & F_3 & F_4 & F_5 & F_6 & F_7 \end{pmatrix} \]

N <- c(20, 40, 32, 45, 80, 50, 10)
Fert <- c( 0,  0,   1,   2,   2,   0,   0)
t(N) %*% (Fert)
     [,1]
[1,]  282

Умножение матриц

Умножать можно только в том случае, если число строк одной из матриц равно числу столбцов у другой матрицы.

Матрица A

     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12

Матрица B

     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12

кзультата матричного умножения \(\textbf{A} \times \textbf{B}\)

A %*% B
     [,1] [,2] [,3] [,4]
[1,]  107  122  137  152
[2,]  122  140  158  176
[3,]  137  158  179  200
[4,]  152  176  200  224

У первой матрицы (A) число строк равно числу столбцов у второй матрицы (B)

Умножение матриц

Порядок имеет смысл!

Если поменять местами матрицы A и B, то изменится и результат

B %*% A
     [,1] [,2] [,3]
[1,]   30   70  110
[2,]   70  174  278
[3,]  110  278  446

Умножение матриц

Удалим из матрицы B последнюю строчку

C <- B[-3, ]
C
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8

Такое произведение возможно

C %*% A
     [,1] [,2] [,3]
[1,]   30   70  110
[2,]   70  174  278

Но невозможно такое произведение

A %*% C
Error in `A %*% C`:
! неподобные аргументы

Умножение матриц

Если матрица (A) не квадратная, то нельзя произвести такое умножение

A %*% A

Схема умножения матриц

M1 <- matrix(c(1, 2, 3, 4, 5, 6), byrow = T, ncol = 3)
M2 <-  matrix(7:12, byrow = T, ncol = 2)
M1 %*% M2  
     [,1] [,2]
[1,]   58   64
[2,]  139  154

\[ \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ \end{pmatrix} \times \begin{pmatrix} 7 & 8 \\ 9 & 10\\ 11 & 12\\ \end{pmatrix} = \begin{pmatrix} 7 \cdot 1 + 9 \cdot 2 + 11 \cdot 3 & 8 \cdot 1 + 10 \cdot 2 + 12 \cdot 3 \\ 7 \cdot 4 + 9 \cdot 5 + 11 \cdot 6 & 8 \cdot 4 + 10 \cdot 5 + 12 \cdot 6 \\ \end{pmatrix} = \\ = \begin{pmatrix} 58 & 64 \\ 139 & 154\\ \end{pmatrix} \]

Схема умножения матриц

\[ \begin{pmatrix} A & B \\ C & D \\ \end{pmatrix} \times \begin{pmatrix} E & F \\ G & H\\ \end{pmatrix} = \begin{pmatrix} (A \cdot E + B \cdot G) & (A \cdot F + B \cdot H ) \\ (C \cdot E + D \cdot G) & (C \cdot F + D \cdot H) \\ \end{pmatrix} \]

Биологическое применение

Простейший пример использования умножения матриц - построение модели динамики демографической структуры популяции Для вычислений необходим начальный демографический вектор и матрица Лесли

\[ \begin{pmatrix} F_1 & F_2 & F_3 & F_4 & F_5 & F_6 & F_7 \\ P_{1-2}& 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & P_{2-3} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & P_{3-4} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & P_{4-5} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & P_{5-6} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & P_{6-7} & 0 \end{pmatrix} \times \begin{pmatrix} N1_t \\ N2_t \\ N3_t \\ N4_t \\ N5_t \\ N6_t \\ N7_t \end{pmatrix} = \begin{pmatrix} N1_{t+1} \\ N2_{t+1} \\ N3_{t+1} \\ N4_{t+1} \\ N5_{t+1} \\ N6_{t+1} \\ N7_{t+1} \end{pmatrix} \]

Простейшая демографическая модель

Демографический вектор в момент времени \(t\)

    Age T1
1     0 20
2  1-10 40
3 11-20 32
4 21-35 45
5 36-45 80
6 46-55 50
7 56-65 10

Матрица Лесли

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]  0.0  0.0  1.0  2.0  2.0  0.0    0
[2,]  0.6  0.0  0.0  0.0  0.0  0.0    0
[3,]  0.0  0.7  0.0  0.0  0.0  0.0    0
[4,]  0.0  0.0  0.8  0.0  0.0  0.0    0
[5,]  0.0  0.0  0.0  0.7  0.0  0.0    0
[6,]  0.0  0.0  0.0  0.0  0.6  0.0    0
[7,]  0.0  0.0  0.0  0.0  0.0  0.2    0

Демографическая струкутра в момент времени \(t+1\)

Демографическая струкутра в момент времени \(t+1\)

Вычисление корреляций через произведение матриц

Коэффициент корреляции \[r_{x,y} = \frac{\sum(x_i-\bar{x})(y_i-\bar{y})} {\sqrt{\sum(x_i-\bar{x})^2}\sqrt{\sum(y_i-\bar{y})^2}} = \frac{cov_{x,y}} {\sigma_x \sigma_y}\] В матричной форме

\[\bf{Cor} = \frac{\bf{{X}'_{stand}}\bf{X_{stand}}}{N-1}\]

Используем известные нам данные по размеру головного мозга

Code
brain <- read.csv("data/IQ_brain.csv", header = TRUE)
br <- brain[complete.cases(brain), -1]
br <- as.matrix(br)
br_scaled <- scale(br) #Стандартизуем значения

Матричное произведение

t(br_scaled) %*% br_scaled / (nrow(br_scaled) - 1) 
              FSIQ     VIQ      PIQ   Weight  Height MRINACount
FSIQ        1.0000  0.9451  0.93443 -0.05148 -0.1184      0.334
VIQ         0.9451  1.0000  0.77602 -0.07609 -0.1190      0.300
PIQ         0.9344  0.7760  1.00000  0.00251 -0.0932      0.378
Weight     -0.0515 -0.0761  0.00251  1.00000  0.6996      0.513
Height     -0.1184 -0.1190 -0.09316  0.69961  1.0000      0.588
MRINACount  0.3337  0.3003  0.37778  0.51338  0.5884      1.000

Важные положения теории матриц

Если существует произведение матриц \(\textbf{BC}\), то не обязательно существует \(\textbf{CB}\)

B <- matrix(1:24, ncol = 4)
C <- matrix(1:12, ncol = 3)

B %*% C
     [,1] [,2] [,3]
[1,]  130  290  450
[2,]  140  316  492
[3,]  150  342  534
[4,]  160  368  576
[5,]  170  394  618
[6,]  180  420  660

НО!

C %*% B
Error in `C %*% B`:
! неподобные аргументы

Такое произведение невозможно

Всегда существует такое произведение матриц \(\textbf{CС'}\) и \(\textbf{C'С}\)

C %*% t(C)
     [,1] [,2] [,3] [,4]
[1,]  107  122  137  152
[2,]  122  140  158  176
[3,]  137  158  179  200
[4,]  152  176  200  224
t(C) %*% C
     [,1] [,2] [,3]
[1,]   30   70  110
[2,]   70  174  278
[3,]  110  278  446

Произведение матриц \(\textbf{BC}\) как правило не равно \(\textbf{CB}\)

B <- matrix(1:9, ncol = 3)
C <- matrix(11:19, ncol = 3)

B %*%  C
     [,1] [,2] [,3]
[1,]  150  186  222
[2,]  186  231  276
[3,]  222  276  330
C %*% B
     [,1] [,2] [,3]
[1,]   90  216  342
[2,]   96  231  366
[3,]  102  246  390

\([\textbf{BC}]' = \textbf{C}'\textbf{B}'\)

t(B %*% C)
     [,1] [,2] [,3]
[1,]  150  186  222
[2,]  186  231  276
[3,]  222  276  330
t(C) %*% t(B)
     [,1] [,2] [,3]
[1,]  150  186  222
[2,]  186  231  276
[3,]  222  276  330

Произведение \(\textbf{BB'}\) и \(\textbf{B'B}\) всегда дает симметричную матрицу

B %*% t(B)
     [,1] [,2] [,3]
[1,]   66   78   90
[2,]   78   93  108
[3,]   90  108  126
t(B) %*% B
     [,1] [,2] [,3]
[1,]   14   32   50
[2,]   32   77  122
[3,]   50  122  194

Определтель матрицы

Определитель матрицы — это некоторое число.

По значению этого числа можно определить есть ли у матрицы некоторые свойства (например, обратима ли матрица).

Определитель бывает только у квадратных матриц.

Определтель матрицы

Матрицы, имеющие определитель равный нулю, называются сингулярными матрицами.

B
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
det(B)
[1] 0



НО!

B[3,3] <- B[3,3] + 1
B
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6   10
det(B)
[1] -3

Обращение (инверсия) матриц

В матричной алгебре нет процедуры деления. Вместо нее используют обращение матриц.

\[ \textbf{X}^{-1}\textbf{X} = \textbf{I} \]

Произведение инверсии матрицы и исходной матрицы дает единичную матрицу

Обращение (инверсия) матриц

Только некоторые квадратные матрицы, могут иметь обратную матрицу.

Поэтому для квадратных матриц справедливо \[\textbf{X} \textbf{X}^{-1} = \textbf{X}^{-1} \textbf{X}\]

Обращение (инверсия) матриц

Для матриц, определитель которых равен нулю, обратной матрицы не существует

Обращение матриц в среде R

Создадим матрицу

X <- matrix(c(seq(1, 8),10), nrow = 3, byrow = T)
X
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8   10

Ее определитель

det(X)
[1] -3

Обращение матриц в среде R

Обратная матрица

solve(X)
       [,1]  [,2] [,3]
[1,] -0.667 -1.33    1
[2,] -0.667  3.67   -2
[3,]  1.000 -2.00    1

По определению, \(\textbf{X}^{-1}\textbf{X} = \textbf{I}\)

round(solve(X) %*% X )
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

Примнение обращенных матриц

Решение систем линейных уравнений

Простейший случай использования обратных матриц - решение систем линейных уравнений

\[ \begin{cases} 1x + 2y + 3z = 2\\ 4x + 5y + 6z = 4 \\ 7x + 8y + 10z = 10 \end{cases} \]

Эту систему можно представить в матричном виде

\[ \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 10 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} 2 \\ 4 \\ 10 \end{pmatrix} \]

Тогда

\[ \begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 10 \end{pmatrix} ^{-1} \begin{pmatrix} 2 \\ 4 \\ 10 \end{pmatrix} \]

Зачем нужен определитель

Если определитель матрицы коэффициентов не равен нулю, то для этой матрицы возможно обращение, а система уравнений имеет решение.

С помощью определителя матрицы коэффициентов можно определить будет ли система уравнений иметь решение.

Если определитель равен нулю, то инверсия матрицы невозможна и система уравнений имеет бесконечное количество решений.

Задание

Решите приведенную систему уравнений с использованием матричной алгебры

\[ \begin{cases} 1x + 2y + 3z = 2\\ 4x + 5y + 6z = 4 \\ 7x + 8y + 10z = 10 \end{cases} \]

Решение

Coef <- matrix(c(1 , 2 , 3 ,
         4 , 5 , 6 ,
         7 , 8 , 10), byrow = T, ncol = 3)
Val <- c(2,4,10)

solve(Coef) %*% Val
      [,1]
[1,]  3.33
[2,] -6.67
[3,]  4.00

Проверим

round(1*3.33 + 2*-6.67 + 3*4.00)
[1] 2
round(4*3.33 + 5*-6.67 + 6*4.00)
[1] 4
round(7*3.33 + 8*-6.67 + 10*4.00)
[1] 10

Задание

Решите систему из двух уравнений

\[ \left\{\begin{array}{c}2x-4y=-6\\ 2x+y=4\end{array}\right. \]

Подбор параметров линейной регрессии методом наименьших квадратов с использованием матричной алгебры

Линейная регрессия в матричном виде

Модель линейной регрессии \(y_i = b_0 + b_1x_i + e_i\) можно записать в матричном виде

\[ \mathbf{y} = \mathbf{X}{\pmb{\beta}} + \mathbf{e} \] Или

\[ \mathbf{\hat{y}} = \mathbf{X}{\pmb{\beta}} \]

Здесь

\(\textbf{y}\) - вектор наблюдаемых значений зависимой переменной

\(\mathbf{\hat{y}}\) - вектор предсказанных значений зависимой переменной

\(\textbf{X}\) - модельная матрица

\(\pmb{\beta}\) - вектор коэффициентов модели \(b_0 + b_1x\): \((b_0, b_1)\)

\(\mathbf{e}\) — вектор остатков.

Модельная матрица

Это матрица, описывающая “поведение” предикторов.

Для модели \(y_i = b_0 + b_1x_i + e_i\)

\[ \mathbf{X} = \begin{pmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots\\ 1 & x_n \\ \end{pmatrix} \]

Решение этого уравнения

Умножим обе части уравнения на транспонированную матрицу \(\textbf{X}'\)

\[ \textbf{X}' \textbf{y} = \textbf{X}'\textbf{X} \pmb{\beta}\]

Матрица \(\textbf{X}'\textbf{X}\) - это всегда квадратная матрица. Ее можно обратить.

Тогда \[ \pmb{\beta} = [\textbf{X}'\textbf{X}]^{-1}[\textbf{X}'\textbf{y}] \]

Подбираем коэффициенты с помощью фунции lm()

brain <- read.csv("data/IQ_brain.csv", header = TRUE)

M_brain <- lm(PIQ ~ MRINACount, data = brain)

Коэффициенты модели

coefficients(M_brain)
(Intercept)  MRINACount 
    1.74376     0.00012 

Графическое отражение, построенное с помощью geom_smooth()

library(ggplot2)
theme_set(theme_bw())

ggplot(brain, aes(x = MRINACount, y = PIQ)) + geom_point() + geom_smooth(method = "lm")

Мы научимся все эти результаты получать вручную!

Вычисление коэффициентов линейной регрессии вручную

Находим вектор коэффициентов на основе матичного уравнения \[\pmb{\beta} = [\textbf{X}'\textbf{X}]^{-1}[\textbf{X}'\textbf{y}]\]

X <- model.matrix(~MRINACount, data = brain)
X <- as.matrix(data.frame(Intercept = 1, MRINACount = brain$MRINACount)) # Аналогичная запись 

Y <- brain$PIQ
betas <- solve(t(X) %*% X) %*% (t(X) %*% Y)
betas
              [,1]
Intercept  1.74376
MRINACount 0.00012

Применение матричной алгебры для построения графиков регрессионных моделей

Шаг 1. Формируем искусственный датасет со всеми возможными значениями предиктора

MyData <- data.frame(MRINACount = seq(min(brain$MRINACount), max(brain$MRINACount), length.out = 100))

head(MyData)
  MRINACount
1     790619
2     793537
3     796456
4     799374
5     802293
6     805211

Применение матричной алгебры для построения графиков регрессионных моделей

Шаг 2. Формируем модельную матрицу для искусственно созданных данных

X2 <- model.matrix( ~ MRINACount, data = MyData)

head(X2)
  (Intercept) MRINACount
1           1     790619
2           1     793537
3           1     796456
4           1     799374
5           1     802293
6           1     805211

Применение матричной алгебры для построения графиков регрессионных моделей

Шаг 3. Вычисляем предсказанные значения для искусственно созданных данных

MyData$Predicted <- X2 %*% betas

head(MyData)
  MRINACount Predicted
1     790619      96.8
2     793537      97.2
3     796456      97.5
4     799374      97.9
5     802293      98.2
6     805211      98.6

Применение матричной алгебры для построения графиков регрессионных моделей

Шаг 4. Рисуем график

library(ggplot2)
Pl <-
ggplot(MyData, aes(x = MRINACount, y = Predicted)) + 
  geom_line(color = "blue", size = 1.5) + 
  geom_point(data = brain, aes(x = MRINACount, y = PIQ))

Pl

Осталось нанести границы доверительной области. Но для этого необходимо ввести понятие вариационно-ковариационной матрицы

Вычисление вариационно-ковариационной матрицы

Подобранные параметры (вектор коэффициентов) – это лишь оценки некоторых параметров, описывающих связь между зависимой переменной и предиктором в генеральной совокупности.

Варьирование параметров описывает вриацонно-ковариационная матрица.

В среде R, если модель задана, например, с помощью функции lm(), эта матрица вычисляется так

vcov(M_brain)
            (Intercept)     MRINACount
(Intercept)  1797.11410 -0.00196543266
MRINACount     -0.00197  0.00000000216

Вычисление вариационно-ковариационной матрицы вручную

\[ \textbf{V}_b = s^2[\textbf{X}'\textbf{X}]^{-1} \] где

\[ s^2 = \frac{\sum\limits_{\substack{i=1}}^n e_i^2}{n-k} \]

\(\sum\limits_{\substack{i=1}}^n e_i^2\) - сумма квадратов остатков

\(n\) - объем

\(k\) - количество параметров модели

Вычисление вариационно-ковариационной матрицы вручную

Вычислем вектор остатков

RES <- brain$PIQ - X %*% betas

Вычисляем \(s^2 = \frac{\sum\limits_{\substack{i=1}}^n e_i^2}{n-k}\)

S2 <- sum(RES^2)/(nrow(brain) - length(betas))

Вычисляем \(\textbf{V}_b = s^2[\textbf{X}'\textbf{X}]^{-1}\)

VarCov <- S2 * solve(t(X) %*% X)
            Intercept     MRINACount
Intercept  1797.11410 -0.00196543266
MRINACount   -0.00197  0.00000000216

Применяем вариационно-ковариационую матрицу для вычисления доверительных интервалов

Для построения доверительного интервала надо вычислить значение стандартных ошибок для всех возможных значений предиктора.

Мы можем взять любую новую модельную матрицу \(\textbf{X}\), связанную с этой моделью.

Тогда \[ SE = \sqrt{diag(\textbf{X} \times \textbf{V} \times \textbf{X}')} \]

MyData$se <- sqrt(diag(X2 %*% VarCov  %*% t(X2)))

MyData$CiUp  <- MyData$Predicted + 1.96 *MyData$se
MyData$CiLow  <- MyData$Predicted - 1.96 *MyData$se

Откладываем на графике доверительные интервалы

Pl +
  geom_line(data = MyData, aes(y = CiLow), linetype = 2) +
  geom_line(data = MyData, aes(y = CiUp), linetype = 2)