Wir brauchen bessere Startwerte. Passen Sie ein Nicht-Plateau-Modell, model0, an und verwenden Sie die Parameter daraus, um alle Datenpunkte anzupassen, die das Modell ergeben. Verwenden Sie dann a und b daraus und ein Wertegitter für clx (aufgrund seiner Problematik), das model.Ab und model.La. (Beachten Sie, dass es nicht möglich ist, Anpassungen von einigen Startwerten des Rasters zu erzeugen, was zu Fehlermeldungen führt. Nls2 verarbeitet jedoch weiterhin weitere Startwerte, sodass diese Fehler ignoriert werden können.)

  Bibliothek (nls2)

# Stellen Sie sicher, dass die Daten zum Plotten sortiert sind
o <- mit (Heilung, Ordnung (Typ, Tage))
h <-Heilung [o,]

Das letzte Argument # gibt an, ob ein Plateau vorhanden ist oder nicht
Quadplat = Funktion (x, a, b, clx, plat = TRUE) {
  if (plat) x <-pmin (x, clx)
  a + b * x + (-0,5 * b / clx) * x * x
}}

# Kein Plateau-Modell mit allen Daten versehen
st <c (a = 1, b = 1, clx = 1)
model0 <nls (Area ~ quadplat (Tage, a, b, clx, FALSE), h, start = st)

# passt alle Datenmodell
Modell <nls (Fläche ~ Quadplat (Tage, a, b, clx), h, Start = Coef (Modell 0))
co <-coef (Modell)
 

Wir können jetzt die Teilmengenmodelle mithilfe der oben in den Startwerten berechneten Werte anpassen und zeichnen.

  if (existiert ("model.Ab")) rm (model.Ab)
model.Ab <nls2 (Fläche ~ Quadplat (Tage, a, b, clx), h, Teilmenge = h $ Typ == "Abrieb",
  start = data.frame (a = co [[1]], b = co [[2]], clx = 0: 140))

if (existiert ("model.La")) rm (model.La)
model.La <nls2 (Fläche ~ Quadplat (Tage, a, b, clx), h, Teilmenge = h $ Typ == "Schnittwunde",
  start = data.frame (a = co [[1]], b = co [[2]], clx = 0: 140))
  
cols <c (Abrieb = "rot", Schnittwunde = "blau")
Diagramm (Fläche ~ Tage, h, col = cols [Typ], pch = 20, cex = 1,5)
Linien (angepasst (Modell.Ab) ~ Tage, Teilmenge (h, Typ == "Abrieb"),
  col = cols ["Abrieb"])
Linien (angepasst (Modell.La) ~ Tage, Teilmenge (h, Typ == "Schnittwunde"),
  col = cols ["Schnittwunde"])
 

(Fortsetzung nach Grafik)

Alternatives Modell

Wenn es in Ordnung ist, andere Modelle in Betracht zu ziehen, hat dieses Modell nur zwei Parameter, ist einfacher anzupassen und hat trotz weniger Parameter geringere Restquadratsummen.

  model.Ab2 <nls (Bereich ~ a * (1 - exp (- b * Tage)), h,
   Teilmenge = Typ == "Abrieb", Start = c (a = 100, b = 0,1))
 
model.La2 <nls (Fläche ~ a * (1 - exp (- b * Tage)), h,
   Teilmenge = Typ == "Schnittwunde", Start = c (a = 100, b = 0,1))

# Handlung
cols <c (Abrieb = "rot", Schnittwunde = "blau")
Diagramm (Fläche ~ Tage, h, col = cols [Typ], pch = 20, cex = 1,5)
Linien (angepasst (Modell.Ab2) ~ Tage, Teilmenge (h, Typ == "Abrieb"),
  col = cols ["Abrieb"])
Linien (angepasst (Modell.La2) ~ Tage, Teilmenge (h, Typ == "Schnittwunde"),
  col = cols ["Schnittwunde"])
 

(Fortsetzung nach Grafik)

Ein Parametermodell

Wenn wir im 2-Parameter-Modell des letzten Abschnitts a = 100 festlegen, erhalten wir ein 1-Parameter-Modell, das statistisch nicht vom 2-Parameter-Modell unterscheidbar ist. Dies geht aus dem in den Anovas gezeigten p-Wert hervor, der größer als 0,05 ist, was darauf hinweist, dass wir die Nullhypothese, dass die 1- und 2-Parametermodelle die Daten für jede der beiden Teilmengen gleich gut beschreiben, nicht ablehnen können.

  model.Ab3 <nls (Fläche ~ 100 * (1 - exp (- b * Tage)), h,
   Teilmenge = Typ == "Abrieb", Start = c (b = .1))
 
model.La3 <nls (Fläche ~ 100 * (1 - exp (- b * Tage)), h,
   Teilmenge = Typ == "Schnittwunde", Start = c (b = 0,1))

anova (model.Ab3, model.Ab2)

anova (model.La3, model.La2)
 

Beachten Sie auch, dass der Punkt, an dem y = 95 erreicht wird, d. h. in der Nähe des Plateaus, -log (1 - 95/100) / b ist (basierend auf der Invertierung der Modellgleichung). Der Zähler ist ungefähr 3, also erreicht er 95 bei ungefähr 3 / b .

Andere

Wenn m <nls (...) ist, gibt die Zusammenfassung (m) Standardfehler von Koeffizienten und anderen Informationen an.

Wenn jemand da draußen gut darin ist, Formeln zu interpretieren, können Sie mir helfen, indem Sie diesen Code in eine lesbare Formel schreiben?

  -Funktion (x, a, b, clx) {
  ifelse (x < clx, a + b * x + (-0,5 * b / clx) * x * x,
         a + b * clx + (-0,5 * b / clx) * clx * clx)}
 

$$ f (x, a, b, x_ {cl}) = \ begin {Fälle} a + bx + (\ frac {-0,5b} {x_ {cl}}) \ times x ^ 2 , & \ text {if} \ x < x_ {cl} \\ a + bx_ {cl} + (\ frac {-0,5b} {x_ {cl}}) \ times {x_ {cl}} ^ 2 , & \ text {sonst} \ end {Fälle} $$ span>

was vereinfacht zu:

$$ f (x, a, b, x_ {cl}) = \ begin {Fälle} a + bx \ left (1 - \ frac {x} {2x_ {cl}} \ right) , & \ text {if} \ x < x_ {cl} \\ a + \ frac {bx_ {cl}} {2} , & \ text {sonst} \ end {Fälle} $$ span>

wobei ich $ x_ {cl} $ span> durch clx ersetzt habe, um die Lesbarkeit zu verbessern.