knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.path = "README-" )
Tools for expanding the non-linear regression method nls and nlsList from nlme.
Please report any issues/suggestions for improvement in the issues link for the repository. Or please email d.padfield@exeter.ac.uk.
This package is licensed under GPL-3.
The mainstay of this package is nlsLoop::nlsLoop()
. If you have a dataset where you want to fit the same model over many levels of a factor using non-linear regression, you may use nlme::nlsList()
. However, nlsList()
only allows for one set of starting values so it is possible that not all of the models will converge if there is large variance in the data or the shape of the relationships are different, and thus the parameter values, are likely to be very different.
nlsLoop()
allows for a range of starting values, looping through and attempting a fit for thousands of unique combination of starting values at each level of the factor, picking the best fit for each model using AIC scores.
It also returns predictions for each level of the factor in tidy format, making it simple to plot and evaluate each fit using ggplot2
.
A more in-depth tutorial and explanation of how to use nlsLoop()
can be found by calling vignette('nlsLoop')
from within R.
# install package devtools::install_github("padpadpadpad/nlsLoop", build_vignettes = TRUE)
# load in nlsLoop and other packages library(nlsLoop) library(ggplot2) # load nlsLoop vignette vignette('nlsLoop') # load in example data set data("Chlorella_TRC") # define the Sharpe-Schoolfield equation schoolfield_high <- function(lnc, E, Eh, Th, temp, Tc) { Tc <- 273.15 + Tc k <- 8.62e-5 boltzmann.term <- lnc + log(exp(E/k*(1/Tc - 1/temp))) inactivation.term <- log(1/(1 + exp(Eh/k*(1/Th - 1/temp)))) return(boltzmann.term + inactivation.term) }
# run nlsLoop fits <- nlsLoop(ln.rate ~ schoolfield_high(lnc, E, Eh, Th, temp = K, Tc = 20), data = Chlorella_TRC, tries = 500, id_col = 'curve_id', param_bds = c(-10, 10, 0.1, 2, 0.5, 5, 285, 330), r2 = 'Y', supp_errors = 'Y', AICc = 'Y', na.action = na.omit, lower = c(lnc = -10, E = 0, Eh = 0, Th = 0))
head(fits$params) head(fits$predictions)
# plot a single curve plot_id_nlsLoop(data = Chlorella_TRC, param_data = fits, id = '1')
# create pdf of each curve plot_all_nlsLoop('path/of/where/you/want/to/save/me.pdf', data = Chlorella_TRC, param_data = fits)
# get distinct values of process, flux and growth.temp for each value of curve_id d_treatment <- Chlorella_TRC[,c('curve_id','process', 'growth.temp', 'flux')] d_treatment <- d_treatment[!duplicated(d_treatment),] # merge with predictions by curve_id fits$predictions <- merge(fits$predictions, d_treatment, by = 'curve_id') # plot every curve ggplot() + geom_point(aes(K - 273.15, ln.rate, col = flux), size = 2, Chlorella_TRC) + geom_line(aes(K - 273.15, ln.rate, col = flux, group = curve_id), alpha = 0.5, fits$predictions) + facet_wrap(~ growth.temp + process, labeller = labeller(.multi_line = F)) + scale_colour_manual(values = c('green4', 'black')) + theme_bw(base_size = 12, base_family = 'Helvetica') + ylab('log Metabolic rate') + xlab('Assay temperature (ÂșC)') + theme(legend.position = c(0.9, 0.15))
# calculate confidence intervals for each fit CIs <- confint_nlsLoop(Chlorella_TRC, fits) # bind with factors dataframe CIs <- merge(CIs, d_treatment, by = 'curve_id') # plot ggplot(CIs, aes(col = flux)) + geom_point(aes(curve_id, mean)) + facet_wrap(~ param, scale = 'free_x', ncol = 4) + geom_linerange(aes(curve_id, ymin = CI_lwr, ymax = CI_upr)) + coord_flip() + scale_color_manual(values = c('green4', 'black')) + theme_bw(base_size = 12, base_family = 'Helvetica') + theme(legend.position = 'top') + xlab('curve') + ylab('parameter estimate')
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.