cbvn {BLA} | R Documentation |
Fitting boundary line using censored bivariate normal model
Description
This function fits a response model to the upper limits of a scatter plot of
of x
and y
to determine the most efficient response of y
as a function of x
(given a measurement error of y
) based on a
censored distribution (Milne et al., 2016). The location of censor in the data
cloud is determined based on the maximum likelihood approach. This is done using
optimization procedure and hence requires some starting guess parameters for the
proposed model. It then compares the results with an uncensored normal bivariate
distribution to access the appropriateness of the censored model.
Usage
cbvn(data,model="lp", equation=NULL, start, sigh, UpLo="U", optim.method="BFGS",
Hessian=FALSE, plot=TRUE, line_smooth=1000, lwd=2, l_col="red",...)
Arguments
data |
A dataframe with two numeric columns, independent ( |
model |
Selects the functional form of the boundary line. It includes
|
equation |
A custom model function writen in the form of an R function. Applies
only when argument |
start |
A numeric vector of initial starting values for optimization in fitting the boundary model. Its length and arrangement depend on the suggested model:
|
sigh |
Standard deviation of the measurement error. |
UpLo |
Selects the type of boundary. |
optim.method |
Describes the method used to optimize the model as in the
|
Hessian |
If |
plot |
If |
line_smooth |
Parameter that describes the smoothness of the boundary line. (default is 1000). The higher the value, the smoother the line. |
lwd |
Determines the thickness of the boundary line on the plot (default is 1). |
l_col |
Selects the color of the boundary line. |
... |
Additional graphical parameters as in the |
Details
Some inbuilt models are available for the cbvn()
function. The suggest model
forms are as follows:
Linear model (
"blm"
)y=\beta_1 + \beta_2x
where
\beta_1
is the intercept and\beta_2
is the slope.Linear plateau model (
"lp"
)y= {\rm min}(\beta_1+\beta_2x, \beta_0)
where
\beta_1
is the intercept ,\beta_2
is the slope and\beta_0
is the maximum response.The logistic (
"logistic"
) and inverse logistic ("inv-logistic"
) modelsy= \frac{\beta_0}{1+e^{\beta_2(\beta_1-x)}}
y= \beta_0 - \frac{\beta_0}{1+e^{\beta_2(\beta_1-x)}}
where
\beta_1
is a scaling parameter ,\beta_2
is a shape parameter and\beta_0
is the maximum response.Logistic model (
"logisticND"
) (Nelder (1961))y= \frac{\beta_0}{1+(\beta_1 \times e^{-\beta_2x})}
where
\beta_1
is a scaling parameter,\beta_2
is a shape parameter and\beta_0
is the maximum response.Double logistic model (
"double-logistic"
)y= \frac{\beta_{0,1}}{1+e^{\beta_2(\beta_1-x)}} - \frac{\beta_{0,2}}{1+e^{\beta_4(\beta_3-x)}}
where
\beta_1
is a scaling parameter one,\beta_2
is shape parameter one,\beta_{0,1}
and\beta_{0,2}
are the maximum response ,\beta_3
is a scaling parameter two and\beta_4
is a shape parameter two.Quadratic model (
"qd"
)y=\beta_1 + \beta_2x + \beta_3x^2
where
\beta_1
is a constant,\beta_2
is a linear coefficient and\beta_3
is the quadratic coefficient.Trapezium model (
"trapezium"
)y={\rm min}(\beta_1+\beta_2x, \beta_0, \beta_3 + \beta_4x)
where
\beta_1
is the intercept one,\beta_2
is the slope one,\beta_0
is the maximum response,\beta_3
is the intercept two and\beta_3
is the slope two.Mitscherlich model (
"mit"
)y= \beta_0 - \beta_1*\beta_2^x
where
\beta_1
is the intercept,\beta_2
is a shape parameter and\beta_0
is the maximum response.Schmidt model (
"schmidt"
)y= \beta_0 + \beta_1(x-\beta_2)^2
where
\beta_1
is a scaling parameter,\beta_2
is a shape parameter (x-value at maximum response ) and\beta_0
is the maximum response .
The function cbvn()
utilities the optimization procedure of the
optim()
function to determine the model parameters. There is a tendency
for optimization algorithms to settle at a local optimum. To remove the risk of
settling for local optimum parameters, it is advised that the function is run using
several starting values and the results with the smallest likelihood (or AIC)
can be taken as a representation of the global optimum.
The common errors encountered due to poor start values
function cannot be evaluated at initial parameters
initial value in 'vmmin' is not finite
Value
A list of length 5 consisting of the fitted model, equation form, parameters of the boundary line, AIC (for boundary line model and a null model) and a hessian matrix. Additionally, a graphical representation of the boundary line on the scatter plot is produced.
Author(s)
Chawezi Miti chawezi.miti@nottingham.ac.uk
Richard Murray Lark murray.lark@nottingham.ac.uk
References
Nelder, J.A. 1961. The fitting of a generalization of the logistic curve. Biometrics 17: 89–110.
Lark, R. M., & Milne, A. E. (2016). Boundary line analysis of the effect of water filled pore space on nitrous oxide emission from cores of arable soil. European Journal of Soil Science, 67 , 148-159.
Lark, R. M., Gillingham, V., Langton, D., & Marchant, B. P. (2020). Boundary line models for soil nutrient concentrations and wheat yield in national-scale datasets. European Journal of Soil Science, 71 , 334-351.
Milne, A. E., Ferguson, R. B., & Lark, R. M. (2006). Estimating a boundary line model for a biological response by maximum likelihood.Annals of Applied Biology, 149, 223–234.
Phillips, B.F. & Campbell, N.A. 1968. A new method of fitting the von Bertelanffy growth curve using data on the whelk. Dicathais, Growth 32: 317–329.
Schmidt, U., Thöni, H., & Kaupenjohann, M. (2000). Using a boundary line approach to analyze N2O flux data from agricultural soils. Nutrient Cycling in Agroecosystems, 57, 119-129.
Examples
x<-evapotranspiration$`ET(mm)`
y<-evapotranspiration$`yield(t/ha)`
data<-data.frame(x,y)
start<-c(0.5,0.02,289.6,2.4,83.7,1.07,0.287)
cbvn(data, start=start, model = "blm", sigh=0.51,
xlab=expression("ET/ mm ha"^-1),
ylab=expression("Yield/ ton ha"^-1),
pch=16, col="grey", line_smooth = 100)