Практическая работа №4

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

Глушков Егор Александрович, гр. 20.М04-мм
Вариант 4

---

1. Промоделировать нелинейную модель \(y = f(x, a, b) + \delta\) с несмещённой нормально распределенной ошибкой, дисперсия которой равна \(\varepsilon\), считая, \(x\) стандартно нормально распределённой случайной величиной.

• \(f(x, a, b) = a \sin x + bx^2, \; a=2, b=1, \varepsilon=0.8\)

2. Оценить параметры нелинейной модели по методу наименьших квадратов (численно). Применить к модельным данным линейную модель и оценить параметры. Построить на двумерной диаграмме основную и линейную модель. Сравнить невязки для обеих моделей.

3. Для линейной модели выполнить дисперсионный анализ, проверить значимость прогноза и коэффициентов регрессии. Сравнить непосредственные вычисления с результатами встроенной функции.

4. Промоделировать данные для множественной регрессии. Применить функцию \(lm\). Ответить на вопросы о значимости коэффициента детерминации, частных коэффициентов регрессии, о коэффициенте корреляции между остатком и независимыми переменными.

---

Задание линейной и нелинейной моделей регрессии

Нелинейная модель \[y_i = a \sin x_i + bx_i^2 + \delta_i, \quad \delta_i \sim \mathcal{N}(0, \sigma)\] Модель одномерной линейной регрессии

\[y_i = \alpha + \beta x_i, \quad \delta_i \sim \mathcal{N}(0, \sigma)\]

Задаем нелинейную и линейную модели и их функции ошибки L

N <- 100
f <- function(x, ab) ab[1] * sin(x) + ab[2] * x^2
L <- function(X, Y, ab) sum((Y - f(X, ab))^2)

f_lin <- function(x, ab_lin) ab_lin[1] + ab_lin[2] * x
L_lin <- function(X, Y, ab_lin) sum((Y - f_lin(X, ab_lin))^2)

Моделируем данные нелинейной модели

ab <- c(2, 1)
eps <- 0.8

X <- rnorm(N)
Y <- f(X, ab) + rnorm(N, 0, eps^2)

Оценка параметров линейной модели и наилучший прогноз

\[ \hat{\beta} = \frac{\sum\nolimits_i x_i y_i - n \bar{x} \bar{y}}{{\sum\nolimits_i x_i^2 - n \bar{x}^2}}, \quad \hat{\alpha} = \bar{y} - \hat{\beta} \bar{x}, \quad \hat{y_i} = \hat{\alpha} - \hat{\beta} x_i\]

b. <- (sum(X * Y) - N * mean(X) * mean(Y)) / (sum(X * X) - N * mean(X)^2)
a. <- mean(Y) - b. * mean(X)
ab_lin_est <- c(a., b.)
Y. <- f_lin(X, ab_lin_est)
ab_lin_est

[1] 0.9219083 1.2337556

Дисперсионный анализ

Источники вариации: общий \(Q_T\), обусловленный регрессией \(Q_R\), невязка \(Q_E\). Коэффициент детерминации \(R^2\)

QT <- sum((Y - mean(Y))^2)
QR <- sum((Y. - mean(Y))^2)
QE <- sum((Y - Y.)^2)
R2 <- QR / QT

c(QT=QT, QR=QR, QE=QE, R2=R2)

         QT          QR          QE          R2 
348.7500256 135.4274537 213.3225718   0.3883224 

38.8% дисперсии объясняется данной моделью – линией регресии.
Проверим равенство источников вариации [посчитаны верно]:

print(paste("(QR + QE == QT):", QR + QE - QT < 1e-9))

[1] "(QR + QE == QT): TRUE"

Дисперсии общая и коэффициентов регрессии

xx <- sum((X - mean(X))^2)
S2 <- QE / (N - 2)
S2.a <- S2 * sum(X^2) / (xx * N)
S2.b <- S2 / xx

c(S2, S2.a, S2.b)

[1] 2.17676094 0.02176906 0.02446596

Статистики для проверки значимости прогноза и коэффициентов регрессии

F_stat <- QR / QE * (N - 2)

Pf <- 1 - pf(F_stat, 1, N - 2)
T.a <- ab_lin_est[1] / sqrt(S2.a)
T.b <- ab_lin_est[2] / sqrt(S2.b)

P.a <- 2*(1 - pt(abs(T.a), N - 2))
P.b <- 2*(1 - pt(abs(T.b), N - 2))

P.model <- 1 - pf(F_stat, 1, N - 2)

Выведем самостоятельно посчитанные величины:

show_stats_coef <- rbind(c(ab_lin_est[1], sqrt(S2.a), T.a, P.a), c(ab_lin_est[2], sqrt(S2.b), T.b, P.b))
colnames(show_stats_coef) <- c("Estimate", "Std.Error","t value", "Pr(>|t|)")
rownames(show_stats_coef) <- c("(Intercept)", "X"); show_stats_coef

             Estimate Std.Error  t value     Pr(>|t|)
(Intercept) 0.9219083 0.1475434 6.248386 1.075902e-08
X           1.2337556 0.1564160 7.887656 4.401812e-12

show_stats_model <- c(sd=sqrt(S2), R2=R2, F_stats=F_stat, p.value.F=P.model); show_stats_model

          sd           R2      F_stats    p.value.F 
1.475385e+00 3.883224e-01 6.221512e+01 4.401812e-12 

Проверим при помощи встроенной функции lm:

summary(lm(Y~X))

Call:
lm(formula = Y ~ X)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.1924 -0.8860 -0.3422  0.5329  6.1879 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.9219     0.1475   6.248 1.08e-08 ***
X             1.2338     0.1564   7.888 4.40e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.475 on 98 degrees of freedom
Multiple R-squared:  0.3883,    Adjusted R-squared:  0.3821 
F-statistic: 62.22 on 1 and 98 DF,  p-value: 4.402e-12

Как видно, все посчитанные величины совпадают. Вычисления были верны. Что касается самой модели, то итоговое p-value (близкое к нулю) говорит о значимости регрессии, также малые Pr говорят о том, что соответствующие этим уровням значимости коэффициенты отличны от нуля и вносят вклад в модель. \(R^2\) сообщает, что 38.8% дисперсии объясняется данной моделью – линией регресии.

Оценка параметров нелинейной модели по МНК

Параметры нелинейной модели оцениваем методом наименьших квадратов (МНК)

NLM <- nlm(function(ab) L(X, Y, ab), ab)
ab. <- NLM$estimate
rbind(ab=ab, ab.=ab.)

        [,1]      [,2]
ab  2.000000 1.0000000
ab. 1.864256 0.9972268

plot(X, Y)
f_ <- function(x) f(x, ab); curve(f_, -3, 3, add=TRUE, col=2)

f_ <- function(x) f(x, ab.); curve(f_, -3, 3, add=TRUE, col=3)
f_ <- function(x) f_lin(x, ab_lin_est); curve(f_, -3, 3, add=TRUE, col=4, lty=2)

legend('bottomright', c('hyp', 'mnk', 'linear'), pch=20, col=2:4)

Невязки

c(Q.linear=L_lin(X, Y, ab_lin_est), Q.hyp=L(X, Y, ab), Q.mnk=L(X, Y, ab.))

 Q.linear     Q.hyp     Q.mnk 
213.32257  42.52355  41.80144 

Заметим, что линейная модель оценивает достаточно грубо, тогда как нелинейная модель достаточно точно описывает данные. Конечно, меньшей невязкой как суммы квадратов отклонений точек от графика обладает модель, построенная по МНК.

Множественная регрессия

N <- 100
a <- c(2, -1, 5)
eps <- 3
X1 <- rnorm(N, 1, 1.5)
X2 <- rnorm(N, 2, 0.5)
Y <- a[1] * X1 + a[2] * X2 + a[3] + rnorm(N, 0, eps^2)

SLM <- summary(lm(Y ~ X1 + X2)); SLM

Call:
lm(formula = Y ~ X1 + X2)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.1076  -6.0017  -0.6972   4.7822  21.5682 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   6.9812     3.3939   2.057  0.04237 * 
X1            1.5445     0.5706   2.707  0.00802 **
X2           -1.3288     1.6065  -0.827  0.41018   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.294 on 97 degrees of freedom
Multiple R-squared:  0.08022,   Adjusted R-squared:  0.06125 
F-statistic:  4.23 on 2 and 97 DF,  p-value: 0.01733

Согласно p-value < 0.05 отвергаем гипотезу о том, что данная модель не зависит от предикторов (т.е. их коэффициенты равны нулю – отвергаем). Отвергается и незначимость коэффициента при X1, а вот X2, как следует, не вносит весомый вклад (\(p.value = 0.41 > 0.05\), то не отвергается факт, коэффициент незначим)

Проверим некоррелированность остатков [практически равны 0]

c(cor(X1, SLM$residuals), cor(X1, SLM$residuals))

[1] -5.386815e-17 -5.386815e-17

op <- par(mfrow = c(2, 2))

plot(X1, SLM$residuals)
title(sub=paste("r", round(cor(SLM$residuals, X1), 3), sep="="))

plot(X2, SLM$residuals)
title(sub=paste("r", round(cor(SLM$residuals, X2), 3), sep="="))

plot(X1, Y)
title(sub=paste("r", round(cor(X1, Y),3), sep="="))

plot(X2, Y)
title(sub=paste("r", round(cor(X2, Y), 3), sep="="))

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

sh_t <- shapiro.test(SLM$residuals); sh_t

    Shapiro-Wilk normality test

data:  SLM$residuals
W = 0.97998, p-value = 0.1325

Согласно критерию Шапиро-Уилка, нормальность остатков не отвергается с \(p.value = 0.133\)

qqnorm(SLM$residuals)