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