Thank you, I'll try it! On Thu, Mar 18, 2021 at 9:46 PM Rui Barradas <ruipbarra...@sapo.pt> wrote: > > Hello, > > Maybe a bit late but there is a contributed package [1] for quantitative > PCR fitting non-linear models with the Levenberg-Marquardt algorithm. > > estim and vector R below are your model and your fitted values vector. > The RMSE of this fit is smaller than your model's. > > > Isn't this simpler? > > > library(qpcR) > > df1 <- data.frame(Cycles = seq_along(high), high) > > fit <- pcrfit( > data = df1, > cyc = 1, > fluo = 2 > ) > summary(fit) > > coef(estim) > coef(fit) > > > sqrt(sum(resid(estim)^2)) > #[1] 1724.768 > sqrt(sum(resid(fit)^2)) > #[1] 1178.318 > > > highpred <- predict(fit, newdata = df1) > > plot(1:45, high, type = "l", col = "red") > points(1:45, R, col = "blue") > points(1:45, highpred$Prediction, col = "cyan", pch = 3) > > > [1] https://CRAN.R-project.org/package=qpcR > > Hope this helps, > > Rui Barradas > > Às 06:51 de 18/03/21, Luigi Marongiu escreveu: > > It worked. I re-written the equation as: > > ``` > > rutledge_param <- function(p, x, y) ( (p$M / ( 1 + exp(-(x-p$m)/p$s)) > > ) + p$B ) - y > > ``` > > and used Desmos to estimate the slope, so: > > ``` > > estim <- nls.lm(par = list(m = halfCycle, s = 2.77, M = MaxFluo, B = > > high[1]), > > fn = rutledge_param, x = 1:45, y = high) > > summary(estim) > > R <- rutledge(list(half_fluorescence = 27.1102, slope = 2.7680, > > max_fluorescence = 11839.7745, back_fluorescence = > > -138.8615) , 1:45) > > points(1:45, R, type="l", col="red") > > ``` > > > > Thanks > > > > On Tue, Mar 16, 2021 at 8:29 AM Luigi Marongiu <marongiu.lu...@gmail.com> > > wrote: > >> > >> Just an update: > >> I tried with desmos and the fitting looks good. Desmos calculated the > >> parameters as: > >> Fmax = 11839.8 > >> Chalf = 27.1102 (with matches with my estimate of 27 cycles) > >> k = 2.76798 > >> Fb = -138.864 > >> I forced R to accept the right parameters using a single named list > >> and re-written the formula (it was a bit unclear in the paper): > >> ``` > >> rutledge <- function(p, x) { > >> m = p$half_fluorescence > >> s = p$slope > >> M = p$max_fluorescence > >> B = p$back_fluorescence > >> y = (M / (1+exp( -((x-m)/s) )) ) + B > >> return(y) > >> } > >> ``` > >> but when I apply it I get a funny graph: > >> ``` > >> desmos <- rutledge(list(half_fluorescence = 27.1102, slope = 2.76798, > >> max_fluorescence = 11839.8, back_fluorescence > >> = -138.864) , high) > >> ``` > >> > >> On Mon, Mar 15, 2021 at 7:39 AM Luigi Marongiu <marongiu.lu...@gmail.com> > >> wrote: > >>> > >>> Hello, > >>> the negative data comes from the machine. Probably I should use raw > >>> data directly, although in the paper this requirement is not reported. > >>> The p$x was a typo. Now I corrected it and I got this error: > >>> ``` > >>> > >>>> rutledge_param <- function(p, x, y) ((p$M / (1 + exp(-1*(x-p$m)/p$s))) + > >>>> p$B) - y > >>>> estim <- nls.lm(par = list(m = halfFluo, s = slopes, M = MaxFluo, B = > >>>> high[1]), > >>> + fn = rutledge_param, x = 1:45, y = high) > >>> Error in dimnames(x) <- dn : > >>> length of 'dimnames' [2] not equal to array extent > >>> ``` > >>> Probably because 'slopes' is a vector instead of a scalar. Since the > >>> slope is changing, I don't think is right to use a scalar, but I tried > >>> and I got: > >>> ``` > >>>> estim <- nls.lm(par = list(m = halfFluo, s = 1, M = MaxFluo, B = > >>>> high[1]), > >>> + fn = rutledge_param, x = 1:45, y = high) > >>>> estim > >>> Nonlinear regression via the Levenberg-Marquardt algorithm > >>> parameter estimates: 6010.94, 1, 12021.88, 4700.49288888889 > >>> residual sum-of-squares: 1.14e+09 > >>> reason terminated: Relative error in the sum of squares is at most `ftol'. > >>> ``` > >>> The values reported are the same I used at the beginning apart from > >>> the last (the background parameter) which is 4700 instead of zero. If > >>> I plug it, I get an L shaped plot that is worse than that at the > >>> beginning: > >>> ``` > >>> after = init = rutledge(halfFluo, 1, MaxFluo, 4700.49288888889, high) > >>> points(1:45, after, type="l", col="blue") > >>> ``` > >>> What did I get wrong here? > >>> Thanks > >>> > >>> On Sun, Mar 14, 2021 at 8:05 PM Bill Dunlap <williamwdun...@gmail.com> > >>> wrote: > >>>> > >>>>> rutledge_param <- function(p, x, y) ((p$M / (1 + > >>>>> exp(-1*(p$x-p$m)/p$s))) + p$B) - y > >>>> > >>>> Did you mean that p$x to be just x? As is, this returns numeric(0) > >>>> for the p that nls.lm gives it because p$x is NULL and NULL-aNumber is > >>>> numeric(). > >>>> > >>>> -Bill > >>>> > >>>> On Sun, Mar 14, 2021 at 9:46 AM Luigi Marongiu > >>>> <marongiu.lu...@gmail.com> wrote: > >>>>> > >>>>> Hello, > >>>>> I would like to use the Rutledge equation > >>>>> (https://pubmed.ncbi.nlm.nih.gov/15601990/) to model PCR data. The > >>>>> equation is: > >>>>> Fc = Fmax / (1+exp(-(C-Chalf)/k)) + Fb > >>>>> I defined the equation and another that subtracts the values from the > >>>>> expectations. I used minpack.lm to get the parameters, but I got an > >>>>> error: > >>>>> ``` > >>>>> > >>>>>> library("minpack.lm") > >>>>>> h <- c(120.64, 66.14, 34.87, 27.11, 8.87, -5.8, 4.52, -7.16, -17.39, > >>>>> + -14.29, -20.26, -14.99, -21.05, -20.64, -8.03, -21.56, -1.28, > >>>>> 15.01, > >>>>> + 75.26, 191.76, 455.09, 985.96, 1825.59, 2908.08, 3993.18, > >>>>> 5059.94, > >>>>> + 6071.93, 6986.32, 7796.01, 8502.25, 9111.46, 9638.01, 10077.19, > >>>>> + 10452.02, 10751.81, 11017.49, 11240.37, 11427.47, 11570.07, > >>>>> 11684.96, > >>>>> + 11781.77, 11863.35, 11927.44, 11980.81, 12021.88, 12058.35, > >>>>> 12100.63, > >>>>> + 12133.57, 12148.89, 12137.09) > >>>>>> high <- h[1:45] > >>>>>> MaxFluo <- max(high) > >>>>>> halfFluo <- MaxFluo/2 > >>>>>> halfCycle = 27 > >>>>>> find_slope <- function(X, Y) { > >>>>> + Slope <- c(0) > >>>>> + for (i in 2:length(X)) { > >>>>> + delta_x <- X[i] - X[i-1] > >>>>> + delta_y <- Y[i] - Y[i-1] > >>>>> + Slope[i] <- delta_y/delta_x > >>>>> + } > >>>>> + return(Slope) > >>>>> + } > >>>>>> slopes <- find_slope(1:45, high) > >>>>>> > >>>>>> rutledge <- function(m, s, M, B, x) { > >>>>> + divisor = 1 + exp(-1* ((x-m)/s) ) > >>>>> + y = (M/divisor) + B > >>>>> + return(y) > >>>>> + } > >>>>>> rutledge_param <- function(p, x, y) ((p$M / (1 + > >>>>>> exp(-1*(p$x-p$m)/p$s))) + p$B) - y > >>>>>> > >>>>>> > >>>>>> init = rutledge(halfFluo, slopes, MaxFluo, 0, high) > >>>>>> points(1:45, init, type="l", col="red") > >>>>>> estim <- nls.lm(par = list(m = halfFluo, s = slopes, M = MaxFluo, B = > >>>>>> high[1]), > >>>>> + fn = rutledge_param, x = 1:45, y = high) > >>>>> Error in nls.lm(par = list(m = halfFluo, s = slopes, M = MaxFluo, B = > >>>>> high[1]), : > >>>>> evaluation of fn function returns non-sensible value! > >>>>> ``` > >>>>> > >>>>> Where could the error be? > >>>>> > >>>>> > >>>>> -- > >>>>> Best regards, > >>>>> Luigi > >>>>> > >>>>> ______________________________________________ > >>>>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > >>>>> https://stat.ethz.ch/mailman/listinfo/r-help > >>>>> PLEASE do read the posting guide > >>>>> http://www.R-project.org/posting-guide.html > >>>>> and provide commented, minimal, self-contained, reproducible code. > >>> > >>> > >>> > >>> -- > >>> Best regards, > >>> Luigi > >> > >> > >> > >> -- > >> Best regards, > >> Luigi > > > > > >
-- Best regards, Luigi ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
