Finding initial values for Weibull curve using nls2

General
1

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
2

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

image

3

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!

4

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.