Показать сообщение отдельно
Старый 04.12.2012, 15:31   #12
Hogfather
Platinum Member
 
Аватар для Hogfather
 
Регистрация: 22.07.2010
Адрес: Санкт-Петербург
Сообщений: 3,281
По умолчанию Краткий ликбез по моделям, или как читать то, что мы получаем в результате подгонки.

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

Рассмотрим пример, который был приведен выше.
Код:
> 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

В данном случае нулевая гипотеза выглядит так: H0: β12=...β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
Обратите внимание на R2
А теперь смотрите на график



Что должно быть на графике:
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
Hogfather вне форума   Ответить с цитированием
Реклама