Initialwert bestimmen für Weibullkurve mit nls2

Modelle zur Korrelations- und Regressionsanalyse

Moderator: EDi

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

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag 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 
> 
eimichae
Beiträge: 36
Registriert: Mo Okt 10, 2016 8:48 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag 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
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...

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 
schubbiaschwilli
Beiträge: 253
Registriert: Di Jun 27, 2017 12:09 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag 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")
eimichae
Beiträge: 36
Registriert: Mo Okt 10, 2016 8:48 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

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








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

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag 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
eimichae
Beiträge: 36
Registriert: Mo Okt 10, 2016 8:48 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag 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!
schubbiaschwilli
Beiträge: 253
Registriert: Di Jun 27, 2017 12:09 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag 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."?
eimichae
Beiträge: 36
Registriert: Mo Okt 10, 2016 8:48 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag 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
schubbiaschwilli
Beiträge: 253
Registriert: Di Jun 27, 2017 12:09 pm

Re: Initialwert bestimmen für Weibullkurve mit nls2

Beitrag 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")
Antworten