Методы и критерии верификации регрессионных моделей

Автор: Владимир Шитиков

В предыдущих сообщениях неоднократно поднимались важные проблемы, связанные с диагностикой регрессионных моделей: проверка статистических допущений, лежащих в основе используемого метода построения модели, оценка адекватности структуры систематической части модели, тестирование чувствительности модели к аномалиям в структуре исходных данных и др. (см. часть 1, часть 2, часть 3). В данном сообщения попробуем в сжатой форме дать сводку методов и критериев верификации регрессионных моделей и подвести некоторые итоги того, что было частично рассмотрено ранее.

Действительно, первое, что следует сделать сразу после построения модели, так это убедиться в том, насколько она полезна и эффективна. Проверка работоспособности полученных моделей (англ. model verification) представляет собой сложный и многостадийный процесс, по итогам которого может быть достигнут приемлемый уровень уверенности исследователя в том, что результаты, полученные при практическом использовании модели, окажутся правильными. Верификация включает оценку качества модели (англ. model evaluation) и ее адекватности выборочным данным (англ. model validation) по совокупности объективных критериев, анализ вклада отдельных предикторов и отбор их оптимальной комбинации (англ. model selection) и, наконец, независимое тестирование модели с целью настройки ее параметров (англ. model tuning).

Тестирование регрессионных моделей ставит своей задачей выявление следующих основных проблем:

• Смещение (англ. bias), или систематическая ошибка модели (например, сдвиг предсказываемых значений на какую-то трудно объяснимую величину);
• Высокая дисперсия прогноза, определяемая чаще всего излишней чувствительностью модели к небольшим изменениям в распределении обучающих данных;
• Неадекватность - свойство модели, которая не отражает основных закономерностей в генеральной совокупности и в силу случайного характера обучающей выборки;
• Переусложнение, или переобучение модели (англ. overfitting), которое «...так же вредно, как и ее недоусложнение» (Ивахненко, 1982). Переобученная модель дает хорошие результаты подгонки на обучающей совокупности, но часто ведет себя непредсказуемо при прогнозировании на "свежих" данных, не участвовавших в построении модели. Построение модели оптимальной сложности, в которой сбалансирована точность и устойчивость прогнозных значений, является фундаментальной задачей совершенствования методов моделирования. Тема переобучения моделей уже поднималась нами ранее, однако, безусловно, заслуживает дополнительного обсуждения.

Данные, используемые в примерах

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

Здесь и далее для построения графиков воспользуемся одним из наиболее популярных графических пакетов для R – ggplot2, который позволяет строить как всевозможные простые диаграммы, так и гораздо более сложные графики, такие как двухмерные диаграммы распределения плотности (см., например, Мастицкий, 2016). Легко показать линию простой регрессии с ее доверительными интервалами, которая в рассматриваемом нами примере не вполне адекватна имеющимся данным:

library(ggplot2)
library(DAAG)
data("fruitohms")
ggplot(fruitohms, aes(x = juice, y = ohms)) + geom_point() +
       stat_smooth(method = "lm") + 
       xlab("Содержание сока, %") + 
       ylab("Сопротивление, Ом")

Сглаживание данных c использованием базовой R-функции loess() лучше подходит для описания этой зависимости (в качестве примера, вместо обычной диаграммы рассеяния данные будут изображены в виде "сотовой диаграммы", где наблюдения распределены по ячейкам и число наблюдений в каждой ячейке представлено цветом соответствующего оттенка):

library(hexbin)
ggplot(fruitohms, aes(x = juice, y = ohms)) +
       geom_hex(binwidth = c(3, 500)) +
       geom_smooth(color = "red", se = FALSE) + 
       xlab("Содержание сока, %") + 
       ylab("Сопротивление, Ом")

Оценка адекватности (валидности) модели

При выполнении поиска оптимальной модели для тех или иных данных полезно иметь набор "эталонов" для сравнения. Если мы зададимся каким-нибудь подходящим критерием качества, то сможем оценить, насколько тестируемая модель отличается по этому критерию от эталона. Такими эталонными моделями могут быть следующие:

• Нулевая модель (англ. null model) – модель, которая не содержит никаких предикторов. В случае с количественной зависимой переменной нулевая модель описывает ни что иное, как среднее значение этой переменной в имеющейся обучающей выборке. В случае же с качественной зависимой переменной (задача классификации) нулевая модель предсказывает наиболее часто встречающуюся категорию. Если тестируемая нами предсказательная модель не превосходит по своей эффективности нулевую, то результат моделирования можно трактовать как неудачу.
• Насыщенная модель (англ. saturated model) – модель, включающая все имеющиеся предикторы. Если качество тестируемой модели по тому или иному критерию значимо превышает таковое у нулевой модели и приближается к качеству максимально насыщенной модели, то процесс отбора оптимальной модели можно считать завершенным.
• Модель с одним наиболее информативным предиктором. Если в процессе отбора более сложные модели по своей эффективности не превосходят модель с одним наиболее информативным предиктором, то включение дополнительных предикторов вряд ли имеет смысл.

Процедуры оценки адекватности классической модели регрессии (обозначим ее как f) основаны на анализе остатков (англ. residuals), т.е. разностей между прогнозируемыми f(x[i,]) и наблюдаемыми y[i] значениями, где x - матрица независимых переменных (предикторов). Базовыми критериями являются сумма квадратов остатков RSS (англ. residual sum of squares), корень из среднеквадратичной ошибки RMSE (англ. root mean square error) и стандартное отклонение остатков RSE (англ. residual standard error): В рассматриваемом примере у нас есть лишь одна зависимая переменная (электропроводность) и один предиктор (содержание сока в мякоти фруктов). Для удобства обозначим их как y и x соответственно:

y <- fruitohms$juice
x <- fruitohms$ohms

#  Строим модель с одним предиктором:
n <- dim(as.matrix(x))[1]
m <- dim(as.matrix(x))[2] 
M_reg <- lm(y ~ x)
pred <- predict(M_reg)
RSS <- sum((y - pred) * (y - pred))
RMSE <- sqrt(RSS/n)
RSE <- sqrt(RSS/(n - m - 1)) 
c(RSS, RMSE, RSE)
[1] 12548.645657     9.901328     9.979600

Другими важными критериями качества аппроксимации данных являются коэффициент детерминации (Rsquared), дисперсионное отношение Фишера (F) и соответствующее ему р-значение:

Rsquared <- 1 - RSS/sum((mean(y) - y)^2)
F <- (sum((mean(y) - pred)^2)/m)/(RSS/(n - m - 1))
p <- pf(q = F, df1 = m, df2 = (n - m - 1), lower.tail = FALSE)
c(Rsquared, F, p)
[1] 6.387089e-01 2.227492e+02 1.234110e-29

Разумеется, все эти величины можно было бы получить и с помощью команды summary(M_reg), но нам показалось интересным показать всю "кухню" их расчетов. Отметим, что F-критерий интерпретируется как мера превышения качества построенной модели над качеством нулевой модели, которая не содержит объясняющих переменных.

summary(M_reg)
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 66.1364158  2.2556380   29.32   <2e-16 ***
x           -0.0071064  0.0004762  -14.93   <2e-16 ***
---
Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1  1 
Residual standard error: 9.98 on 126 degrees of freedom
Multiple R-squared: 0.6387,  Adjusted R-squared: 0.6358 
F-statistic: 222.7 on 1 and 126 DF,  p-value: < 2.2e-16

Чтобы выяснить, насколько объясненная нашей моделью сумма квадратов остатков превысила общую сумму квадратов остатков (т.е. сумму квадратов остатков у нулевой модели) можно выполнить дисперсионный анализ (см. также здесь):

anova(M_reg)
Analysis of Variance Table
Response: y
           Df Sum Sq Mean Sq F value    Pr(>F)    
x           1  22184 22184.1  222.75 < 2.2e-16 ***
Residuals 126  12549    99.6

В случае обобщенных линейных моделей (Generalized Linear Models, GLM) вместо минимизации суммы квадратов отклонений выполняется поиск экстремума логарифма функции максимального правдоподобия (англ. maximum likelihood), вид которой зависит от характера распределения данных. В общем случае значение функции правдоподобия численно равно вероятности того, что модель правильно предсказывает любое предъявленное ей наблюдение из заданной выборки. Логарифм функции правдоподобия (LL) всегда будет отрицательным, и поскольку \(\log(1) = 0\), то для наиболее оптимальной модели, прогнозирующей значение отклика с наименьшей ошибкой, значение LL будет стремиться к 0. Часто для оценки расхождений между наблюдаемыми и прогнозируемыми значениями вместо LL используют девианс (англ. deviance; в русскоязычной литературе используются также термины "девиация" и "аномалия"), который определен как \(D = -2(LL - S)\), где S – правдоподобие модели с минимально возможным уровнем ошибки (такую ошибку называют также неустранимой, или (в задачах классификации) байесовской). Обычно принимают S = 0, поскольку истинно оптимальная модель, как правило, выполняет верные предсказания.

Чем меньше остаточный девианс D некоторой модели, тем лучше эта модель. Аналогично можно рассчитать логарифм правдоподобия и девианс для нулевой модели (ниже обозначен как D.null). Тогда адекватность тестируемой модели можно описать при помощи псевдо-коэффициента детерминации (Prsquar), который рассчитывается с учетом того, во сколько раз D превышает D.null. Подобный анализ является обобщением техники дисперсионного анализа. При этом статистическую значимость разности девиансов можно оценить по критерию \(\chi^2\):

M_glm <- glm(y ~ x)
lgLik <- logLik(M_glm)
D.null <- M_glm$null.deviance 
D <- M_glm$deviance
df <- with(M_glm, df.null - df.residual)
p <- pchisq(D.null - D, df, lower.tail = FALSE)
PRsquar = 1 - D/D.null
c(lgLik, D , D.null, p, PRsquar)
[1]  -475.08574  12548.64565 34732.77929 0.00000 0.63871

То же можно получить и с использованием функции summary():

summary(M_glm)
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 66.1364158  2.2556380   29.32   <2e-16 ***
x           -0.0071064  0.0004762  -14.93   <2e-16 ***
---
Signif. codes:  0 С***Т 0.001 С**Т 0.01 С*Т 0.05 С.Т 0.1 С Т 1 
(Dispersion parameter for gaussian family 99.59243)
    Null deviance: 34733  on 127  degrees of freedom
Residual deviance: 12549  on 126  degrees of freedom
AIC: 956.17

with(M_glm, null.deviance - deviance)
[1] 22184.13

anova(M_glm, glm(x ~ 1) , test = "Chisq")
Analysis of Deviance Table
     Df Deviance Resid. Df Resid. Dev P(>|Chi|)    
NULL                   127      34733              
x     1    22184       126      12549 < 2.2e-16 ***

Нетрудно заметить, что классическая (т.е. рассчитываемая по методу наименьших квадратов) и обобщенная линейные модели в случае нормального распределения дают идентичные результаты.


Оптимизация параметров моделей с использованием информационных критериев

Поиск оптимальной модели часто реализуется с использованием тех или иных критериев качества, налагающих "штраф" на добавление новых параметров. Если D - выборочный остаточный девианс, k – число свободных параметров модели, а n - объем обучающей выборки, то можно определить семейство следующих информационных критериев (Burnham, Anderson, 2002):

• классический критерий Акаике: AIC = D + 2 * k
• байесовский критерий Шварца: BIC = D + k * ln(n)
• скорректированный критерий Акаике: AICс = AIC + 2 * k * (k+1)/(n-k-1)

При построении параметрических моделей с использованием R рассчитать эти критерии можно различными способами:

k <- extractAIC(M_glm)[1]
AIC <- extractAIC(M_glm)[2]
AIC <- AIC(M_glm)
AICc <- AIC(M_glm) + 2 * k * (k + 1)/(n - k - 1)
BIC <- BIC(M_glm)
BIC <- AIC(M_glm, k = log(n))
c(AIC, AICc, BIC)
[1] 956.1715 956.2675 964.7276

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

# Построение моделей со степенью полинома от 1 до 7:
max.poly <- 7

# Создание пустой таблицы для хранения значений AIC и BIC,
# рассчитанных для всех моделей и ее заполнение:
AIC.BIC <- data.frame(criterion = c(rep("AIC", max.poly),
                                    rep("BIC", max.poly)),
                      value = numeric(max.poly * 2),
                      degree = rep(1:max.poly, times = 2))
for(i in 1:max.poly) {
     AIC.BIC[i, 2] <- AIC(lm(y ~ poly(x, i)))
     AIC.BIC[i + max.poly, 2] <- BIC(lm(y ~ poly(x, i)))
}

# График AIC и BIC для разных степеней полинома:
qplot(degree, value, data = AIC.BIC, 
      geom = "line", linetype = criterion) +
xlab("Степень полинома") + ylab("Значение критерия")

Минимальные значения информационных критериев приводят к выводу о том, что оптимальным является полином 5-й степени.


Использование алгоритмов ресэмплинга для тестирования и оптимизации параметров моделей

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

Перекрестная проверка или "кросс-проверка" (англ. cross-validation, CV) является общей процедурой эмпирического оценивания параметров моделей по выборочным данным. В ходе многократного случайного разбиения исходной выборки на обучающее и проверочное подмножества оценка параметров модели выполняется только на обучающей выборке, тогда как оценка погрешности прогноза выполняется для наблюдений из проверочной совокупности.

Одной из часто применяемых является методика r×k-кратной перекрестной проверки (англ. r×k-fold cross-validation), когда исходная выборка случайным образом r раз разбивается на k блоков (англ. folds) равной (или почти равной) длины. При реализации каждой повторности r (англ. replication) один из блоков по очереди становится проверочной последовательностью, а все остальные блоки – обучающей выборкой. При этом генерируется r×n модельных значений отклика (n – общее количество наблюдений), по которым оценивается ошибка перекрестной проверки.

Выполнить разбиение исходной выборки y на k блоков можно с использованием следующих функций из пакета caret (Kuhn, Johnson, 2013), о которых мы упоминали ранее при рассмотрении функции createDataPartition():

createFolds (y, k = 10)
createMultiFolds (y, k = 10, times = r)

Частным случаем является "скользящий контроль", или перекрестная проверка с последовательным исключением одного наблюдения (англ. leave-one-out CV), т.е. k = n. При этом строится n моделей по (n – 1) выборочным значениям, а исключенное наблюдение каждый раз используется для расчета ошибки прогноза. В. Н. Вапник (1984) теоретически обосновал применение скользящего контроля и показал, что если исходные выборки независимы, то средняя ошибка перекрестной проверки дает несмещенную оценку ошибки модели. Это выгодно отличает ее от средней ошибки на обучающей выборке, которая при переусложнении модели может оказаться оптимистично заниженной.

Выше мы рассмотрели пример выбора оптимального числа параметров полиномиальной регрессии с использованием информационных критериев. Попробуем достичь той же цели, выполнив перекрестную проверку моделей в режиме "leave-one-out" с использованием самостоятельно оформленной функции:

# Функция для выполнения скользящего контроля для модели y ~ poly(x, degree):
crossvalidate <- function(x, y, degree) {
   
   preds <- numeric(length(x))

   for (i in 1:length(x)) {
      x.in <- x[-i]
      x.out <- x[i]
      y.in <- y[-i]
      y.out <- x[i]
      m <- lm(y.in ~ poly(x.in, degree = degree) )
      new <- data.frame(x.in = seq(-3, 3, by = 0.1))
      preds[i] <- predict(m, newdata = data.frame(x.in = x.out))
   }

  # Тестовая статистика - сумма квадратов отклонений:
  return(sum((y - preds)^2))
}

# Заполнение таблицы результатами перекрестной проверки 
# и сохранение квадрата ошибки в таблице "a":
a <- data.frame(cross = numeric(max.poly))
for (i in 1:max.poly)
{
  a[i, 1] <- crossvalidate(x, y, degree = i)
}

# График суммы квадратов ошибки при перекрестной проверке:
qplot(1:max.poly, cross, data = a, geom = c("line")) +
      xlab("Степень полинома") + ylab("Квадратичная ошибка")

Из представленного выше графика следует, что в нашем конкретном случае оба метода в качестве оптимальной выбрали модель с полиномом 5-й степени. Проверим, имеется ли значимое превышение ее качества над моделью с одним предиктором:

M_poly <- glm(y ~ poly(x, 5))
summary(M_poly)
Coefficients:
(Intercept)  poly(x, 5)1  poly(x, 5)2  poly(x, 5)3  poly(x, 5)4  poly(x, 5)5  
      35.15      -148.94        -5.54        51.08       -13.90       -23.53  
Degrees of Freedom: 127 Total (i.e. Null);  122 Residual
Null Deviance:      34730 
Residual Deviance: 9162         AIC: 923.9 

anova(M_glm, M_poly, test = "Chisq")
Analysis of Deviance Table
Model 1: y ~ x
Model 2: y ~ poly(x, 5)
  Resid. Df Resid. Dev Df Deviance P(>|Chi|)    
1       126      12549                          
2       122       9162  4   3386.6 3.798e-09 ***

Как видим, имеет место статистически значимое снижение ошибки модели.

В R перекрестную проверку модели можно выполнить с использованием не только самостоятельно оформленных процедур, но и функций из нескольких общедоступных пакетов:

• CVlm(df, form.lm, m = 3) из пакета DAAG, где df – исходная таблица с данными, form.lm – формула линейной модели, m – число блоков (сюда подставляется значение k);
• cv.glm(df, glmfit, K = n) из пакета boot, где df – исходная таблица с данными, glmfit – объект обобщенной линейной модели класса glm, K – число блоков (по умолчанию выполняется скользящий контроль);
• cvLm(lmfit, folds = cvFolds(n, K, R), cost) из пакета cvTools, где lmfit – объект обобщенной линейной модели, K – число блоков, R – число репликаций, cost – наименование одного из пяти разновидностей выводимых критериев перекрестной проверки;
• train(form, df, method, trainControl, ...) из пакета caret, где df – исходная таблица с данными, form, method – формула и метод построения модели, trainControl – объект, в котором прописываются все необходимые параметры перекрестной проверки.

Использование функции train(...) будет предметом наших дальнейших сообщений.


Литературные источники

• Вапник В. Н. (1984) Алгоритмы и программы восстановления зависимостей. М.: Наука, 816 с.
• Ивахненко А. Г. (1982) Индуктивный метод самоорганизации моделей сложных систем. Киев: Наук. думка, 296 с.
• Мастицкий С. Э. (2016) Визуализация данных с помощью ggplot2. М.: ДМК Пресс, 222 с.
• Burnham K. P., Anderson D. R. (2002) Model selection and multimodel inference: a practical information-theoretic approach. N.Y.: Springer-Verlag, 496 p.
• Kuhn M., Johnson K. (2013) Applied Predictive Modeling. Springer, 600 p.