So passen Sie ein stückweise konstantes Modell (oder ein Schrittfunktionsmodell) an und vergleichen es mit dem logistischen Modell in R.

7

Ich habe x, y-Daten, wobei x die Position (entlang eines Transekts) und y eine kontinuierliche Variable (z. B. Temperatur) ist.

x<-c(115,116,117,118,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151)
y<-c(68.3864344929207,69.3317985238957,70.5644430507252,68.7245413735345,68.2444929220624,69.2335442981487,68.9683874048675,69.7104166614154,70.6047681836632,71.1807115091657,70.3333333333333,70,69.735105818315,69.6487585375593,70,69.4943374527374,69,70.3590104484011,70.283981899468,68.9087146344906,68.6666666666667,68.5232663969658,68.088410889283,67.567883666713,66.9865494087022,66,66.3333333333333,66.0165138638852,65.6666666666667,65.7465975732667,66.3333333333333,66.6579079386334)

Ich möchte ein lineares Modell, ein 4-Parameter-Logistikmodell und eine stückweise konstante Funktion (auch bekannt als Schrittfunktion) mit einem Haltepunkt an diese Daten anpassen und Modelle vergleichen, um festzustellen, ob sich y allmählich, abrupter als linear oder abrupt in einem Schritt ändert. wie Mode.

(HINWEIS: Das Zeichnen der Daten zeigt, dass das Logistikmodell in diesem Fall am besten geeignet ist, aber ich muss letztendlich viele verschiedene abhängige Variablen ohne diese Informationen im Voraus durchlaufen.)

Ich verwende nls () mit manueller Formel (y ~ int + Steigung) für das lineare Modell und SSfpl für die 4-Parameter-Logistik.

l<-nls(y~int+x*slope,data=data.frame(cbind(x,y)),start=list(int=0,slope=1)) # linear model 
fl<-nls(y~SSfpl(x,A,B,xmid,scal),data=data.frame(cbind(x,y))) # four parameter logistic

Mein Problem ist die Anpassung eines stückweise konstanten Modells. Es scheinen Funktionen für ein stückweise lineares Modell vorhanden zu sein, aber nicht viele für eine Schrittfunktion. Ich fand:

library(cumSeg)
pc<-jumpoints(y,k=1,output="1")

Das gibt mir den Index eines optimalen Haltepunkts in x:

breakpt<-pc$psi

Ich bin mir jedoch nicht sicher, wie ich die Ausgabe anders mit den linearen oder logistischen Modellen vergleichen soll. Ich denke zum Beispiel, ich muss einen Modellvergleich mit AIC durchführen, bin mir aber nicht sicher, wie ich dies basierend auf der Klasse des von jumpoints () zurückgegebenen Objekts berechnen soll. Daher habe ich manuell versucht, das Modell mit dem optimalen Haltepunkt zu generieren:

pcm<-lm(y~I(x<=x[breakpt])+I(x>x[breakpt]))

Jetzt kann ich AIC-Werte für alle Modelle erhalten. Aber ich habe festgestellt, dass, wenn ich die Anpassung meiner manuellen Schrittfunktion (pcm) an die in cumSeg (pc) erzeugte Darstellung zeichne, die Anpassung meines manuellen Modells nicht so gut aussieht ...

par(mfrow=c(1,2))
plot(x,y)
lines(x,predict(cpm) # versus
plot(pc)

Sie fragen sich, warum dies geschieht?

Und auch nach alternativen Vorschlägen suchen, um das Gesamtproblem anzugehen. Auf der Modellanpassungsseite: Das Hauptproblem besteht darin, wie ich feststellen kann, ob eine Variable gestuft ist oder nicht, während sie sich ansonsten mehr oder weniger linear ändert. (Ich kann ein Logistikmodell haben, das besser passt als ein lineares Modell, aber bei xmid ähnliche Steigungen aufweist ... das ist also keine wirkliche Antwort).

Zum Ende des Modellvergleichs: Ich verwende die AICc () - Funktion aus dem MuMIn-Paket, um die AIC-Werte für die verschiedenen Modelle abzurufen, und die akaike.weights () -Funktion aus dem qpcR-Paket, um den Delta-AIC- "Test" durchzuführen. ... aber vielleicht gibt es noch andere Ideen?

Vielen Dank! Sehr geschätzt!

user2414840
quelle
@ user241480 Eine Sache, bei der mir der Thread in Ihrer Beschreibung fehlt, ist die Verwendung eines Logistikmodells. Die von Ihnen geposteten Daten haben keine Binärwerte y! Könnten Sie ein wenig Ihre Verwendung von Logistikmodellen beschreiben, damit ich möglicherweise einige Ratschläge zum Vergleich Ihrer Modelle mit meiner Antwort ändern kann?
Matthew Drury
Ich habe diese Referenz durch Suchen nach "Vier-Parameter-Logistikmodell" gefunden. Ist das richtig? psg.hitachi-solutions.com/masterplex/blog/…
Matthew Drury
Warum verwenden Sie nls anstelle von lm für das lineare Modell? Sie können einen Klassifizierungs- und Regressionsbaum (CART) für den stückweise konstanten Teil verwenden. Der Paketbaum bietet eine Implementierung in R.
Shane

Antworten:

13

Ich denke, eine der Hauptschwierigkeiten besteht darin, dass stückweise konstante Regressionen normalerweise nicht als "stückweise konstante Regressionen" bezeichnet werden. Sie werden normalerweise Regressionsbäume genannt , was ein schöner visueller Name ist, aber nicht besonders gut, wenn Sie nicht bereits wissen, wie die Leute sie nennen! Sie können in Rdas integrierte rpart- Paket integriert werden (ich glaube, es rpartsteht für "rekursive Partitionierung", was unser dritter Name für dasselbe Konzept ist).

Hier ist rpartin Aktion auf Ihre Daten:

df <- data.frame(x=x, y=y)
tree <- rpart(y ~ x, data=df)

Ich habe eine kleine plot_treeFunktion geschrieben, die die Vorhersagen zeigt

plot_tree <- function(tree, x, y) {
  s <- seq(110, 155, by=.5)
  plot(x, y)
  lines(s, predict(tree, data.frame(x=s)))
}

Wenn es auf den Standardbaum Ihrer Daten angewendet wird, sieht es so aus

plot_tree(tree, x, y)

Geben Sie hier die Bildbeschreibung ein

Sie können die Granularität der Anpassung an Ihre Daten mithilfe von steuern rpart.controll

tree <- rpart(y ~ x, data=df, control=rpart.control(minsplit=5, cp=.0001))
plot_tree(tree, x, y)

Geben Sie hier die Bildbeschreibung ein

Matthew Drury
quelle
Vielen Dank für den Vorschlag Matthew Drury. Ich weiß nicht viel über Regressionsbäume, aber dies scheint eine ziemlich andere Herangehensweise an das Problem zu sein. Was ich jedoch anpassen möchte, ist eine einstufige Funktion: von einer konstanten Steigung zur nächsten und in einem Rahmen, der es mir ermöglicht, die Anpassung mit linearen und nichtlinearen Modellen zu vergleichen. Ich sehe nicht wirklich, wie man erzwingt ein einzelner Schritt in rpart (ohne möglicherweise den cp-Parameter in rpart.control zu optimieren) oder wie man schließlich den AIC berechnet ...
user2414840
Es sollte einen Parameter geben, den Sie für die maximale Tiefe des Baums oder die maximale Anzahl von Teilungen festlegen können. Beide sollten in der Lage sein, einen Baum mit einer Diskontinuität zu erhalten.
Matthew Drury
Nur um nachzufolgen ... Ich habe das noch nicht gelöst. Es scheint keine Möglichkeit zu geben, eine maximale Anzahl von Teilungen in den Einstellungen von rpart.control festzulegen ...
rpart.control