Показать сообщение отдельно
Старый 23.12.2012, 12:37   #21
Hogfather
Platinum Member
 
Аватар для Hogfather
 
Регистрация: 22.07.2010
Адрес: Санкт-Петербург
Сообщений: 3,298
По умолчанию gnu r регрессионный анализ. Ответы на вопросы.

Как выяснилось в переписке, вопрос касался также 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

Последний раз редактировалось Hogfather; 23.12.2012 в 21:08.
---------
DNF is not an option
Hogfather вне форума   Ответить с цитированием
Реклама