Сложность тестирования линейности в регрессии

В статистическом моделировании: две культуры Лев Брейман пишет

В настоящее время применяется практика проверки соответствия модели данных с помощью тестов соответствия и анализа остаточных данных. Однажды, несколько лет назад, я поставил задачу симулированной регрессии в семи измерениях с контролируемой степенью нелинейности. Стандартные тесты на пригодность к приему не отклоняли линейность, пока нелинейность не была экстремальной.

Брейман не дает подробностей своего моделирования. Он ссылается на статью, которая, по его словам, дает теоретическое обоснование его наблюдения, но статья не опубликована.

Кто-нибудь видел опубликованный результат моделирования или теоретическую статью, подтверждающую утверждение Бримана?

Я создал симуляцию, которая отвечала бы описанию Бреймана, и нашел только очевидное: результат зависит от контекста и того, что подразумевается под «экстремальным».

Можно сказать очень много, но позвольте мне ограничиться одним примером, который проводится с помощью легко модифицируемого Rкода, который заинтересованные читатели могут использовать в своих собственных исследованиях. Этот код начинается с настройки матрицы проектирования, состоящей из приблизительно равномерно распределенных независимых значений, которые приблизительно ортогональны (чтобы мы не сталкивались с проблемами мультиколлинеарности). Он вычисляет одиночное квадратичное (то есть нелинейное) взаимодействие между первыми двумя переменными: это только один из многих видов «нелинейностей», которые можно изучить, но, по крайней мере, это распространенная, хорошо понятая. Затем он стандартизирует все так, чтобы коэффициенты были сопоставимы:

set.seed(41)
p <- 7                                            # Dimensions
n <- 2^p                                          # Observations
x <- as.matrix(do.call(expand.grid, lapply(as.list(1:p), function(i) c(-1,1))))
x <- x + runif(n*p, min=-1, max=1)
x <- cbind(x, x.12 = x[,1]*x[,2])                 # The nonlinear part
x <- apply(x, 2, function(y) (y - mean(y))/sd(y)) # Standardization

Для базовой модели OLS (без нелинейности) необходимо указать некоторые коэффициенты и стандартное отклонение остаточной ошибки. Вот набор единичных коэффициентов и сопоставимый SD:

beta <- rep(c(1,-1), p)[1:p]
sd <- 1

$1/4$$1$$-1$

gamma = 1/4          # The standardized interaction term
df <- data.frame(x)
df$y <- x %*% c(beta, gamma) + rnorm(n, sd=sd)
summary(df)
cor(df)*100
plot(df, lower.panel=function(x,y) lines(lowess(y~x)), 
     upper.panel=function(x,y) points(x,y, pch=".", cex=4))
summary(lm(df$y ~ x))

Вместо того, чтобы просмотреть все выходные данные, давайте посмотрим на эти данные, используя выходные данные plotкоманды:

Трассы низких значений в нижнем треугольнике показывают, по существу, отсутствие линейных отношений между взаимодействие ( x.12) и зависимой переменной ( y) и скромные линейные отношения между другими переменными и y. Результаты OLS подтверждают это; взаимодействие едва ли значимо:

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.0263     0.0828    0.32    0.751    
xVar1         0.9947     0.0833   11.94   <2e-16 ***
xVar2        -0.8713     0.0842  -10.35   <2e-16 ***
xVar3         1.0709     0.0836   12.81   <2e-16 ***
xVar4        -1.0007     0.0840  -11.92   <2e-16 ***
xVar5         1.0233     0.0836   12.24   <2e-16 ***
xVar6        -0.9514     0.0835  -11.40   <2e-16 ***
xVar7         1.0482     0.0835   12.56   <2e-16 ***
xx.12         0.1902     0.0836    2.27    0.025 *  

Я возьму p-значение члена взаимодействия как критерий нелинейности: когда это p-значение достаточно низкое (вы можете выбрать, насколько низко), мы обнаружим нелинейность.

(Здесь есть тонкость в том, что именно мы ищем. На практике нам может потребоваться изучить все 7 * 6/2 = 21 возможных таких квадратичных взаимодействий, а также, возможно, еще 7 квадратичных терминов, вместо того, чтобы сосредоточиться на одном термине как это сделано здесь. Мы хотели бы сделать исправление для этих 28 взаимосвязанных тестов. Я не делаю здесь этого явного исправления, потому что вместо этого я отображаю смоделированное распределение р-значений. Вы можете прочитать показатели обнаружения непосредственно из гистограммы в конце на основе ваших порогов значимости.)

Но давайте не будем делать этот анализ только один раз; давайте делать это много раз, генерируя новые значения yв каждой итерации в соответствии с той же моделью и той же самой матрицей проекта. Для этого мы используем функцию для выполнения одной итерации и возвращаем p-значение члена взаимодействия:

test <- function(gamma, sd=1) {
  y <- x %*% c(beta, gamma) + rnorm(n, sd=sd)
  fit <- summary(lm(y ~ x))
  m <- coef(fit)
  n <- dim(m)[1]
  m[n, 4]
}

Я предпочитаю представлять результаты моделирования в виде гистограмм значений p, варьирующих стандартизированный коэффициент gammaчлена взаимодействия. Сначала гистограммы:

h <- function(g, n.trials=1000) {
  hist(replicate(n.trials, test(g, sd)), xlim=c(0,1), 
       main=toString(g), xlab="x1:x2 p-value")
}
par(mfrow=c(2,2)) # Draw a 2 by 2 panel of results

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

temp <- sapply(2^(-3:0) * gamma, h)

Результаты:

xsdbeta$1/4$$1/8$$1/16$gamma$1/2$

$1/32$$1/4$xsdbetasd

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

Не уверен, что это дает окончательный ответ на вопрос, но я бы посмотрел на это . Особенно пункт 2. См. Также обсуждение в приложении А2 документа .