Deutsches R-Forum

Deutschsprachiges Forum zur Statistikumgebung R
http://forum.r-statistik.de/

Initialwert bestimmen für Weibullkurve mit nls2

http://forum.r-statistik.de/viewtopic.php?t=3058

Seite 3 von 4

Re: Initialwert bestimmen für Weibullkurve mit nls2

Verfasst: Di Mär 16, 2021 2:56 pm
von schubbiaschwilli
Gude!

Bei mir nicht:

Code: Alles auswählen

> fitdip3
Nonlinear regression model
  model: Photo ~ m * exp(-1 * (PDLWP/b)^c)
   data: dip3dta
      m       b       c 
42.8571  1.9048  0.4857 
 residual sum-of-squares: 58.18

Number of iterations to convergence: 10648 
Achieved convergence tolerance: NA
> finalfitdip3 <- nls(Photo ~ m*exp(-1*(PDLWP/b)^c), data=dip3dta, start=as.list(coef(fitdip3)))
Fehler in numericDeriv(form[[3L]], names(ind), env) : 
  Fehlender Wert oder etwas Unendliches durch das Modell erzeugt
> 
> # 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
Warnmeldung:
Computation failed in `stat_smooth()`:
Objekt 'm' nicht gefunden 
>

Re: Initialwert bestimmen für Weibullkurve mit nls2

Verfasst: Di Mär 16, 2021 3:31 pm
von eimichae
Habs soeben ausprobiert. Komisch. Es sollte alles da sein. Nur copy past in die R-console nötig.
Habe es unten nochmals reinkopiert

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

Re: Initialwert bestimmen für Weibullkurve mit nls2

Verfasst: Di Mär 16, 2021 3:41 pm
von schubbiaschwilli
Gude!

Hm...

Code: Alles auswählen

> finalfitdip3
Nonlinear regression model
  model: Photo ~ m * exp(-1 * (PDLWP/b)^c)
   data: dip3dta
        m         b         c 
3.697e+02 8.163e-03 2.069e-01 
 residual sum-of-squares: 56.11

Algorithm "port", convergence message: relative convergence (4)
> # 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
Warnmeldung:
Computation failed in `stat_smooth()`:
Fehlender Wert oder etwas Unendliches durch das Modell erzeugt

Re: Initialwert bestimmen für Weibullkurve mit nls2

Verfasst: Di Mär 16, 2021 3:44 pm
von schubbiaschwilli
Gude!

Du solltest mal

Code: Alles auswählen

m*exp(-1*(x/b)^c)
definieren...
Nachtrag: Hast du deine Umgebung/Variablen gesäubert?

Code: Alles auswählen

rm(list = ls())

Code: Alles auswählen

# 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")

Re: Initialwert bestimmen für Weibullkurve mit nls2

Verfasst: Di Mär 16, 2021 4:14 pm
von eimichae
das ist ja sehr komisch. Ja habe alles gesäubert.
Dann mach ichs Schritt für Schritt.
Input data
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")








Re: Initialwert bestimmen für Weibullkurve mit nls2

Verfasst: Di Mär 16, 2021 4:16 pm
von eimichae
Code:
#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

Re: Initialwert bestimmen für Weibullkurve mit nls2

Verfasst: Di Mär 16, 2021 4:19 pm
von eimichae
Klappt es jetzt?

Ich hätte 2 Fragen:

1:
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?

2: 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. Warum?


Vielen Dank!

Re: Initialwert bestimmen für Weibullkurve mit nls2

Verfasst: Di Mär 16, 2021 4:26 pm
von schubbiaschwilli
Gude!
Warum?
Naja, wie ich bereits schrieb: Du solltest mal

Code: Alles auswählen

m*exp(-1*(x/b)^c)
definieren...

Nachtrag:

Code: Alles auswählen

Warnmeldung:
Computation failed in `stat_smooth()`:
Fehlender Wert oder etwas Unendliches durch das Modell erzeugt
Das verstehst du unter "Mein Beispiel sollten doch reproduzierbar sein."?

Re: Initialwert bestimmen für Weibullkurve mit nls2

Verfasst: Di Mär 16, 2021 4:52 pm
von eimichae
Hallo schubbiaschwilli,
Bitte entschuldige. Ich glaube wir reden/schreiben aneinander vorbei.
Das verstehst du unter "Mein Beispiel sollten doch reproduzierbar sein."?
Ich habe das so verstanden, dass ich ein Beispilecode bereitstellen soll, mit dem jeder nachvollziehen kann wo mein Problem liegt. Das habe ich (meines Wissens nach) gemacht.

Bezüglich deiner Aussage ich solle
y=m*exp(-1*(x/b)^c)
definieren.
y= "Photo" in meine Beispieldatensatz
x= PDLWP in meine Beispieldatensatz
m, b, c = Kurvenparameter.

Weiss nicht ob du das meinst mit "definieren"
Gruss und Dank

Re: Initialwert bestimmen für Weibullkurve mit nls2

Verfasst: Di Mär 16, 2021 5:02 pm
von schubbiaschwilli
Gude!

OK, was hast du so für Werte?

Code: Alles auswählen

> m
Fehler: Objekt 'm' nicht gefunden
> x
Fehler: Objekt 'x' nicht gefunden
> b
Fehler: Objekt 'b' nicht gefunden
> c
function (...)  .Primitive("c")
Alle Zeiten sind UTC+02:00
Seite 3 von 4

Powered by phpBB® Forum Software © phpBB Limited

Deutsche Übersetzung durch phpBB.de