Dificuldade de testar a linearidade na regressão

Em Modelagem Estatística: As Duas Culturas Leo Breiman escreve

A prática atual aplicada é verificar o ajuste do modelo de dados usando testes de ajuste de qualidade e análise residual. Em um ponto, alguns anos atrás, montei um problema de regressão simulado em sete dimensões com uma quantidade controlada de não linearidade. Testes padrão de qualidade do ajuste não rejeitaram a linearidade até que a não linearidade fosse extrema.

Breiman não fornece os detalhes de sua simulação. Ele faz referência a um trabalho que, segundo ele, fornece justificativa teórica para sua observação, mas o trabalho não é publicado.

Alguém viu um resultado de simulação publicado ou artigo teórico para apoiar a alegação de Brieman?

Criei uma simulação que responderia à descrição de Breiman e encontrei apenas o óbvio: o resultado depende do contexto e do significado de "extremo".

Muita coisa pode ser dita, mas deixe-me limitar a apenas um exemplo conduzido por meio de Rcódigo facilmente modificado para que os leitores interessados ​​usem em suas próprias investigações. Esse código começa configurando uma matriz de design que consiste em valores independentes aproximadamente uniformemente distribuídos que são aproximadamente ortogonais (para que não entremos em problemas de multicolinearidade). Ele calcula uma única interação quadrática (isto é, não linear) entre as duas primeiras variáveis: esse é apenas um dos muitos tipos de "não linearidades" que poderiam ser estudados, mas pelo menos é comum e bem compreendido. Depois, padroniza tudo para que os coeficientes sejam comparáveis:

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

Para o modelo OLS base (sem não linearidade), devemos especificar alguns coeficientes e o desvio padrão do erro residual. Aqui está um conjunto de coeficientes de unidade e um SD comparável:

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))

Em vez de percorrer toda a saída aqui, vamos ver esses dados usando a saída do plotcomando:

Os traços inferiores no triângulo inferior mostram essencialmente nenhuma relação linear entre a interação ( x.12) e a variável dependente ( y) e modestas relações lineares entre as outras variáveis ​​e y. Os resultados do OLS confirmam isso; a interação é pouco significativa:

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 *  

Tomarei o valor p do termo de interação como um teste de não linearidade: quando esse valor p for suficientemente baixo (você pode escolher o quão baixo), teremos detectado a não linearidade.

(Há uma sutileza aqui sobre o que exatamente procuramos. Na prática, talvez seja necessário examinar todas as 7 * 6/2 = 21 possíveis interações quadráticas possíveis, bem como talvez mais 7 termos quadráticos, em vez de focar em um único termo como é feito aqui. Gostaríamos de fazer uma correção para esses 28 testes inter-relacionados. Não faço explicitamente essa correção aqui, porque em vez disso, mostro a distribuição simulada dos valores-p. Você pode ler as taxas de detecção diretamente em os histogramas no final, com base nos seus limites de significância.)

Mas não vamos fazer essa análise apenas uma vez; vamos fazer isso muitas vezes, gerando novos valores yem cada iteração de acordo com o mesmo modelo e a mesma matriz de design. Para fazer isso, usamos uma função para realizar uma iteração e retornar o valor p do termo de interação:

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]
}

Escolho apresentar os resultados da simulação como histogramas dos valores p, variando o coeficiente padronizado gammado termo de interação. Primeiro, os histogramas:

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

Agora, faça o trabalho. Demora alguns segundos para 1000 tentativas por simulação (e quatro simulações independentes, começando com o valor determinado do termo de interação e diminuindo pela metade sucessivamente):

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

Os resultados:

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

$1/32$$1/4$xsdbetasd

Em resumo, uma simulação como essa pode provar o que você quiser se você apenas a configurar e interpretar da maneira certa. Isso sugere que o estatístico individual deve conduzir suas próprias explorações, adequadas aos problemas específicos que enfrentam, a fim de obter uma compreensão pessoal e profunda das capacidades e fraquezas dos procedimentos que estão usando.

Não tenho certeza de que dê uma resposta final à pergunta, mas eu daria uma olhada nisso . Especialmente o ponto 2. Ver também a discussão no apêndice A2 do artigo .