Initialwert bestimmen für Weibullkurve mit nls2

Modelle zur Korrelations- und Regressionsanalyse

Moderator: EDi

eimichae
Beiträge: 36
Registriert: Mo Okt 10, 2016 8:48 pm

Initialwert bestimmen für Weibullkurve mit nls2

Beitrag von eimichae »

Hallo zuammen,
Ich habe mehrere Datensätze, die den Zusammenhang zwischen Trockenstress ( X=PDLWP) und Photosyntheseleistung (Y=Photo) einzelner Pflanzen beschreiben. Aus experimentellen Gründen möchte ich den Zusammenhang zwischen den X und Y Variblen mit 3-Parameter Weibullkurven beschreiben ( y ~ m*exp(-1*(x/b)^c) ).
Ich benutze das R-package nls2 um die Initialwerte meiner Kurven mit dem Algorythmus "brute-force" zu bestimmen. Die Kurven werden dann mit dem Package nls gefittet und schlussendlich in ggplot2 geplottet.
In einigen Fällen jedoch, kann ich keine geeigneten Initialwerte bestimmen. Wenn ich die Kurven in nls fitten will erhalte ich die Fehlermeldung: "Fehlender Wert oder Unendlichkeit durch Modell erzeugt". In den meisten Fällen lässt sich das Problem lösen indem ich den Initialwertbereich in nls2 leicht anpasse. Für einige wenige Kurven, allerdings klappt das nicht.

Ein Kollege benutzte für die Initalwertbestimmung und Kurvenplotten die "solver" Funktion in Excel. Er hatte damit keine Problemen. Ich weiss daher, dass auch meine problematischen Datnesätze mit der 3-Paramerterweibullkurve beschrieben werden können. Aber auch wenn ich die Initialwerte von Excel benutze (für das untenstehende Beispiel m=17.84, b=4.87, c=0.53) klappt es nicht.
Untenstehend habe ich einen Beispielcode für einen problematischen Datensatz gepostet.
Könnte mir jemand einen Tipp geben was ich ändern muss oder wie ich besser Initialwerte bestimmen kann? Vielen herzlichen Dank!

Code: Alles auswählen

#Load packages
library(ggplot2)
library(nlme)
library(nls2)
library(proto)

#model structure: 3 parameter Weibull
#y ~ m*exp(-1*(x/b)^c)

dip3dta<-structure(list(ploidy = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                    1L, 1L, 1L, 1L, 1L, 1L), .Label = c("dip", "trp"), class = "factor"), 
               geno = structure(c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                  3L, 3L, 3L), .Label = c("dip1", "dip2", "dip3", "dip4", "dip5", 
                                                          "dip6", "dip7", "dip8", "trp1", "trp2", "trp3", "trp4", "trp5", 
                                                          "trp6", "trp7", "trp8"), class = "factor"), Photo = c(10.03907124, 
                                                                                                                16.04016877, 5.799933798, 6.256058037, 1.34916505, 9.609508391, 
                                                                                                                12.84023945, 8.436093321, 7.732332332, 15.38729611, 2.157795186, 
                                                                                                                5.93553951, 3.37322132), WBPhoto = c(11.77970983, 13.52705488, 
                                                                                                                                                     7.585920181, 6.118582453, 2.570461685, 10.80358492, 9.445462376, 
                                                                                                                                                     5.386306724, 5.840252952, 15.84494637, 3.60398487, 9.32456564, 
                                                                                                                                                     3.437440219), PDLWP = c(3.1, 2.6, 5.8, 7.7, 19, 3.5, 4.25, 
                                                                                                                                                                             9, 8.16, 2.25, 13.92, 4.33, 14.58), Treatment = structure(c(1L, 
                                                                                                                                                                                                                                         1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "DC", class = "factor")), row.names = 27:39, class = "data.frame")

#Define the outer boundaries to search for initial values
grddip3 <- data.frame(m=c(0,45),
                      b=c(0,8),
                      c=c(0,0.6))

#Brute-force initial values
fitdip3 <- nls2(Photo ~ m*exp(-1*(PDLWP/b)^c),
                data=dip3dta,
                start = grddip3,
                algorithm = "brute-force",
                control=list(maxiter=10000))
fitdip3
finalfitdip3 <- nls(Photo ~ m*exp(-1*(PDLWP/b)^c), data=dip3dta, start=as.list(coef(fitdip3)))

# Plot in ggplot
exampledip3<-ggplot(dip3dta, aes(x=PDLWP,y=Photo)) + geom_point(size=2)+
  theme_classic()+
  theme(legend.position="none")+
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=17,),axis.title.y=element_text(margin=margin(0,20,0,0)))+
  stat_smooth(method = "lm", formula = y ~ m*exp(-1*(x/b)^c), size = 0.9, se = FALSE, colour = "black")

exampledip3
schubbiaschwilli
Beiträge: 253
Registriert: Di Jun 27, 2017 12:09 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag von schubbiaschwilli »

Gude!

Hm, hast du mal andere Methoden zum Fitten benutzt?
Levenberg-Marquardt lässt sich recht schnell umsetzen:
https://de.wikipedia.org/wiki/Levenberg ... lgorithmus
https://www.r-bloggers.com/2012/07/a-better-nls/

Auch Nelder-Mead oder BFGS können funktionieren:
https://stat.ethz.ch/R-manual/R-devel/l ... optim.html

Und es ist auch nix ungewöhnliches, wenn man verschiedene Parameterkombinationen für die Startwerte durchprobieren muss - Der Wert RSS ist hier das Maß.

Dank&Gruß
Schubbiaschwilli
eimichae
Beiträge: 36
Registriert: Mo Okt 10, 2016 8:48 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag von eimichae »

Guten Abend,
Erstmal vielen Dank für deine Antwort Schubbiaschwilli.

Ich habe erstmal eine grundsätzliche Frage:
In meinem Beispiel "Brute-force" ich die Initialwerte mit der nls2 Funktion. Dadurch erhalte ich das Objekt(?) "fitdip3"
Anschliessend benutze ich noch die "nls" Funktion um mittels den gefundene Initialwerten (fitdip3) das Objekt "finalfitdip3" zu erhalten. Dieser zweite Schritt ist komplett unnötig oder? Ich könnte gleich meinen Initialwerte die ich über "nls2" erhalten habe für ggplot2 verwenden. Oder würde das finalfitdip3 Objekt mir idealerweise meine Finalen Kurvenparameter (m, b, c) geben?

Ich bin nun auch etwas weiter mit meiner Kurve. Ich habe nun die untere Grenze für den Parameter "b" etwas herabgesetze und den verwende den Algoritmus "port". Nun schaffe ich es allerdings noch nicht die Kurve in GGplot mit der stat_smooth option zu plotten (ich möchte nicht die "loess" Methodik von geom_smooth benutzen). Es kommt die Warnmeldung "Warnmeldung:
Fehlender Wert oder etwas Unendliches durch das Modell erzeugt

Kannst du mir da weiterhelfen? Unten der aktualisierte Code:
Vielen Dank!

Code: Alles auswählen

#Load packages
library(ggplot2)
library(nlme)
library(nls2)
library(proto)

#model structure: 3 parameter Weibull
#y ~ m*exp(-1*(x/b)^c)

dip3dta<-structure(list(ploidy = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                             1L, 1L, 1L, 1L, 1L, 1L), .Label = c("dip", "trp"), class = "factor"), 
                        geno = structure(c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                           3L, 3L, 3L), .Label = c("dip1", "dip2", "dip3", "dip4", "dip5", 
                                                                   "dip6", "dip7", "dip8", "trp1", "trp2", "trp3", "trp4", "trp5", 
                                                                   "trp6", "trp7", "trp8"), class = "factor"), Photo = c(10.03907124, 
                                                                                                                         16.04016877, 5.799933798, 6.256058037, 1.34916505, 9.609508391, 
                                                                                                                         12.84023945, 8.436093321, 7.732332332, 15.38729611, 2.157795186, 
                                                                                                                         5.93553951, 3.37322132), WBPhoto = c(11.77970983, 13.52705488, 
                                                                                                                                                              7.585920181, 6.118582453, 2.570461685, 10.80358492, 9.445462376, 
                                                                                                                                                              5.386306724, 5.840252952, 15.84494637, 3.60398487, 9.32456564, 
                                                                                                                                                              3.437440219), PDLWP = c(3.1, 2.6, 5.8, 7.7, 19, 3.5, 4.25, 
                                                                                                                                                                                      9, 8.16, 2.25, 13.92, 4.33, 14.58), Treatment = structure(c(1L, 
                                                                                                                                                                                                                                                  1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "DC", class = "factor")), row.names = 27:39, class = "data.frame")

#Define the outer boundaries to search for initial values
grddip3 <- data.frame(m=c(0,40),
                      b=c(0,12),
                      c=c(0,8))

#Brute-force initial values
fitdip3 <- nls2(Photo ~ m*exp(-1*(PDLWP/b)^c),
                data=dip3dta,
                start = grddip3,
                algorithm = "brute-force",
                control=list(maxiter=10000))
fitdip3

finalfitdip3 <- nls(Photo ~ m*exp(-1*(PDLWP/b)^c),
                    data=dip3dta,
                    start=as.list(coef(fitdip3)),
                    algorithm = "port",
                    lower = coef(fitdip3)/700,
                    control=list(maxiter=1000))
finalfitdip3
# Plot in ggplot
exampledip3<-ggplot(dip3dta, aes(x=PDLWP,y=Photo)) + geom_point(size=2)+
  theme_classic()+
  theme(legend.position="none")+
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=17,),axis.title.y=element_text(margin=margin(0,20,0,0)))+
  stat_smooth(method = "nls", formula = y ~ m*exp(-1*(x/b)^c), method.args=list(start=as.list(coef(finalfitdip3))), size = 0.9, se = FALSE, colour = "black")

exampledip3
schubbiaschwilli
Beiträge: 253
Registriert: Di Jun 27, 2017 12:09 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag von schubbiaschwilli »

Gude!

Ich habe sowas bisher nur mit dem normalen plot, und nicht mit ggplot2 gebaut - Falls das interessant ist, suche ich das morgen mal.

Dank&Gruß
Schubbiaschwilli
Athomas
Beiträge: 768
Registriert: Mo Feb 26, 2018 8:19 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag von Athomas »

Kannst Du bitte mal erklären, wie Du auf "3-Parameter Weibullkurven", insbesondere ihre Beschreibung durch "y ~ m*exp(-1*(x/b)^c)", kommst?
schubbiaschwilli
Beiträge: 253
Registriert: Di Jun 27, 2017 12:09 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag von schubbiaschwilli »

Gude!

Hm, ich (bzw. Minitab) kenne die (und in einer Vorlesung über Statistische Methoden im Qualitätsmanagement wurde die auch mal behandelt): https://support.minitab.com/de-de/minit ... tribution/
Ob die Formel stimmt, hab' ich nicht geschaut.

Dank&Gruß
Schubbiaschwilli
Athomas
Beiträge: 768
Registriert: Mo Feb 26, 2018 8:19 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag von Athomas »

schubbiaschwilli hat geschrieben: Mo Mär 15, 2021 11:30 am Gude!

Hm, ich (bzw. Minitab) kenne die (und in einer Vorlesung über Statistische Methoden im Qualitätsmanagement wurde die auch mal behandelt): https://support.minitab.com/de-de/minit ... tribution/
Ob die Formel stimmt, hab' ich nicht geschaut.

Dank&Gruß
Schubbiaschwilli
Die "Weibull-Verteilung" ist mir selbstredend auch ein Begriff :D , es geht mir hier eher um die konkrete Verwendung.

Was ich von der Verwendung von Verteilungs- (oder Dichte-) Funktionen zum Fitten von "nicht-Verteilungs-Daten" halte, hatte ich Dir - glaube ich - schon an anderem Orte gesagt...
schubbiaschwilli
Beiträge: 253
Registriert: Di Jun 27, 2017 12:09 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag von schubbiaschwilli »

Gude!

Code: Alles auswählen

Die "Weibull-Verteilung" ist mir selbstredend auch ein Begriff (...)
Ich dachte, du meinst die 3-Parameter - Ich weiß noch, dass ich damals (in der guten alten Zeit in der Vorlesung - hach) in der Literatur auch erst mal nur 1- und 2-Parameter fand - Ja, die Begriffsverwirrungen in der Statistik sind schon legendär.

Code: Alles auswählen

Was ich von der Verwendung von Verteilungs- (oder Dichte-) Funktionen zum Fitten von "nicht-Verteilungs-Daten" halte, hatte ich Dir - glaube ich - schon an anderem Orte gesagt...
Tue dir an dieser Stelle mal Modelle aus der Wirtschafts- und Finanzmathematik an, und wie da gebogen wird, bis es halbwegs passt.

Dank&Gruß
Schubbiaschwilli
Athomas
Beiträge: 768
Registriert: Mo Feb 26, 2018 8:19 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag von Athomas »

Tue dir an dieser Stelle mal Modelle aus der Wirtschafts- und Finanzmathematik an, und wie da gebogen wird, bis es halbwegs passt.
Kannst Du dafür ein Beispiel geben?
schubbiaschwilli
Beiträge: 253
Registriert: Di Jun 27, 2017 12:09 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag von schubbiaschwilli »

Gude!
Kannst Du dafür ein Beispiel geben?
Schau dir mal an, wie viele verschiedene Modelle es zur Bewertung von Derivate gibt - Also wie Entwickeln sich Preise, Zinsen, Korrelationen usw.usf. in der Zukunft... In meinem Blog spiel ich ja verschiedene Modelle zur Bewertung von Optionen auf Aktien/Aktienindizes durch, und bin nicht viel weiter als - Ich steh' am Anfang - Und das bei Aktien; Zinsen sind noch komplexer; ich verweise mal auf Brigo & Mercurio "Interest Rate Models - Theory and Practice" (https://www.springer.com/de/book/9783540221494) - Das sind knapp 1000 Seiten Zinsmodelle (und ein wenig Kredite) - Und keines davon ist 'richtig' - Und die Diskussion ist ja nicht akademisch, wenn's da falsch läuft, gibt's 'ne Weltwirtschaftskrise - So wie hier: https://en.wikipedia.org/wiki/Long-Term ... Management

Ich nehme an, du verstehst jetzt meine gelegentlich kritischen Anmerkungen, wenn der nächste BWLer hier aufschlägt, und irgendwas mit Aktien analysieren will...

Dank&Gruß
Schubbiaschwilli
Antworten