Schwierigkeiten beim Testen der Linearität in der Regression

In der statistischen Modellierung: Die zwei Kulturen schreibt Leo Breiman

Die derzeitige angewandte Praxis besteht darin, die Anpassung des Datenmodells mithilfe von Anpassungstests und Restanalyse zu überprüfen. Vor einigen Jahren habe ich einmal ein simuliertes Regressionsproblem in sieben Dimensionen mit einem kontrollierten Maß an Nichtlinearität erstellt. Standardtests der Anpassungsgüte lehnten die Linearität erst ab, wenn die Nichtlinearität extrem war.

Breiman gibt keine Details seiner Simulation an. Er bezieht sich auf einen Artikel, der seiner Ansicht nach theoretisch gerechtfertigt ist, aber der Artikel ist unveröffentlicht.

Hat jemand ein veröffentlichtes Simulationsergebnis oder eine theoretische Abhandlung gesehen, die Briemans Behauptung stützt?

Ich erstellte eine Simulation, die Breimans Beschreibung entsprach, und fand nur das Offensichtliche: Das Ergebnis hängt vom Kontext ab und davon, was mit "extrem" gemeint ist.

Man könnte eine Menge sagen, aber lassen Sie mich dies auf ein Beispiel beschränken, das mit Hilfe eines leicht zu modifizierenden RCodes erstellt wurde, den interessierte Leser für ihre eigenen Untersuchungen verwenden können. Dieser Code beginnt mit der Erstellung einer Entwurfsmatrix, die aus ungefähr gleichmäßig verteilten unabhängigen Werten besteht, die ungefähr orthogonal sind (damit wir nicht auf Multikollinearitätsprobleme eingehen). Es berechnet eine einzige quadratische (dh nichtlineare) Wechselwirkung zwischen den ersten beiden Variablen: Dies ist nur eine von vielen Arten von "Nichtlinearitäten", die untersucht werden könnten, aber zumindest eine allgemein bekannte. Dann standardisiert es alles so, dass die Koeffizienten vergleichbar sind:

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

Für das Basis-OLS-Modell (ohne Nichtlinearität) müssen einige Koeffizienten und die Standardabweichung des Restfehlers angegeben werden. Hier ist eine Reihe von Einheitskoeffizienten und eine vergleichbare 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))

Anstatt die gesamte Ausgabe hier durchzugehen, schauen wir uns diese Daten mit der Ausgabe des plotBefehls an:

Die geringen Spuren im unteren Dreieck zeigen im Wesentlichen keine lineare Beziehung zwischen der Wechselwirkung ( x.12) und der abhängigen Variablen ( y) und bescheidene lineare Beziehungen zwischen den anderen Variablen und y. Die OLS-Ergebnisse bestätigen dies; Die Wechselwirkung ist kaum signifikant:

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 *  

Ich nehme den p-Wert des Interaktionsterms als Test für die Nichtlinearität: Wenn dieser p-Wert ausreichend niedrig ist (Sie können wählen, wie niedrig er sein soll), haben wir die Nichtlinearität erkannt.

(Hier gibt es eine subtile Frage, wonach genau wir suchen. In der Praxis müssen wir möglicherweise alle 7 * 6/2 = 21 möglichen quadratischen Wechselwirkungen sowie vielleicht 7 weitere quadratische Terme untersuchen, anstatt uns auf einen einzelnen Term zu konzentrieren Wie hier beschrieben, möchten wir eine Korrektur für diese 28 miteinander in Zusammenhang stehenden Tests vornehmen. Diese Korrektur nehme ich hier nicht explizit vor, da stattdessen die simulierte Verteilung der p-Werte angezeigt wird. Sie können die Erkennungsraten direkt ablesen die Histogramme am Ende basieren auf Ihren Signifikanzschwellen.)

Aber lassen Sie uns diese Analyse nicht nur einmal durchführen. Lass es uns viele Male tun, indem wir yin jeder Iteration neue Werte nach demselben Modell und derselben Entwurfsmatrix generieren . Um dies zu erreichen, verwenden wir eine Funktion, um eine Iteration durchzuführen und den p-Wert des Interaktionsterms zurückzugeben:

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

Ich beschließe, die Simulationsergebnisse als Histogramme der p-Werte darzustellen, wobei der standardisierte Koeffizient gammades Interaktionsterms variiert wird . Zunächst die Histogramme:

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

Nun, um die Arbeit zu erledigen. Für 1000 Versuche pro Simulation (und vier unabhängige Simulationen, die mit dem angegebenen Wert des Interaktionsbegriffs beginnen und ihn jedes Mal nacheinander halbieren) sind einige Sekunden erforderlich:

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

Die Ergebnisse:

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

$1/32$$1/4$xsdbetasd

Kurz gesagt, eine Simulation wie diese kann beweisen, was Sie wollen, wenn Sie sie nur einrichten und richtig interpretieren. Das bedeutet, dass der einzelne Statistiker seine eigenen Untersuchungen durchführen sollte, die auf die jeweiligen Probleme zugeschnitten sind, um ein persönliches und tiefes Verständnis für die Fähigkeiten und Schwächen der von ihm verwendeten Verfahren zu erlangen.

Nicht sicher , es gibt eine endgültige Antwort auf die Frage, aber ich würde einen Blick auf gibt diese . Insbesondere Punkt 2. Siehe auch die Diskussion in Anhang A2 des Papiers .