| Title: | Nonlinear Regression for Agricultural Applications |
|---|---|
| Description: | Additional nonlinear regression functions using self-start (SS) algorithms. One of the functions is the Beta growth function proposed by Yin et al. (2003) <doi:10.1093/aob/mcg029>. There are several other functions with breakpoints (e.g. linear-plateau, plateau-linear, exponential-plateau, plateau-exponential, quadratic-plateau, plateau-quadratic and bilinear), a non-rectangular hyperbola and a bell-shaped curve. Twenty eight (28) new self-start (SS) functions in total. This package also supports the publication 'Nonlinear regression Models and applications in agricultural research' by Archontoulis and Miguez (2015) <doi:10.2134/agronj2012.0506>, a book chapter with similar material <doi:10.2134/appliedstatistics.2016.0003.c15> and a publication by Oddi et. al. (2019) in Ecology and Evolution <doi:10.1002/ece3.5543>. The function 'nlsLMList' uses 'nlsLM' for fitting, but it is otherwise almost identical to 'nlme::nlsList'.In addition, this release of the package provides functions for conducting simulations for 'nlme' and 'gnls' objects as well as bootstrapping. These functions are intended to work with the modeling framework of the 'nlme' package. It also provides four vignettes with extended examples. |
| Authors: | Fernando Miguez [aut, cre]  ,
Caio dos Santos [ctb] (author of SSscard),
José Pinheiro [ctb, cph] (author of nlme::nlsList, nlme::predict.gnls,
nlme::predict.nlme),
Douglas Bates [ctb, cph] (author of nlme::nlsList, nlme::predict.gnls,
nlme::predict.nlme),
R-core [ctb, cph] |
| Maintainer: | Fernando Miguez <[email protected]> |
| License: | GPL-3 |
| Version: | 1.9.9 |
| Built: | 2024-11-13 05:16:48 UTC |
| Source: | https://github.com/femiguez/nlraa |
Data from a paper by Arild Vold on response of barley to nitrogen fertilizer
barleybarley
A data frame with 76 rows and 3 columns
Year when the trial was conducted (1970-1988).
Nitrogen fertilizer (g/m^2).
Grain yield of barley (g/m^2).
Aril Vold (1998). A generalization of ordinary yield response functions. Ecological Applications. 108:227-236.
Bootstraping for linear models
boot_lm( object, f = NULL, R = 999, psim = 2, resid.type = c("resample", "normal", "wild"), data = NULL, verbose = TRUE, ... )boot_lm( object, f = NULL, R = 999, psim = 2, resid.type = c("resample", "normal", "wild"), data = NULL, verbose = TRUE, ... )
object |
object of class |
f |
function to be applied (and bootstrapped), default coef |
R |
number of bootstrap samples, default 999 |
psim |
simulation level for |
resid.type |
either “resample”, “normal” or “wild”. |
data |
optional data argument (useful/needed when data are not in an available environment). |
verbose |
logical (default TRUE) whether to print message if model does not converge. (rare for linear models). |
... |
additional arguments to be passed to function |
The residuals can either be generated by resampling with replacement (default), from a normal distribution (parameteric) or by changing their signs (wild). This last one is called “wild bootstrap”.
at the moment, when the argument data is used, it is not possible to check that it matches the original data used to fit the model. It will also override the fetching of data.
require(car) data(barley, package = "nlraa") ## Fit a linear model (quadratic) fit.lm <- lm(yield ~ NF + I(NF^2), data = barley) ## Bootstrap coefficients by default fit.lm.bt <- boot_lm(fit.lm) ## Compute confidence intervals confint(fit.lm.bt, type = "perc") ## Visualize hist(fit.lm.bt, 1, ci = "perc", main = "Intercept") hist(fit.lm.bt, 2, ci = "perc", main = "NF term") hist(fit.lm.bt, 3, ci = "perc", main = "I(NF^2) term")require(car) data(barley, package = "nlraa") ## Fit a linear model (quadratic) fit.lm <- lm(yield ~ NF + I(NF^2), data = barley) ## Bootstrap coefficients by default fit.lm.bt <- boot_lm(fit.lm) ## Compute confidence intervals confint(fit.lm.bt, type = "perc") ## Visualize hist(fit.lm.bt, 1, ci = "perc", main = "Intercept") hist(fit.lm.bt, 2, ci = "perc", main = "NF term") hist(fit.lm.bt, 3, ci = "perc", main = "I(NF^2) term")
Bootstraping tools for linear mixed-models using a consistent interface
bootstrap function for objects of class gls
boot_lme( object, f = NULL, R = 999, psim = 1, cores = 1L, data = NULL, verbose = TRUE, ... ) boot_gls( object, f = NULL, R = 999, psim = 1, cores = 1L, data = NULL, verbose = TRUE, ... )boot_lme( object, f = NULL, R = 999, psim = 1, cores = 1L, data = NULL, verbose = TRUE, ... ) boot_gls( object, f = NULL, R = 999, psim = 1, cores = 1L, data = NULL, verbose = TRUE, ... )
object |
|
f |
function to be applied (and bootstrapped), default coef (gls) or fixef (lme) |
R |
number of bootstrap samples, default 999 |
psim |
simulation level for vector of fixed parameters either for |
cores |
number of cores to use for parallel computation |
data |
optional data argument (useful/needed when data are not in an available environment). |
verbose |
logical (default TRUE) whether to print a message if model does not converge. |
... |
additional arguments to be passed to function |
This function is inspired by Boot, which does not
seem to work with ‘gls’ or ‘lme’ objects. This function makes multiple copies
of the original data, so it can be very hungry in terms of memory use, but
I do not believe this to be a big problem given the models we typically fit.
require(nlme) require(car) data(Orange) fm1 <- lme(circumference ~ age, random = ~ 1 | Tree, data = Orange) fm1.bt <- boot_lme(fm1, R = 50) hist(fm1.bt)require(nlme) require(car) data(Orange) fm1 <- lme(circumference ~ age, random = ~ 1 | Tree, data = Orange) fm1.bt <- boot_lme(fm1, R = 50) hist(fm1.bt)
Bootstraping tools for nonlinear models using a consistent interface
bootstrap function for objects of class gnls
boot_nlme( object, f = NULL, R = 999, psim = 1, cores = 1L, data = NULL, verbose = TRUE, ... ) boot_gnls( object, f = NULL, R = 999, psim = 1, cores = 1L, data = NULL, verbose = TRUE, ... )boot_nlme( object, f = NULL, R = 999, psim = 1, cores = 1L, data = NULL, verbose = TRUE, ... ) boot_gnls( object, f = NULL, R = 999, psim = 1, cores = 1L, data = NULL, verbose = TRUE, ... )
object |
|
f |
function to be applied (and bootstrapped), default coef (gnls) or fixef (nlme) |
R |
number of bootstrap samples, default 999 |
psim |
simulation level for vector of fixed parameters either for |
cores |
number of cores to use for parallel computation |
data |
optional data argument (useful/needed when data are not in an available environment). |
verbose |
logical (default TRUE) whether to print a message if model does not converge. |
... |
additional arguments to be passed to function |
This function is inspired by Boot, which does not
seem to work with 'gnls' or 'nlme' objects. This function makes multiple copies
of the original data, so it can be very hungry in terms of memory use, but
I do not believe this to be a big problem given the models we typically fit.
require(car) require(nlme) data(barley, package = "nlraa") barley2 <- subset(barley, year < 1974) fit.lp.gnls2 <- gnls(yield ~ SSlinp(NF, a, b, xs), data = barley2) barley2$year.f <- as.factor(barley2$year) cfs <- coef(fit.lp.gnls2) fit.lp.gnls3 <- update(fit.lp.gnls2, params = list(a + b + xs ~ year.f), start = c(cfs[1], 0, 0, 0, cfs[2], 0, 0, 0, cfs[3], 0, 0, 0)) ## This will take a few seconds fit.lp.gnls.Bt3 <- boot_nlme(fit.lp.gnls3, R = 300) confint(fit.lp.gnls.Bt3, type = "perc")require(car) require(nlme) data(barley, package = "nlraa") barley2 <- subset(barley, year < 1974) fit.lp.gnls2 <- gnls(yield ~ SSlinp(NF, a, b, xs), data = barley2) barley2$year.f <- as.factor(barley2$year) cfs <- coef(fit.lp.gnls2) fit.lp.gnls3 <- update(fit.lp.gnls2, params = list(a + b + xs ~ year.f), start = c(cfs[1], 0, 0, 0, cfs[2], 0, 0, 0, cfs[3], 0, 0, 0)) ## This will take a few seconds fit.lp.gnls.Bt3 <- boot_nlme(fit.lp.gnls3, R = 300) confint(fit.lp.gnls.Bt3, type = "perc")
Bootstraping for nonlinear models
boot_nls( object, f = NULL, R = 999, psim = 2, resid.type = c("resample", "normal", "wild"), data = NULL, verbose = TRUE, ... )boot_nls( object, f = NULL, R = 999, psim = 2, resid.type = c("resample", "normal", "wild"), data = NULL, verbose = TRUE, ... )
object |
object of class |
f |
function to be applied (and bootstrapped), default coef |
R |
number of bootstrap samples, default 999 |
psim |
simulation level for |
resid.type |
either “resample”, “normal” or “wild”. |
data |
optional data argument (useful/needed when data are not in an available environment). |
verbose |
logical (default TRUE) whether to print a message if model does not converge. |
... |
additional arguments to be passed to function |
The residuals can either be generated by resampling with replacement
(default or non-parametric), from a normal distribution (parameteric) or by changing
their signs (wild). This last one is called “wild bootstrap”.
There is more information in boot_lm.
at the moment, when the argument data is used, it is not possible to check that it matches the original data used to fit the model. It will also override the fetching of data.
require(car) data(barley, package = "nlraa") ## Fit a linear-plateau fit.nls <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley) ## Bootstrap coefficients by default ## Keeping R small for simplicity, increase R for a more realistic use fit.nls.bt <- boot_nls(fit.nls, R = 1e2) ## Compute confidence intervals confint(fit.nls.bt, type = "perc") ## Visualize hist(fit.nls.bt, 1, ci = "perc", main = "Intercept") hist(fit.nls.bt, 2, ci = "perc", main = "linear term") hist(fit.nls.bt, 3, ci = "perc", main = "xs break-point term")require(car) data(barley, package = "nlraa") ## Fit a linear-plateau fit.nls <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley) ## Bootstrap coefficients by default ## Keeping R small for simplicity, increase R for a more realistic use fit.nls.bt <- boot_nls(fit.nls, R = 1e2) ## Compute confidence intervals confint(fit.nls.bt, type = "perc") ## Visualize hist(fit.nls.bt, 1, ci = "perc", main = "Intercept") hist(fit.nls.bt, 2, ci = "perc", main = "linear term") hist(fit.nls.bt, 3, ci = "perc", main = "xs break-point term")
Confidence interval methods for (non)-linear models.
Method ‘wald’ for objects of class ‘nls’ uses
confint.default
confidence_intervals( x, method = c("wald", "profile", "bootstrap", "simple-bayes", "all"), parm, level = 0.95, verbose = FALSE, ... )confidence_intervals( x, method = c("wald", "profile", "bootstrap", "simple-bayes", "all"), parm, level = 0.95, verbose = FALSE, ... )
x |
|
method |
method or methods to use. Possible options are: ‘wald’, ‘profile’, ‘bootstrap’, ‘all’ |
parm |
optional character string to select the parameter |
level |
probability level |
verbose |
logical (default FALSE) whether to print messages |
... |
additional arguments to be passed to function |
require(car) data(barley, package = "nlraa") ## Fit a linear model (quadratic) fit.lm <- lm(yield ~ NF + I(NF^2), data = barley) cfs.int <- confidence_intervals(fit.lm, method = c("wald", "bootstrap")) fit.nls <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley) cfs.int2 <- confidence_intervals(fit.nls, method = c("wald", "profile", "bootstrap"))require(car) data(barley, package = "nlraa") ## Fit a linear model (quadratic) fit.lm <- lm(yield ~ NF + I(NF^2), data = barley) cfs.int <- confidence_intervals(fit.lm, method = c("wald", "bootstrap")) fit.nls <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley) cfs.int2 <- confidence_intervals(fit.nls, method = c("wald", "profile", "bootstrap"))
object for confidence bands vignette fm1.P.at.x.0.4
fm1.P.at.x.0.4fm1.P.at.x.0.4
An object of class ‘boot’
object created in the vignette in chunk ‘Puromycin-6’
this package vignette
object for confidence bands vignette fm1.P.bt
fm1.P.btfm1.P.bt
An object of class ‘boot’
object created in the vignette in chunk ‘Puromycin-2’
this package vignette
object for confidence bands vignette fm1.P.bt.ft
fm1.P.bt.ftfm1.P.bt.ft
An object of class ‘boot’
object created in the vignette in chunk ‘Puromycin-4’
this package vignette
object for confidence bands vignette fm2.Lob.bt
fm2.Lob.btfm2.Lob.bt
An object of class ‘boot’
object created in the vignette in chunk ‘Loblolly-methods-2’
this package vignette
object for confidence bands vignette fmm1.bt
fmm1.btfmm1.bt
An object of class ‘boot’
object created in the vignette in chunk ‘maizeleafext-2’
this package vignette
Indexes of agreement
printing function for IA_tab
plotting function for a IA_tab, it requires ‘ggplot2’
IA_tab(obs, sim, object, null.object) ## S3 method for class 'IA_tab' print(x, ..., digits = 2) ## S3 method for class 'IA_tab' plot(x, y, ..., type = c("OvsS", "RvsS"))IA_tab(obs, sim, object, null.object) ## S3 method for class 'IA_tab' print(x, ..., digits = 2) ## S3 method for class 'IA_tab' plot(x, y, ..., type = c("OvsS", "RvsS"))
obs |
vector with observed data |
sim |
vector with simulated data (should be the same length as observed) |
object |
alternative to the previous two arguments. An object of class ‘lm’, ‘nls’, ‘lme’ or ‘nlme’ |
null.object |
optional object which represents the ‘null’ model. It is an intercept-only model by default. (Not used at the moment). |
x |
object of class ‘IA_tab’. |
... |
additional plotting arguments (none use at the moment). |
digits |
number of digits for rounding (default is 2) |
y |
not used at the moment |
type |
either “OvsS” (observed vs. simulated) or “RvsS” (residuals vs. simulated). |
This function returns several indexes that might be useful for interpretation
For objects of class ‘lm’ and ‘nls’
bias: mean(obs - sim)
intercept: intercept of the model obs ~ beta_0 + beta_1 * sim + error
slope: slope of the model obs ~ beta_0 + beta_1 * sim + error
RSS (deviance): residual sum of squares of the previous model
MSE (RSS / n): mean squared error; where n is the number of observations
RMSE: squared root of the previous index
R2.1: R-squared extracted from an ‘lm’ object
R2.2: R-squared computed as the correlation between observed and simulated to the power of 2.
ME: model efficiency
NME: Normalized model efficiency
Corr: correlation between observed and simulated
ConCorr: concordance correlation
For objects of class ‘gls’, ‘gnls’, ‘lme’ or ‘nlme’ there are additional metrics such as:
https://en.wikipedia.org/wiki/Coefficient_of_determination
https://en.wikipedia.org/wiki/Nash-Sutcliffe_model_efficiency_coefficient
https://en.wikipedia.org/wiki/Concordance_correlation_coefficient
require(nlme) require(ggplot2) ## Fit a simple model and then compute IAs data(swpg) #' ## Linear model fit0 <- lm(lfgr ~ ftsw + I(ftsw^2), data = swpg) ias0 <- IA_tab(object = fit0) ias0 ## Nonlinear model fit1 <- nls(lfgr ~ SSblin(ftsw, a, b, xs, c), data = swpg) ias1 <- IA_tab(object = fit1) ias1 plot(ias1) ## Linear Mixed Models data(barley, package = "nlraa") fit2 <- lme(yield ~ NF + I(NF^2), random = ~ 1 | year, data = barley) ias2 <- IA_tab(object = fit2) ias2 ## Nonlinear Mixed Model barleyG <- groupedData(yield ~ NF | year, data = barley) fit3L <- nlsLMList(yield ~ SSquadp3(NF, a, b, c), data = barleyG) fit3 <- nlme(fit3L, random = pdDiag(a + b ~ 1)) ias3 <- IA_tab(object = fit3) ias3 plot(ias3) ## Plotting model prds <- predict_nlme(fit3, interval = "conf", plevel = 0) barleyGA <- cbind(barleyG, prds) ggplot(data = barleyGA, aes(x = NF, y = yield)) + geom_point() + geom_line(aes(y = Estimate)) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.2) ## R2M for model 2 R2M(fit2) ## R2M for model 3 R2M(fit3) ## Using IA_tab without a model IA_tab(obs = swpg$lfgr, sim = fitted(fit0))require(nlme) require(ggplot2) ## Fit a simple model and then compute IAs data(swpg) #' ## Linear model fit0 <- lm(lfgr ~ ftsw + I(ftsw^2), data = swpg) ias0 <- IA_tab(object = fit0) ias0 ## Nonlinear model fit1 <- nls(lfgr ~ SSblin(ftsw, a, b, xs, c), data = swpg) ias1 <- IA_tab(object = fit1) ias1 plot(ias1) ## Linear Mixed Models data(barley, package = "nlraa") fit2 <- lme(yield ~ NF + I(NF^2), random = ~ 1 | year, data = barley) ias2 <- IA_tab(object = fit2) ias2 ## Nonlinear Mixed Model barleyG <- groupedData(yield ~ NF | year, data = barley) fit3L <- nlsLMList(yield ~ SSquadp3(NF, a, b, c), data = barleyG) fit3 <- nlme(fit3L, random = pdDiag(a + b ~ 1)) ias3 <- IA_tab(object = fit3) ias3 plot(ias3) ## Plotting model prds <- predict_nlme(fit3, interval = "conf", plevel = 0) barleyGA <- cbind(barleyG, prds) ggplot(data = barleyGA, aes(x = NF, y = yield)) + geom_point() + geom_line(aes(y = Estimate)) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.2) ## R2M for model 2 R2M(fit2) ## R2M for model 3 R2M(fit3) ## Using IA_tab without a model IA_tab(obs = swpg$lfgr, sim = fitted(fit0))
Information criteria table with weights
IC_tab(..., criteria = c("AIC", "AICc", "BIC"), sort = TRUE)IC_tab(..., criteria = c("AIC", "AICc", "BIC"), sort = TRUE)
... |
model fit objects fitted to the same data |
criteria |
either ‘AIC’, ‘AICc’ or ‘BIC’. |
sort |
whether to sort by weights (default to TRUE) |
The delta and weights are calculated based on the ‘criteria’
Live fuel moisture content
lfmclfmc
A data frame with 247 rows and 5 variables:
-factor- Species for which data was recorded ("Grass E", "Grass W", "M. spinosum", "S. bracteolactus")
-integer- time in days 1-80
-factor- plot with levels 1-6 (discrete)
-factor- either P ("East") or SR ("West")
-numeric- Live fuel moisture content (percent)
grouping for regression
A dataset containing the leaf.type, time, plot, site and lfmc (live fuel mass concentration)
object for confidence bands vignette Lob.bt.pe
Lob.bt.peLob.bt.pe
An object of class ‘boot’
object created in the vignette in chunk ‘Loblolly-bootstrap-estimates-1’
this package vignette
Data on leaf extension rate as a response to meristem temperature in maize. The data are re-created liberally from Walls, W.R., 1971. Role of temperature in the regulation of leaf extension in Zea mays. Nature, 229: 46-47. The data points are not the same as in the original paper. Some additional points were inserted to fill in the blanks and allow for reasonable parameter estimates
maizeleafextmaizeleafext
A data frame with 10 rows and 2 columns
Meristem temperature (in Celsius).
Leaf extension rate (relative to 25 degrees).
Walls, W.R., 1971. Role of temperature in the regulation of leaf extension in Zea mays. Nature, 229: 46-47.
Environment which stores indecies and data for bootstraping mostly
nlraa.envnlraa.env
An object of class environment of length 2.
Create an nlraa environment for bootstrapping
This function is a copy of 'nlsList' from the 'nlme' package modified to use the 'nlsLM' function in addition to (optionally) 'nls'. By changing the algorithm argument it is possible to use 'nls' as well
nlsLMList( model, data, start, control, level, subset, na.action = na.fail, algorithm = c("LM", "default", "port", "plinear"), lower = NULL, upper = NULL, pool = TRUE, warn.nls = NA ) ## S3 method for class 'selfStart' nlsLMList( model, data, start, control, level, subset, na.action = na.fail, algorithm = c("LM", "default", "port", "plinear"), lower = NULL, upper = NULL, pool = TRUE, warn.nls = NA )nlsLMList( model, data, start, control, level, subset, na.action = na.fail, algorithm = c("LM", "default", "port", "plinear"), lower = NULL, upper = NULL, pool = TRUE, warn.nls = NA ) ## S3 method for class 'selfStart' nlsLMList( model, data, start, control, level, subset, na.action = na.fail, algorithm = c("LM", "default", "port", "plinear"), lower = NULL, upper = NULL, pool = TRUE, warn.nls = NA )
model |
either a nonlinear model formula, with the response on the left of a ~ operator and an expression involving parameters, covariates, and a grouping factor separated by the | operator on the right, or a selfStart function. |
data |
a data frame |
start |
list with starting values |
control |
control list, see |
level |
an optional integer specifying the level of grouping to be used when multiple nested levels of grouping are present. |
subset |
subset of rows to use |
na.action |
a function that indicates what should happen when the data contain NAs. The default action (na.fail) causes nlsList to print an error message and terminate if there are any incomplete observations. |
algorithm |
choice of algorithm. Default is ‘LM’ which uses ‘nlsLM’ from the minpack.lm package. Other options are: “default”, “port” and “plinear” (nls). |
lower |
vectors of lower and upper bounds, replicated to be as long as start. If unspecified, all parameters are assumed to be unconstrained. Bounds can only be used with the “port” algorithm. They are ignored, with a warning, if given for other algorithms. |
upper |
see ‘lower’ |
pool |
an optional logical value that is preserved as an attribute of the returned value. This will be used as the default for pool in calculations of standard deviations or standard errors for summaries. |
warn.nls |
logical indicating if nls errors (all of which are caught by tryCatch) should be signalled as a “summarizing” warning. |
See function nlsList and nlsLM. This function is a copy of nlsList but with minor changes to use LM instead as the default algorithm. The authors of the original function are Pinheiro and Bates.
an object of class ‘nlsList’
an object of class ‘nlsList’
Jose C. Pinheiro and Douglas M. Bates [email protected] wrote the original nlsList. Fernando E. Miguez made minor changes to use nlsLM in addition to (optionally) nls. R-Core maintains copyright after 2006.
formula method for nlsLMList
## S3 method for class 'formula' nlsLMList( model, data, start = NULL, control, level, subset, na.action = na.fail, algorithm = c("LM", "default", "port", "plinear"), lower = NULL, upper = NULL, pool = TRUE, warn.nls = NA )## S3 method for class 'formula' nlsLMList( model, data, start = NULL, control, level, subset, na.action = na.fail, algorithm = c("LM", "default", "port", "plinear"), lower = NULL, upper = NULL, pool = TRUE, warn.nls = NA )
model |
see |
data |
see |
start |
see |
control |
see |
level |
see |
subset |
see |
na.action |
see |
algorithm |
choice of algorithm default is ‘LM’ which uses ‘nlsLM’ from the minpack.lm package. |
lower |
vectors of lower and upper bounds, replicated to be as long as start. If unspecified, all parameters are assumed to be unconstrained. Bounds can only be used with the “port” algorithm. They are ignored, with a warning, if given for other algorithms. |
upper |
see ‘lower’ |
pool |
see |
warn.nls |
see |
an object of class ‘nlsList’
Largely based on predict.gam, but with some minor modifications to make it compatible
with predict_nls
predict_gam( object, newdata = NULL, type = "link", se.fit = TRUE, terms = NULL, exclude = NULL, block.size = NULL, newdata.guaranteed = FALSE, na.action = na.pass, unconditional = FALSE, iterms.type = NULL, interval = c("none", "confidence", "prediction"), level = 0.95, tvalue = NULL, ... )predict_gam( object, newdata = NULL, type = "link", se.fit = TRUE, terms = NULL, exclude = NULL, block.size = NULL, newdata.guaranteed = FALSE, na.action = na.pass, unconditional = FALSE, iterms.type = NULL, interval = c("none", "confidence", "prediction"), level = 0.95, tvalue = NULL, ... )
object |
object of class ‘gam’ or as returned by function ‘gamm’ |
newdata |
see |
type |
see |
se.fit |
see |
terms |
see |
exclude |
see |
block.size |
see |
newdata.guaranteed |
see |
na.action |
see |
unconditional |
see |
iterms.type |
see |
interval |
either ‘none’, ‘confidence’ or ‘prediction’. |
level |
probability level for the interval (default 0.95) |
tvalue |
t-value statistic used for constructing the intervals |
... |
additional arguments to be passed to |
numeric vector of the same length as the fitted object when interval is equal to ‘none’. Otherwise, a data.frame with columns named (for a 0.95 level) ‘Estimate’, ‘Est.Error’, ‘Q2.5’ and ‘Q97.5’
this is a very simple wrapper for predict.gam.
predict.lm, predict.nls, predict.gam, simulate_nls, simulate_gam
require(ggplot2) require(mgcv) data(barley) fm.G <- gam(yield ~ s(NF, k = 6), data = barley) ## confidence and prediction intervals cis <- predict_gam(fm.G, interval = "conf") pis <- predict_gam(fm.G, interval = "pred") barleyA.ci <- cbind(barley, cis) barleyA.pi <- cbind(barley, pis) ggplot() + geom_point(data = barleyA.ci, aes(x = NF, y = yield)) + geom_line(data = barleyA.ci, aes(x = NF, y = Estimate)) + geom_ribbon(data = barleyA.ci, aes(x = NF, ymin = Q2.5, ymax = Q97.5), color = "red", alpha = 0.3) + geom_ribbon(data = barleyA.pi, aes(x = NF, ymin = Q2.5, ymax = Q97.5), color = "blue", alpha = 0.3) + ggtitle("95% confidence and prediction bands")require(ggplot2) require(mgcv) data(barley) fm.G <- gam(yield ~ s(NF, k = 6), data = barley) ## confidence and prediction intervals cis <- predict_gam(fm.G, interval = "conf") pis <- predict_gam(fm.G, interval = "pred") barleyA.ci <- cbind(barley, cis) barleyA.pi <- cbind(barley, pis) ggplot() + geom_point(data = barleyA.ci, aes(x = NF, y = yield)) + geom_line(data = barleyA.ci, aes(x = NF, y = Estimate)) + geom_ribbon(data = barleyA.ci, aes(x = NF, ymin = Q2.5, ymax = Q97.5), color = "red", alpha = 0.3) + geom_ribbon(data = barleyA.pi, aes(x = NF, ymin = Q2.5, ymax = Q97.5), color = "blue", alpha = 0.3) + ggtitle("95% confidence and prediction bands")
Computes weights based on AIC, AICc, or BIC and it generates weighted predictions by the relative value of the IC values
predict function for objects of class lme
predict function for objects of class gnls
predict function for objects of class gls
predict_nlme( ..., criteria = c("AIC", "AICc", "BIC"), interval = c("none", "confidence", "prediction", "new-prediction"), level = 0.95, nsim = 1000, plevel = 0, newdata = NULL, weights ) predict_lme( ..., criteria = c("AIC", "AICc", "BIC"), interval = c("none", "confidence", "prediction", "new-prediction"), level = 0.95, nsim = 1000, plevel = 0, newdata = NULL, weights ) predict_gnls( ..., criteria = c("AIC", "AICc", "BIC"), interval = c("none", "confidence", "prediction", "new-prediction"), level = 0.95, nsim = 1000, plevel = 0, newdata = NULL, weights ) predict_gls( ..., criteria = c("AIC", "AICc", "BIC"), interval = c("none", "confidence", "prediction", "new-prediction"), level = 0.95, nsim = 1000, plevel = 0, newdata = NULL, weights )predict_nlme( ..., criteria = c("AIC", "AICc", "BIC"), interval = c("none", "confidence", "prediction", "new-prediction"), level = 0.95, nsim = 1000, plevel = 0, newdata = NULL, weights ) predict_lme( ..., criteria = c("AIC", "AICc", "BIC"), interval = c("none", "confidence", "prediction", "new-prediction"), level = 0.95, nsim = 1000, plevel = 0, newdata = NULL, weights ) predict_gnls( ..., criteria = c("AIC", "AICc", "BIC"), interval = c("none", "confidence", "prediction", "new-prediction"), level = 0.95, nsim = 1000, plevel = 0, newdata = NULL, weights ) predict_gls( ..., criteria = c("AIC", "AICc", "BIC"), interval = c("none", "confidence", "prediction", "new-prediction"), level = 0.95, nsim = 1000, plevel = 0, newdata = NULL, weights )
... |
nlme, lme, gls or gnls objects. |
criteria |
either ‘AIC’, ‘AICc’ or ‘BIC’. |
interval |
either ‘none’, ‘confidence’ or ‘prediction’. It is also possible to choose ‘new-prediction’, which is a prediction that resamples the random effects (only relevant for ‘lme’ or ‘nlme’ objects.) |
level |
probability level for the interval (default 0.95) |
nsim |
number of simulations to perform for intervals. Default 1000. |
plevel |
parameter level prediction to be passed to prediciton functions. |
newdata |
new data frame for predictions |
weights |
vector of weights of the same length as the number of models. It should sum up to one and it will override the information-criteria based weights. The weights should match the order of the models. |
numeric vector of the same length as the fitted object.
all the objects should be fitted to the same data. The weights are based on the IC value.
predict.nlme predict.lme predict.gnls
## Example require(ggplot2) require(nlme) data(Orange) ## All models should be fitted using Maximum Likelihood fm.L <- nlme(circumference ~ SSlogis(age, Asym, xmid, scal), random = pdDiag(Asym + xmid + scal ~ 1), method = "ML", data = Orange) fm.G <- nlme(circumference ~ SSgompertz(age, Asym, b2, b3), random = pdDiag(Asym + b2 + b3 ~ 1), method = "ML", data = Orange) fm.F <- nlme(circumference ~ SSfpl(age, A, B, xmid, scal), random = pdDiag(A + B + xmid + scal ~ 1), method = "ML", data = Orange) fm.B <- nlme(circumference ~ SSbg4rp(age, w.max, lt.e, ldtm, ldtb), random = pdDiag(w.max + lt.e + ldtm + ldtb ~ 1), method = "ML", data = Orange) ## Print the table with weights IC_tab(fm.L, fm.G, fm.F, fm.B) ## Each model prediction is weighted according to their AIC values prd <- predict_nlme(fm.L, fm.G, fm.F, fm.B) ggplot(data = Orange, aes(x = age, y = circumference)) + geom_point() + geom_line(aes(y = predict(fm.L, level = 0), color = "Logistic")) + geom_line(aes(y = predict(fm.G, level = 0), color = "Gompertz")) + geom_line(aes(y = predict(fm.F, level = 0), color = "4P-Logistic")) + geom_line(aes(y = predict(fm.B, level = 0), color = "Beta")) + geom_line(aes(y = prd, color = "Avg. Model"), linewidth = 1.2)## Example require(ggplot2) require(nlme) data(Orange) ## All models should be fitted using Maximum Likelihood fm.L <- nlme(circumference ~ SSlogis(age, Asym, xmid, scal), random = pdDiag(Asym + xmid + scal ~ 1), method = "ML", data = Orange) fm.G <- nlme(circumference ~ SSgompertz(age, Asym, b2, b3), random = pdDiag(Asym + b2 + b3 ~ 1), method = "ML", data = Orange) fm.F <- nlme(circumference ~ SSfpl(age, A, B, xmid, scal), random = pdDiag(A + B + xmid + scal ~ 1), method = "ML", data = Orange) fm.B <- nlme(circumference ~ SSbg4rp(age, w.max, lt.e, ldtm, ldtb), random = pdDiag(w.max + lt.e + ldtm + ldtb ~ 1), method = "ML", data = Orange) ## Print the table with weights IC_tab(fm.L, fm.G, fm.F, fm.B) ## Each model prediction is weighted according to their AIC values prd <- predict_nlme(fm.L, fm.G, fm.F, fm.B) ggplot(data = Orange, aes(x = age, y = circumference)) + geom_point() + geom_line(aes(y = predict(fm.L, level = 0), color = "Logistic")) + geom_line(aes(y = predict(fm.G, level = 0), color = "Gompertz")) + geom_line(aes(y = predict(fm.F, level = 0), color = "4P-Logistic")) + geom_line(aes(y = predict(fm.B, level = 0), color = "Beta")) + geom_line(aes(y = prd, color = "Avg. Model"), linewidth = 1.2)
Computes weights based on AIC, AICc, or BIC and it generates weighted predictions by the relative value of the IC values
predict function for objects of class gam
predict_nls( ..., criteria = c("AIC", "AICc", "BIC"), interval = c("none", "confidence", "prediction"), level = 0.95, nsim = 1000, resid.type = c("none", "resample", "normal", "wild"), newdata = NULL, weights ) predict2_gam( ..., criteria = c("AIC", "AICc", "BIC"), interval = c("none", "confidence", "prediction"), level = 0.95, nsim = 1000, resid.type = c("none", "resample", "normal", "wild"), newdata = NULL, weights )predict_nls( ..., criteria = c("AIC", "AICc", "BIC"), interval = c("none", "confidence", "prediction"), level = 0.95, nsim = 1000, resid.type = c("none", "resample", "normal", "wild"), newdata = NULL, weights ) predict2_gam( ..., criteria = c("AIC", "AICc", "BIC"), interval = c("none", "confidence", "prediction"), level = 0.95, nsim = 1000, resid.type = c("none", "resample", "normal", "wild"), newdata = NULL, weights )
... |
‘nls’ or ‘lm’ objects (‘glm’ and ‘gam’ objects inherit ‘lm’). |
criteria |
either ‘AIC’, ‘AICc’ or ‘BIC’. |
interval |
either ‘none’, ‘confidence’ or ‘prediction’. |
level |
probability level for the interval (default 0.95) |
nsim |
number of simulations to perform for intervals. Default 1000. |
resid.type |
either ‘none’, “resample”, “normal” or “wild”. |
newdata |
new data frame for predictions |
weights |
vector of weights of the same length as the number of models. It should sum up to one and it will override the information-criteria based weights. The weights should match the order of the models. |
numeric vector of the same length as the fitted object when interval is equal to ‘none’. Otherwise, a data.frame with columns named (for a 0.95 level) ‘Estimate’, ‘Est.Error’, ‘Q2.5’ and ‘Q97.5’
all the objects should be fitted to the same data. Weights are
based on the chosen IC value (exp(-0.5 * delta IC)).
For models of class gam there is very limited support.
predict.lm, predict.nls, predict.gam, simulate_nls, simulate_gam
## Example require(ggplot2) require(mgcv) data(barley, package = "nlraa") fm.L <- lm(yield ~ NF, data = barley) fm.Q <- lm(yield ~ NF + I(NF^2), data = barley) fm.A <- nls(yield ~ SSasymp(NF, Asym, R0, lrc), data = barley) fm.LP <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley) fm.QP <- nls(yield ~ SSquadp3(NF, a, b, c), data = barley) fm.BL <- nls(yield ~ SSblin(NF, a, b, xs, c), data = barley) fm.G <- gam(yield ~ s(NF, k = 6), data = barley) ## Print the table with weights IC_tab(fm.L, fm.Q, fm.A, fm.LP, fm.QP, fm.BL, fm.G) ## Each model prediction is weighted according to their AIC values prd <- predict_nls(fm.L, fm.Q, fm.A, fm.LP, fm.QP, fm.BL, fm.G) ggplot(data = barley, aes(x = NF, y = yield)) + geom_point() + geom_line(aes(y = fitted(fm.L), color = "Linear")) + geom_line(aes(y = fitted(fm.Q), color = "Quadratic")) + geom_line(aes(y = fitted(fm.A), color = "Asymptotic")) + geom_line(aes(y = fitted(fm.LP), color = "Linear-plateau")) + geom_line(aes(y = fitted(fm.QP), color = "Quadratic-plateau")) + geom_line(aes(y = fitted(fm.BL), color = "Bi-linear")) + geom_line(aes(y = fitted(fm.G), color = "GAM")) + geom_line(aes(y = prd, color = "Avg. Model"), linewidth = 1.2)## Example require(ggplot2) require(mgcv) data(barley, package = "nlraa") fm.L <- lm(yield ~ NF, data = barley) fm.Q <- lm(yield ~ NF + I(NF^2), data = barley) fm.A <- nls(yield ~ SSasymp(NF, Asym, R0, lrc), data = barley) fm.LP <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley) fm.QP <- nls(yield ~ SSquadp3(NF, a, b, c), data = barley) fm.BL <- nls(yield ~ SSblin(NF, a, b, xs, c), data = barley) fm.G <- gam(yield ~ s(NF, k = 6), data = barley) ## Print the table with weights IC_tab(fm.L, fm.Q, fm.A, fm.LP, fm.QP, fm.BL, fm.G) ## Each model prediction is weighted according to their AIC values prd <- predict_nls(fm.L, fm.Q, fm.A, fm.LP, fm.QP, fm.BL, fm.G) ggplot(data = barley, aes(x = NF, y = yield)) + geom_point() + geom_line(aes(y = fitted(fm.L), color = "Linear")) + geom_line(aes(y = fitted(fm.Q), color = "Quadratic")) + geom_line(aes(y = fitted(fm.A), color = "Asymptotic")) + geom_line(aes(y = fitted(fm.LP), color = "Linear-plateau")) + geom_line(aes(y = fitted(fm.QP), color = "Quadratic-plateau")) + geom_line(aes(y = fitted(fm.BL), color = "Bi-linear")) + geom_line(aes(y = fitted(fm.G), color = "GAM")) + geom_line(aes(y = prd, color = "Avg. Model"), linewidth = 1.2)
The method used in this function is described in Battes and Watts (2007) Nonlinear Regression Analysis and Its Applications (see pages 58-59). It is known as the Delta method.
predict2_nls( object, newdata = NULL, interval = c("none", "confidence", "prediction"), level = 0.95 )predict2_nls( object, newdata = NULL, interval = c("none", "confidence", "prediction"), level = 0.95 )
object |
object of class ‘nls’ |
newdata |
data frame with values for the predictor |
interval |
either ‘none’, ‘confidence’ or ‘prediction’ |
level |
probability level (default is 0.95) |
This method is approximate and it works better when the distribution of the parameter estimates are normally distributed. This assumption can be evaluated by using bootstrap.
The method currently works well for any nonlinear function, but if predictions are needed on new data, then it is required that a selfStart function is used.
a data frame with Estimate, Est.Error, lower interval bound and upper interval bound. For example, if the level = 0.95, the lower bound would be named Q2.5 and the upper bound would be name Q97.5
require(ggplot2) require(nlme) data(Soybean) SoyF <- subset(Soybean, Variety == "F" & Year == 1988) fm1 <- nls(weight ~ SSlogis(Time, Asym, xmid, scal), data = SoyF) ## The SSlogis also supplies analytical derivatives ## therefore the predict function returns the gradient too prd1 <- predict(fm1, newdata = SoyF) ## Gradient head(attr(prd1, "gradient")) ## Prediction method using gradient prds <- predict2_nls(fm1, interval = "conf") SoyFA <- cbind(SoyF, prds) ggplot(data = SoyFA, aes(x = Time, y = weight)) + geom_point() + geom_line(aes(y = Estimate)) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) + ggtitle("95% Confidence Bands") ## This is equivalent fm2 <- nls(weight ~ Asym/(1 + exp((xmid - Time)/scal)), data = SoyF, start = c(Asym = 20, xmid = 56, scal = 8)) ## Prediction interval prdi <- predict2_nls(fm1, interval = "pred") SoyFA.PI <- cbind(SoyF, prdi) ## Make prediction interval plot ggplot(data = SoyFA.PI, aes(x = Time, y = weight)) + geom_point() + geom_line(aes(y = Estimate)) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) + ggtitle("95% Prediction Band") ## For these data we should be using gnls instead with an increasing variance fmg1 <- gnls(weight ~ SSlogis(Time, Asym, xmid, scal), data = SoyF, weights = varPower()) IC_tab(fm1, fmg1) prdg <- predict_gnls(fmg1, interval = "pred") SoyFA.GPI <- cbind(SoyF, prdg) ## These prediction bands are not perfect, but they could be smoothed ## to eliminate the ragged appearance ggplot(data = SoyFA.GPI, aes(x = Time, y = weight)) + geom_point() + geom_line(aes(y = Estimate)) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) + ggtitle("95% Prediction Band. NLS model which \n accomodates an increasing variance")require(ggplot2) require(nlme) data(Soybean) SoyF <- subset(Soybean, Variety == "F" & Year == 1988) fm1 <- nls(weight ~ SSlogis(Time, Asym, xmid, scal), data = SoyF) ## The SSlogis also supplies analytical derivatives ## therefore the predict function returns the gradient too prd1 <- predict(fm1, newdata = SoyF) ## Gradient head(attr(prd1, "gradient")) ## Prediction method using gradient prds <- predict2_nls(fm1, interval = "conf") SoyFA <- cbind(SoyF, prds) ggplot(data = SoyFA, aes(x = Time, y = weight)) + geom_point() + geom_line(aes(y = Estimate)) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) + ggtitle("95% Confidence Bands") ## This is equivalent fm2 <- nls(weight ~ Asym/(1 + exp((xmid - Time)/scal)), data = SoyF, start = c(Asym = 20, xmid = 56, scal = 8)) ## Prediction interval prdi <- predict2_nls(fm1, interval = "pred") SoyFA.PI <- cbind(SoyF, prdi) ## Make prediction interval plot ggplot(data = SoyFA.PI, aes(x = Time, y = weight)) + geom_point() + geom_line(aes(y = Estimate)) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) + ggtitle("95% Prediction Band") ## For these data we should be using gnls instead with an increasing variance fmg1 <- gnls(weight ~ SSlogis(Time, Asym, xmid, scal), data = SoyF, weights = varPower()) IC_tab(fm1, fmg1) prdg <- predict_gnls(fmg1, interval = "pred") SoyFA.GPI <- cbind(SoyF, prdg) ## These prediction bands are not perfect, but they could be smoothed ## to eliminate the ragged appearance ggplot(data = SoyFA.GPI, aes(x = Time, y = weight)) + geom_point() + geom_line(aes(y = Estimate)) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) + ggtitle("95% Prediction Band. NLS model which \n accomodates an increasing variance")
R-squared ‘modified’ for nonlinear (mixed) models
R2M(x, ...) ## S3 method for class 'nls' R2M(x, ...) ## S3 method for class 'lm' R2M(x, ...) ## S3 method for class 'gls' R2M(x, ...) ## S3 method for class 'gnls' R2M(x, ...) ## S3 method for class 'lme' R2M(x, ...) ## S3 method for class 'nlme' R2M(x, ...)R2M(x, ...) ## S3 method for class 'nls' R2M(x, ...) ## S3 method for class 'lm' R2M(x, ...) ## S3 method for class 'gls' R2M(x, ...) ## S3 method for class 'gnls' R2M(x, ...) ## S3 method for class 'lme' R2M(x, ...) ## S3 method for class 'nlme' R2M(x, ...)
x |
object of class ‘lm’, ‘nls’, ‘gls’, ‘gnls’, ‘lme’ or ‘nlme’ . |
... |
additional arguments (none use at the moment). |
I have read some papers about computing an R-squared for (non)linear (mixed) models
and I am not sure that it makes sense at all. However, here they are and
the method is general enough that it extends to nonlinear mixed models. What do
these numbers mean and why would you want to compute them are good questions to
ponder...
Recommended reading:
Nakagawa and Schielzeth Methods in Ecology and Evolution doi:10.1111/j.2041-210x.2012.00261.x
https://stats.oarc.ucla.edu/other/mult-pkg/faq/general/faq-what-are-pseudo-r-squareds/
Spiess, AN., Neumeyer, N. An evaluation of R2 as an inadequate measure for nonlinear models in
pharmacological and biochemical research: a Monte Carlo approach. BMC Pharmacol 10, 6 (2010).
doi:10.1186/1471-2210-10-6
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q1/003363.html
Other R pacakges which calculate some version of an R-squared: performance, rcompanion, MuMIn
it returns a list with the following structure:
for an object of class ‘lm’, ‘nls’, ‘gls’ or ‘gnls’,
R2: R-squared
var.comp: variance components with var.fixed and var.resid
var.perc: variance components percents (should add up to 100)
for an object of class ‘lme’ or ‘nlme’ in addition it also returns:
R2.marginal: marginal R2 which only includes the fixed effects
R2.conditional: conditional R2 which includes both the fixed and random effects
var.random: the variance contribution of the random effects
The references here strongly discourage the use of R-squared in anything but linear models.
require(nlme) data(barley, package = "nlraa") fit2 <- lme(yield ~ NF + I(NF^2), random = ~ 1 | year, data = barley) R2M(fit2) ## Nonlinear Mixed Model barleyG <- groupedData(yield ~ NF | year, data = barley) fit3L <- nlsLMList(yield ~ SSquadp3(NF, a, b, c), data = barleyG) fit3 <- nlme(fit3L, random = pdDiag(a + b ~ 1)) R2M(fit3)require(nlme) data(barley, package = "nlraa") fit2 <- lme(yield ~ NF + I(NF^2), random = ~ 1 | year, data = barley) R2M(fit2) ## Nonlinear Mixed Model barleyG <- groupedData(yield ~ NF | year, data = barley) fit3L <- nlsLMList(yield ~ SSquadp3(NF, a, b, c), data = barleyG) fit3 <- nlme(fit3L, random = pdDiag(a + b ~ 1)) R2M(fit3)
gam
By sampling from the vector of coefficients it is possible to simulate data from a ‘gam’ model.
simulate_gam( object, nsim = 1, psim = 1, resid.type = c("none", "resample", "normal", "wild"), value = c("matrix", "data.frame"), ... )simulate_gam( object, nsim = 1, psim = 1, resid.type = c("none", "resample", "normal", "wild"), value = c("matrix", "data.frame"), ... )
object |
|
nsim |
number of simulations to perform |
psim |
parameter simulation level (an integer, 0, 1, 2 or 3). |
resid.type |
type of residual to use. ‘none’, ‘resample’, ‘normal’ or ‘wild’. |
value |
either ‘matrix’ or ‘data.frame’ |
... |
additional arguments (none used at the moment) |
This function is probably redundant. Simply using simulate
generates data from the correct distribution for objects which inherit
class lm. The difference is that I'm trying to add the
uncertainty in the parameter estimates.
These are the options that control the parameter simulation level
returns the fitted values
simulates from a beta vector (mean response)
simulates observations according to the residual type (similar to observed data)
simulates a beta vector, considers uncertainty in the variance covariance matrix of beta and adds residuals (prediction)
The residual type (resid.type) controls how the residuals are generated. They are either resampled, simulated from a normal distribution or ‘wild’ where the Rademacher distribution is used (https://en.wikipedia.org/wiki/Rademacher_distribution). Resampled and normal both assume iid, but ‘normal’ makes the stronger assumption of normality. ‘wild’ does not assume constant variance, but it assumes symmetry.
matrix or data.frame with responses
psim = 3 is not implemented at the moment.
The purpose of this function is to make it compatible with other functions in this package. It has some limitations compared to the functions in the ‘see also’ section.
Generalized Additive Models. An Introduction with R. Second Edition. (2017) Simon N. Wood. CRC Press.
predict, predict.gam, simulate and simulate_lm.
require(ggplot2) require(mgcv) ## These count data are from GAM book by Simon Wood (pg. 132) - see reference y <- c(12, 14, 33, 50, 67, 74, 123, 141, 165, 204, 253, 246, 240) t <- 1:13 dat <- data.frame(y = y, t = t) fit <- gam(y ~ t + I(t^2), family = poisson, data = dat) sims <- simulate_gam(fit, nsim = 100, value = "data.frame") ggplot(data = sims) + geom_line(aes(x = t, y = sim.y, group = ii), color = "gray", alpha = 0.5) + geom_point(aes(x = t, y = y))require(ggplot2) require(mgcv) ## These count data are from GAM book by Simon Wood (pg. 132) - see reference y <- c(12, 14, 33, 50, 67, 74, 123, 141, 165, 204, 253, 246, 240) t <- 1:13 dat <- data.frame(y = y, t = t) fit <- gam(y ~ t + I(t^2), family = poisson, data = dat) sims <- simulate_gam(fit, nsim = 100, value = "data.frame") ggplot(data = sims) + geom_line(aes(x = t, y = sim.y, group = ii), color = "gray", alpha = 0.5) + geom_point(aes(x = t, y = y))
gls
Simulate values from an object of class gls. Unequal variances,
as modeled using the ‘weights’ option are supported, and there is experimental
code for dealing with the ‘correlation’ structure. This generates just one simulation
from these type of models. To generate multiple simulations use simulate_lme
simulate_gls( object, psim = 1, na.action = na.fail, naPattern = NULL, data = NULL, ... )simulate_gls( object, psim = 1, na.action = na.fail, naPattern = NULL, data = NULL, ... )
object |
object of class |
psim |
parameter simulation level, 0: for fitted values, 1: for simulation from fixed parameters (assuming a fixed vcov matrix), 2: for simulation considering the uncertainty in the residual standard error (sigma), this returns data which will appear similar to the observed values |
na.action |
default ‘na.fail’. See |
naPattern |
missing value pattern. See |
data |
the data argument is needed when using this function inside user defined functions. It should be identical to the data used to fit the model. |
... |
additional arguments (it is possible to supply a newdata this way) |
This function is based on predict.gls function
It uses function mvrnorm to generate new values for the coefficients
of the model using the Variance-Covariance matrix vcov. This variance-covariance matrix
refers to the one for the parameters ‘beta’, not the one for the residuals.
It returns a vector with simulated values with length equal to the number of rows in the original data
require(nlme) data(Orange) fit.gls <- gls(circumference ~ age, data = Orange, weights = varPower()) ## Visualize covariance matrix fit.gls.vc <- var_cov(fit.gls) image(log(fit.gls.vc[,ncol(fit.gls.vc):1])) sim <- simulate_gls(fit.gls)require(nlme) data(Orange) fit.gls <- gls(circumference ~ age, data = Orange, weights = varPower()) ## Visualize covariance matrix fit.gls.vc <- var_cov(fit.gls) image(log(fit.gls.vc[,ncol(fit.gls.vc):1])) sim <- simulate_gls(fit.gls)
gnls
Simulate values from an object of class gnls. Unequal variances, as modeled using the ‘weights’ option are supported, and there is experimental code for dealing with the ‘correlation’ structure.
simulate_gnls( object, psim = 1, na.action = na.fail, naPattern = NULL, data = NULL, ... )simulate_gnls( object, psim = 1, na.action = na.fail, naPattern = NULL, data = NULL, ... )
object |
object of class |
psim |
parameter simulation level, 0: for fitted values, 1: for simulation from fixed parameters (assuming a fixed vcov matrix), 2: for simulation considering the uncertainty in the residual standard error (sigma), this returns data which will appear similar to the observed values |
na.action |
default ‘na.fail’. See |
naPattern |
missing value pattern. See |
data |
the data argument is needed when using this function inside user defined functions. It should be identical to the data used to fit the model. |
... |
additional arguments (it is possible to supply a newdata this way) |
This function is based on predict.gnls function
It uses function mvrnorm to generate new values for the coefficients
of the model using the Variance-Covariance matrix vcov. This variance-covariance matrix
refers to the one for the parameters ‘beta’, not the one for the residuals.
It returns a vector with simulated values with length equal to the number of rows in the original data
require(nlme) data(barley, package = "nlraa") fit.gnls <- gnls(yield ~ SSlinp(NF, a, b, xs), data = barley) sim <- simulate_gnls(fit.gnls)require(nlme) data(barley, package = "nlraa") fit.gnls <- gnls(yield ~ SSlinp(NF, a, b, xs), data = barley) sim <- simulate_gnls(fit.gnls)
lm
The function simulate does not consider the
uncertainty in the estimation of the model parameters. This function will attempt
to do this.
simulate_lm( object, psim = 1, nsim = 1, resid.type = c("none", "resample", "normal", "wild"), value = c("matrix", "data.frame"), data = NULL, ... )simulate_lm( object, psim = 1, nsim = 1, resid.type = c("none", "resample", "normal", "wild"), value = c("matrix", "data.frame"), data = NULL, ... )
object |
object of class |
psim |
parameter simulation level (an integer, 0, 1, 2, 3 or 4). |
nsim |
number of simulations to perform |
resid.type |
type of residual to include (none, resample, normal or wild) |
value |
either ‘matrix’ or ‘data.frame’ |
data |
the data argument might be needed when using this function inside user defined functions. At least it is expected to be safer. |
... |
additional arguments (none used at the moment) |
Simulate responses from a linear model lm
These are the options that control the parameter simulation level
returns the fitted values
simulates a beta vector (mean response)
simulates a beta vector and adds resampled residuals (similar to observed data)
simulates a beta vector, considers uncertainty in the variance covariance matrix of beta and adds residuals (prediction)
only adds residuals according to resid.type (similar to simulate.lm)
The residual type (resid.type) controls how the residuals are generated.
They are either resampled, simulated from a normal distribution or ‘wild’ where the
Rademacher distribution is used (https://en.wikipedia.org/wiki/Rademacher_distribution).
Resampled and normal both assume iid, but ‘normal’ makes the stronger assumption of normality.
When psim = 2 and resid.type = none, simulate is used instead.
‘wild’ does not assume constant variance, but it assumes symmetry.
matrix or data.frame with responses
See “Inference Based on the Wild Bootstrap” James G. MacKinnon https://www.math.kth.se/matstat/gru/sf2930/papers/wild.bootstrap.pdf “Bootstrap in Nonstationary Autoregression” Zuzana Praskova https://dml.cz/bitstream/handle/10338.dmlcz/135473/Kybernetika_38-2002-4_1.pdf “Jackknife, Bootstrap and other Resampling Methods in Regression Analysis” C. F. J. Wu. The Annals of Statistics. 1986. Vol 14. 1261-1295.
require(ggplot2) data(Orange) fit <- lm(circumference ~ age, data = Orange) sims <- simulate_lm(fit, nsim = 100, value = "data.frame") ggplot(data = sims) + geom_line(aes(x = age, y = sim.y, group = ii), color = "gray", alpha = 0.5) + geom_point(aes(x = age, y = circumference))require(ggplot2) data(Orange) fit <- lm(circumference ~ age, data = Orange) sims <- simulate_lm(fit, nsim = 100, value = "data.frame") ggplot(data = sims) + geom_line(aes(x = age, y = sim.y, group = ii), color = "gray", alpha = 0.5) + geom_point(aes(x = age, y = circumference))
lme
Simulate values from an object of class lme. Unequal variances, as modeled using the ‘weights’ option are supported, and there is experimental code for considering the ‘correlation’ structure.
simulate_lme( object, nsim = 1, psim = 1, value = c("matrix", "data.frame"), data = NULL, ... )simulate_lme( object, nsim = 1, psim = 1, value = c("matrix", "data.frame"), data = NULL, ... )
object |
|
nsim |
number of samples, default 1 |
psim |
parameter simulation level, 0: for fitted values, 1: for simulation from fixed parameters (assuming a fixed vcov matrix), 2: for simulation considering the uncertainty in the residual standard error (sigma), this returns data which will appear similar to the observed values. 3: in addition samples a new set of random effects. |
value |
whether to return a matrix (default) or an augmented data frame |
data |
the data argument is needed when using this function inside user defined functions. |
... |
additional arguments (it is possible to supply a newdata this way) |
This function is based on predict.lme function
It uses function mvrnorm to generate new values for the coefficients
of the model using the Variance-Covariance matrix vcov. This variance-covariance matrix
refers to the one for the parameters 'beta', not the one for the residuals.
It returns a vector with simulated values with length equal to the number of rows in the original data
I find the simulate.merMod in the lme4 pacakge confusing. There is use.u and several versions of re.form. From the documentation it seems that if use.u = TRUE, then the current values of the random effects are used. This would mean that it is equivalent to psim = 2 in this function. Then use.u = FALSE, would be equivalent to psim = 3. re.form allows for specifying the formula of the random effects.
predict.lme and ‘simulate.merMod’ in the ‘lme4’ package.
require(nlme) data(Orange) fm1 <- lme(circumference ~ age, random = ~ 1 | Tree, data = Orange) sims <- simulate_lme(fm1, nsim = 10)require(nlme) data(Orange) fm1 <- lme(circumference ~ age, random = ~ 1 | Tree, data = Orange) sims <- simulate_lme(fm1, nsim = 10)
Simulate multiple samples from a nonlinear model
simulate_nlme( object, nsim = 1, psim = 1, value = c("matrix", "data.frame"), data = NULL, ... )simulate_nlme( object, nsim = 1, psim = 1, value = c("matrix", "data.frame"), data = NULL, ... )
object |
|
nsim |
number of samples, default 1 |
psim |
simulation level for fixed and random parameters see |
value |
whether to return a matrix (default) or an augmented data frame |
data |
the data argument is needed when using this function inside user defined functions. |
... |
additional arguments to be passed to either |
The details can be found in either simulate_gnls or simulate_nlme_one.
This function is very simple and it only sets up a matrix and a loop in order to simulate several instances of
model outputs.
It returns a matrix with simulated values from the original object
with number of rows equal to the number of rows of fitted and number
of columns equal to the number of simulated samples (‘nsim’). In the case of ‘data.frame’
it returns an augmented data.frame, which can potentially be a very large object, but which
makes furhter plotting more convenient.
require(nlme) data(barley, package = "nlraa") barley2 <- subset(barley, year < 1974) fit.lp.gnls2 <- gnls(yield ~ SSlinp(NF, a, b, xs), data = barley2) barley2$year.f <- as.factor(barley2$year) cfs <- coef(fit.lp.gnls2) fit.lp.gnls3 <- update(fit.lp.gnls2, params = list(a + b + xs ~ year.f), start = c(cfs[1], 0, 0, 0, cfs[2], 0, 0, 0, cfs[3], 0, 0, 0)) sims <- simulate_nlme(fit.lp.gnls3, nsim = 3)require(nlme) data(barley, package = "nlraa") barley2 <- subset(barley, year < 1974) fit.lp.gnls2 <- gnls(yield ~ SSlinp(NF, a, b, xs), data = barley2) barley2$year.f <- as.factor(barley2$year) cfs <- coef(fit.lp.gnls2) fit.lp.gnls3 <- update(fit.lp.gnls2, params = list(a + b + xs ~ year.f), start = c(cfs[1], 0, 0, 0, cfs[2], 0, 0, 0, cfs[3], 0, 0, 0)) sims <- simulate_nlme(fit.lp.gnls3, nsim = 3)
nlme
This function is based on predict.nlme function
simulate_nlme_one( object, psim = 1, level = Q, asList = FALSE, na.action = na.fail, naPattern = NULL, data = NULL, ... )simulate_nlme_one( object, psim = 1, level = Q, asList = FALSE, na.action = na.fail, naPattern = NULL, data = NULL, ... )
object |
object of class |
psim |
parameter simulation level, 0: for fitted values, 1: for simulation from fixed parameters (assuming a fixed vcov matrix), 2: for simulation considering the residual error (sigma), this returns data which will appear similar to the observed values. Currently, working on psim = 3, which will simulate a new set of random effects. This can be useful when computing prediction intervals at the subject-level. |
level |
level at which simulations are performed. See |
asList |
optional logical value. See |
na.action |
missing value action. See |
naPattern |
missing value pattern. See |
data |
the data argument is needed when using this function inside user defined functions. |
... |
additional arguments to be passed (possible to pass newdata this way) |
It uses function mvrnorm to generate new values for the coefficients
of the model using the Variance-Covariance matrix vcov
This function should return a vector with the same dimensions as the original data, unless newdata is provided.
nls
Simulate values from an object of class nls.
simulate_nls( object, nsim = 1, psim = 1, resid.type = c("none", "resample", "normal", "wild"), value = c("matrix", "data.frame"), data = NULL, ... )simulate_nls( object, nsim = 1, psim = 1, resid.type = c("none", "resample", "normal", "wild"), value = c("matrix", "data.frame"), data = NULL, ... )
object |
object of class |
nsim |
number of simulations to perform |
psim |
parameter simulation level, 0: for fitted values, 1: for simulation from fixed parameters (assuming a fixed vcov matrix), 2: simulation from sampling both from the parameters and the residuals, 3: for simulation considering the uncertainty in the residual standard error only (sigma) and fixing the parameter estimates at their original value; this will result in simulations similar to the observed values. |
resid.type |
either ‘none’, “resample”, “normal” or “wild”. |
value |
either ‘matrix’ or ‘data.frame’ |
data |
the data argument is needed when using this function inside user defined functions. |
... |
additional arguments (it is possible to supply a newdata this way) |
This function is based on predict.gnls function
It uses function mvrnorm to generate new values for the coefficients
of the model using the Variance-Covariance matrix vcov. This variance-covariance matrix
refers to the one for the parameters ‘beta’, not the one for the residuals.
It returns a vector with simulated values with length equal to the number of rows in the original data
The default behavior is that simulations are perfomed for the mean function only.
When ‘psim = 2’ this function will silently choose ‘resample’ as the
‘resid.type’. This is not ideal design for this function, but I made this choice for
compatibility with other types of simulation originating from glm and
gam.
data(barley, package = "nlraa") fit <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley) sim <- simulate_nls(fit, nsim = 100)data(barley, package = "nlraa") fit <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley) sim <- simulate_nls(fit, nsim = 100)
Sorghum and Maize growth in Greece
smsm
A data frame with 235 rows and 5 columns
-integer- Day of the year 141-303
-integer- Block in the experimental design 1-4
-integer- Input level 1 (Low) or 2 (High)
-factor- either F (Fiber Sorghum), M (Maize), S (Sweet Sorghum)
-numeric- Biomass yield in Mg/ha
A dataset containing growth data for sorghum and maize in Greece.
Danalatos, N.G., S.V. Archontoulis, and K. Tsiboukas. 2009. Comparative analysis of sorghum versus corn growing under optimum and under water/nitrogen limited conditions in central Greece. In: From research to industry and markets: Proceedings of the 17th European Biomass Conference, Hamburg, Germany. 29 June–3 July 2009. ETA–Renewable Energies, Florence, Italy. p. 538–544.
See above reference. (Currently available on ResearchGate).
Self starter for a type of bell-shaped curve
agauss(x, eta, beta, delta, sigma1, sigma2) SSagauss(x, eta, beta, delta, sigma1, sigma2)agauss(x, eta, beta, delta, sigma1, sigma2) SSagauss(x, eta, beta, delta, sigma1, sigma2)
x |
input vector |
eta |
maximum value of y |
beta |
parameter controlling the lower values |
delta |
break-point separating the first and second half-bell curve |
sigma1 |
scale for the first half |
sigma2 |
scale for the second half |
This function is described in (doi:10.3390/rs12050827).
a numeric vector of the same length as x containing parameter estimates for equation specified
agauss: vector of the same length as x using an asymmetric bell-shaped Gaussian curve
require(ggplot2) set.seed(1234) x <- 1:30 y <- agauss(x, 10, 2, 10, 2, 6) + rnorm(length(x), 0, 0.02) dat <- data.frame(x = x, y = y) fit <- minpack.lm::nlsLM(y ~ SSagauss(x, eta, beta, delta, sigma1, sigma2), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) set.seed(1234) x <- 1:30 y <- agauss(x, 10, 2, 10, 2, 6) + rnorm(length(x), 0, 0.02) dat <- data.frame(x = x, y = y) fit <- minpack.lm::nlsLM(y ~ SSagauss(x, eta, beta, delta, sigma1, sigma2), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for a type of bell-shaped curve
bell(x, ymax, a, b, xc) SSbell(x, ymax, a, b, xc)bell(x, ymax, a, b, xc) SSbell(x, ymax, a, b, xc)
x |
input vector |
ymax |
maximum value of y |
a |
parameter controlling the spread (associated with a quadratic term) |
b |
parameter controlling the spread (associated with a cubic term) |
xc |
centering parameter |
This function is described in Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506). One example application is Hammer et al. (2009) (doi:10.2135/cropsci2008.03.0152).
a numeric vector of the same length as x containing parameter estimates for equation specified
bell: vector of the same length as x using a bell-shaped curve
require(ggplot2) set.seed(1234) x <- 1:20 y <- bell(x, 8, -0.0314, 0.000317, 13) + rnorm(length(x), 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSbell(x, ymax, a, b, xc), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) set.seed(1234) x <- 1:20 y <- bell(x, 8, -0.0314, 0.000317, 13) + rnorm(length(x), 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSbell(x, ymax, a, b, xc), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for Beta 5-parameter function
beta5(temp, mu, tb, a, tc, b) SSbeta5(temp, mu, tb, a, tc, b)beta5(temp, mu, tb, a, tc, b) SSbeta5(temp, mu, tb, a, tc, b)
temp |
input vector which is normally ‘temperature’ |
mu |
mu parameter (see equation) |
tb |
base (low) temperature at which no expansion accurs |
a |
parameter describing the initial increasing phase |
tc |
critical (high) temperature at which no expasion occurs |
b |
parameter describing the decreasing phase |
For details see the publication by Yin et al. (1995) “A nonlinear model for crop development as a function of temperature ”. Agricultural and Forest Meteorology 77 (1995) 1-16
The form of the equation is:
.
beta5: vector of the same length as x (temp) using the beta5 function
require(minpack.lm) require(ggplot2) ## Temperature response example data(maizeleafext) ## Fit model fit <- nlsLM(rate ~ SSbeta5(temp, mu, tb, a, tc, b), data = maizeleafext) ## Visualize ndat <- data.frame(temp = 0:45) ndat$rate <- predict(fit, newdata = ndat) ggplot() + geom_point(data = maizeleafext, aes(temp, rate)) + geom_line(data = ndat, aes(x = temp, y = rate))require(minpack.lm) require(ggplot2) ## Temperature response example data(maizeleafext) ## Fit model fit <- nlsLM(rate ~ SSbeta5(temp, mu, tb, a, tc, b), data = maizeleafext) ## Visualize ndat <- data.frame(temp = 0:45) ndat$rate <- predict(fit, newdata = ndat) ggplot() + geom_point(data = maizeleafext, aes(temp, rate)) + geom_line(data = ndat, aes(x = temp, y = rate))
Self starter for Beta Growth function with parameters w.max, lt.e, ldtm, ldtb
bg4rp(time, w.max, lt.e, ldtm, ldtb) SSbg4rp(time, w.max, lt.e, ldtm, ldtb)bg4rp(time, w.max, lt.e, ldtm, ldtb) SSbg4rp(time, w.max, lt.e, ldtm, ldtb)
time |
input vector (x) which is normally ‘time’, the smallest value should be close to zero. |
w.max |
value of weight or mass at its peak |
lt.e |
log of the time at which the maximum weight or mass has been reached. |
ldtm |
log of the difference between time at which the weight or mass reaches its peak and half its peak. |
ldtb |
log of the difference between time at which the weight or mass reaches its peak and when it starts growing |
For details see the publication by Yin et al. (2003) “A Flexible Sigmoid Function of Determinate Growth”.
This is a reparameterization of the beta growth function (4 parameters) with guaranteed constraints, so it is expected to
behave numerically better than SSbgf4.
Reparameterizing the four parameter beta growth
ldtm = log(t.e - t.m)
ldtb = log(t.m - t.b)
t.e = exp(lt.e)
t.m = exp(lt.e) - exp(ldtm)
t.b = (exp(lt.e) - exp(ldtm)) - exp(ldtb)
The form of the equation is:
This is a reparameterized version of the Beta-Growth function in which the parameters are unconstrained, but they are expressed in the log-scale.
bg4rp: vector of the same length as x (time) using the beta growth function with four parameters
require(ggplot2) set.seed(1234) x <- 1:100 y <- bg4rp(x, 20, log(70), log(30), log(20)) + rnorm(100, 0, 1) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSbg4rp(x, w.max, lt.e, ldtm, ldtb), data = dat) ## We are able to recover the original values exp(coef(fit)[2:4]) ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) set.seed(1234) x <- 1:100 y <- bg4rp(x, 20, log(70), log(30), log(20)) + rnorm(100, 0, 1) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSbg4rp(x, w.max, lt.e, ldtm, ldtb), data = dat) ## We are able to recover the original values exp(coef(fit)[2:4]) ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for Beta Growth function with parameters w.max, t.m and t.e
bgf(time, w.max, t.e, t.m) SSbgf(time, w.max, t.e, t.m) bgf2(time, w.max, w.b, t.e, t.m, t.b)bgf(time, w.max, t.e, t.m) SSbgf(time, w.max, t.e, t.m) bgf2(time, w.max, w.b, t.e, t.m, t.b)
time |
input vector (x) which is normally ‘time’, the smallest value should be close to zero. |
w.max |
value of weight or mass at its peak |
t.e |
time at which the weight or mass reaches its peak. |
t.m |
time at which half of the maximum weight or mass has been reached. |
w.b |
weight or biomass at initial time |
t.b |
initial time offset |
For details see the publication by Yin et al. (2003) “A Flexible Sigmoid Function of Determinate Growth”.
The form of the equation is:
.
Given this function weight is expected to decay and reach zero again at .
bgf: vector of the same length as x (time) using the beta growth function
bgf2: a numeric vector of the same length as x (time) containing parameter estimates for equation specified
## See extended example in vignette 'nlraa-AgronJ-paper' x <- seq(0, 17, by = 0.25) y <- bgf(x, 5, 15, 7) plot(x, y)## See extended example in vignette 'nlraa-AgronJ-paper' x <- seq(0, 17, by = 0.25) y <- bgf(x, 5, 15, 7) plot(x, y)
Self starter for Beta Growth function with parameters w.max, t.e, t.m and t.b
bgf4(time, w.max, t.e, t.m, t.b) SSbgf4(time, w.max, t.e, t.m, t.b)bgf4(time, w.max, t.e, t.m, t.b) SSbgf4(time, w.max, t.e, t.m, t.b)
time |
input vector (x) which is normally ‘time’. |
w.max |
value of weight or mass at its peak. |
t.e |
time at which the weight or mass reaches its peak. |
t.m |
time at which half of the maximum weight or mass has been reached. |
t.b |
time at which growth starts. |
For details see the publication by Yin et al. (2003) “A Flexible Sigmoid Function of Determinate Growth”.
This is a difficult function to fit because the linear constrains are not explicitly introduced
in the optimization process. See function SSbgrp for a reparameterized version.
This is equation 11 (pg. 368) in the publication by Yin, but with correction (https://doi.org/10.1093/aob/mcg091) and with ‘w.b’ equal to zero.
a numeric vector of the same length as x (time) containing parameter estimates for equation specified
bgf4: vector of the same length as x (time) using the beta growth function with four parameters
data(sm) ## Let's just pick one crop sm2 <- subset(sm, Crop == "M") fit <- nls(Yield ~ SSbgf4(DOY, w.max, t.e, t.m, t.b), data = sm2) plot(Yield ~ DOY, data = sm2) lines(sm2$DOY,fitted(fit)) ## For this particular problem it could be better to 'fix' t.b and w.b fit0 <- nls(Yield ~ bgf2(DOY, w.max, w.b = 0, t.e, t.m, t.b = 141), data = sm2, start = list(w.max = 16, t.e= 240, t.m = 200)) x <- seq(0, 17, by = 0.25) y <- bgf4(x, 20, 15, 10, 2) plot(x, y)data(sm) ## Let's just pick one crop sm2 <- subset(sm, Crop == "M") fit <- nls(Yield ~ SSbgf4(DOY, w.max, t.e, t.m, t.b), data = sm2) plot(Yield ~ DOY, data = sm2) lines(sm2$DOY,fitted(fit)) ## For this particular problem it could be better to 'fix' t.b and w.b fit0 <- nls(Yield ~ bgf2(DOY, w.max, w.b = 0, t.e, t.m, t.b = 141), data = sm2, start = list(w.max = 16, t.e= 240, t.m = 200)) x <- seq(0, 17, by = 0.25) y <- bgf4(x, 20, 15, 10, 2) plot(x, y)
Self starter for Beta Growth function with parameters w.max, lt.m and ldt
bgrp(time, w.max, lt.e, ldt) SSbgrp(time, w.max, lt.e, ldt)bgrp(time, w.max, lt.e, ldt) SSbgrp(time, w.max, lt.e, ldt)
time |
input vector (x) which is normally ‘time’, the smallest value should be close to zero. |
w.max |
value of weight or mass at its peak |
lt.e |
log of the time at which the maximum weight or mass has been reached. |
ldt |
log of the difference between time at which the weight or mass reaches its peak and half its peak ( |
For details see the publication by Yin et al. (2003) “A Flexible Sigmoid Function of Determinate Growth”.
This is a reparameterization of the beta growth function with guaranteed constraints, so it is expected to
behave numerically better than SSbgf.
The form of the equation is:
.
Given this function weight is expected to decay and reach zero again at . This is a reparameterized version
of the Beta-Growth function in which the parameters are unconstrained, but they are expressed in the log-scale.
bgrp: vector of the same length as x (time) using the beta growth function (reparameterized).
In a few tests it seems that zero values of ‘time’ can cause the error message ‘NA/NaN/Inf in foreign function call (arg 1)’, so it might be better to remove them before running this function.
require(ggplot2) x <- 1:30 y <- bgrp(x, 20, log(25), log(5)) + rnorm(30, 0, 1) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSbgrp(x, w.max, lt.e, ldt), data = dat) ## We are able to recover the original values exp(coef(fit)[2:3]) ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) x <- 1:30 y <- bgrp(x, 20, log(25), log(5)) + rnorm(30, 0, 1) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSbgrp(x, w.max, lt.e, ldt), data = dat) ## We are able to recover the original values exp(coef(fit)[2:3]) ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for a bilinear function with parameters a (intercept), b (first slope), xs (break-point), c (second slope)
blin(x, a, b, xs, c) SSblin(x, a, b, xs, c)blin(x, a, b, xs, c) SSblin(x, a, b, xs, c)
x |
input vector |
a |
the intercept |
b |
the first-phase slope |
xs |
break-point of transition between first-phase linear and second-phase linear |
c |
the second-phase slope |
This is a special case with just two parts but a more general approach is to consider a segmented function with several breakpoints and linear segments. Splines would be even more general. Also this model assumes that there is a break-point that needs to be estimated.
a numeric vector of the same length as x containing parameter estimates for equation specified
blin: vector of the same length as x using the bilinear function
package segmented.
require(ggplot2) set.seed(1234) x <- 1:30 y <- blin(x, 0, 0.75, 15, 1.75) + rnorm(30, 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSblin(x, a, b, xs, c), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## Minimal example ## This is probably about the smallest dataset you ## should use with this function dat2 <- data.frame(x = 1:5, y = c(1.1, 1.9, 3.1, 2, 0.9)) fit2 <- nls(y ~ SSblin(x, a, b, xs, c), data = dat2) ggplot(data = dat2, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit2)))require(ggplot2) set.seed(1234) x <- 1:30 y <- blin(x, 0, 0.75, 15, 1.75) + rnorm(30, 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSblin(x, a, b, xs, c), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## Minimal example ## This is probably about the smallest dataset you ## should use with this function dat2 <- data.frame(x = 1:5, y = c(1.1, 1.9, 3.1, 2, 0.9)) fit2 <- nls(y ~ SSblin(x, a, b, xs, c), data = dat2) ggplot(data = dat2, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit2)))
Self starter for cardinal temperature response function
card3(x, tb, to, tm) SScard3(x, tb, to, tm)card3(x, tb, to, tm) SScard3(x, tb, to, tm)
x |
input vector (x) which is normally ‘temperature’. |
tb |
base temperature |
to |
optimum temperature |
tm |
maximum temperature |
This function is described in Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506)
card3: vector of the same length as x using a card3 function
Caio dos Santos and Fernando Miguez
## A temperature response function require(ggplot2) set.seed(1234) x <- 1:50 y <- card3(x, 13, 25, 36) + rnorm(length(x), sd = 0.05) dat1 <- data.frame(x = x, y = y) fit1 <- nls(y ~ SScard3(x, tb, to, tm), data = dat1) ggplot(data = dat1, aes(x, y)) + geom_point() + geom_line(aes(y = fitted(fit1)))## A temperature response function require(ggplot2) set.seed(1234) x <- 1:50 y <- card3(x, 13, 25, 36) + rnorm(length(x), sd = 0.05) dat1 <- data.frame(x = x, y = y) fit1 <- nls(y ~ SScard3(x, tb, to, tm), data = dat1) ggplot(data = dat1, aes(x, y)) + geom_point() + geom_line(aes(y = fitted(fit1)))
Self starter for declining logistic function with parameters asym, a2, xmid and scal
dlf(time, asym, a2, xmid, scal) SSdlf(time, asym, a2, xmid, scal)dlf(time, asym, a2, xmid, scal) SSdlf(time, asym, a2, xmid, scal)
time |
input vector (x) which is normally ‘time’, the smalles value should be close to zero. |
asym |
value of weight or mass at its peak (maximum) |
a2 |
value of weight or mass at its trough (minimum) |
xmid |
time at which half of the maximum weight or mass has bean reached. |
scal |
scale parameter which controls the spread also interpreted in terms of time to go from xmid to approx. 0.63 asym |
Response function:
.
asym: upper asymptote
xmid: time when y is midway between w and a
scal: controls the spread
a2: lower asymptote
The four parameter logistic SSfpl is essentially equivalent to this function,
but here the interpretation of the parameters is different and this is the form used in
Oddi et. al. (2019) (see vignette).
a numeric vector of the same length as x (time) containing parameter estimates for equation specified
dlf: vector of the same length as x (time) using the declining logistic function
## Extended example in the vignette 'nlraa-Oddi-LFMC' x <- seq(0, 17, by = 0.25) y <- dlf(x, 2, 10, 8, 1) plot(x, y, type = "l")## Extended example in the vignette 'nlraa-Oddi-LFMC' x <- seq(0, 17, by = 0.25) y <- dlf(x, 2, 10, 8, 1) plot(x, y, type = "l")
Self starter for a simple exponential function
expf(x, a, c) SSexpf(x, a, c)expf(x, a, c) SSexpf(x, a, c)
x |
input vector (x) |
a |
represents the value at x = 0 |
c |
represents the exponential rate |
This is the exponential function
For more details see: Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506).
a numeric vector of the same length as x containing parameter estimates for equation specified
expf: vector of the same length as x using the profd function
require(ggplot2) set.seed(1234) x <- 1:15 y <- expf(x, 10, -0.3) + rnorm(15, 0, 0.2) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSexpf(x, a, c), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) set.seed(1234) x <- 1:15 y <- expf(x, 10, -0.3) + rnorm(15, 0, 0.2) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSexpf(x, a, c), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for an exponential-plateau function
expfp(x, a, c, xs) SSexpfp(x, a, c, xs)expfp(x, a, c, xs) SSexpfp(x, a, c, xs)
x |
input vector (x) |
a |
represents the value at x = 0 |
c |
represents the exponential rate |
xs |
represents the breakpoint at which the plateau starts |
This is the exponential-plateua function, where ‘xs’ is the break-point
For more details see: Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506).
a numeric vector of the same length as x containing parameter estimates for equation specified
expfp: vector of the same length as x using the expfp function
require(ggplot2) set.seed(12345) x <- 1:30 y <- expfp(x, 10, 0.1, 15) + rnorm(30, 0, 1.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSexpfp(x, a, c, xs), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) set.seed(12345) x <- 1:30 y <- expfp(x, 10, 0.1, 15) + rnorm(30, 0, 1.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSexpfp(x, a, c, xs), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for an exponential-linear growth equation
explin(t, cm, rm, tb) SSexplin(t, cm, rm, tb)explin(t, cm, rm, tb) SSexplin(t, cm, rm, tb)
t |
input vector (time) |
cm |
parameter related to the maximum growth during the linear phase |
rm |
parameter related to the maximum growth during the exponential phase |
tb |
time at which switch happens |
J. GOUDRIAAN, J. L. MONTEITH, A Mathematical Function for Crop Growth Based on Light Interception and Leaf Area Expansion, Annals of Botany, Volume 66, Issue 6, December 1990, Pages 695–701, doi:10.1093/oxfordjournals.aob.a088084
The equation is:
This function is described in Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506).
a numeric vector of the same length as x containing parameter estimates for equation specified
explin: vector of the same length as x using a explin function
require(ggplot2) set.seed(12345) x <- seq(1,100, by = 5) y <- explin(x, 20, 0.14, 30) + rnorm(length(x), 0, 5) y <- abs(y) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSexplin(x, cm, rm, tb), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) set.seed(12345) x <- seq(1,100, by = 5) y <- explin(x, 20, 0.14, 30) + rnorm(length(x), 0, 5) y <- abs(y) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSexplin(x, cm, rm, tb), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for a harmonic regression
harm1(x, b0, b1, cos1, sin1) SSharm1(x, b0, b1, cos1, sin1)harm1(x, b0, b1, cos1, sin1) SSharm1(x, b0, b1, cos1, sin1)
x |
input vector |
b0 |
intercept of the harmonic regression |
b1 |
slope of the harmonic regression |
cos1 |
coefficient associated with the cosine of the harmonic regression |
sin1 |
coefficient associated with the sine of the harmonic regression |
Harmonic regression is actually a type of linear regression. Just adding it for convenience.
a numeric vector of the same length as x containing parameter estimates for equation specified
harm1: vector of the same length as x using a harmonic regression
require(ggplot2) set.seed(1234) x <- seq(0, 3, length.out = 100) y <- harm1(x, 0, 0, 0.05, 0) + rnorm(length(x), 0, 0.002) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSharm1(x, b0, b1, cos1, sin1), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) set.seed(1234) x <- seq(0, 3, length.out = 100) y <- harm1(x, 0, 0, 0.05, 0) + rnorm(length(x), 0, 0.002) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSharm1(x, b0, b1, cos1, sin1), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for Hill function with parameters Ka, n and a
hill1(x, Ka) SShill1(x, Ka) hill2(x, Ka, n) SShill2(x, Ka, n) hill3(x, Ka, n, a) SShill3(x, Ka, n, a)hill1(x, Ka) SShill1(x, Ka) hill2(x, Ka, n) SShill2(x, Ka, n) hill3(x, Ka, n, a) SShill3(x, Ka, n, a)
x |
input vector (x). Concentration of substrate in the original Hill model. |
Ka |
parameter representing the concentration at which half of maximum y is attained |
n |
parameter which controls the curvature |
a |
parameter which controls the maximum value of the response (asymptote) |
For details see https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)
The form of the equations are:
hill1:
.
hill2:
.
hill3:
.
hill1: vector of the same length as x (time) using the Hill 1 function
hill2: vector of the same length as x (time) using the Hill 2 function
hill3: vector of the same length as x (time) using the Hill 3 function
Zero values are not allowed.
require(ggplot2) ## Example for hill1 set.seed(1234) x <- 1:20 y <- hill1(x, 10) + rnorm(20, sd = 0.03) dat1 <- data.frame(x = x, y = y) fit1 <- nls(y ~ SShill1(x, Ka), data = dat1) ## Example for hill2 y <- hill2(x, 10, 1.5) + rnorm(20, sd = 0.03) dat2 <- data.frame(x = x, y = y) fit2 <- nls(y ~ SShill2(x, Ka, n), data = dat2) ## Example for hill3 y <- hill3(x, 10, 1.5, 5) + rnorm(20, sd = 0.03) dat3 <- data.frame(x = x, y = y) fit3 <- nls(y ~ SShill3(x, Ka, n, a), data = dat3) ggplot(data = dat3, aes(x, y)) + geom_point() + geom_line(aes(y = fitted(fit3)))require(ggplot2) ## Example for hill1 set.seed(1234) x <- 1:20 y <- hill1(x, 10) + rnorm(20, sd = 0.03) dat1 <- data.frame(x = x, y = y) fit1 <- nls(y ~ SShill1(x, Ka), data = dat1) ## Example for hill2 y <- hill2(x, 10, 1.5) + rnorm(20, sd = 0.03) dat2 <- data.frame(x = x, y = y) fit2 <- nls(y ~ SShill2(x, Ka, n), data = dat2) ## Example for hill3 y <- hill3(x, 10, 1.5, 5) + rnorm(20, sd = 0.03) dat3 <- data.frame(x = x, y = y) fit3 <- nls(y ~ SShill3(x, Ka, n, a), data = dat3) ggplot(data = dat3, aes(x, y)) + geom_point() + geom_line(aes(y = fitted(fit3)))
Self starter for linear-plateau function with parameters a (intercept), b (slope), xs (break-point)
linp(x, a, b, xs) SSlinp(x, a, b, xs)linp(x, a, b, xs) SSlinp(x, a, b, xs)
x |
input vector |
a |
the intercept |
b |
the slope |
xs |
break-point of transition between linear and plateau |
This function is linear when and flat () when .
a numeric vector of the same length as x containing parameter estimates for equation specified
linp: vector of the same length as x using the linear-plateau function
package segmented.
require(ggplot2) set.seed(123) x <- 1:30 y <- linp(x, 0, 1, 20) + rnorm(30, 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSlinp(x, a, b, xs), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## Confidence intervals confint(fit)require(ggplot2) set.seed(123) x <- 1:30 y <- linp(x, 0, 1, 20) + rnorm(30, 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSlinp(x, a, b, xs), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## Confidence intervals confint(fit)
Self starter for a five-parameter logistic function.
logis5(x, asym1, asym2, xmid, iscal, theta) SSlogis5(x, asym1, asym2, xmid, iscal, theta)logis5(x, asym1, asym2, xmid, iscal, theta) SSlogis5(x, asym1, asym2, xmid, iscal, theta)
x |
input vector (x) |
asym1 |
asymptotic value for low values of x |
asym2 |
asymptotic value for high values of x |
xmid |
value of x at which y = (asym1 + asym2)/2 (only when theta = 1) |
iscal |
steepness of transition from asym1 to asym2 (inverse of the scale) |
theta |
asymmetry parameter, if it is equal to 1, this is the four parameter logistic |
The equation for this function is:
This is known as the Richards' function or the log-logistic and it is described in Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506).
a numeric vector of the same length as x (time) containing parameter estimates for equation specified
logis5: vector of the same length as x (time) using the 5-parameter logistic
require(ggplot2) set.seed(1234) x <- seq(0, 2000, 100) y <- logis5(x, 35, 10, 800, 5, 2) + rnorm(length(x), 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSlogis5(x, asym1, asym2, xmid, iscal, theta), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) x <- seq(0, 2000) y <- logis5(x, 30, 10, 800, 5, 2) plot(x, y)require(ggplot2) set.seed(1234) x <- seq(0, 2000, 100) y <- logis5(x, 35, 10, 800, 5, 2) + rnorm(length(x), 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSlogis5(x, asym1, asym2, xmid, iscal, theta), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) x <- seq(0, 2000) y <- logis5(x, 30, 10, 800, 5, 2) plot(x, y)
Self starter for modified Hyperbola with parameters: asymp, xmin and k
moh(x, asym, xmin, k) SSmoh(x, asym, xmin, k)moh(x, asym, xmin, k) SSmoh(x, asym, xmin, k)
x |
input vector (x) which is normally a controlling variable such as nitrogen |
asym |
asymptotic value when x tends to infinity |
xmin |
value of x for which y equals zero |
k |
curvature parameter |
This function is described in Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506). See Table S3 (Eq. 3.8)
a numeric vector of the same length as x containing parameter estimates for equation specified
moh: vector of the same length as x (time) using the modified hyperbola
require(ggplot2) set.seed(1234) x <- seq(3, 30) y <- moh(x, 35, 3, 0.83) + rnorm(length(x), 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSmoh(x, asym, xmin, k), data = dat) ## Visualize observed and simulated ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## Testing predict function prd <- predict_nls(fit, interval = "confidence") datA <- cbind(dat, prd) ## Plotting ggplot(data = datA, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) x <- seq(0, 20) y <- moh(x, 30, 3, 0.9) plot(x, y)require(ggplot2) set.seed(1234) x <- seq(3, 30) y <- moh(x, 35, 3, 0.83) + rnorm(length(x), 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSmoh(x, asym, xmin, k), data = dat) ## Visualize observed and simulated ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## Testing predict function prd <- predict_nls(fit, interval = "confidence") datA <- cbind(dat, prd) ## Plotting ggplot(data = datA, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) x <- seq(0, 20) y <- moh(x, 30, 3, 0.9) plot(x, y)
Self starter for Non-rectangular Hyperbola with parameters: asymptote, quantum efficiency, curvature and dark respiration
nrh(x, asym, phi, theta, rd) SSnrh(x, asym, phi, theta, rd)nrh(x, asym, phi, theta, rd) SSnrh(x, asym, phi, theta, rd)
x |
input vector (x) which is normally light intensity (PPFD, Photosynthetic Photon Flux Density). |
asym |
asymptotic value for photosynthesis |
phi |
quantum efficiency (mol CO2 per mol of photons) or initial slope of the light response |
theta |
curvature parameter for smooth transition between limitations |
rd |
dark respiration or value of CO2 uptake at zero light levels |
This function is described in Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506).
a numeric vector of the same length as x (time) containing parameter estimates for equation specified
nrh: vector of the same length as x (time) using the non-rectangular hyperbola
require(ggplot2) set.seed(1234) x <- seq(0, 2000, 100) y <- nrh(x, 35, 0.04, 0.83, 2) + rnorm(length(x), 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSnrh(x, asym, phi, theta, rd), data = dat) ## Visualize observed and simulated ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## Testing predict function prd <- predict_nls(fit, interval = "confidence") datA <- cbind(dat, prd) ## Plotting ggplot(data = datA, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) x <- seq(0, 2000) y <- nrh(x, 30, 0.04, 0.85, 2) plot(x, y)require(ggplot2) set.seed(1234) x <- seq(0, 2000, 100) y <- nrh(x, 35, 0.04, 0.83, 2) + rnorm(length(x), 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSnrh(x, asym, phi, theta, rd), data = dat) ## Visualize observed and simulated ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## Testing predict function prd <- predict_nls(fit, interval = "confidence") datA <- cbind(dat, prd) ## Plotting ggplot(data = datA, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) x <- seq(0, 2000) y <- nrh(x, 30, 0.04, 0.85, 2) plot(x, y)
Self starter for an plateau-exponential function
pexpf(x, a, xs, c) SSpexpf(x, a, xs, c)pexpf(x, a, xs, c) SSpexpf(x, a, xs, c)
x |
input vector (x) |
a |
represents the value for the plateau |
xs |
represents the breakpoint at which the plateau ends |
c |
represents the exponential rate |
The equation is: .
a numeric vector of the same length as x containing parameter estimates for equation specified
pexpf: vector of the same length as x using the pexpf function
require(ggplot2) set.seed(1234) x <- 1:30 y <- pexpf(x, 20, 15, -0.2) + rnorm(30, 0, 1) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSpexpf(x, a, xs, c), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) set.seed(1234) x <- 1:30 y <- pexpf(x, 20, 15, -0.2) + rnorm(30, 0, 1) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSpexpf(x, a, xs, c), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for plateau-linear function with parameters a (plateau), xs (break-point), b (slope)
plin(x, a, xs, b) SSplin(x, a, xs, b)plin(x, a, xs, b) SSplin(x, a, xs, b)
x |
input vector |
a |
the initial plateau |
xs |
break-point of transition between plateau and linear |
b |
the slope |
Initial plateau with a second linear phase. When and when .
a numeric vector of the same length as x containing parameter estimates for equation specified
plin: vector of the same length as x using the plateau-linear function
require(ggplot2) set.seed(123) x <- 1:30 y <- plin(x, 10, 20, 1) + rnorm(30, 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSplin(x, a, xs, b), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## Confidence intervals confint(fit)require(ggplot2) set.seed(123) x <- 1:30 y <- plin(x, 10, 20, 1) + rnorm(30, 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSplin(x, a, xs, b), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## Confidence intervals confint(fit)
Self starter for plateau-quadratic function with parameters a (plateau), xs (break-point), b (slope), c (quadratic)
pquad(x, a, xs, b, c) SSpquad(x, a, xs, b, c)pquad(x, a, xs, b, c) SSpquad(x, a, xs, b, c)
x |
input vector |
a |
the plateau value |
xs |
break-point of transition between plateau and quadratic |
b |
the slope (linear term) |
c |
quadratic term |
Reference for nonlinear regression Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506).
a numeric vector of the same length as x containing parameter estimates for equation specified
pquad: vector of the same length as x using the plateau-quadratic function
require(ggplot2) set.seed(12345) x <- 1:40 y <- pquad(x, 5, 20, 1.7, -0.04) + rnorm(40, 0, 0.6) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSpquad(x, a, xs, b, c), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) confint(fit)require(ggplot2) set.seed(12345) x <- 1:40 y <- pquad(x, 5, 20, 1.7, -0.04) + rnorm(40, 0, 0.6) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSpquad(x, a, xs, b, c), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) confint(fit)
Self starter for plateau-quadratic function with (three) parameters a (intercept), b (slope), c (quadratic)
pquad3(x, a, b, c) SSpquad3(x, a, b, c)pquad3(x, a, b, c) SSpquad3(x, a, b, c)
x |
input vector |
a |
the intercept |
b |
the slope |
c |
quadratic term |
Reference for nonlinear regression Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506).
a numeric vector of the same length as x containing parameter estimates for equation specified
quadp: vector of the same length as x using the quadratic-plateau function
require(ggplot2) require(minpack.lm) set.seed(123) x <- 1:30 y <- pquad3(x, 20.5, 0.36, -0.012) + rnorm(30, 0, 0.3) dat <- data.frame(x = x, y = y) fit <- nlsLM(y ~ SSpquad3(x, a, b, c), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) require(minpack.lm) set.seed(123) x <- 1:30 y <- pquad3(x, 20.5, 0.36, -0.012) + rnorm(30, 0, 0.3) dat <- data.frame(x = x, y = y) fit <- nlsLM(y ~ SSpquad3(x, a, b, c), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for a decay of a variable within a canopy (e.g. nitrogen concentration).
profd(x, a, b, c, d) SSprofd(x, a, b, c, d)profd(x, a, b, c, d) SSprofd(x, a, b, c, d)
x |
input vector (x) |
a |
represents the maximum value |
b |
represents the minimum value |
c |
represents the rate of decay |
d |
represents an empirical coefficient which provides flexibility |
This function is described in Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506) and originally in Johnson et al. (2010) Annals of Botany 106: 735–749, 2010. (doi:10.1093/aob/mcq183).
a numeric vector of the same length as x containing parameter estimates for equation specified
profd: vector of the same length as x using the profd function
require(ggplot2) set.seed(1234) x <- 1:10 y <- profd(x, 0.3, 0.05, 0.5, 4) + rnorm(10, 0, 0.01) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSprofd(x, a, b, c, d), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## profiling ## It does not work at a lower alphamax level ## Use this if you want to look at all four parameters ## par(mfrow=c(2,2)) plot(profile(fit, alphamax = 0.016)) ## Reset graphical parameter as appropriate: par(mfrow=c(1,1)) ## parameter 'd' is not well constrained confint(fit, level = 0.9)require(ggplot2) set.seed(1234) x <- 1:10 y <- profd(x, 0.3, 0.05, 0.5, 4) + rnorm(10, 0, 0.01) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSprofd(x, a, b, c, d), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## profiling ## It does not work at a lower alphamax level ## Use this if you want to look at all four parameters ## par(mfrow=c(2,2)) plot(profile(fit, alphamax = 0.016)) ## Reset graphical parameter as appropriate: par(mfrow=c(1,1)) ## parameter 'd' is not well constrained confint(fit, level = 0.9)
Self starter for quadratic plateau function with parameters a (intercept), b (slope), c (quadratic), xs (break-point)
quadp(x, a, b, c, xs) SSquadp(x, a, b, c, xs)quadp(x, a, b, c, xs) SSquadp(x, a, b, c, xs)
x |
input vector |
a |
the intercept |
b |
the slope |
c |
quadratic term |
xs |
break point of transition between quadratic and plateau |
Reference for nonlinear regression Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506).
a numeric vector of the same length as x containing parameter estimates for equation specified
quadp: vector of the same length as x using the quadratic-plateau function
require(ggplot2) set.seed(123) x <- 1:30 y <- quadp(x, 5, 1.7, -0.04, 20) + rnorm(30, 0, 0.6) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSquadp(x, a, b, c, xs), data = dat, algorithm = "port") ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) set.seed(123) x <- 1:30 y <- quadp(x, 5, 1.7, -0.04, 20) + rnorm(30, 0, 0.6) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSquadp(x, a, b, c, xs), data = dat, algorithm = "port") ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for quadratic plateau function with (three) parameters a (intercept), b (slope), c (quadratic)
quadp3(x, a, b, c) SSquadp3(x, a, b, c)quadp3(x, a, b, c) SSquadp3(x, a, b, c)
x |
input vector |
a |
the intercept |
b |
the slope |
c |
quadratic term |
The equation is, for a response (y) and a predictor (x):
where the break-point (xs) is -0.5*b/c
and the asymptote is (a + (-b^2)/(4 * c))
In this model the parameter ‘xs’ is not directly estimated. If this is required, the model ‘SSquadp3xs’ should be used instead.
a numeric vector of the same length as x containing parameter estimates for equation specified
quadp: vector of the same length as x using the quadratic-plateau function
require(ggplot2) set.seed(123) x <- 1:30 y <- quadp3(x, 5, 1.7, -0.04) + rnorm(30, 0, 0.6) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSquadp3(x, a, b, c), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) set.seed(123) x <- 1:30 y <- quadp3(x, 5, 1.7, -0.04) + rnorm(30, 0, 0.6) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSquadp3(x, a, b, c), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for quadratic plateau function with (three) parameters a (intercept), b (slope), xs (break-point)
quadp3xs(x, a, b, xs) SSquadp3xs(x, a, b, xs)quadp3xs(x, a, b, xs) SSquadp3xs(x, a, b, xs)
x |
input vector |
a |
the intercept |
b |
the slope |
xs |
break-point |
The equation is, for a response (y) and a predictor (x):
where the quadratic term is (c) is -0.5*b/xs
and the asymptote is (a + (b^2)/(4 * c)).
This model does not estimate the quadratic parameter ‘c’ directly. If this is required, the model ‘SSquadp3’ should be used instead.
a numeric vector of the same length as x containing parameter estimates for equation specified
quadp3xs: vector of the same length as x using the quadratic-plateau function
require(ggplot2) set.seed(123) x <- 1:30 y <- quadp3xs(x, 5, 1.7, 20) + rnorm(30, 0, 0.6) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSquadp3xs(x, a, b, xs), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) set.seed(123) x <- 1:30 y <- quadp3xs(x, 5, 1.7, 20) + rnorm(30, 0, 0.6) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSquadp3xs(x, a, b, xs), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for quadratic plateau function with (four) parameters a (intercept), b (slope), xs (break-point), dxs (difference between breakpoint and second break-point)
quadpq(x, a, b, xs, ldxs) SSquadpq(x, a, b, xs, ldxs)quadpq(x, a, b, xs, ldxs) SSquadpq(x, a, b, xs, ldxs)
x |
input vector |
a |
the intercept |
b |
the slope |
xs |
break-point |
ldxs |
log of the difference between break-point and second break-point |
The equation is, for a response (y) and a predictor (x):
where the quadratic term is (c) is -0.5*b/xs
and the asymptote is (a + (b^2)/(4 * c)).
This model does not estimate the quadratic parameter ‘c’ directly. If this is required, the model ‘SSquadp3’ should be used instead.
a numeric vector of the same length as x containing parameter estimates for equation specified
quadpq: vector of the same length as x using the quadratic-plateau-quadratic function
require(ggplot2) set.seed(123) x <- 0:25 y <- quadpq(x, 1, 0.5, 10, 1.5) + rnorm(length(x), 0, 0.3) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSquadpq(x, a, b, xs, dxs), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) + geom_vline(aes(xintercept = coef(fit)[3]), linetype = 2) + geom_vline(aes(xintercept = coef(fit)[3] + exp(coef(fit)[4])), linetype = 3)require(ggplot2) set.seed(123) x <- 0:25 y <- quadpq(x, 1, 0.5, 10, 1.5) + rnorm(length(x), 0, 0.3) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSquadpq(x, a, b, xs, dxs), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) + geom_vline(aes(xintercept = coef(fit)[3]), linetype = 2) + geom_vline(aes(xintercept = coef(fit)[3] + exp(coef(fit)[4])), linetype = 3)
Self starter for a rational curve
ratio(x, a, b, c, d) SSratio(x, a, b, c, d)ratio(x, a, b, c, d) SSratio(x, a, b, c, d)
x |
input vector |
a |
parameter related to the maximum value of the response (numerator) |
b |
power exponent for numerator |
c |
parameter related to the maximum value of the response (denominator) |
d |
power exponent for denominator |
The equation is:
This function is described in Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506). One example application is in Bril et al. (1994) https://edepot.wur.nl/333930 - pages 19 and 21. The parameters are difficult to interpret, but the function is very flexible. I have not tested this, but it might be beneficial to re-scale x and y to the (0,1) range if this function is hard to fit. https://en.wikipedia.org/wiki/Rational_function.
a numeric vector of the same length as x containing parameter estimates for equation specified
ratio: vector of the same length as x using a rational function
require(ggplot2) require(minpack.lm) set.seed(1234) x <- 1:100 y <- ratio(x, 1, 0.5, 1, 1.5) + rnorm(length(x), 0, 0.025) dat <- data.frame(x = x, y = y) fit <- nlsLM(y ~ SSratio(x, a, b, c, d), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) require(minpack.lm) set.seed(1234) x <- 1:100 y <- ratio(x, 1, 0.5, 1, 1.5) + rnorm(length(x), 0, 0.025) dat <- data.frame(x = x, y = y) fit <- nlsLM(y ~ SSratio(x, a, b, c, d), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for Ricker function with parameters a and b
ricker(time, a, b) SSricker(time, a, b)ricker(time, a, b) SSricker(time, a, b)
time |
input vector (x) which is normally ‘time’, the smallest value should be close to zero. |
a |
which is related to the initial growth slope |
b |
which is related to the slowing down or decline |
This function is described in Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506) and originally in Ricker, W. E. (1954) Stock and Recruitment Journal of the Fisheries Research Board of Canada, 11(5): 559–623. (doi:10.1139/f54-039).
The equation is: .
a numeric vector of the same length as x (time) containing parameter estimates for equation specified
ricker: vector of the same length as x (time) using the ricker function
require(ggplot2) set.seed(123) x <- 1:30 y <- 30 * x * exp(-0.3 * x) + rnorm(30, 0, 0.25) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSricker(x, a, b), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) set.seed(123) x <- 1:30 y <- 30 * x * exp(-0.3 * x) + rnorm(30, 0, 0.25) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SSricker(x, a, b), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for smooth cardinal temperature response function
scard3(x, tb, to, tm, curve = 2) SSscard3(x, tb, to, tm, curve = 2)scard3(x, tb, to, tm, curve = 2) SSscard3(x, tb, to, tm, curve = 2)
x |
input vector (x) which is normally ‘temperature’. |
tb |
base temperature |
to |
optimum temperature |
tm |
maximum temperature |
curve |
curvature (default is 2) |
An example application can be found in (doi:10.1016/j.envsoft.2014.04.009)
This function is described in Archontoulis and Miguez (2015) - (doi:10.2134/agronj2012.0506) - Equation 5.1 in Table 1.
scard3: vector of the same length as x using a scard3 function
Caio dos Santos and Fernando Miguez
## A temperature response function require(ggplot2) set.seed(1234) x <- 1:50 y <- scard3(x, 13, 25, 36) + rnorm(length(x), sd = 0.05) dat1 <- data.frame(x = x, y = y) fit1 <- nls(y ~ SSscard3(x, tb, to, tm), data = dat1) ggplot(data = dat1, aes(x, y)) + geom_point() + geom_line(aes(y = fitted(fit1)))## A temperature response function require(ggplot2) set.seed(1234) x <- 1:50 y <- scard3(x, 13, 25, 36) + rnorm(length(x), sd = 0.05) dat1 <- data.frame(x = x, y = y) fit1 <- nls(y ~ SSscard3(x, tb, to, tm), data = dat1) ggplot(data = dat1, aes(x, y)) + geom_point() + geom_line(aes(y = fitted(fit1)))
Self starter for temperature response function
sharp(temp, r_tref, e, el, tl, eh, th, tref = 25) SSsharp(temp, r_tref, e, el, tl, eh, th, tref = 25)sharp(temp, r_tref, e, el, tl, eh, th, tref = 25) SSsharp(temp, r_tref, e, el, tl, eh, th, tref = 25)
temp |
input vector (x) which is normally ‘temperature’. |
r_tref |
rate at the standardised temperature, tref |
e |
activation energy (eV) |
el |
low temperature de-activation energy (eV) |
tl |
temperature at which the enzyme is half active and half suppressed due to low temperatures |
eh |
high temperature de-activation energy (eV) |
th |
temperature at which enzyme is half active and half suppressed dut to high temperatures |
tref |
standardisation temperature in degrees centigrade. Temperature at which rates are not inactivated by either high or low temperatures. Typically, 25 degrees. |
For details see Schoolfield, R. M., Sharpe, P. J. & Magnuson, C. E. Non-linear regression of biological temperature-dependent rate models based on absolute reaction-rate theory. Journal of Theoretical Biology 88, 719-731 (1981)
sharp: vector of the same length as x using a sharp function
I do not recommend using this function.
require(ggplot2) require(minpack.lm) temp <- 0:45 rate <- sharp(temp, 1, 0.03, 1.44, 28, 19, 44) + rnorm(length(temp), 0, 0.05) dat <- data.frame(temp = temp, rate = rate) ## Fit model fit <- nlsLM(rate ~ SSsharp(temp, r_tref, e, el, tl, eh, th, tref = 20), data = dat) ## Visualize ggplot(data = dat, aes(temp, rate)) + geom_point() + geom_line(aes(y = fitted(fit)))require(ggplot2) require(minpack.lm) temp <- 0:45 rate <- sharp(temp, 1, 0.03, 1.44, 28, 19, 44) + rnorm(length(temp), 0, 0.05) dat <- data.frame(temp = temp, rate = rate) ## Fit model fit <- nlsLM(rate ~ SSsharp(temp, r_tref, e, el, tl, eh, th, tref = 20), data = dat) ## Visualize ggplot(data = dat, aes(temp, rate)) + geom_point() + geom_line(aes(y = fitted(fit)))
Self starter for Collatz temperature response function
temp3(x, t.m, t.l, t.h) SStemp3(x, t.m, t.l, t.h)temp3(x, t.m, t.l, t.h) SStemp3(x, t.m, t.l, t.h)
x |
input vector (x) which is normally ‘temperature’. |
t.m |
medium temperature |
t.l |
low temperature |
t.h |
high temperature |
Collatz GJ , Ribas-Carbo M Berry JA (1992) Coupled Photosynthesis-Stomatal Conductance Model for Leaves of C4 Plants. Functional Plant Biology 19, 519-538. https://doi.org/10.1071/PP9920519
temp3: vector of the same length as x using a temp function
## A temperature response function require(ggplot2) set.seed(1234) x <- 1:50 y <- temp3(x, 25, 13, 36) + rnorm(length(x), sd = 0.05) dat1 <- data.frame(x = x, y = y) fit1 <- nls(y ~ SStemp3(x, t.m, t.l, t.h), data = dat1) ggplot(data = dat1, aes(x, y)) + geom_point() + geom_line(aes(y = fitted(fit1)))## A temperature response function require(ggplot2) set.seed(1234) x <- 1:50 y <- temp3(x, 25, 13, 36) + rnorm(length(x), sd = 0.05) dat1 <- data.frame(x = x, y = y) fit1 <- nls(y ~ SStemp3(x, t.m, t.l, t.h), data = dat1) ggplot(data = dat1, aes(x, y)) + geom_point() + geom_line(aes(y = fitted(fit1)))
Self starter for a tri-linear function with parameters a (intercept), b (first slope), xs1 (first break-point), c (second slope), xs2 (second break-point) and d (third slope)
trlin(x, a, b, xs1, c, xs2, d) SStrlin(x, a, b, xs1, c, xs2, d)trlin(x, a, b, xs1, c, xs2, d) SStrlin(x, a, b, xs1, c, xs2, d)
x |
input vector |
a |
the intercept |
b |
the first-phase slope |
xs1 |
first break-point of transition between first-phase linear and second-phase linear |
c |
the second-phase slope |
xs2 |
second break-point of transition between second-phase linear and third-phase linear |
d |
the third-phase slope |
This is a special case with just three parts (and two break points) but a more general approach is to consider a segmented function with several breakpoints and linear segments. Splines would be even more general. Also this model assumes that there are two break-points that needs to be estimated. The guess for the initial values splits the dataset in half, so it this will work when one break-point is in the first half and the second is in the second half.
a numeric vector of the same length as x containing parameter estimates for equation specified
trlin: vector of the same length as x using the tri-linear function
package segmented.
require(ggplot2) set.seed(1234) x <- 1:30 y <- trlin(x, 0.5, 2, 10, 0.1, 20, 1.75) + rnorm(30, 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SStrlin(x, a, b, xs1, c, xs2, d), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## Minimal example ## This is probably about the smallest dataset you ## should use with this function dat2 <- data.frame(x = 1:8, y = c(1.1, 1.9, 3.1, 2.5, 1.4, 0.9, 2.2, 2.9)) fit2 <- nls(y ~ SStrlin(x, a, b, xs1, c, xs2, d), data = dat2) ## expangin for plotting ndat <- data.frame(x = seq(1, 8, by = 0.1)) ndat$prd <- predict(fit2, newdata = ndat) ggplot() + geom_point(data = dat2, aes(x = x, y = y)) + geom_line(data = ndat, aes(x = x, y = prd))require(ggplot2) set.seed(1234) x <- 1:30 y <- trlin(x, 0.5, 2, 10, 0.1, 20, 1.75) + rnorm(30, 0, 0.5) dat <- data.frame(x = x, y = y) fit <- nls(y ~ SStrlin(x, a, b, xs1, c, xs2, d), data = dat) ## plot ggplot(data = dat, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit))) ## Minimal example ## This is probably about the smallest dataset you ## should use with this function dat2 <- data.frame(x = 1:8, y = c(1.1, 1.9, 3.1, 2.5, 1.4, 0.9, 2.2, 2.9)) fit2 <- nls(y ~ SStrlin(x, a, b, xs1, c, xs2, d), data = dat2) ## expangin for plotting ndat <- data.frame(x = seq(1, 8, by = 0.1)) ndat$prd <- predict(fit2, newdata = ndat) ggplot() + geom_point(data = dat2, aes(x = x, y = y)) + geom_line(data = ndat, aes(x = x, y = prd))
Utility function to summarize the output from ‘simulate’ functions in this package
summary_simulate( object, probs = c(0.025, 0.975), robust = FALSE, data, by, ... )summary_simulate( object, probs = c(0.025, 0.975), robust = FALSE, data, by, ... )
object |
nobs x nsim matrix where nobs are the number of observations in the dataset and nsim are the number of simulations |
probs |
the percentiles to be computed by the quantile function |
robust |
If FALSE (the default) the mean is used as the measure of central tendency and the standard deviation as the measure of variability. If TRUE, the median and the median absolute deviation (MAD) are applied instead. |
data |
the original data.frame used to fit the model. A data.frame will be returned instead of a matrix in this case. |
by |
opionally aggregate the results by some factor in the data.frame. It will be coarced to a formula. This should either be a character or a formula (starting with ‘~’). The aggregation follows the ‘robust’ argument above. |
... |
additional arguments to be passed. (none used at the moment) |
By default it returns a matrix unless the ‘data’ argument is present and then it will return a data.frame
data(barley, package = "nlraa") fit <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley) sim <- simulate_nls(fit, nsim = 100) sims <- summary_simulate(sim) ## If we want to combine the data.frame simd <- summary_simulate(sim, data = barley) ## If we also want to aggregate by nitrogen rate simda <- summary_simulate(sim, data = barley, by = "NF") ## The robust option uses the median instead simdar <- summary_simulate(sim, data = barley, by = "NF", robust = TRUE)data(barley, package = "nlraa") fit <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley) sim <- simulate_nls(fit, nsim = 100) sims <- summary_simulate(sim) ## If we want to combine the data.frame simd <- summary_simulate(sim, data = barley) ## If we also want to aggregate by nitrogen rate simda <- summary_simulate(sim, data = barley, by = "NF") ## The robust option uses the median instead simdar <- summary_simulate(sim, data = barley, by = "NF", robust = TRUE)
Simulated data based on obseved data presented in Sinclair (1986) - Fig. 1A
swpgswpg
A data frame with 20 rows and 3 columns
Fraction of Transpirable Soil Water (0-1)
Relative Leaf Growth scaled from 0 to 1
Sinclair, T.R. Water and Nitrogen Limitations in Soybean Grain Production I. Model Development. Field Crops Research. 125-141.
Simluated data (much cleaner than original) based on the above publication
Extracts the variance covariance matrix (residuals, random or all)
var_cov( object, type = c("residual", "random", "all", "conditional", "marginal"), aug = FALSE, sparse = FALSE, data = NULL )var_cov( object, type = c("residual", "random", "all", "conditional", "marginal"), aug = FALSE, sparse = FALSE, data = NULL )
object |
|
type |
“residual” for the variance-covariance for the residuals, “random” for the variance-covariance of the random effects or “all” for the sum of both. |
aug |
whether to augment the matrix of the random effects to the dimensions of the data |
sparse |
whether to return a sparse matrix (default is FALSE) |
data |
optional passing of ‘data’, probably needed when using this function inside other functions. |
Variance Covariance matrix for (non)linear mixed models
It returns a matrix or a sparse matrix Matrix.
See Chapter 5 of Pinheiro and Bates. This returns potentially a very large
matrix of N x N, where N is the number of rows in the data.frame.
The function getVarCov only works well for
lme objects.
The equivalence is more or less:
getVarCov type = “random.effects” equivalent to var_cov type = “random”.
getVarCov type = “conditional” equivalent to var_cov type = “residual”.
getVarCov type = “marginal” equivalent to var_cov type = “all”.
The difference is that getVarCov has an argument that specifies the individual
for which the matrix is being retrieved and var_cov returns the full matrix only.
require(graphics) require(nlme) data(ChickWeight) ## First a linear model flm <- lm(weight ~ Time, data = ChickWeight) vlm <- var_cov(flm) ## First model with no modeling of the Variance-Covariance fit0 <- gls(weight ~ Time, data = ChickWeight) v0 <- var_cov(fit0) ## Only modeling the diagonal (weights) fit1 <- gls(weight ~ Time, data = ChickWeight, weights = varPower()) v1 <- var_cov(fit1) ## Only the correlation structure is defined and there are no groups fit2 <- gls(weight ~ Time, data = ChickWeight, correlation = corAR1()) v2 <- var_cov(fit2) ## The correlation structure is defined and there are groups present fit3 <- gls(weight ~ Time, data = ChickWeight, correlation = corCAR1(form = ~ Time | Chick)) v3 <- var_cov(fit3) ## There are both weights and correlations fit4 <- gls(weight ~ Time, data = ChickWeight, weights = varPower(), correlation = corCAR1(form = ~ Time | Chick)) v4 <- var_cov(fit4) ## Tip: you can visualize these matrices using image(log(v4[,ncol(v4):1]))require(graphics) require(nlme) data(ChickWeight) ## First a linear model flm <- lm(weight ~ Time, data = ChickWeight) vlm <- var_cov(flm) ## First model with no modeling of the Variance-Covariance fit0 <- gls(weight ~ Time, data = ChickWeight) v0 <- var_cov(fit0) ## Only modeling the diagonal (weights) fit1 <- gls(weight ~ Time, data = ChickWeight, weights = varPower()) v1 <- var_cov(fit1) ## Only the correlation structure is defined and there are no groups fit2 <- gls(weight ~ Time, data = ChickWeight, correlation = corAR1()) v2 <- var_cov(fit2) ## The correlation structure is defined and there are groups present fit3 <- gls(weight ~ Time, data = ChickWeight, correlation = corCAR1(form = ~ Time | Chick)) v3 <- var_cov(fit3) ## There are both weights and correlations fit4 <- gls(weight ~ Time, data = ChickWeight, weights = varPower(), correlation = corCAR1(form = ~ Time | Chick)) v4 <- var_cov(fit4) ## Tip: you can visualize these matrices using image(log(v4[,ncol(v4):1]))