GNU R: Вопросы и ответы - Страница 2

Hogfather · 02.12.2012, 16:14

Попробуем резюмировать.
Для подгонки модели МНК, используем код

Код:

> fit <- lm(lnZ ~ x + y, data=myData)
> summary(fit)

Call:
lm(formula = lnZ ~ x + y, data = myData)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.61700 -0.30351 -0.03218  0.19850  0.81873 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.5939071  0.3067107  -1.936   0.0678 .  
x            0.0016541  0.0012529   1.320   0.2025    
y            0.0029021  0.0002526  11.488  5.4e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.4538 on 19 degrees of freedom
Multiple R-squared: 0.889,	Adjusted R-squared: 0.8773 
F-statistic: 76.08 on 2 and 19 DF,  p-value: 8.53e-10

Также можно вызвать следующие функции

Код:

> coefficients(fit) # коэффициенты модели
 (Intercept)            x            y 
-0.593907127  0.001654092  0.002902143 
> confint(fit, level=0.95) # доверительные интервалы для параметров модели 
                    2.5 %      97.5 %
(Intercept) -1.2358599183 0.048045664
x           -0.0009683542 0.004276538
y            0.0023733723 0.003430914
> fitted(fit) [COLOR="rgb(46, 139, 87)"]# вычисленные значения[/COLOR]
        1         2         3         4         5         6         7         8 
0.6184535 1.3439892 2.4322929 3.8833645 1.0225736 0.7011581 2.8777654 1.8308140 
        9        10        11        12        13        14        15        16 
3.2818856 2.0898271 4.0487737 0.8252150 2.6390544 3.0431746 1.9962232 3.4472948 
       17        18        19        20        21        22 
4.1728306 3.8514150 0.9906242 2.2965886 4.2555352 3.5299994 
> residuals(fit) # остатки модели
          1           2           3           4           5           6           7 
-0.24551154 -0.03430598  0.50459127  0.03730105 -0.35167208  0.11950244  0.81873492 
          8           9          10          11          12          13          14 
-0.61649494  0.58275224 -0.59378185 -0.03005108  0.13451833 -0.43331068 -0.31506039 
         15          16          17          18          19          20          21 
 0.02651567 -0.14198445 -0.07559024 -0.26886877  0.21982761  0.81652665 -0.61700212 
         22 
 0.46336392 
> anova(fit) [COLOR="rgb(46, 139, 87)"]# дисперсионный анализ[/COLOR]
Analysis of Variance Table

Response: lnZ
          Df  Sum Sq Mean Sq F value    Pr(>F)    
x          1  4.1589  4.1589  20.192 0.0002488 ***
y          1 27.1795 27.1795 131.963 5.395e-10 ***
Residuals 19  3.9133  0.2060                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> vcov(fit) # ковариационная матрица для параметров модели
              (Intercept)             x             y
(Intercept)  9.407143e-02 -2.162456e-04 -4.568268e-05
x           -2.162456e-04  1.569874e-06 -8.310037e-08
y           -4.568268e-05 -8.310037e-08  6.382438e-08
> influence(fit) # анализ регрессионной модели
$hat
         1          2          3          4          5          6          7 
0.19665640 0.14809675 0.14788589 0.28317818 0.13829845 0.15275498 0.09302145 
         8          9         10         11         12         13         14 
0.06408668 0.09164807 0.05113364 0.14270474 0.15836003 0.04654710 0.05902938 
        15         16         17         18         19         20         21 
0.09846547 0.08567971 0.13738973 0.12649809 0.29922020 0.15799923 0.18148452 
        22 
0.13986129 

$coefficients
    (Intercept)             x             y
1  -0.090100375  1.924032e-04  4.151880e-05
2  -0.009639360  2.941447e-05  2.351114e-06
3   0.092492681 -5.221326e-04  3.423996e-05
4   0.002356970 -5.638034e-05  1.107148e-05
5  -0.098292644  1.997504e-04  4.375226e-05
6   0.034179275 -3.504510e-05 -2.200750e-05
7   0.068581503 -4.974500e-04  6.895179e-05
8  -0.087542609  1.378093e-04  3.287562e-05
9   0.014114197 -2.636425e-04  6.738312e-05
10 -0.056330941  3.061248e-05  2.300041e-05
11  0.002092611  1.126132e-05 -6.043748e-06
12  0.026144651  5.165601e-05 -2.977428e-05
13 -0.011340740 -3.227886e-05 -3.357501e-06
14  0.009716284 -7.069650e-05 -1.206585e-05
15  0.000820820  1.626464e-05 -2.654591e-06
16  0.012887806 -5.454754e-05 -1.004487e-05
17  0.012131630 -2.194206e-05 -1.245702e-05
18  0.042158631 -1.512501e-04 -2.872849e-05
19  0.018378098  3.404810e-04 -7.109334e-05
20 -0.039976009  8.764986e-04 -8.455088e-05
21  0.143930409 -4.760290e-04 -9.195059e-05
22 -0.072988620  3.945329e-04  2.397855e-05

$sigma
        1         2         3         4         5         6         7         8 
0.4617764 0.4661857 0.4481136 0.4661523 0.4576378 0.4652627 0.4199357 0.4414127 
        9        10        11        12        13        14        15        16 
0.4434361 0.4435794 0.4662052 0.4649853 0.4543848 0.4599405 0.4662215 0.4649526 
       17        18        19        20        21        22 
0.4658732 0.4613113 0.4621416 0.4164320 0.4376836 0.4511520 

$wt.res
          1           2           3           4           5           6           7 
-0.24551154 -0.03430598  0.50459127  0.03730105 -0.35167208  0.11950244  0.81873492 
          8           9          10          11          12          13          14 
-0.61649494  0.58275224 -0.59378185 -0.03005108  0.13451833 -0.43331068 -0.31506039 
         15          16          17          18          19          20          21 
 0.02651567 -0.14198445 -0.07559024 -0.26886877  0.21982761  0.81652665 -0.61700212 
         22 
 0.46336392

Сравнение двух моделей

Код:

> fit1 <- lm(lnZ ~ y, data=myData)
> anova(fit,fit1)
Analysis of Variance Table

Model 1: lnZ ~ x + y
Model 2: lnZ ~ y
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     19 3.9133                           
2     20 4.2723 -1  -0.35896 1.7428 0.2025

Остановлюсь на этом месте поподробнее. Когда мы добавляем или убавляем регрессоры, мы должны быть уверены, что принципиально это на дисперсию остатков не влияет. В данном случае используется нулевая гипотеза, что коэффициенты регрессоров равны между собой. Примерно то же самое мы делаем, когда проверяем саму модель на значимость(см. ниже), но в таком случаем сравниваем её с уравнением, где все регрессоры имеют коэффициенты 0. Аналогично описанному ниже, если Pr(>F) превышает 0.05, то гипотеза принимается. Т.е. в приведенном выше примере Pr(>F)=0.2025, отсюда разница между Model 1: lnZ ~ x + y и Model 2: lnZ ~ y является статистически незначимой.

Автоматическая подгонка модели посредством исключения незначимых членов.

Код:

> library(MASS)
> step <- stepAIC(fit, direction="both")
Start:  AIC=-31.99
lnZ ~ x + y

       Df Sum of Sq     RSS     AIC
- x     1     0.359  4.2723 -32.056
<none>               3.9133 -31.987
- y     1    27.180 31.0928  11.611

Step:  AIC=-32.06
lnZ ~ y

       Df Sum of Sq    RSS     AIC
<none>               4.272 -32.056
+ x     1     0.359  3.913 -31.987
- y     1    30.979 35.252  12.372
> step$anova # отображаем результаты
Stepwise Model Path 
Analysis of Deviance Table

Initial Model:
lnZ ~ x + y

Final Model:
lnZ ~ y


  Step Df  Deviance Resid. Df Resid. Dev       AIC
1                          19   3.913305 -31.98653
2  - x  1 0.3589588        20   4.272264 -32.05577

Прочую красоту смотреть вот здесь. А подробно проверка регрессионной модели описана вот тут.

Hogfather · 04.12.2012, 15:31

Хотя меня особо не спрашивают, но имея опыт работы со студентами постараюсь кратко и без лишнего занудства изложить, что же Пытливый Ум может увидеть в сборище букв на языке эвентуального противника, когда строит модель.

Рассмотрим пример, который был приведен выше.

Код:

> summary(fit)

Call:
lm(formula = lnZ ~ x + y, data = myData)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.61700 -0.30351 -0.03218  0.19850  0.81873 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.5939071  0.3067107  -1.936   0.0678 .  
x            0.0016541  0.0012529   1.320   0.2025    
y            0.0029021  0.0002526  11.488  5.4e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.4538 on 19 degrees of freedom
Multiple R-squared: 0.889,	Adjusted R-squared: 0.8773 
F-statistic: 76.08 on 2 and 19 DF,  p-value: 8.53e-10

Итак, первое на что мы смотрим, это не на разные Эр квадрат и прочие коэффициенты, нас интересует только вот эта строчка:
F-statistic: 76.08 on 2 and 19 DF, p-value: 8.53e-10

В данном случае нулевая гипотеза выглядит так: H₀: β₁=β₂=...β_n=0, т.е. все коэффициенты кроме свободного члена равны 0.

Это мы проводим проверку значимости модели. Если p-value у нас превышает 0.05, то дальнейшая работа с моделью просто бессмысленна, поскольку модель незначимая.

Пример незначимой модели

Код:

> lm0<-lm(y~x,data=test)
> summary(lm0)

Call:
lm(formula = y ~ x, data = test)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3245 -1.4632  0.6214  1.3782  2.6417 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   3.4424     0.8937   3.852  0.00117 **
x             0.1102     0.0746   1.478  0.15681   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.924 on 18 degrees of freedom
Multiple R-squared: 0.1082,	Adjusted R-squared: 0.05862 
F-statistic: 2.183 on 1 and 18 DF,  p-value: 0.1568

Второе, на что мы смотрим, это на остатки модели (Residuals).

Код:

Residuals:
     Min       1Q   Median       3Q      Max 
-0.61700 -0.30351 -0.03218  0.19850  0.81873

Алгоритм МНК подразумевает, что матожидание остатков должно быть равно 0 и мы надеемся, что остатки имеют нормальное распределение. Напоминаю, что у нормального распределения матожидание и медиана совпадают. В данном случае -- не очень. Поскольку мы видим отрицательную медиану (-0.03218), то распределение остатков незначительно смещено влево. Квартили, по идее, должны быть симметричны относительно 0, здесь же они значительно смещены влево. Максимум и минимум позволяют оценить выбросы (outliers) модели.

Третье, на что мы смотрим -- коэффициенты.

Код:

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.5939071  0.3067107  -1.936   0.0678 .  
x            0.0016541  0.0012529   1.320   0.2025    
y            0.0029021  0.0002526  11.488  5.4e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Для каждого коэффициента регрессора мы проверяем его значимость, вычисляя соответствующую t-статистику и проверяя нулевую гипотезу о незначимости его коэффициента. В дальнейшем мы будем по очереди убирать регрессоры исходя из наибольшего значения в последней колонке. Впрочем, об этом уже писалось.

Дальше, Residual standard error: 1.924 on 18 degrees of freedom говорит нам о средней квадратической ошибке остатков (σ).

Multiple R-squared: 0.889, Adjusted R-squared: 0.8773 говорит о коэффициентах детерминации и качестве модели. Чем они ближе к единице, тем лучше.

В то же время, проверку адекватности модели следует произвести графически. Именно для этого все графики и строились..

Пример неадекватной модели

Код:

> test<-data.frame(x<-1:20)
> test$y<-sqrt(test$x)
> lm0<-lm(y~x,data=test)
> summary(lm0)

Call:
lm(formula = y ~ x, data = test)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46467 -0.08467  0.04310  0.13672  0.17227 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.294290   0.081968   15.79 5.44e-12 ***
x           0.170382   0.006843   24.90 2.13e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.1765 on 18 degrees of freedom
Multiple R-squared: 0.9718,	Adjusted R-squared: 0.9702 
F-statistic:   620 on 1 and 18 DF,  p-value: 2.127e-15

Обратите внимание на R²
А теперь смотрите на график

Что должно быть на графике:
1. Точки на графике Residuals vs Fitted должны быть случайно разбросаны, без явно выраженного тренда. Тут этого нет.
2. Точки на графике Normal Q-Q должны более-менее лежать на линии, показывая что остатки имеют нормальное распредение
3. Точки на двух остальных графиках должны быть в группах не слишком далеко от центра.

Т.е. модель у нас неадекватна, не смотря на милый Эр Квадрат. Такие дела.
Для сравнения посмотрите график на предыдущей странице

Добавлено через 4 часа 3 минуты

Использование метода Бокса-Кокса для преобразования регрессионной модели

В предыдущем примере мы рассмотрели случай, когда модель выбрана неадекватной. Судя по параболе, которую мы наблюдаем, здесь необходимо какое-то степенное преобразование.
Попробуем использовать преобразование Бокса-Кокса, которое входит в пакет MASS. Т.е. найдем такую λ, чтобы y^λ сделал нам аж боже мой.

Исходные данные у нас те же и та же модель.

Код:

> test<-data.frame(x<-1:20)
> test$y<-sqrt(test$x)
> lm0<-lm(y~x,data=test)

А дальше начинается колдунство. В данном случае lambda взята в интервале 0 до 4 с шагом 0.01 ( по умолчанию -2 до 2 с шагом 0.1)

Код:

> library(MASS)
> bc<-boxcox(lm0,lambda = seq(0,4,0.01))

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

К сожалению, конкретного результата, который нужно подставить у нас нигде нет. Но мы можем его получить. Функция which.max(bc$y) вернёт нам номер элемента, у которого Игрек наибольший, а bc$x[which.max(bc$y)]) -- собственно нужную лямбду.

Код:

> (lambda<-bc$x[which.max(bc$y)])
[1] 2

Лямбда получилась равной 2. Ах, удивительно!
Делаем преобразование и строим новую модель.

Код:

> lm1<-lm(I(y^lambda)~x,data=test)
> summary(lm1)

Call:
lm(formula = I(y^lambda) ~ x, data = test)

Residuals:
       Min         1Q     Median         3Q        Max 
-7.049e-15 -1.023e-15 -1.436e-16  1.235e-15  7.585e-15 

Coefficients:
             Estimate Std. Error   t value Pr(>|t|)    
(Intercept) 0.000e+00  1.324e-15 0.000e+00        1    
x           1.000e+00  1.105e-16 9.047e+15   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 2.85e-15 on 18 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1 
F-statistic: 8.185e+31 on 1 and 18 DF,  p-value: < 2.2e-16

Нарисуем первый график из четырех диагностических графиков модели.

Код:

> plot(lm1,which=1)

Воот! Теперь совсем другое дело.

Следует отметить особые случаи в преобразовании Бокса-Кокса. Если λ=1, то преобразование не требуется. Если λ=0, то следует использовать логарифмическое преобразование log(y). Пример.

Код:

> test<-data.frame(x<-1:20)
> test$y<-exp(-test$x)
> lm0<-lm(y~x,data=test)
> summary(lm0)

Call:
lm(formula = y ~ x, data = test)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.06529 -0.04258 -0.01349  0.02358  0.26464 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)  0.111047   0.034464   3.222  0.00473 **
x           -0.007805   0.002877  -2.713  0.01426 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.07419 on 18 degrees of freedom
Multiple R-squared: 0.2902,	Adjusted R-squared: 0.2508 
F-statistic: 7.359 on 1 and 18 DF,  p-value: 0.01426 

> bc<-boxcox(lm0,lambda = seq(-4,4,0.01))
> (lambda<-bc$x[which.max(bc$y)])
[1] 0
> lm1<-lm(log(y)~x,data=test)
> summary(lm1)

Call:
lm(formula = log(y) ~ x, data = test)

Residuals:
       Min         1Q     Median         3Q        Max 
-2.038e-15 -6.074e-16 -3.526e-16 -8.400e-18  6.208e-15 

Coefficients:
              Estimate Std. Error    t value Pr(>|t|)    
(Intercept)  4.766e-15  7.696e-16  6.193e+00  7.6e-06 ***
x           -1.000e+00  6.425e-17 -1.557e+16  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.657e-15 on 18 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1 
F-statistic: 2.423e+32 on 1 and 18 DF,  p-value: < 2.2e-16

И для закрепления материала два графика. Слева -- плохой, негодный, а справа -- хороший, годный.

Скрипт для графика

Код:

> par(mfrow=c(1,2))
> plot(lm0,which=1)
> plot(lm1,which=1)
> par(mfrow=c(1,1))

Добавлено через 18 часов 12 минут

Кореллограммы

В рассмотренной задаче (#10) рисовалась корреляционная матрица в виде соответствующих диаграмм рассеяния. Также можно получить корреляционную матрицу в виде коэффициентов корреляции в текстовом виде.

Код:

> # Инициализация данных
> x <- c(75, 75, 75, 75, 100, 125, 125, 150, 150, 175, 175, 200, 200, 225, 250, 250, 250, 275, 300, 300, 300, 300)
> y <- c(375, 625, 1000, 1500, 500, 375, 1125, 750, 1250, 825, 1500, 375, 1000, 1125, 750, 1250, 1500, 1375, 375, 825, 1500, 1250)
> z <- c(1.452, 3.705, 18.857, 50.434, 1.956, 2.272, 40.306, 3.368, 47.686, 4.464, 55.63, 2.611, 9.077, 15.304, 7.559, 27.257, 60.174, 35.965, 3.355, 22.491, 38.036, 54.237)
> myData<-data.frame(x,y,z)
> myData$lnZ<-log(myData$z)
>
> # Собственно, вывод корреляционной матрицы
> cor(myData[,c("x","y","z")])
          x         y         z
x 1.0000000 0.2625290 0.2394832
y 0.2625290 1.0000000 0.8857714
z 0.2394832 0.8857714 1.0000000

Обратите внимание, мы взяли не 4 столбца таблицы, а только три. Столбец lnZ опущен. Хотя данная матрица неплохо читается, её можно сделать более удобной для восприятия с помощью функции symnum, заменяющей числовые значения на символы.

Код:

> symnum(cor(myData[,c("x","y","z")]))
  x y z
x 1    
y   1  
z   + 1
attr(,"legend")
[1] 0 ‘ ’ 0.3 ‘.’ 0.6 ‘,’ 0.8 ‘+’ 0.9 ‘*’ 0.95 ‘B’ 1

Для сравнения

Код:

> symnum(cor(myData[,c("x","y","lnZ")]))
    x y l
x   1    
y     1  
lnZ . * 1
attr(,"legend")
[1] 0 ‘ ’ 0.3 ‘.’ 0.6 ‘,’ 0.8 ‘+’ 0.9 ‘*’ 0.95 ‘B’ 1

Однако, есть и другой способ сделать это, да и к тому же красиво, нарисовав коррелограмму. Такой рисунок несказанно обогатит статью и поразит всех читателей. Для этого нужна установленная библиотека corrgram.
Примеры вызовов функции приведены ниже.

Существует 4 возможных отображения корреляции:
panel.pts -- диаграмма разброса
panel.pie -- круговая диаграмма, меняет цвет и заполненность в зависимость от коэффициента корреляции (отрицательная -- красная, положительная -- синяя).
panel.ellipse -- рисует красоту, как показана на рисунке 3. В общем, границы и аппроксимацию.
panel.shade -- закрашивает прямоугольник в зависимости от коэффициента корреляции.

Их присваивают, соответственно lower.panel или upper.panel, в зависимости от того, где нужен рисунок: сверху или снизу.
Также можно выводить название переменных и минимум/максимум значений, как это показано ниже.

Код:

> 
> library(corrgram)
> 
> 
> corrgram(myData[,c("x","y","z")],  lower.panel=panel.pie,
+          upper.panel=panel.pts, text.panel=panel.txt,
+          diag.panel=panel.minmax, main="Кореллограмма 1")
>

Код:

> corrgram(myData[,c("x","y","lnZ")], lower.panel=panel.pts,
+          upper.panel=panel.pie, text.panel=panel.txt,
+          diag.panel=panel.minmax, main="Кореллограмма 2")
>

Код:

> corrgram(log(myData[,c("x","y","z")]), order=TRUE, lower.panel=panel.shade,
+          upper.panel=panel.ellipse, text.panel=panel.txt,
+          main="Кореллограмма 3\n Логарифмы величин")

Понятно, что на трех переменных это не сильно показательно, но если их будет десяток, то картинка будет гораздо симпатичнее.

Hogfather · 05.12.2012, 13:18

Чем хорош R, что в нём можно много всякого забавного сделать. например, захотелось поиграть в регрессионные деревья -- никто слова плохого не скажет.

Рассмотрим ещё раз нашу задачу.

Код:

> library(rpart)
> (fit <- rpart(lnZ ~ x + y, data = myData))
n= 22 

node), split, n, deviance, yval
      * denotes terminal node

1) root 22 35.251650 2.508130  
  2) y< 912.5 10  5.460743 1.319059 *
  3) y>=912.5 12  3.869613 3.499022 *
> summary(fit)
Call:
rpart(formula = lnZ ~ x + y, data = myData)
  n= 22 

         CP nsplit rel error    xerror      xstd
1 0.7353214      0 1.0000000 1.0863277 0.1579386
2 0.0100000      1 0.2646786 0.7911578 0.1943923

Variable importance
 y  x 
91  9 

Node number 1: 22 observations,    complexity param=0.7353214
  mean=2.50813, MSE=1.602348 
  left son=2 (10 obs) right son=3 (12 obs)
  Primary splits:
      y < 912.5 to the left,  improve=0.7353214, (0 missing)
      x < 212.5 to the left,  improve=0.1397332, (0 missing)
  Surrogate splits:
      x < 112.5 to the left,  agree=0.591, adj=0.1, (0 split)

Node number 2: 10 observations
  mean=1.319059, MSE=0.5460743 

Node number 3: 12 observations
  mean=3.499022, MSE=0.3224678 

> par(xpd = NA) # сильное колдунство, чтобы текст не обрезался
> plot(fit)
> text(fit, use.n = TRUE)

Результат, как видите следующий. Если y<912.5, то ln(z) равен 1.319, а если y>=912.5, то 3.499. Смысла для данной задачи никакого, разве что лишний раз убедились, что Икс у нас не при делах.

Как это должно работать

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

В качестве исходных данных берется таблица данных по кифозу у детей из 81 строк and 4 столбцов, описывающая детей у которых были корректирующие операции на позвоночнике

Таблица содержит следующие колонки:

Kyphosis - качественные данные (factor) с вариантами "absent" и "present", показывающая был ли кифоз после операции.
Age- возраст в месяцах
Number - число затронутых позвонков
Start - номер верхнего позвонка, на котором происходила операция

Код:

> # Первые строки таблицы
> head(kyphosis)
  Kyphosis Age Number Start
1   absent  71      3     5
2   absent 158      3    14
3  present 128      4     5
4   absent   2      5     1
5   absent   1      4    15
6   absent   1      2    16
> (fit <- rpart(Kyphosis ~ Age + Number + Start, data = kyphosis))
n= 81 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

 1) root 81 17 absent (0.79012346 0.20987654)  
   2) Start>=8.5 62  6 absent (0.90322581 0.09677419)  
     4) Start>=14.5 29  0 absent (1.00000000 0.00000000) *
     5) Start< 14.5 33  6 absent (0.81818182 0.18181818)  
      10) Age< 55 12  0 absent (1.00000000 0.00000000) *
      11) Age>=55 21  6 absent (0.71428571 0.28571429)  
        22) Age>=111 14  2 absent (0.85714286 0.14285714) *
        23) Age< 111 7  3 present (0.42857143 0.57142857) *
   3) Start< 8.5 19  8 present (0.42105263 0.57894737) *
> (fit2 <- rpart(Kyphosis ~ Age + Number + Start, data = kyphosis,
+               parms = list(prior = c(.65,.35), split = "information")))
n= 81 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

1) root 81 28.350000 absent (0.65000000 0.35000000)  
  2) Start>=12.5 46  3.335294 absent (0.91563089 0.08436911) *
  3) Start< 12.5 35 16.453120 present (0.39676840 0.60323160)  
    6) Age< 34.5 10  1.667647 absent (0.81616742 0.18383258) *
    7) Age>=34.5 25  9.049219 present (0.27932897 0.72067103) *
> par(mfrow = c(1,2), xpd = NA) # otherwise on some devices the text is clipped
> plot(fit)
> text(fit, use.n = TRUE)
> plot(fit2)
> text(fit2, use.n = TRUE)
> par(mfrow = c(1,1))
> ############################################
> summary(fit)
Call:
rpart(formula = Kyphosis ~ Age + Number + Start, data = kyphosis)
  n= 81 

          CP nsplit rel error    xerror      xstd
1 0.17647059      0 1.0000000 1.0000000 0.2155872
2 0.01960784      1 0.8235294 0.9411765 0.2107780
3 0.01000000      4 0.7647059 1.0000000 0.2155872

Variable importance
 Start    Age Number 
    64     24     12 

Node number 1: 81 observations,    complexity param=0.1764706
  predicted class=absent   expected loss=0.2098765  P(node) =1
    class counts:    64    17
   probabilities: 0.790 0.210 
  left son=2 (62 obs) right son=3 (19 obs)
  Primary splits:
      Start  < 8.5  to the right, improve=6.762330, (0 missing)
      Number < 5.5  to the left,  improve=2.866795, (0 missing)
      Age    < 39.5 to the left,  improve=2.250212, (0 missing)
  Surrogate splits:
      Number < 6.5  to the left,  agree=0.802, adj=0.158, (0 split)

Node number 2: 62 observations,    complexity param=0.01960784
  predicted class=absent   expected loss=0.09677419  P(node) =0.7654321
    class counts:    56     6
   probabilities: 0.903 0.097 
  left son=4 (29 obs) right son=5 (33 obs)
  Primary splits:
      Start  < 14.5 to the right, improve=1.0205280, (0 missing)
      Age    < 55   to the left,  improve=0.6848635, (0 missing)
      Number < 4.5  to the left,  improve=0.2975332, (0 missing)
  Surrogate splits:
      Number < 3.5  to the left,  agree=0.645, adj=0.241, (0 split)
      Age    < 16   to the left,  agree=0.597, adj=0.138, (0 split)

Node number 3: 19 observations
  predicted class=present  expected loss=0.4210526  P(node) =0.2345679
    class counts:     8    11
   probabilities: 0.421 0.579 

Node number 4: 29 observations
  predicted class=absent   expected loss=0  P(node) =0.3580247
    class counts:    29     0
   probabilities: 1.000 0.000 

Node number 5: 33 observations,    complexity param=0.01960784
  predicted class=absent   expected loss=0.1818182  P(node) =0.4074074
    class counts:    27     6
   probabilities: 0.818 0.182 
  left son=10 (12 obs) right son=11 (21 obs)
  Primary splits:
      Age    < 55   to the left,  improve=1.2467530, (0 missing)
      Start  < 12.5 to the right, improve=0.2887701, (0 missing)
      Number < 3.5  to the right, improve=0.1753247, (0 missing)
  Surrogate splits:
      Start  < 9.5  to the left,  agree=0.758, adj=0.333, (0 split)
      Number < 5.5  to the right, agree=0.697, adj=0.167, (0 split)

Node number 10: 12 observations
  predicted class=absent   expected loss=0  P(node) =0.1481481
    class counts:    12     0
   probabilities: 1.000 0.000 

Node number 11: 21 observations,    complexity param=0.01960784
  predicted class=absent   expected loss=0.2857143  P(node) =0.2592593
    class counts:    15     6
   probabilities: 0.714 0.286 
  left son=22 (14 obs) right son=23 (7 obs)
  Primary splits:
      Age    < 111  to the right, improve=1.71428600, (0 missing)
      Start  < 12.5 to the right, improve=0.79365080, (0 missing)
      Number < 4.5  to the left,  improve=0.07142857, (0 missing)

Node number 22: 14 observations
  predicted class=absent   expected loss=0.1428571  P(node) =0.1728395
    class counts:    12     2
   probabilities: 0.857 0.143 

Node number 23: 7 observations
  predicted class=present  expected loss=0.4285714  P(node) =0.08641975
    class counts:     3     4
   probabilities: 0.429 0.571

Olafson · 11.12.2012, 10:36

Hogfather,

мне бывает нужен МНК на очень подробных сетках (1000\times 1000, например). Потянет это оболочка? (без ограничения производительности железа)

Hogfather · 11.12.2012, 12:05

Цитата:

Сообщение от Olafson

мне бывает нужен МНК на очень подробных сетках (1000\times 1000, например). Потянет это оболочка? (без ограничения производительности железа)

Специально для Вас прогнал простой тест. Помучался с подобром матрицы, поскольку сингулярные матрицы решаются на раз.
Обратите внимание, я специально выводил время. 13 секунд на модель. Дольше таблицу инициализировал...
Процесcор старый, Intel Core Duo T7250 2 ГГц. Озу 3 Гб

Скрытый текст

Код:

> # Инициализируем матрицу 1000х1000 для x1..x1000, с матожидаением i
> # y присваеваем вектор в 1000 элементов
> Sys.time()
[1] "2012-12-11 12:01:09 MSK"
> MyData<-data.frame(y=rep(0,1000))
> for(i in 1:1000) { 
+ MyData[paste("x",i,sep="")]<-(1:1000)^(1/i)
+ MyData["y"]<-MyData["y"]+MyData[paste("x",i,sep="")]*i+rnorm(10)
+ }
> 
> # Понеслась
> Sys.time()
[1] "2012-12-11 12:01:49 MSK"
> test<-lm(y~.,data=MyData)
> Sys.time()
[1] "2012-12-11 12:02:03 MSK"
> summary(test)

Call:
lm(formula = y ~ ., data = MyData)

Residuals:
    Min      1Q  Median      3Q     Max 
-55.804 -16.049   4.162  14.695  61.563 

Coefficients: (993 not defined because of singularities)
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.290e+06  2.890e+06   0.446    0.655
x1          -2.337e+00  1.785e+01  -0.131    0.896
x2           1.271e+04  5.447e+04   0.233    0.815
x3          -6.552e+05  2.619e+06  -0.250    0.803
x4           5.647e+06  2.185e+07   0.258    0.796
x5          -1.424e+07  5.407e+07  -0.263    0.792
x6           1.048e+07  3.930e+07   0.267    0.790
x7                  NA         NA      NA       NA
x8                  NA         NA      NA       NA
x9                  NA         NA      NA       NA
x10                 NA         NA      NA       NA
x11                 NA         NA      NA       NA
x12                 NA         NA      NA       NA
x13                 NA         NA      NA       NA
x14                 NA         NA      NA       NA
x15                 NA         NA      NA       NA
x16                 NA         NA      NA       NA
x17                 NA         NA      NA       NA
x18                 NA         NA      NA       NA
x19                 NA         NA      NA       NA
x20                 NA         NA      NA       NA
x21         -2.029e+06  7.406e+06  -0.274    0.784
x22                 NA         NA      NA       NA
x23                 NA         NA      NA       NA
x24                 NA         NA      NA       NA
x25                 NA         NA      NA       NA
x26                 NA         NA      NA       NA
x27                 NA         NA      NA       NA
x28                 NA         NA      NA       NA
x29                 NA         NA      NA       NA
.....
x993                NA         NA      NA       NA
x994                NA         NA      NA       NA
x995                NA         NA      NA       NA
x996                NA         NA      NA       NA
x997                NA         NA      NA       NA
x998                NA         NA      NA       NA
x999                NA         NA      NA       NA
x1000               NA         NA      NA       NA

Residual standard error: 30.74 on 992 degrees of freedom
Multiple R-squared: 0.9994,     Adjusted R-squared: 0.9994 
F-statistic: 2.538e+05 on 7 and 992 DF,  p-value: < 2.2e-16

Olafson · 11.12.2012, 12:10

Hogfather, спасибо!

Синтаксис не точно понимаю, но я понял, что R справился. Как он отдает результаты?

Hogfather · 11.12.2012, 12:11

Цитата:

Сообщение от Olafson

Как он отдает результаты?

Можно конкретизировать вопрос: а как хотите?

Olafson · 11.12.2012, 12:16

Цитата:

Сообщение от Hogfather

а как хотите?

Скорее всего, он визуализирует и препарирует любые агрегаты из данных сам.
Но может понадобиться что-нибудь сугубо универсальное, типа текстового файла, для демонстрации результатов через всякие Excel'и с Maple'ами. (ведь R, как я понял, -- GNU, а те далеко нет)

Hogfather · 11.12.2012, 13:14

Цитата:

Сообщение от Olafson

Но может понадобиться что-нибудь сугубо универсальное, типа текстового файла

Функция write.table() позволяет сохранить нужные данные в формате csv. А дальше уже дело техники.

Дублирую табличку из темы про GRETL.

Откуда берем данные	Формат обмена	Куда перегоняем	Примечание
Excel	xls, xlsx или csv	Gretl	Читается прямо в формате Excel
Excel	csv	GNU R	Можно и напрямую, но с доп. библиотекой, а также через clipboard (см. read.table("clipboard",header=T) )
Gretl	csv	Excel	Используйте правильный разделитель данных. Лучше ";". Например, сохраняем таблицу MyData: write.table(MyData,"myData.csv", sep = ";", dec = "."). Для небольшого объема данных можно через clipboard, заменив имя файла на "clipboard" и использовав sep = "\t" (табулятор).
Gretl	GNU R (.RData), csv	GNU R	Лучше использовать внутренний формат .RData, потом в R воспользоваться командой load("имя файла")
GNU R	csv,xls	Excel	Используйте правильный разделитель данных. Лучше "\t". Для экспорта в Excel можно использовать библиотеку xlsReadWrite
GNU R	csv	gretl	В случае, если запущена сессия из gretl, можно также воспользоваться командой gretl.export(). Подробности см. в руководстве пользователя раздел "Gretl and R"

Olafson · 12.12.2012, 00:08

Hogfather, merci за консультацию.

02.12.2012, 16:14	#11
Hogfather Platinum Member Регистрация: 22.07.2010 Адрес: Санкт-Петербург Сообщений: 3,304	Множественная регрессия. Краткий справочник. Попробуем резюмировать. Для подгонки модели МНК, используем код Код: > fit <- lm(lnZ ~ x + y, data=myData) > summary(fit) Call: lm(formula = lnZ ~ x + y, data = myData) Residuals: Min 1Q Median 3Q Max -0.61700 -0.30351 -0.03218 0.19850 0.81873 Coefficients: Estimate Std. Error t value Pr(>\|t\|) (Intercept) -0.5939071 0.3067107 -1.936 0.0678 . x 0.0016541 0.0012529 1.320 0.2025 y 0.0029021 0.0002526 11.488 5.4e-10 * --- Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.4538 on 19 degrees of freedom Multiple R-squared: 0.889, Adjusted R-squared: 0.8773 F-statistic: 76.08 on 2 and 19 DF, p-value: 8.53e-10 Также можно вызвать следующие функции Код: > coefficients(fit) # коэффициенты модели (Intercept) x y -0.593907127 0.001654092 0.002902143 > confint(fit, level=0.95) # доверительные интервалы для параметров модели 2.5 % 97.5 % (Intercept) -1.2358599183 0.048045664 x -0.0009683542 0.004276538 y 0.0023733723 0.003430914 > fitted(fit) [COLOR="rgb(46, 139, 87)"]# вычисленные значения[/COLOR] 1 2 3 4 5 6 7 8 0.6184535 1.3439892 2.4322929 3.8833645 1.0225736 0.7011581 2.8777654 1.8308140 9 10 11 12 13 14 15 16 3.2818856 2.0898271 4.0487737 0.8252150 2.6390544 3.0431746 1.9962232 3.4472948 17 18 19 20 21 22 4.1728306 3.8514150 0.9906242 2.2965886 4.2555352 3.5299994 > residuals(fit) # остатки модели 1 2 3 4 5 6 7 -0.24551154 -0.03430598 0.50459127 0.03730105 -0.35167208 0.11950244 0.81873492 8 9 10 11 12 13 14 -0.61649494 0.58275224 -0.59378185 -0.03005108 0.13451833 -0.43331068 -0.31506039 15 16 17 18 19 20 21 0.02651567 -0.14198445 -0.07559024 -0.26886877 0.21982761 0.81652665 -0.61700212 22 0.46336392 > anova(fit) [COLOR="rgb(46, 139, 87)"]# дисперсионный анализ[/COLOR] Analysis of Variance Table Response: lnZ Df Sum Sq Mean Sq F value Pr(>F) x 1 4.1589 4.1589 20.192 0.0002488 * y 1 27.1795 27.1795 131.963 5.395e-10 * Residuals 19 3.9133 0.2060 --- Signif. codes: 0 ‘*’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 > vcov(fit) # ковариационная матрица для параметров модели (Intercept) x y (Intercept) 9.407143e-02 -2.162456e-04 -4.568268e-05 x -2.162456e-04 1.569874e-06 -8.310037e-08 y -4.568268e-05 -8.310037e-08 6.382438e-08 > influence(fit) # анализ регрессионной модели $hat 1 2 3 4 5 6 7 0.19665640 0.14809675 0.14788589 0.28317818 0.13829845 0.15275498 0.09302145 8 9 10 11 12 13 14 0.06408668 0.09164807 0.05113364 0.14270474 0.15836003 0.04654710 0.05902938 15 16 17 18 19 20 21 0.09846547 0.08567971 0.13738973 0.12649809 0.29922020 0.15799923 0.18148452 22 0.13986129 $coefficients (Intercept) x y 1 -0.090100375 1.924032e-04 4.151880e-05 2 -0.009639360 2.941447e-05 2.351114e-06 3 0.092492681 -5.221326e-04 3.423996e-05 4 0.002356970 -5.638034e-05 1.107148e-05 5 -0.098292644 1.997504e-04 4.375226e-05 6 0.034179275 -3.504510e-05 -2.200750e-05 7 0.068581503 -4.974500e-04 6.895179e-05 8 -0.087542609 1.378093e-04 3.287562e-05 9 0.014114197 -2.636425e-04 6.738312e-05 10 -0.056330941 3.061248e-05 2.300041e-05 11 0.002092611 1.126132e-05 -6.043748e-06 12 0.026144651 5.165601e-05 -2.977428e-05 13 -0.011340740 -3.227886e-05 -3.357501e-06 14 0.009716284 -7.069650e-05 -1.206585e-05 15 0.000820820 1.626464e-05 -2.654591e-06 16 0.012887806 -5.454754e-05 -1.004487e-05 17 0.012131630 -2.194206e-05 -1.245702e-05 18 0.042158631 -1.512501e-04 -2.872849e-05 19 0.018378098 3.404810e-04 -7.109334e-05 20 -0.039976009 8.764986e-04 -8.455088e-05 21 0.143930409 -4.760290e-04 -9.195059e-05 22 -0.072988620 3.945329e-04 2.397855e-05 $sigma 1 2 3 4 5 6 7 8 0.4617764 0.4661857 0.4481136 0.4661523 0.4576378 0.4652627 0.4199357 0.4414127 9 10 11 12 13 14 15 16 0.4434361 0.4435794 0.4662052 0.4649853 0.4543848 0.4599405 0.4662215 0.4649526 17 18 19 20 21 22 0.4658732 0.4613113 0.4621416 0.4164320 0.4376836 0.4511520 $wt.res 1 2 3 4 5 6 7 -0.24551154 -0.03430598 0.50459127 0.03730105 -0.35167208 0.11950244 0.81873492 8 9 10 11 12 13 14 -0.61649494 0.58275224 -0.59378185 -0.03005108 0.13451833 -0.43331068 -0.31506039 15 16 17 18 19 20 21 0.02651567 -0.14198445 -0.07559024 -0.26886877 0.21982761 0.81652665 -0.61700212 22 0.46336392 Сравнение двух моделей Код: > fit1 <- lm(lnZ ~ y, data=myData) > anova(fit,fit1) Analysis of Variance Table Model 1: lnZ ~ x + y Model 2: lnZ ~ y Res.Df RSS Df Sum of Sq F Pr(>F) 1 19 3.9133 2 20 4.2723 -1 -0.35896 1.7428 0.2025 Остановлюсь на этом месте поподробнее. Когда мы добавляем или убавляем регрессоры, мы должны быть уверены, что принципиально это на дисперсию остатков не влияет. В данном случае используется нулевая гипотеза, что коэффициенты регрессоров равны между собой. Примерно то же самое мы делаем, когда проверяем саму модель на значимость(см. ниже), но в таком случаем сравниваем её с уравнением, где все регрессоры имеют коэффициенты 0. Аналогично описанному ниже, если Pr(>F) превышает 0.05, то гипотеза принимается. Т.е. в приведенном выше примере Pr(>F)=0.2025, отсюда разница между Model 1: lnZ ~ x + y и Model 2: lnZ ~ y является статистически незначимой. Автоматическая подгонка модели посредством исключения незначимых членов. Код: > library(MASS) > step <- stepAIC(fit, direction="both") Start: AIC=-31.99 lnZ ~ x + y Df Sum of Sq RSS AIC - x 1 0.359 4.2723 -32.056 <none> 3.9133 -31.987 - y 1 27.180 31.0928 11.611 Step: AIC=-32.06 lnZ ~ y Df Sum of Sq RSS AIC <none> 4.272 -32.056 + x 1 0.359 3.913 -31.987 - y 1 30.979 35.252 12.372 > step$anova # отображаем результаты Stepwise Model Path Analysis of Deviance Table Initial Model: lnZ ~ x + y Final Model: lnZ ~ y Step Df Deviance Resid. Df Resid. Dev AIC 1 19 3.913305 -31.98653 2 - x 1 0.3589588 20 4.272264 -32.05577 Прочую красоту смотреть вот здесь. А подробно проверка регрессионной модели описана вот тут. Последний раз редактировалось Hogfather; 05.12.2012 в 14:11.*
	--------- DNF is not an option

04.12.2012, 15:31	#12
Hogfather Platinum Member Регистрация: 22.07.2010 Адрес: Санкт-Петербург Сообщений: 3,304	Краткий ликбез по моделям, или как читать то, что мы получаем в результате подгонки. Хотя меня особо не спрашивают, но имея опыт работы со студентами постараюсь кратко и без лишнего занудства изложить, что же Пытливый Ум может увидеть в сборище букв на языке эвентуального противника, когда строит модель. Рассмотрим пример, который был приведен выше. Код: > summary(fit) Call: lm(formula = lnZ ~ x + y, data = myData) Residuals: Min 1Q Median 3Q Max -0.61700 -0.30351 -0.03218 0.19850 0.81873 Coefficients: Estimate Std. Error t value Pr(>\|t\|) (Intercept) -0.5939071 0.3067107 -1.936 0.0678 . x 0.0016541 0.0012529 1.320 0.2025 y 0.0029021 0.0002526 11.488 5.4e-10 * --- Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.4538 on 19 degrees of freedom Multiple R-squared: 0.889, Adjusted R-squared: 0.8773 F-statistic: 76.08 on 2 and 19 DF, p-value: 8.53e-10 Итак, первое на что мы смотрим, это не на разные Эр квадрат и прочие коэффициенты, нас интересует только вот эта строчка: F-statistic: 76.08 on 2 and 19 DF, p-value: 8.53e-10 В данном случае нулевая гипотеза выглядит так: H₀: β₁=β₂=...β_n=0, т.е. все коэффициенты кроме свободного члена равны 0. Это мы проводим проверку значимости модели. Если p-value у нас превышает 0.05, то дальнейшая работа с моделью просто бессмысленна, поскольку модель незначимая. Пример незначимой модели Код: > lm0<-lm(y~x,data=test) > summary(lm0) Call: lm(formula = y ~ x, data = test) Residuals: Min 1Q Median 3Q Max -3.3245 -1.4632 0.6214 1.3782 2.6417 Coefficients: Estimate Std. Error t value Pr(>\|t\|) (Intercept) 3.4424 0.8937 3.852 0.00117 x 0.1102 0.0746 1.478 0.15681 --- Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.924 on 18 degrees of freedom Multiple R-squared: 0.1082, Adjusted R-squared: 0.05862 F-statistic: 2.183 on 1 and 18 DF, p-value: 0.1568 Второе, на что мы смотрим, это на остатки модели (Residuals). Код: Residuals: Min 1Q Median 3Q Max -0.61700 -0.30351 -0.03218 0.19850 0.81873 Алгоритм МНК подразумевает, что матожидание остатков должно быть равно 0 и мы надеемся, что остатки имеют нормальное распределение. Напоминаю, что у нормального распределения матожидание и медиана совпадают. В данном случае -- не очень. Поскольку мы видим отрицательную медиану (-0.03218), то распределение остатков незначительно смещено влево. Квартили, по идее, должны быть симметричны относительно 0, здесь же они значительно смещены влево. Максимум и минимум позволяют оценить выбросы (outliers) модели. Третье, на что мы смотрим -- коэффициенты. Код: Coefficients: Estimate Std. Error t value Pr(>\|t\|) (Intercept) -0.5939071 0.3067107 -1.936 0.0678 . x 0.0016541 0.0012529 1.320 0.2025 y 0.0029021 0.0002526 11.488 5.4e-10 * --- Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 Для каждого коэффициента регрессора мы проверяем его значимость, вычисляя соответствующую t-статистику и проверяя нулевую гипотезу о незначимости его коэффициента. В дальнейшем мы будем по очереди убирать регрессоры исходя из наибольшего значения в последней колонке. Впрочем, об этом уже писалось. Дальше, Residual standard error: 1.924 on 18 degrees of freedom говорит нам о средней квадратической ошибке остатков (σ). Multiple R-squared: 0.889, Adjusted R-squared: 0.8773 говорит о коэффициентах детерминации и качестве модели. Чем они ближе к единице, тем лучше. В то же время, проверку адекватности модели следует произвести графически. Именно для этого все графики и строились.. Пример неадекватной модели Код: > test<-data.frame(x<-1:20) > test$y<-sqrt(test$x) > lm0<-lm(y~x,data=test) > summary(lm0) Call: lm(formula = y ~ x, data = test) Residuals: Min 1Q Median 3Q Max -0.46467 -0.08467 0.04310 0.13672 0.17227 Coefficients: Estimate Std. Error t value Pr(>\|t\|) (Intercept) 1.294290 0.081968 15.79 5.44e-12 * x 0.170382 0.006843 24.90 2.13e-15 * --- Signif. codes: 0 ‘*’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1765 on 18 degrees of freedom Multiple R-squared: 0.9718, Adjusted R-squared: 0.9702 F-statistic: 620 on 1 and 18 DF, p-value: 2.127e-15 Обратите внимание на R² А теперь смотрите на график Что должно быть на графике: 1. Точки на графике Residuals vs Fitted должны быть случайно разбросаны, без явно выраженного тренда. Тут этого нет. 2. Точки на графике Normal Q-Q должны более-менее лежать на линии, показывая что остатки имеют нормальное распредение 3. Точки на двух остальных графиках должны быть в группах не слишком далеко от центра. Т.е. модель у нас неадекватна, не смотря на милый Эр Квадрат. Такие дела. Для сравнения посмотрите график на предыдущей странице Добавлено через 4 часа 3 минуты* Использование метода Бокса-Кокса для преобразования регрессионной модели В предыдущем примере мы рассмотрели случай, когда модель выбрана неадекватной. Судя по параболе, которую мы наблюдаем, здесь необходимо какое-то степенное преобразование. Попробуем использовать преобразование Бокса-Кокса, которое входит в пакет MASS. Т.е. найдем такую λ, чтобы y^λ сделал нам аж боже мой. Исходные данные у нас те же и та же модель. Код: > test<-data.frame(x<-1:20) > test$y<-sqrt(test$x) > lm0<-lm(y~x,data=test) А дальше начинается колдунство. В данном случае lambda взята в интервале 0 до 4 с шагом 0.01 ( по умолчанию -2 до 2 с шагом 0.1) Код: > library(MASS) > bc<-boxcox(lm0,lambda = seq(0,4,0.01)) Получаем следующий график зависимости логарифмической функции правдоподобия от лямбды. К сожалению, конкретного результата, который нужно подставить у нас нигде нет. Но мы можем его получить. Функция which.max(bc$y) вернёт нам номер элемента, у которого Игрек наибольший, а bc$x[which.max(bc$y)]) -- собственно нужную лямбду. Код: > (lambda<-bc$x[which.max(bc$y)]) [1] 2 Лямбда получилась равной 2. Ах, удивительно! Делаем преобразование и строим новую модель. Код: > lm1<-lm(I(y^lambda)~x,data=test) > summary(lm1) Call: lm(formula = I(y^lambda) ~ x, data = test) Residuals: Min 1Q Median 3Q Max -7.049e-15 -1.023e-15 -1.436e-16 1.235e-15 7.585e-15 Coefficients: Estimate Std. Error t value Pr(>\|t\|) (Intercept) 0.000e+00 1.324e-15 0.000e+00 1 x 1.000e+00 1.105e-16 9.047e+15 <2e-16 * --- Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.85e-15 on 18 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 8.185e+31 on 1 and 18 DF, p-value: < 2.2e-16 Нарисуем первый график из четырех диагностических графиков модели. Код: > plot(lm1,which=1) Воот! Теперь совсем другое дело. Следует отметить особые случаи в преобразовании Бокса-Кокса. Если λ=1, то преобразование не требуется. Если λ=0, то следует использовать логарифмическое преобразование log(y). Пример. Код: > test<-data.frame(x<-1:20) > test$y<-exp(-test$x) > lm0<-lm(y~x,data=test) > summary(lm0) Call: lm(formula = y ~ x, data = test) Residuals: Min 1Q Median 3Q Max -0.06529 -0.04258 -0.01349 0.02358 0.26464 Coefficients: Estimate Std. Error t value Pr(>\|t\|) (Intercept) 0.111047 0.034464 3.222 0.00473 ** x -0.007805 0.002877 -2.713 0.01426 * --- Signif. codes: 0 ‘*’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.07419 on 18 degrees of freedom Multiple R-squared: 0.2902, Adjusted R-squared: 0.2508 F-statistic: 7.359 on 1 and 18 DF, p-value: 0.01426 > bc<-boxcox(lm0,lambda = seq(-4,4,0.01)) > (lambda<-bc$x[which.max(bc$y)]) [1] 0 > lm1<-lm(log(y)~x,data=test) > summary(lm1) Call: lm(formula = log(y) ~ x, data = test) Residuals: Min 1Q Median 3Q Max -2.038e-15 -6.074e-16 -3.526e-16 -8.400e-18 6.208e-15 Coefficients: Estimate Std. Error t value Pr(>\|t\|) (Intercept) 4.766e-15 7.696e-16 6.193e+00 7.6e-06 x -1.000e+00 6.425e-17 -1.557e+16 < 2e-16 * --- Signif. codes: 0 ‘*’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.657e-15 on 18 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 2.423e+32 on 1 and 18 DF, p-value: < 2.2e-16 И для закрепления материала два графика. Слева -- плохой, негодный, а справа -- хороший, годный. Скрипт для графика Код: > par(mfrow=c(1,2)) > plot(lm0,which=1) > plot(lm1,which=1) > par(mfrow=c(1,1)) Добавлено через 18 часов 12 минут* Кореллограммы В рассмотренной задаче (#10) рисовалась корреляционная матрица в виде соответствующих диаграмм рассеяния. Также можно получить корреляционную матрицу в виде коэффициентов корреляции в текстовом виде. Код: > # Инициализация данных > x <- c(75, 75, 75, 75, 100, 125, 125, 150, 150, 175, 175, 200, 200, 225, 250, 250, 250, 275, 300, 300, 300, 300) > y <- c(375, 625, 1000, 1500, 500, 375, 1125, 750, 1250, 825, 1500, 375, 1000, 1125, 750, 1250, 1500, 1375, 375, 825, 1500, 1250) > z <- c(1.452, 3.705, 18.857, 50.434, 1.956, 2.272, 40.306, 3.368, 47.686, 4.464, 55.63, 2.611, 9.077, 15.304, 7.559, 27.257, 60.174, 35.965, 3.355, 22.491, 38.036, 54.237) > myData<-data.frame(x,y,z) > myData$lnZ<-log(myData$z) > > # Собственно, вывод корреляционной матрицы > cor(myData[,c("x","y","z")]) x y z x 1.0000000 0.2625290 0.2394832 y 0.2625290 1.0000000 0.8857714 z 0.2394832 0.8857714 1.0000000 Обратите внимание, мы взяли не 4 столбца таблицы, а только три. Столбец lnZ опущен. Хотя данная матрица неплохо читается, её можно сделать более удобной для восприятия с помощью функции symnum, заменяющей числовые значения на символы. Код: > symnum(cor(myData[,c("x","y","z")])) x y z x 1 y 1 z + 1 attr(,"legend") [1] 0 ‘ ’ 0.3 ‘.’ 0.6 ‘,’ 0.8 ‘+’ 0.9 ‘’ 0.95 ‘B’ 1 Для сравнения Код: > symnum(cor(myData[,c("x","y","lnZ")])) x y l x 1 y 1 lnZ . 1 attr(,"legend") [1] 0 ‘ ’ 0.3 ‘.’ 0.6 ‘,’ 0.8 ‘+’ 0.9 ‘’ 0.95 ‘B’ 1 Однако, есть и другой способ сделать это, да и к тому же красиво, нарисовав коррелограмму. Такой рисунок несказанно обогатит статью и поразит всех читателей. Для этого нужна установленная библиотека corrgram. Примеры вызовов функции приведены ниже. Существует 4 возможных отображения корреляции: panel.pts -- диаграмма разброса panel.pie -- круговая диаграмма, меняет цвет и заполненность в зависимость от коэффициента корреляции (отрицательная -- красная, положительная -- синяя). panel.ellipse -- рисует красоту, как показана на рисунке 3. В общем, границы и аппроксимацию. panel.shade -- закрашивает прямоугольник в зависимости от коэффициента корреляции. Их присваивают, соответственно lower.panel или upper.panel, в зависимости от того, где нужен рисунок: сверху или снизу. Также можно выводить название переменных и минимум/максимум значений, как это показано ниже. Код: > > library(corrgram) > > > corrgram(myData[,c("x","y","z")], lower.panel=panel.pie, + upper.panel=panel.pts, text.panel=panel.txt, + diag.panel=panel.minmax, main="Кореллограмма 1") > Код: > corrgram(myData[,c("x","y","lnZ")], lower.panel=panel.pts, + upper.panel=panel.pie, text.panel=panel.txt, + diag.panel=panel.minmax, main="Кореллограмма 2") > Код: > corrgram(log(myData[,c("x","y","z")]), order=TRUE, lower.panel=panel.shade, + upper.panel=panel.ellipse, text.panel=panel.txt, + main="Кореллограмма 3\n Логарифмы величин") Понятно, что на трех переменных это не сильно показательно, но если их будет десяток, то картинка будет гораздо симпатичнее. Последний раз редактировалось Hogfather; 05.12.2012 в 14:07.*
	--------- DNF is not an option

05.12.2012, 13:18	#13
Hogfather Platinum Member Регистрация: 22.07.2010 Адрес: Санкт-Петербург Сообщений: 3,304	Регресионные деревья Чем хорош R, что в нём можно много всякого забавного сделать. например, захотелось поиграть в регрессионные деревья -- никто слова плохого не скажет. Рассмотрим ещё раз нашу задачу. Код: > library(rpart) > (fit <- rpart(lnZ ~ x + y, data = myData)) n= 22 node), split, n, deviance, yval * denotes terminal node 1) root 22 35.251650 2.508130 2) y< 912.5 10 5.460743 1.319059 * 3) y>=912.5 12 3.869613 3.499022 * > summary(fit) Call: rpart(formula = lnZ ~ x + y, data = myData) n= 22 CP nsplit rel error xerror xstd 1 0.7353214 0 1.0000000 1.0863277 0.1579386 2 0.0100000 1 0.2646786 0.7911578 0.1943923 Variable importance y x 91 9 Node number 1: 22 observations, complexity param=0.7353214 mean=2.50813, MSE=1.602348 left son=2 (10 obs) right son=3 (12 obs) Primary splits: y < 912.5 to the left, improve=0.7353214, (0 missing) x < 212.5 to the left, improve=0.1397332, (0 missing) Surrogate splits: x < 112.5 to the left, agree=0.591, adj=0.1, (0 split) Node number 2: 10 observations mean=1.319059, MSE=0.5460743 Node number 3: 12 observations mean=3.499022, MSE=0.3224678 > par(xpd = NA) # сильное колдунство, чтобы текст не обрезался > plot(fit) > text(fit, use.n = TRUE) Результат, как видите следующий. Если y<912.5, то ln(z) равен 1.319, а если y>=912.5, то 3.499. Смысла для данной задачи никакого, разве что лишний раз убедились, что Икс у нас не при делах. Как это должно работать В данном случае я немного модифицировал штатный пример из документации В качестве исходных данных берется таблица данных по кифозу у детей из 81 строк and 4 столбцов, описывающая детей у которых были корректирующие операции на позвоночнике Таблица содержит следующие колонки: Kyphosis - качественные данные (factor) с вариантами "absent" и "present", показывающая был ли кифоз после операции. Age- возраст в месяцах Number - число затронутых позвонков Start - номер верхнего позвонка, на котором происходила операция Код: > # Первые строки таблицы > head(kyphosis) Kyphosis Age Number Start 1 absent 71 3 5 2 absent 158 3 14 3 present 128 4 5 4 absent 2 5 1 5 absent 1 4 15 6 absent 1 2 16 > (fit <- rpart(Kyphosis ~ Age + Number + Start, data = kyphosis)) n= 81 node), split, n, loss, yval, (yprob) * denotes terminal node 1) root 81 17 absent (0.79012346 0.20987654) 2) Start>=8.5 62 6 absent (0.90322581 0.09677419) 4) Start>=14.5 29 0 absent (1.00000000 0.00000000) * 5) Start< 14.5 33 6 absent (0.81818182 0.18181818) 10) Age< 55 12 0 absent (1.00000000 0.00000000) * 11) Age>=55 21 6 absent (0.71428571 0.28571429) 22) Age>=111 14 2 absent (0.85714286 0.14285714) * 23) Age< 111 7 3 present (0.42857143 0.57142857) * 3) Start< 8.5 19 8 present (0.42105263 0.57894737) * > (fit2 <- rpart(Kyphosis ~ Age + Number + Start, data = kyphosis, + parms = list(prior = c(.65,.35), split = "information"))) n= 81 node), split, n, loss, yval, (yprob) * denotes terminal node 1) root 81 28.350000 absent (0.65000000 0.35000000) 2) Start>=12.5 46 3.335294 absent (0.91563089 0.08436911) * 3) Start< 12.5 35 16.453120 present (0.39676840 0.60323160) 6) Age< 34.5 10 1.667647 absent (0.81616742 0.18383258) * 7) Age>=34.5 25 9.049219 present (0.27932897 0.72067103) * > par(mfrow = c(1,2), xpd = NA) # otherwise on some devices the text is clipped > plot(fit) > text(fit, use.n = TRUE) > plot(fit2) > text(fit2, use.n = TRUE) > par(mfrow = c(1,1)) > ############################################ > summary(fit) Call: rpart(formula = Kyphosis ~ Age + Number + Start, data = kyphosis) n= 81 CP nsplit rel error xerror xstd 1 0.17647059 0 1.0000000 1.0000000 0.2155872 2 0.01960784 1 0.8235294 0.9411765 0.2107780 3 0.01000000 4 0.7647059 1.0000000 0.2155872 Variable importance Start Age Number 64 24 12 Node number 1: 81 observations, complexity param=0.1764706 predicted class=absent expected loss=0.2098765 P(node) =1 class counts: 64 17 probabilities: 0.790 0.210 left son=2 (62 obs) right son=3 (19 obs) Primary splits: Start < 8.5 to the right, improve=6.762330, (0 missing) Number < 5.5 to the left, improve=2.866795, (0 missing) Age < 39.5 to the left, improve=2.250212, (0 missing) Surrogate splits: Number < 6.5 to the left, agree=0.802, adj=0.158, (0 split) Node number 2: 62 observations, complexity param=0.01960784 predicted class=absent expected loss=0.09677419 P(node) =0.7654321 class counts: 56 6 probabilities: 0.903 0.097 left son=4 (29 obs) right son=5 (33 obs) Primary splits: Start < 14.5 to the right, improve=1.0205280, (0 missing) Age < 55 to the left, improve=0.6848635, (0 missing) Number < 4.5 to the left, improve=0.2975332, (0 missing) Surrogate splits: Number < 3.5 to the left, agree=0.645, adj=0.241, (0 split) Age < 16 to the left, agree=0.597, adj=0.138, (0 split) Node number 3: 19 observations predicted class=present expected loss=0.4210526 P(node) =0.2345679 class counts: 8 11 probabilities: 0.421 0.579 Node number 4: 29 observations predicted class=absent expected loss=0 P(node) =0.3580247 class counts: 29 0 probabilities: 1.000 0.000 Node number 5: 33 observations, complexity param=0.01960784 predicted class=absent expected loss=0.1818182 P(node) =0.4074074 class counts: 27 6 probabilities: 0.818 0.182 left son=10 (12 obs) right son=11 (21 obs) Primary splits: Age < 55 to the left, improve=1.2467530, (0 missing) Start < 12.5 to the right, improve=0.2887701, (0 missing) Number < 3.5 to the right, improve=0.1753247, (0 missing) Surrogate splits: Start < 9.5 to the left, agree=0.758, adj=0.333, (0 split) Number < 5.5 to the right, agree=0.697, adj=0.167, (0 split) Node number 10: 12 observations predicted class=absent expected loss=0 P(node) =0.1481481 class counts: 12 0 probabilities: 1.000 0.000 Node number 11: 21 observations, complexity param=0.01960784 predicted class=absent expected loss=0.2857143 P(node) =0.2592593 class counts: 15 6 probabilities: 0.714 0.286 left son=22 (14 obs) right son=23 (7 obs) Primary splits: Age < 111 to the right, improve=1.71428600, (0 missing) Start < 12.5 to the right, improve=0.79365080, (0 missing) Number < 4.5 to the left, improve=0.07142857, (0 missing) Node number 22: 14 observations predicted class=absent expected loss=0.1428571 P(node) =0.1728395 class counts: 12 2 probabilities: 0.857 0.143 Node number 23: 7 observations predicted class=present expected loss=0.4285714 P(node) =0.08641975 class counts: 3 4 probabilities: 0.429 0.571 Последний раз редактировалось Hogfather; 06.12.2012 в 08:28.
	--------- DNF is not an option

11.12.2012, 10:36	#14
Olafson Gold Member Регистрация: 08.02.2009 Сообщений: 1,408	Hogfather, мне бывает нужен МНК на очень подробных сетках (1000\times 1000, например). Потянет это оболочка? (без ограничения производительности железа)

11.12.2012, 12:10	#16
Olafson Gold Member Регистрация: 08.02.2009 Сообщений: 1,408	Hogfather, спасибо! Синтаксис не точно понимаю, но я понял, что R справился. Как он отдает результаты?

12.12.2012, 00:08	#20
Olafson Gold Member Регистрация: 08.02.2009 Сообщений: 1,408	Hogfather, merci за консультацию.

Реклама