install.packages("BLA")
Introducing BLA: A Comprehensive R Package for Boundary Line Analysis
We are excited to announce the release of a new R
package, BLA
, on CRAN. This package is designed to provide a robust set of tools for fitting boundary line models to data sets. Whether you’re dealing with crop yield predictions or analyzing ecological data, BLA
offers an array of methods to help you derive meaningful insights from your data. The boundary line approach is suitable for analyzing biological response datasets from surveys i.e. cases in which multiple potential limiting factors to biological response occur but are not controlled as done in experimental conditions. In the case, one can determine the largest expected biological response (also called the boundary value) for a given value of a factor interest. Additionally, the most limiting factor to the biological response can be modeled from the observed potential limiting factors.
Overview of BLA
The BLA
package builds on the foundational work by Webb (1972) and incorporates a variety of heuristic and statistical boundary line fitting methods developed since then. We give an overview of these methods in a recent paper (Miti et al., 2024a). These includes the (1) Binning Methods that involves selecting of boundary points by segmenting the independent variable into bins of equal width from which the boundary point is chosen based on a given criteria e.g. 95\(^{\rm th}\) percentile of the response variable, (2) Bolides Algorithm, an advanced technique for selecting boundary points introduced by Schnug et al. (1995), (3) Quantile Regression, which models the boundary line as a given quantile (Cade and Noon, 2003), and (4) Censored Bivariate Normal Model, a statistical method developed by Milne et al. (2006) that is based on the maximum likelihood approach assuming the data takes a censored bivariate normal distribution.
Additionally, BLA
provides data exploration tools to identify outliers and checking for evidence of boundary presence in your data before diving into the boundary line analysis. Post boundary line fit functions are also available for a comprehensive analysis. The Table below shows a summary of the available functions in BLA
.
Exploratory | Model fitting | Post-hoc |
---|---|---|
summastat() |
blbin() |
predictBL() |
exp_bounary() |
bolides() |
limfactor() |
startValues() |
blqr() |
|
ble_profile() |
cbvn() |
Installation and loading
Getting started with BLA
is straightforward. You can install the current version of the package directly from CRAN with the following command:
The package can now be loaded using the library()
function. Additionally, the aplpack
package is loaded which contains a function for the identification of outliers.
library(BLA)
library(aplpack)
Use case for BLA
BLA
is particularly useful in scenarios where you need to model the relationship between a response variable,y
, and a potentially limiting factor, x
. For example, in agricultural studies, you might want to analyze the relationship between crop yield and soil P concentration data from non-experimental settings (Lark et al., 2020). BLA
helps you determine the most efficient biological response (yield) for a given value of soil P concentration. Below is an example of how you can use the censored bivariate normal model (the cbvn()
function) to carryout boundary line analysis. This process includes initial data exploration, fitting the boundary line model and prediction of boundary values.
Data exploration
There are three important exploratory steps required in boundary line analysis. These are (1) checking the distribution of the potential limiting and response factors to fulfill the assumption of normality or normality with a censor, (2) detection of outliers and (3) assessing the evidence that the dataset exhibits a boundary which limits the response variable (Miti et al., 2024b).
1. Checking variable distribution indices and plots
The function summastat()
can be used for this purpose. We Start with soil P concentration which we will denote as x
.
<- soil$P
x
summastat(x)
Mean Median Quartile.1 Quartile.3 Variance SD Skewness
[1,] 25.9647 22 16 32 207.0066 14.38772 1.840844
Octile skewness Kurtosis No. outliers
[1,] 0.3571429 5.765138 43
From the plots and the summary statistics, the x
variable can not be assumed to be from a normal distribution and hence requires transformation to fulfill the assumption of normality. You can do this by taking the natural log of x
.
summastat(log(x))
Mean Median Quartile.1 Quartile.3 Variance SD Skewness
[1,] 3.126361 3.091042 2.772589 3.465736 0.2556936 0.5056615 0.1297406
Octile skewness Kurtosis No. outliers
[1,] 0.08395839 -0.05372586 0
Normality can now be assumed after transformation. Next, you check the yield variable which we denote as y
.
<- soil$yield
y
summastat(y)
Mean Median Quartile.1 Quartile.3 Variance SD Skewness
[1,] 9.254813 9.36468 8.203703 10.39477 3.456026 1.859039 -0.4819805
Octile skewness Kurtosis No. outliers
[1,] -0.05793291 1.292635 7
From the plots and the summary statistics, the variable y
can be assumed to be from a normal distribution.
2. Outlier identification and removal
This is done using the bagplot()
function from the aplpack
package. Its input is a matrix
and hence we assign the x
and y
variables to a matrix data_ur
.
<- length(soil$P)
nobs <- matrix(NA, nobs, 2) # Create a matrix
data_ur 1] <- log(soil$P)
data_ur[,2] <- soil$yield
data_ur[,
<- bagplot(data_ur, show.whiskers = F,create.plot = TRUE) # Identify outliers bag
All the points outside the bag are considered outliers and can be removed from the analysis.
<- rbind(bag$pxy.bag, bag$pxy.outer) # Exclude bivariate outliers data
3. Evidence of boundary presence
This is done using the function expl_boundary()
to check if the are grounds for us to fit a boundary line model to the dataset. there should be evidence of a limiting relationship of x
and y
at the upper bounds of the data scatter.
<- data[,1]
x <- data[,2]
y
expl_boundary(x,y,10,1000)
Note: This function might take longer to execute when running a large number of simulations.
Index Section value
1 sd Left 1.045711
2 sd Right 1.115379
3 Mean sd Left 1.131335
4 Mean sd Right 1.203422
5 p_value Left 0.032000
6 p_value Right 0.025000
The p-values in the left and right sections are less than 0.05, indicating evidence of an upper boundary existence as opposed to a random scatter of points. This justifies the fitting of a boundary line model to the dataset.
Boundary line fitting
Based on the structure of the data at the upper edge, a trapezium model can be fitted. Take note that any other model of your choice that is biologically plausible can be fitted.
It arguments include (1) a dataframe,data
, containing the x
and y
variables, (2) a vector of initial start
values for the optimization that includes the parameters of the boundary model and for the distribution (i.e. mean(x)
, mean(y)
, sd(x)
, sd(y)
and cor(x,y)
), and (3) the measurement error value for the response variable, sigh
. The rest of the inputs are related to the plot features as in the function plot()
from base R
.
The start values for the preferred model can be determined using the startValues()
function. Set the argument model
to the desired model e.g. model="trapezium"
, run the function and click on the plot of y
against x
, the points that make up the structure of the model at the upper edges of the data scatter.
startValues("trapezium")
Using the obtained values for the model, create a vector of start values.
<-data.frame(x,y)
data<-c(4,3,13.6, 35,-5,3,9,0.50,1.9,0.05) # initial start values for optimization
start
<- cbvn(data, start = start, model = "trapezium", sigh=0.7,
model xlab = expression("Soil P/ log(mg l"^-1*")"),
ylab = expression("Yield/ t ha"^-1),
pch = 16, plot = TRUE, col = "grey40", cex = 0.8)
model
$Model
[1] "trapezium"
$Equation
[1] y = min (β₁ + β₂x, β₀, β₃ + β₄x)
$Parameters
Estimate Standard error
β₁ 4.29795522 1.035391840
β₂ 3.23397375 0.460850826
β₀ 13.15187257 0.206656034
β₃ 33.17267393 1.693458789
β₄ -5.22857503 0.393758941
mux 3.12597270 0.006451086
muy 9.30482617 0.022879293
sdx 0.50053107 0.004561474
sdy 1.60754420 0.018448808
rcorr 0.04150832 0.014806076
$AIC
mvn 32429.55
BL 32391.07
The output gives the name of model fitted which is a trapezium in the case, its equation form, its parameters (the first 5 rows) and corresponding standard errors, and the AIC values for the boundary model, BL, and a corresponding multivariate normal model, mvn. Since the AIC for the BL is smaller than mvn, the boundary line model is appropriate. We interpret the first rising section of the boundary model as showing the P-limited response of yield, and the soil P concentration at which this intersects the horizontal section as a critical value (it is consistent with advisory index values). We think that the falling section of the line is likely to reflect the accumulation of P in the soil in areas with consistent large limitations on crop yield (e.g. soil compaction in headlands).
Predicting boundary yield
Using the fitted model, you can predict the boundary yield for each data point based on soil P values:
<- log(soil$P)
P_values is.na(x)] <- mean(x, na.rm = TRUE) # replace missing values with mean pH
P_values[<- predictBL(model, P_values)
predicted_yield
head(predicted_yield) # predicted yield for the first six farms
[1] 11.74445 11.40372 11.74445 11.40372 11.40372 12.33408
Determining critical soil P concentration value
The critical soil P concentration value is the point beyond which yield increase is not expected. You can calculate it from the model parameters as follows:
<- model$Parameters[1]
intercept <- model$Parameters[2]
slope <- model$Parameters[3]
plateau
<- (plateau - intercept) / slope
critical_P
print(exp(critical_P))# results are in mg/l
[1] 15.45268
Closing remarks
The BLA
package is a powerful tool for researchers and analysts working with data sets where analyzing boundary response is crucial. With methods ranging from the heuristic binning to advanced statistical algorithms like censored bivariate normal model, BLA
provides all the necessary tools for a comprehensive boundary line analysis. For further details and advanced usage, visit the website https://chawezimiti.github.io/BLA/.
We hope you find BLA
useful and look forward to your feedback and contributions (https://github.com/chawezimiti/BLA/issues).