Как выяснилось в переписке,
вопрос касался также GNU R.
В R, на мой взгляд, это все делается проще.
Пусть у нас есть некие данные (данные немного другие, я не сохранил оригинал того примера), которые представлены в таблице (фрейме)
Код:
> head(MyD)
A B C Y
1 1 1 0.1 6.6
2 2 1 0.2 9.2
3 3 2 0.3 32.3
4 4 2 0.4 40.4
5 5 3 0.5 88.0
6 6 3 0.6 105.6
Тогда линейна модель МНК Y от всех без исключения параметров выглядит так:
Код:
> summary(mdl<-lm(Y~.,data=MyD))
Call:
lm(formula = Y ~ ., data = MyD)
Residuals:
Min 1Q Median 3Q Max
-56.588 -26.540 -0.223 20.860 79.223
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -21.69 52.98 -0.409 0.691
A -76.39 44.15 -1.730 0.114
B 64.29 77.05 0.834 0.424
C 548.60 761.95 0.720 0.488
Residual standard error: 46.28 on 10 degrees of freedom
Multiple R-squared: 0.987, Adjusted R-squared: 0.9831
F-statistic: 252.6 on 3 and 10 DF, p-value: 1.01e-09
Если нам нужны еще произведения параметров, то делаем вот так:
Код:
> summary(mdl<-lm(Y~(.)^2,data=MyD))
Call:
lm(formula = Y ~ (.)^2, data = MyD)
Residuals:
Min 1Q Median 3Q Max
-0.68607 -0.04158 -0.00208 0.11227 0.68607
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.6861 1.4453 6.010 0.000537 ***
A 6.2334 0.6792 9.177 3.76e-05 ***
B -4.1133 2.1816 -1.885 0.101355
C -121.0582 9.2645 -13.067 3.58e-06 ***
A:B 0.4813 0.5852 0.822 0.437927
A:C -2.2817 2.3619 -0.966 0.366206
B:C 79.8960 0.2492 320.614 7.58e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4427 on 7 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 1.398e+06 on 6 and 7 DF, p-value: < 2.2e-16
Знаки "*" и "^" имеют в описании модели сакральный смысл.
Результат, аналогичный вышеописанному можно также получить
Код:
> summary(mdl<-lm(Y~(.)*(.),data=MyD))
Теперь рассмотрим именно нашу задачу с созданием без создания фиктивных переменных.
Создание фиктивных переменных делается весьма просто
Код:
> MyD$A2=MyD$A^2
> MyD$AB=MyD$A*MyD$B
> MyD$AC=MyD$A*MyD$C
> MyD$BC=MyD$B*MyD$C
> MyD$B2=MyD$B^2
> MyD$C2=MyD$C^2
> MyD$ABC=MyD$A*MyD$B*MyD$C
> head(MyD)
A B C Y A2 AB AC BC B2 C2 ABC
1 1 1 0.1 6.6 1 1 0.1 0.1 1 0.01 0.1
2 2 1 0.2 9.2 4 2 0.4 0.2 1 0.04 0.4
3 3 2 0.3 32.3 9 6 0.9 0.6 4 0.09 1.8
4 4 2 0.4 40.4 16 8 1.6 0.8 4 0.16 3.2
5 5 3 0.5 88.0 25 15 2.5 1.5 9 0.25 7.5
6 6 3 0.6 105.6 36 18 3.6 1.8 9 0.36 10.8
Теперь строим модель, способом показанным выше
Код:
> summary(mdl<-lm(Y~.,data=MyD))
Call:
lm(formula = Y ~ ., data = MyD)
Residuals:
Min 1Q Median 3Q Max
-7.240e-14 -1.505e-14 0.000e+00 2.387e-14 6.347e-14
Coefficients: (3 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.533e-13 5.138e-13 -1.271e+00 0.251
A 1.725e+01 5.911e-13 2.918e+13 <2e-16 ***
B 2.500e+00 4.752e-13 5.261e+12 <2e-16 ***
C -1.740e+02 3.122e-12 -5.573e+13 <2e-16 ***
A2 -1.250e+00 6.484e-14 -1.928e+13 <2e-16 ***
AB -7.663e-14 8.777e-14 -8.730e-01 0.416
AC NA NA NA NA
BC 5.000e+01 1.619e-12 3.089e+13 <2e-16 ***
B2 NA NA NA NA
C2 NA NA NA NA
ABC 5.000e+00 2.707e-13 1.847e+13 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6.341e-14 on 6 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 5.843e+31 on 7 and 6 DF, p-value: < 2.2e-16
Чисто чтоб поржать, сделаем автоматический подбор наилучшей модели.
Код:
> (newmod<-step(mdl))
Start: AIC=-846.76
Y ~ A + B + C + A2 + AB + AC + BC + B2 + C2 + ABC
Step: AIC=-846.76
Y ~ A + B + C + A2 + AB + AC + BC + B2 + ABC
Step: AIC=-846.76
Y ~ A + B + C + A2 + AB + AC + BC + ABC
Step: AIC=-846.76
Y ~ A + B + C + A2 + AB + BC + ABC
Df Sum of Sq RSS AIC
- AB 1 0.0000 0.0000 -852.25
<none> 0.0000 -846.76
- B 1 0.1113 0.1113 -53.68
- ABC 1 1.3721 1.3721 -18.52
- A2 1 1.4944 1.4944 -17.32
- A 1 3.4249 3.4249 -5.71
- BC 1 3.8362 3.8362 -4.12
- C 1 12.4892 12.4892 12.40
Step: AIC=-852.25
Y ~ A + B + C + A2 + BC + ABC
Df Sum of Sq RSS AIC
<none> 0.0000 -852.25
- B 1 0.3263 0.3263 -40.63
- ABC 1 1.5047 1.5047 -19.23
- A2 1 1.5742 1.5742 -18.59
- A 1 3.5163 3.5163 -7.34
- BC 1 4.1983 4.1983 -4.86
- C 1 15.0305 15.0305 12.99
Call:
lm(formula = Y ~ A + B + C + A2 + BC + ABC, data = MyD)
Coefficients:
(Intercept) A B C A2 BC
-3.950e-13 1.725e+01 2.500e+00 -1.740e+02 -1.250e+00 5.000e+01
ABC
5.000e+00
Предупреждения
1: попытка выбора модели для в целом хорошей подгонки -- глупость
2: попытка выбора модели для в целом хорошей подгонки -- глупость
3: попытка выбора модели для в целом хорошей подгонки -- глупость
4: попытка выбора модели для в целом хорошей подгонки -- глупость
5: попытка выбора модели для в целом хорошей подгонки -- глупость
А теперь "на бис", то же самое, но без создания фиктивных переменных. Поскольку символы "^" и "*" имеют особый смысл, посмотрите как это рекомендуется делать (можно было начать с "." через + добавить недостающие параметры, но я сделал в явном виде)
Код:
> (summary(mdl<-lm(Y~A+B+C+A:B+B:C+A:C+A:B:C+I(A^2)+I(B^2)+I(C^2),data=MyD)))
Call:
lm(formula = Y ~ A + B + C + A:B + B:C + A:C + A:B:C + I(A^2) +
I(B^2) + I(C^2), data = MyD)
Residuals:
Min 1Q Median 3Q Max
-7.500e-14 -1.737e-15 0.000e+00 2.102e-15 6.022e-14
Coefficients: (3 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.194e-13 3.875e-13 1.340e+00 0.229
A 1.000e+00 1.183e-12 8.454e+11 <2e-16 ***
B -6.000e-13 4.321e-13 -1.389e+00 0.214
C 1.000e+00 1.052e-11 9.503e+10 <2e-16 ***
I(A^2) 8.061e-14 1.113e-13 7.240e-01 0.496
I(B^2) 5.000e+00 1.551e-13 3.225e+13 <2e-16 ***
I(C^2) -6.361e-12 6.620e-12 -9.610e-01 0.374
A:B NA NA NA NA
B:C NA NA NA NA
A:C NA NA NA NA
A:B:C 5.000e+00 2.042e-13 2.449e+13 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.783e-14 on 6 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 1.027e+32 on 7 and 6 DF, p-value: < 2.2e-16
Обратите внимание, получилось несколько иное решение, также с Эр Квадрат равное 1. А что Вы хотели, кто обещал, что будет легко?
Теперь о
тесте Рамсея RESET.
Естественно, что он есть в GNU R и находится в пакете lmtest.
Код:
> library(lmtest)
> x <- c(1:30)
> y1 <- 1 + x + x^2 + rnorm(30)
> y2 <- 1 + x + rnorm(30)
> resettest(y1 ~ x , power=2, type="regressor")
RESET test
data: y1 ~ x
RESET = 108553.5, df1 = 1, df2 = 27, p-value < 2.2e-16
> resettest(y2 ~ x , power=2, type="regressor")
RESET test
data: y2 ~ x
RESET = 2.976, df1 = 1, df2 = 27, p-value = 0.09594
Ну, а в нашем случае, это выглядит так
Код:
> mdl<-lm(Y~A+B+C,data=MyD)
> resettest(mdl)
RESET test
data: mdl
RESET = 204.7904, df1 = 2, df2 = 8, p-value = 1.347e-07
> mdl<-lm(Y~A+B+C+A:B+B:C+A:C+A:B:C+I(A^2)+I(B^2)+I(C^2),data=MyD)
> resettest(mdl)
RESET test
data: mdl
RESET = 0.925, df1 = 2, df2 = 1, p-value = 0.5923