 # Finding initial values for Weibull curve using nls2

Dear all,
I would like to plot several 3 parameter Weibull curves with the function` y ~ m*exp(-1*(x/b)^c)`.
In some few cases when fitting the final curve using my "brute-forced" initial values I run into an error ("missing or infinitive value generated by model"). I usually could solve the problem by feeding my nls2 model with sligthly different initial values. But for a few models this is not working.

A friend of mine used the "solver" function in excel and found without much problem initial values, so I do know that the data can be described with a 3 parameter Weibull curve. I also used the initial values calculated by excel (m=17.84, b=4.87, c=0.53) in my model, but without much success.
Could anyone help me? do I have to change something in my nls2 settings or my approach to generate the initial model values?
below I posted the code of one example model that caused some trouble.
Help would be very much appreciated.
In my example y = PDLWP and x=Photo

``````#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 = "nls", formula = y ~ m*exp(-1*(x/b)^c), size = 0.9, se = FALSE, colour = "black")

exampledip3
``````

I think that param b going negative is a problem constraint, I tried using lower param on it.

``````#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),
cc=c(0,0.6))

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

# 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)))+
geom_smooth(aes(y=predict(fitdip3,newdata=dip3dta)),color="blue")+
geom_smooth(aes(y=predict(finalfitdip3,newdata=dip3dta)),color="green")

exampledip3
`````` Dear nirgrahamuk,
Thank you for your reply. This solution works great for most of my problematic cases. Thanks a lot!
Allow me to ask a follwoup question.

I already plotted multiple of these weibull curves and combined them 1 plot. For plotting the curves I used the"stat_smooth" command with method="nls" (not loess). Using this approach I was able to combine the curves in a final plot using the approach as shown in the example code below.

``````plotall <- ggplot(dta, aes(x = PDLWP, y = Photo, color = geno)) +
theme_classic() +
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(data = dip1dta, mapping = aes(x = PDLWP, y = Photo), method = "nls", formula = y ~ m * exp(-1 * (x / b)^c), method.args = list(start = as.list(coef(finalfitdip1))), size = 0.9, se = FALSE, colour = "#d6c46750") +
stat_smooth(data = dip2dta, mapping = aes(x = PDLWP, y = Photo), method = "nls", formula = y ~ m * exp(-1 * (x / b)^c), method.args = list(start = as.list(coef(finalfitdip2))), size = 0.9, se = FALSE, colour = "#d6c46750")+
.....etc.
plotall
``````

Hence, I would find it great, if I could use the "stat_smooht" and method="nls" rather than the "geom_smooth" , method="loess" command (like you did ) for plotting the curves. However, when plotting the curves using method="nls" ggplot is complaining saying "computation failed". Do I have to "lower the parameter boundaries" there as well?

Much obliged!

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.