Trudność testowania liniowości w regresji

W Modelowaniu statystycznym: The Two Cultures pisze Leo Breiman

Obecnie stosowaną praktyką jest sprawdzanie dopasowania modelu danych za pomocą testów dopasowania i analizy resztkowej. W pewnym momencie, kilka lat temu, stworzyłem symulowany problem regresji w siedmiu wymiarach z kontrolowaną nieliniowością. Standardowe testy dobroci dopasowania nie odrzucały liniowości, dopóki nieliniowość nie była ekstremalna.

Breiman nie podaje szczegółów swojej symulacji. Odwołuje się do artykułu, który, jak twierdzi, podaje teoretyczne uzasadnienie jego obserwacji, ale artykuł nie został opublikowany.

Czy ktoś widział opublikowany wynik symulacji lub artykuł teoretyczny na poparcie twierdzenia Briemana?

Stworzyłem symulację, która byłaby odpowiedzią na opis Breimana i znalazłem tylko to, co oczywiste: wynik zależy od kontekstu i tego, co należy rozumieć przez „ekstremalność”.

Można powiedzieć bardzo dużo, ale ograniczę się do jednego przykładu przeprowadzonego za pomocą łatwo modyfikowalnego Rkodu, który zainteresowani czytelnicy mogą wykorzystać we własnych badaniach. Ten kod zaczyna się od ustawienia macierzy projektowej składającej się z w przybliżeniu równomiernie rozłożonych niezależnych wartości, które są w przybliżeniu ortogonalne (abyśmy nie wpadli w problemy z wielokoliniowością). Oblicza pojedynczą interakcję kwadratową (tj. Nieliniową) między dwiema pierwszymi zmiennymi: jest to tylko jeden z wielu rodzajów „nieliniowości”, które można badać, ale przynajmniej jest to powszechna, dobrze rozumiana. Następnie standaryzuje wszystko, aby współczynniki były porównywalne:

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

W przypadku podstawowego modelu OLS (bez nieliniowości) musimy określić niektóre współczynniki i odchylenie standardowe błędu resztkowego. Oto zestaw współczynników jednostkowych i porównywalnej 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))

Zamiast przedzierać się przez cały wyjście tutaj, niech spojrzeć na te dane z wykorzystaniem wyjście plotpolecenia:

Ślady lowess w dolnym trójkącie zasadniczo nie wykazują liniowej zależności między interakcją ( x.12) a zmienną zależną ( y) oraz skromne relacje liniowe między innymi zmiennymi i y. Wyniki OLS to potwierdzają; interakcja ma niewielkie znaczenie:

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 *  

Przyjmę wartość p terminu interakcji za test nieliniowości: gdy ta wartość p jest wystarczająco niska (możesz wybrać, jak niska), wykryjemy nieliniowość.

(Jest tu subtelność dotycząca tego, czego dokładnie szukamy. W praktyce może być konieczne zbadanie wszystkich 7 * 6/2 = 21 możliwych takich kwadratowych interakcji, a także być może 7 bardziej kwadratowych warunków, zamiast skupiania się na pojedynczym terminie jak to tutaj zrobiono. Chcielibyśmy wprowadzić korektę dla tych 28 powiązanych ze sobą testów. Nie dokonuję tutaj tej korekty wprost, ponieważ zamiast tego wyświetlam symulowany rozkład wartości p. Możesz odczytać wskaźniki wykrywalności bezpośrednio z histogramy w końcu na bazie swoich progów istotności).

Ale nie róbmy tej analizy tylko raz; zróbmy to wiele razy, generując nowe wartości yw każdej iteracji według tego samego modelu i tej samej matrycy projektowej. Aby to osiągnąć, używamy funkcji do przeprowadzenia jednej iteracji i zwrócenia wartości p terminu interakcji:

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

Zdecydowałem się przedstawić wyniki symulacji jako histogramy wartości p, zmieniając znormalizowany współczynnik gammaterminu interakcji. Po pierwsze, histogramy:

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

Teraz do pracy. 1000 prób na symulację zajmuje kilka sekund (i cztery niezależne symulacje, zaczynając od podanej wartości terminu interakcji i sukcesywnie zmniejszając go o połowę):

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

Wyniki:

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

$1/32$$1/4$xsdbetasd

Krótko mówiąc, taka symulacja może udowodnić, co chcesz, jeśli tylko ją skonfigurujesz i zinterpretujesz we właściwy sposób. Sugeruje to, że indywidualny statystyk powinien przeprowadzić własne eksploracje, odpowiednie do konkretnych problemów, z którymi się borykają, aby uzyskać osobiste i głębokie zrozumienie możliwości i słabości stosowanych procedur.

Nie wiesz, że daje ostatecznej odpowiedzi na pytanie, ale chciałbym dać do obejrzenia tego . Zwłaszcza punkt 2. Zobacz także dyskusję w załączniku A2 do artykułu .