Wiederherstellung von Rohkoeffizienten und Varianzen aus der orthogonalen Polynomregression

Es scheint, als hätte ich ein Regressionsmodell wie $y_i \sim \beta_0 + \beta_1 x_i+\beta_2 x_i^2 +\beta_3 x_i^3$ Ich kann entweder ein Rohpolynom anpassen und unzuverlässige Ergebnisse erhalten oder ein orthogonales Polynom anpassen und Koeffizienten erhalten, die keine direkte physikalische Interpretation haben (z. B. kann ich sie nicht verwenden, um die Orte der Extrema auf der ursprünglichen Skala zu finden). Es scheint, als ob ich das Beste aus beiden Welten haben und die angepassten orthogonalen Koeffizienten und ihre Varianzen zurück in die Rohskala transformieren könnte. Ich habe einen Aufbaustudiengang in angewandter linearer Regression (mit Kutner, 5ed) absolviert und das Kapitel über polynomielle Regression in Draper (3ed, von Kutner referiert) durchgesehen, aber keine Diskussion darüber gefunden, wie dies zu tun ist. Der Hilfetext für diepoly()Funktion in R nicht. Ich habe auch nichts in meiner Websuche gefunden, auch nicht hier. Rekonstruiert rohe Koeffizienten (und erhält ihre Varianzen) aus Koeffizienten, die an ein orthogonales Polynom angepasst sind?

unmöglich zu tun, und ich verschwende meine Zeit.
Möglicherweise möglich, aber nicht bekannt, wie im allgemeinen Fall.
möglich, aber nicht besprochen, weil "wer würde wollen?"
möglich, aber nicht diskutiert, weil "es ist offensichtlich".

Wenn die Antwort 3 oder 4 ist, wäre ich sehr dankbar, wenn jemand die Geduld hätte, dies zu erklären oder auf eine Quelle zu verweisen, die dies tut. Wenn es 1 oder 2 ist, wäre ich immer noch neugierig, was das Hindernis ist. Vielen Dank für das Lesen und ich entschuldige mich im Voraus, wenn ich etwas Offensichtliches übersehen habe.

regression linear-model regression-coefficients polynomial f1r3br4nd
quelle

Ich verstehe deine Punkte nicht. x, x

und x

sind nicht orthogonal. Daher sind sie korreliert und die Regressionsparameter könnten instabil sein, aber es ist nicht automatisch der Fall, dass sie unzuverlässig sind. Die Umwandlung in orthognonale Polynome kann zuverlässiger sein. Aber was macht den Koeffizienten der ursprünglichen Potenzen von x deutbarer als die Koeffizienten der orthogonalen Polynome? Wenn x die einzige Variable ist, wie im Modell y = a + bx, dann ist ∆y = yi-yi-1 = b∆x und b als Änderung in y pro Änderungseinheit in x interpretierbar. Aber mit den involvierten Kräften geht eine solche Interpretation verloren.

^{2}

$^2$

^{3}

$^3$

Michael R. Chernick

Ich habe der Einfachheit halber ein Modell mit nur x als Variable verwendet, aber in Wirklichkeit vergleiche ich Kurven zwischen Behandlungsgruppen. Je nachdem, welche Terme signifikant sind und wie groß sie sind, kann ich sie interpretieren - zum Beispiel als Aufwärts- / Abwärtsverschiebung insgesamt oder als größere / kleinere Anfangssteigung. Wie meine Frage sagt, ist ein natürlicher Vergleich zwischen Kurven die Position der Maxima / Minima, die leichter zu interpretieren sind, wenn sie auf der ursprünglichen Skala liegen. Also, Ihre Stimme ist für Wahl 3, nehme ich es?

f1r3br4nd

Nein, ich habe noch nicht herausgefunden, ob es möglich ist oder nicht. Ich habe gerade verstanden, warum Sie es tun wollen.

Michael R. Chernick

Beachten Sie, dass die Modellanpassung mit orthogonalen Polynomen genau dieselbe Anpassung (dh dasselbe

, dieselben Anpassungswerte usw.) wie die Modellanpassung mit den rohen Polynomtermen aufweist. Wenn Sie dies also auf die ursprünglichen Daten zurückführen möchten, können Sie die Koeffizienten für die Rohterme betrachten, aber die orthogonalen Polynome verwenden, um die Inferenz für die einzelnen Terme auf eine Weise zu erstellen, die die Abhängigkeit zwischen ihnen "berücksichtigt" .

R^{2}

$R^2$

Makro

Wie sich herausstellt, sind kubische Splines und B-Splines eine Klasse für sich und das Beste aus zwei Welten.

Carl

Antworten:

Ja es ist möglich.

Sei der nicht konstante Teil der aus berechneten orthogonalen Polynome . (Jeder ist ein Spaltenvektor.) Wenn man diese gegen das muss eine perfekte Anpassung geben. Sie können dies mit der Software ausführen, auch wenn die Prozeduren zur Berechnung orthogonaler Polynome nicht dokumentiert sind. Die Regression von $z_1, z_2, z_3$ $x_i$ $x_i$ $z_j$ ergibt Koeffizienten für die $\gamma_{ij}$

z_{i j} = γ_{j 0} + x_{i} γ_{j 1} + x_{i}^{2} γ_{j 2} + x_{i}^{3} γ_{j 3} .

$z_{ij} = \gamma_{j0} + x_i\gamma_{j1} + x_i^2\gamma_{j2} + x_i^3\gamma_{j3}.$

Das Ergebnis ist eine Matrix , die bei rechter Multiplikation die Entwurfsmatrix in $4\times 4$ $\Gamma$ $X=\pmatrix{1;&x;&x^2;&x^3}$

\begin{matrix} (1) & Z = (\begin{matrix} 1; & z_{1}; & z_{2}; & z_{3} \end{matrix}) = X Γ . \end{matrix}

$Z=\pmatrix{1;&z_1;&z_2;&z_3} = X\Gamma.\tag{1}$

Nach dem Einbau des Modells

E (Y) = Z β

$\mathbb{E}(Y) = Z\beta$

$\hat\beta$ $(1)$

\hat{Y} = Z \hat{β} = (X Γ) \hat{β} = X (Γ \hat{β}) .

$\hat Y = Z\hat\beta = (X\Gamma)\hat\beta = X(\Gamma\hat\beta).$

$\Gamma\hat\beta$ $x$ .

Der folgende RCode veranschaulicht diese Verfahren und testet sie mit synthetischen Daten.

n <- 10        # Number of observations
d <- 3         # Degree
#
# Synthesize a regressor, its powers, and orthogonal polynomials thereof.
#
x <- rnorm(n)
x.p <- outer(x, 0:d, `^`); colnames(x.p) <- c("Intercept", paste0("x.", 1:d))
z <- poly(x, d)
#
# Compute the orthogonal polynomials in terms of the powers via OLS.
#
xform <- lm(cbind(1, z) ~ x.p-1)
gamma <- coef(xform)
#
# Verify the transformation: all components should be tiny, certainly
# infinitesimal compared to 1.
#
if (!all.equal(as.vector(1 + crossprod(x.p %*% gamma - cbind(1,z)) - 1), 
    rep(0, (d+1)^2)))
  warning("Transformation is inaccurate.")
#
# Fit the model with orthogonal polynomials.
#
y <- x + rnorm(n)
fit <- lm(y ~ z)
#summary(fit)
#
# As a check, fit the model with raw powers.
#
fit.p <- lm(y ~ .-1, data.frame(x.p))
#summary(fit.p)
#
# Compare the results.
#
(rbind(Computed=as.vector(gamma %*% coef(fit)), Fit=coef(fit.p)))

if (!all.equal(as.vector(gamma %*% coef(fit)), as.vector(coef(fit.p))))
  warning("Results were not the same.")

whuber
quelle

Γ

$\Gamma$

1

$1$

10^{- 16}

$10^{-16}$

1

$1$

Zwei Jahre später ... @whuber, ist es möglich, dies auch auf die 95% CIs der Koeffizienten auszudehnen?

user2602640

@ user2602640 Ja. Sie müssen die Varianz-Kovarianz-Matrix der Koeffizienten extrahieren (verwenden vcovin R), um auf einer Basis berechnete Varianzen in Varianzen auf der neuen Basis umzuwandeln, und dann die CIs auf die übliche Weise manuell berechnen.

whuber

@whuber Ich bin deinem Kommentar ungefähr auf halbem Weg gefolgt und habe dich dann völlig verloren ... Gibt es eine Chance, dass du Mitleid mit einem mathematisch herausgeforderten Biologen hast und ihn in Code schreibst?

user2602640