Title: | Biological Geometries |
---|---|
Description: | Is used to simulate and fit biological geometries. 'biogeom' incorporates several novel universal parametric equations that can generate the profiles of bird eggs, flowers, linear and lanceolate leaves, seeds, starfish, and tree-rings (Gielis (2003) <doi:10.3732/ajb.90.3.333>; Shi et al. (2020) <doi:10.3390/sym12040645>), three growth-rate curves representing the ontogenetic growth trajectories of animals and plants against time, and the axially symmetrical and integral forms of all these functions (Shi et al. (2017) <doi:10.1016/j.ecolmodel.2017.01.012>; Shi et al. (2021) <doi:10.3390/sym13081524>). The optimization method proposed by Nelder and Mead (1965) <doi:10.1093/comjnl/7.4.308> was used to estimate model parameters. 'biogeom' includes several real data sets of the boundary coordinates of natural shapes, including avian eggs, fruit, lanceolate and ovate leaves, tree rings, seeds, and sea stars,and can be potentially applied to other natural shapes. 'biogeom' can quantify the conspecific or interspecific similarity of natural outlines, and provides information with important ecological and evolutionary implications for the growth and form of living organisms. Please see Shi et al. (2022) <doi:10.1111/nyas.14862> for details. |
Authors: | Peijian Shi [aut, cre], Johan Gielis [aut], Brady K. Quinn [aut] |
Maintainer: | Peijian Shi <[email protected]> |
License: | GPL (>= 2) |
Version: | 1.4.3 |
Built: | 2024-11-26 06:25:34 UTC |
Source: | CRAN |
adjdata
adjusts the data points in counterclockwise order based on the shortest distance method.
adjdata(x, y, ub.np = 2000, times = 1.2, len.pro = 1/20, index.sp = 1)
adjdata(x, y, ub.np = 2000, times = 1.2, len.pro = 1/20, index.sp = 1)
x |
the |
y |
the |
ub.np |
the upper bound of the number of points eventually retained on the polygon's boundary. |
times |
the number of times |
len.pro |
the proportion of the distance between any two points to the maximum distance between the points on the polygon's boundary, which is used to determine whether the second point needs to be deleted. |
index.sp |
the index of the starting point of a group of indices that regularly
divide the number of points on the polygon's boundary into |
When ub.np
> length(x)
, length(x)
points
on the polygon's boundary are retained.
The quantile
function in package stats
is used to carry out the regular division of data points. From the starting point, the second point is the
one that has the shortest distance from the former. When the distance between the two points
is larger than len.pro
multiplied by the maximum distance between points on the polygon's boundary,
the second point is deleted from the coordinates.
Then, the third point that has the shortest distance from the first point is defined as the second point.
If the distance between the first point and the second point is no more than len.pro
multiplied
by the maximum distance, the first and second points are recorded in a new matrix for the coordinates of the polygon,
and the second point is defined as the first point in the old matrix for the coordinates of the polygon.
The shortest distance method is then used to look for a third point that meets the requirement.
x |
the |
y |
the |
The initial boundary data of a polygon can be obtained by running the M-file based on Matlab (version >= 2009a) developed by Shi et al. (2018) and Su et al. (2019) for a .bmp black and white image of the polygon. See references below.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Ratkowsky, D.A., Li, Y., Zhang, L., Lin, S., Gielis, J. (2018) General leaf-area geometric formula exists for plants - Evidence from the simplified Gielis equation. Forests 9, 714. doi:10.3390/f9110714
Su, J., Niklas, K.J., Huang, W., Yu, X., Yang, Y., Shi, P. (2019) Lamina shape does not correlate with lamina surface area: An analysis based on the simplified Gielis equation. Global Ecology and Conservation 19, e00666. doi:10.1016/j.gecco.2019.e00666
data(eggs) uni.C1 <- sort( unique(eggs$Code) ) ind1 <- 2 Data1 <- eggs[eggs$Code==uni.C1[ind1], ] x0 <- Data1$x y0 <- Data1$y Res1 <- adjdata(x0, y0, ub.np=2000, times=1.2, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( x1, y1, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")), pch=1, col=1 ) Res2 <- adjdata(x0, y0, ub.np=40, times=1, len.pro=1/2, index.sp=20) x2 <- Res2$x y2 <- Res2$y Res3 <- adjdata(x0, y0, ub.np=100, times=1, len.pro=1/2, index.sp=100) x3 <- Res3$x y3 <- Res3$y dev.new() plot( x2, y2, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")), pch=1, col=4 ) points( x3, y3, col=2) data(starfish) uni.C2 <- sort( unique(starfish$Code) ) ind2 <- 2 Data2 <- starfish[starfish$Code==uni.C2[ind2], ] x4 <- Data2$x y4 <- Data2$y dev.new() plot( x4, y4, asp=1, type="l", cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) Res4 <- adjdata(x4, y4, ub.np=500, times=1.2, len.pro=1/20) x5 <- Res4$x y5 <- Res4$y dev.new() plot( x5, y5, asp=1, type="l", cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) graphics.off()
data(eggs) uni.C1 <- sort( unique(eggs$Code) ) ind1 <- 2 Data1 <- eggs[eggs$Code==uni.C1[ind1], ] x0 <- Data1$x y0 <- Data1$y Res1 <- adjdata(x0, y0, ub.np=2000, times=1.2, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( x1, y1, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")), pch=1, col=1 ) Res2 <- adjdata(x0, y0, ub.np=40, times=1, len.pro=1/2, index.sp=20) x2 <- Res2$x y2 <- Res2$y Res3 <- adjdata(x0, y0, ub.np=100, times=1, len.pro=1/2, index.sp=100) x3 <- Res3$x y3 <- Res3$y dev.new() plot( x2, y2, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")), pch=1, col=4 ) points( x3, y3, col=2) data(starfish) uni.C2 <- sort( unique(starfish$Code) ) ind2 <- 2 Data2 <- starfish[starfish$Code==uni.C2[ind2], ] x4 <- Data2$x y4 <- Data2$y dev.new() plot( x4, y4, asp=1, type="l", cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) Res4 <- adjdata(x4, y4, ub.np=500, times=1.2, len.pro=1/20) x5 <- Res4$x y5 <- Res4$y dev.new() plot( x5, y5, asp=1, type="l", cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) graphics.off()
areaGE
is used to calculate the area of the polygon
generated by the Gielis curve within .
areaGE(expr, P, m = 1, simpver = NULL, nval = 1, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
areaGE(expr, P, m = 1, simpver = NULL, nval = 1, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
expr |
the original (or twin) Gielis equation or one of its simplified versions. |
P |
the parameters of the original (or twin) Gielis equation or one of its simplified versions. |
m |
the given |
simpver |
an optional argument to use the simplified version of the original (or twin) Gielis equation. |
nval |
the specified value for |
subdivisions |
please see the arguments for the |
rel.tol |
please see the arguments for the |
abs.tol |
please see the arguments for the |
stop.on.error |
please see the arguments for the |
keep.xy |
please see the arguments for the |
aux |
please see the arguments for the |
The arguments of P
, m
, simpver
, and nval
should correspond
to expr
(i.e., GE
or TGE
). Please note the differences in the simplified
version number and the number of parameters between GE
and TGE
.
The area of the polygon within generated by the original (or twin) Gielis equation
or one of its simplified versions.
simpver
in GE
is different from that in TGE
.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Gielis, J. (2003) A generic geometric transformation that unifies a wide range of natural
and abstract shapes. American Journal of Botany 90, 333338. doi:10.3732/ajb.90.3.333
Li, Y., Quinn, B.K., Gielis, J., Li, Y., Shi, P. (2022) Evidence that supertriangles exist in nature from the vertical projections of Koelreuteria paniculata fruit. Symmetry 14, 23. doi:10.3390/sym14010023
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Ratkowsky, D.A., Gielis, J. (2020) The generalized Gielis geometric equation and its application. Symmetry 12, 645. doi:10.3390/sym12040645
Shi, P., Xu, Q., Sandhu, H.S., Gielis, J., Ding, Y., Li, H., Dong, X. (2015) Comparison
of dwarf bamboos (Indocalamus sp.) leaf parameters to determine relationship
between spatial density of plants and total leaf area per plant. Ecology and Evolution 5,
45784589. doi:10.1002/ece3.1728
Para1 <- c(1.7170, 5.2258, 7.9802) areaGE(GE, P = Para1, m=5, simpver=1) Para2 <- c(2.1066, 3.5449, 0.4619, 10.5697) areaGE(TGE, P = Para2, m=5, simpver=1)
Para1 <- c(1.7170, 5.2258, 7.9802) areaGE(GE, P = Para1, m=5, simpver=1) Para2 <- c(2.1066, 3.5449, 0.4619, 10.5697) areaGE(TGE, P = Para2, m=5, simpver=1)
areaovate
is used to calculate the area of an ovate polygon made
from combing two symmetrical curves generated by a performance equation (e.g., MLRFE
).
areaovate(expr, P, simpver = NULL, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
areaovate(expr, P, simpver = NULL, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
expr |
a performance equation or one of its simplified versions. |
P |
the parameters of the performance equation or one of its simplified versions. |
simpver |
an optional argument to use the simplfied version of the performance equation. |
subdivisions |
please see the arguments for the |
rel.tol |
please see the arguments for the |
abs.tol |
please see the arguments for the |
stop.on.error |
please see the arguments for the |
keep.xy |
please see the arguments for the |
aux |
please see the arguments for the |
The performance equations denote MbetaE
, MBriereE
,
MLRFE
, and their simplified versions.
The arguments of P
and simpver
should correspond
to expr
(i.e., MbetaE
or MBriereE
or MLRFE
).
The area of two symmetrical curves along the -axis
generated by a performance equation or one of its simplified versions.
Here, the user can define other performance equations, but new equations or
their simplified versions should include the lower and upper thresholds in
the -axis corresponding to
, whose indices should
be the same as those in
MbetaE
or MBriereE
or MLRFE
.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Jin, J., Quinn, B.K., Shi, P. (2022) The modified Brière equation and its applications. Plants 11, 1769. doi:10.3390/plants11131769
Shi, P., Fan, M., Ratkowsky, D.A., Huang, J., Wu, H., Chen, L., Fang, S.,
Zhang, C. (2017) Comparison of two ontogenetic growth equations for animals and plants.
Ecological Modelling 349, 110. doi:10.1016/j.ecolmodel.2017.01.012
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Yu, K., Niklas, K.J., Schrader, J., Song, Y., Zhu, R., Li, Y., Wei, H., Ratkowsky, D.A. (2021) A general model for describing the ovate leaf shape. Symmetry 13, 1524. doi:10.3390/sym13081524
curveovate
, fitovate
, MbetaE
, MBriereE
,
MLRFE
, MPerformanceE
, sigmoid
Par1 <- c(1.8175, 2.7795, 7.1557, 1.6030) areaovate(MbetaE, P = Par1, simpver = 1) Par2 <- c(0.0550, 0.3192, 7.1965, 0.5226) areaovate(MBriereE, P = Par2, simpver = 1) Par3 <- c(1.8168, 2.7967, 7.2623, 0.9662) areaovate(MLRFE, P = Par3, simpver = 1) Par4 <- c(2.4, 0.96, 0.64, 7.75, 1.76, 3.68) areaovate(MPerformanceE, P = Par4, simpver = 1)
Par1 <- c(1.8175, 2.7795, 7.1557, 1.6030) areaovate(MbetaE, P = Par1, simpver = 1) Par2 <- c(0.0550, 0.3192, 7.1965, 0.5226) areaovate(MBriereE, P = Par2, simpver = 1) Par3 <- c(1.8168, 2.7967, 7.2623, 0.9662) areaovate(MLRFE, P = Par3, simpver = 1) Par4 <- c(2.4, 0.96, 0.64, 7.75, 1.76, 3.68) areaovate(MPerformanceE, P = Par4, simpver = 1)
The data consist of the boundary data of six leaves of P. incarnata sampled at Nanjing Forestry University campus in early December 2016.
data(bambooleaves)
data(bambooleaves)
In the data set, there are three columns of variables: Code
, x
, and y
.
Code
saves the codes of individual leaves;
x
saves the coordinates of the leaf boundary in the Cartesian coordinate system (cm);
and
y
saves the coordinates of the leaf boundary in the Cartesian coordinate system (cm).
Lin, S., Shao, L., Hui, C., Song, Y., Reddy, G.V.P., Gielis, J., Li, F., Ding, Y., Wei, Q., Shi, P. (2018) Why does not the leaf weight-area allometry of bamboos follow the 3/2-power law? Frontiers in Plant Science 9, 583. doi:10.3389/fpls.2018.00583
Shi, P., Ratkowsky, D.A., Li, Y., Zhang, L., Lin, S., Gielis, J. (2018) General leaf-area geometric formula exists for plants - Evidence from the simplified Gielis equation. Forests 9, 714. doi:10.3390/f9110714
data(bambooleaves) uni.C <- sort( unique(bambooleaves$Code) ) ind <- 1 Data <- bambooleaves[bambooleaves$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y dev.new() plot( x0, y0, asp=1, type="l", cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) length(x0) Res1 <- adjdata(x0, y0, ub.np=600, len.pro=1/20) dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic("x")), ylab=expression(italic("y")) ) graphics.off()
data(bambooleaves) uni.C <- sort( unique(bambooleaves$Code) ) ind <- 1 Data <- bambooleaves[bambooleaves$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y dev.new() plot( x0, y0, asp=1, type="l", cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) length(x0) Res1 <- adjdata(x0, y0, ub.np=600, len.pro=1/20) dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic("x")), ylab=expression(italic("y")) ) graphics.off()
bilat
is used to measure the extent of bilateral (a)symmetry and other measures for a polygon (e.g., a leaf).
bilat(x, y, strip.num = 200, peri.np = NULL, n.loop = 60, auto.search = TRUE, animation.fig = TRUE, time.interval = 0.001, unit = "cm", main = NULL, diff.fig = TRUE, angle = NULL, ratiox = 0.02, ratioy = 0.08, fd.opt = TRUE, frac.fig = TRUE, denomi.range = seq(8, 30, by = 1))
bilat(x, y, strip.num = 200, peri.np = NULL, n.loop = 60, auto.search = TRUE, animation.fig = TRUE, time.interval = 0.001, unit = "cm", main = NULL, diff.fig = TRUE, angle = NULL, ratiox = 0.02, ratioy = 0.08, fd.opt = TRUE, frac.fig = TRUE, denomi.range = seq(8, 30, by = 1))
x |
the |
y |
the |
strip.num |
the number of equidistant strips intersecting with the polygon that are horizontally placed. See Shi et al. (2018, 2020) for details. |
peri.np |
the number of data points on the boundary retained for calculating the perimeter of the polygon. |
n.loop |
the number of data points to randomly sample for calculating the mean perimeter of the polygon. |
auto.search |
an optional argument to automatically search the maximum distance between two points on the polygon's boundary. |
animation.fig |
the option of showing the data points on the polygon's boundary in an animation. |
time.interval |
the time interval at which to suspend execution, in seconds. |
unit |
the units of the |
main |
the main title of the figure. |
diff.fig |
an optional argument to draw the differences in areas between the intersections of the strips
with the upper part of the polygon and the intersections of the strips with the lower part of the polygon.
The polygon is divided into the upper and lower parts by the |
angle |
the angle between the major axis (i.e., the leaf length axis) and the |
fd.opt |
An optional argument to use the box-counting method to calculate the fractal dimension of the polygon's boundary on a log-log scale. |
ratiox |
the |
ratioy |
the |
frac.fig |
an optional argument to draw the results of the linear fitting using the box-counting method to calculate the fractal dimension of the polygon's boundary on a log-log scale. |
denomi.range |
the number of equidistant segments of the maximum range
between the range of the |
The data of x
and y
should be the coordinates adjusted using the adjdata
function.
If peri.np = NULL
, the number of length(x)
is used to calculate the perimeter of the polygon;
if peri.np
is a positive integer, the number of data points retained on the polygon's boundary
is equal to peri.np
and random sampling for retaining peri.np
data points is carried out
n.loop
times for calculating the mean perimeter of the polygon. That is to say, the final output for
the perimeter is the mean of the n.loop
perimeters (i.e., replicates). If the user wants to get a consistent result
for the mean perimeter, the set.seed
function can be used. In addition, if length(x) < peri.np
,
peri.np
then becomes length(x)
rather than the specified value in Arguments.
If the polygon apparently has a major axis (e.g., the leaf length axis for an ovate leaf), auto.search
is appropriate. If the major axis of the polygon is not the straight line through two points on the polygon's
boundary having the maximum distance, the user can define the major axis using the locator
function in graphic by clicking two points on or near the polygon's boundary. The location of the first click
should be northeast of the location of the second click. This means that the angle between the straight line
through the locations of the two clicks and the -axis should range from 0 to
/2. The locations of the
clicks can be on the boundary or be approximate to the boundary. The function will automatically find the nearest
data point on the boundary to the location of each click.
When
angle = NULL
, the observed polygon will be shown at its initial angle in the scanned image;
when angle
is a numerical value (e.g., ) defined by the user, it indicates that the major axis
is rotated
counterclockwise from the
-axis.
x |
the |
y |
the |
phi |
the angle between the length axis (i.e., the major axis) of the polygon and the |
n1 |
the number of data points on the upper boundary of the polygon. |
n2 |
the number of data points on the lower boundary of the polygon. |
n |
the number of data points on the whole polygon's boundary. |
total.poly |
an object of class "ppp" representing a point pattern dataset in the two-dimensional plane, representing the polygon's boundary. |
upper.poly |
an object of class "ppp" representing a point pattern dataset
in the two-dimensional plane, representing the upper boundary of the polygon along the |
lower.poly |
an object of class "ppp" representing a point pattern dataset
in the two-dimensional plane, representing the lower boundary of the polygon along the |
D |
the differences in areas between the upper and lower boundaries of the polygon. |
par.upper.area |
the area of the upper boundary of the polygon along the |
par.lower.area |
the area of the lower boundary of the polygon along the |
SI |
the standardized index for bilateral (a)symmetry for the polygon. |
AR |
the ratio of the areas of the upper to the lower parts of the polygon. |
scan.length |
the length of the polygon. The default is the maximum distance between two points on the polygon's boundary. |
scan.width |
the maximum width of the polygon. |
scan.area |
the area of the polygon. |
scan.perimeter |
the perimeter of the polygon
based on all data points or a mean of |
x.width |
distance from the base to a point on the major axis associated with the maximum width of the polygon. |
width.1e |
the width associated with 1/8 of |
width.2e |
the width associated with 2/8 of |
width.4e |
the width associated with 4/8 of |
width.6e |
the width associated with 6/8 of |
width.7e |
the width associated with 7/8 of |
bi.test |
the testing results for |
a |
the estimate of the intercept obtained using the box-counting method to calculate the fractal dimension of the polygon's boundary. |
sd.a |
the standard deviation of the estimated intercept. |
lci.a |
the lower bound of the 95% confidence interval of the estimated intercept. |
uci.a |
the upper bound of the 95% confidence interval of the estimated intercept. |
b |
the estimate of the slope obtained using the box-counting method to calculate the fractal dimension of the polygon's boundary. |
sd.b |
the standard deviation of the estimated slope. |
lci.a |
the lower bound of the 95% confidence interval of the estimated slope. |
uci.a |
the upper bound of the 95% confidence interval of the estimated slope. |
r.sq |
the coefficient of determination obtained when using the box-counting method to calculate the fractal dimension of the polygon's boundary. |
delta |
the vector of box sizes used in the box-counting method to calculate the fractal dimension of the polygon's boundary. |
N |
the number of boxes that include at least one pixel of the polygon's boundary. |
The polygon is expected to have an apparent major axis (e.g., the straight line through two points
on the polygon's boundary having the maximum distance or one that can be clearly defined to pass by
two landmarks on the polygon's boundary [i.e., the leaf length axis, the egg length axis, etc.]).
The polygon is placed with its major axis overlapping the -axis; the base of the polygon is
located at the origin; the apex of the polygon is located to the right of the base.
phi
is equal to angle
when angle
is not null.
In theory, n1 + n2 = n
, but in most cases n1 + n2
is slightly smaller than n
.
The reason is that very few boundary points fall outside the the lower and upper boundaries of the polygon
when using the intersect.owin
function in spatstat.geom. However, this does not considerably affect the results.
The log-transformed SI
and
the log-transformed AR
are demontrated to have a more symmetrical frequency distribution than their original forms.
This is important when performing an analysis of variance between (or among) groups to compared
their extents of bilateral (a)symmetry. See Shi et al. (2020) for details.
The box-counting approach uses a group of boxes (squares for simplicity) with different
sizes () to divide the leaf vein image into different parts. Let
represent the number
of boxes that include at least one pixel of the polygon's boundary.
The maximum of the range of the
coordinates and the range of the
coordinates
for the pixels of the polygon's boundary is defined as
. Let
represent the vector of
/
denomi.range
. We then used the following equation to calculate the fractal
dimension of the polygon's boundary:
where is the theoretical value of the fractal dimension. We can use its estimate as the
numerical value of the fractal dimension for the polygon's boundary.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Niinemets, Ü., Hui, C., Niklas, K.J., Yu, X., Hölscher, D. (2020) Leaf bilateral symmetry and the scaling of the perimeter vs. the surface area in 15 vine species. Forests 11, 246. doi:10.3390/f11020246
Shi, P., Zheng, X., Ratkowsky, D.A., Li, Y., Wang, P., Cheng, L. (2018) A simple method for measuring the bilateral symmetry of leaves. Symmetry 10, 118. doi:10.3390/sym10040118
data(bambooleaves) uni.C <- sort( unique(bambooleaves$Code) ) ind <- 3 Data <- bambooleaves[bambooleaves$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y dev.new() plot( x0, y0, asp=1, type="l", cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y)) ) Res1 <- adjdata(x0, y0, ub.np=2000, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y Res2 <- bilat( x=x1, y=y1, time.interval=0.00045, peri.np=NULL, auto.search=TRUE, fd.opt=TRUE ) Res2$scan.perimeter set.seed(123) Res3 <- bilat( x=x1, y=y1, time.interval=0.00045, peri.np=500, n.loop=30, auto.search=TRUE, fd.opt=FALSE ) Res3$scan.perimeter set.seed(123) Res4 <- bilat( x=x1, y=y1, time.interval=0.00045, peri.np=500, n.loop=30, auto.search=TRUE, fd.opt=FALSE, angle=pi/4 ) Res4$scan.perimeter set.seed(123) Res5 <- bilat( x=x1, y=y1, time.interval=0.00045, peri.np=500, n.loop=30, auto.search=TRUE, fd.opt=FALSE, angle=0 ) Res5$scan.perimeter if(interactive()){ # The angle between the leaf length axis (namely the straight # line through the leaf apex and base) and the horizontal axis # should be between 0 and pi/2 for a scanned leaf's profile. # Here, the user needs to first click the leaf apex, # and then click the leaf base. set.seed(123) Res6 <- bilat( x=x1, y=y1, time.interval=0.00045, peri.np=500, n.loop=30, auto.search=FALSE, fd.opt=FALSE, angle=NULL ) Res6$scan.perimeter } set.seed(NULL) graphics.off()
data(bambooleaves) uni.C <- sort( unique(bambooleaves$Code) ) ind <- 3 Data <- bambooleaves[bambooleaves$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y dev.new() plot( x0, y0, asp=1, type="l", cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y)) ) Res1 <- adjdata(x0, y0, ub.np=2000, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y Res2 <- bilat( x=x1, y=y1, time.interval=0.00045, peri.np=NULL, auto.search=TRUE, fd.opt=TRUE ) Res2$scan.perimeter set.seed(123) Res3 <- bilat( x=x1, y=y1, time.interval=0.00045, peri.np=500, n.loop=30, auto.search=TRUE, fd.opt=FALSE ) Res3$scan.perimeter set.seed(123) Res4 <- bilat( x=x1, y=y1, time.interval=0.00045, peri.np=500, n.loop=30, auto.search=TRUE, fd.opt=FALSE, angle=pi/4 ) Res4$scan.perimeter set.seed(123) Res5 <- bilat( x=x1, y=y1, time.interval=0.00045, peri.np=500, n.loop=30, auto.search=TRUE, fd.opt=FALSE, angle=0 ) Res5$scan.perimeter if(interactive()){ # The angle between the leaf length axis (namely the straight # line through the leaf apex and base) and the horizontal axis # should be between 0 and pi/2 for a scanned leaf's profile. # Here, the user needs to first click the leaf apex, # and then click the leaf base. set.seed(123) Res6 <- bilat( x=x1, y=y1, time.interval=0.00045, peri.np=500, n.loop=30, auto.search=FALSE, fd.opt=FALSE, angle=NULL ) Res6$scan.perimeter } set.seed(NULL) graphics.off()
Is used to simulate and fit biological geometries. 'biogeom' incorporates several novel universal parametric equations that can generate the profiles of bird eggs, flowers, linear and lanceolate leaves, seeds, starfish, and tree-rings (Gielis, 2003; Shi et al., 2020), three growth-rate curves representing the ontogenetic growth trajectories of animals and plants against time, and the axially symmetrical and integral forms of all these functions (Shi et al., 2017, 2021). The optimization method proposed by Nelder and Mead (1965) was used to estimate model parameters. 'biogeom' includes several real data sets of the boundary coordinates of natural shapes, including avian eggs, fruit, lanceolate and ovate leaves, tree rings, seeds, and sea stars,and can be potentially applied to other natural shapes. 'biogeom' can quantify the conspecific or interspecific similarity of natural outlines, and provides information with important ecological and evolutionary implications for the growth and form of living organisms. Please see Shi et al. (2022) for details.
The DESCRIPTION file:
Package: | biogeom |
Type: | Package |
Title: | Biological Geometries |
Version: | 1.4.3 |
Date: | 2024-03-21 |
Authors@R: | c(person(given="Peijian", family="Shi", email="[email protected]", role=c("aut", "cre")), person(given=c("Johan"), family="Gielis", email="[email protected]", role=c("aut")), person(given=c("Brady K."), family="Quinn", email="[email protected]", role=c("aut"))) |
Author: | Peijian Shi [aut, cre], Johan Gielis [aut], Brady K. Quinn [aut] |
Maintainer: | Peijian Shi <[email protected]> |
Imports: | spatstat.geom (>= 2.4-0) |
Description: | Is used to simulate and fit biological geometries. 'biogeom' incorporates several novel universal parametric equations that can generate the profiles of bird eggs, flowers, linear and lanceolate leaves, seeds, starfish, and tree-rings (Gielis (2003) <doi:10.3732/ajb.90.3.333>; Shi et al. (2020) <doi:10.3390/sym12040645>), three growth-rate curves representing the ontogenetic growth trajectories of animals and plants against time, and the axially symmetrical and integral forms of all these functions (Shi et al. (2017) <doi:10.1016/j.ecolmodel.2017.01.012>; Shi et al. (2021) <doi:10.3390/sym13081524>). The optimization method proposed by Nelder and Mead (1965) <doi:10.1093/comjnl/7.4.308> was used to estimate model parameters. 'biogeom' includes several real data sets of the boundary coordinates of natural shapes, including avian eggs, fruit, lanceolate and ovate leaves, tree rings, seeds, and sea stars,and can be potentially applied to other natural shapes. 'biogeom' can quantify the conspecific or interspecific similarity of natural outlines, and provides information with important ecological and evolutionary implications for the growth and form of living organisms. Please see Shi et al. (2022) <doi:10.1111/nyas.14862> for details. |
Depends: | R (>= 4.3.0) |
License: | GPL (>= 2) |
NeedsCompilation: | no |
Packaged: | 2024-03-21 14:14:58 UTC; PEIJIAN SHI |
Repository: | CRAN |
Date/Publication: | 2024-03-29 06:20:02 UTC |
Index of help topics:
DEPE Calculation of the First-Order Derivative of the Explicit Preston Equation DETE Calculation of the First-Order Derivative of the Explicit Troscianko Equation DNRGE Calculation of the First-Order Derivative of the Narushin-Romanov-Griffin Equation DSGE Calculation of the First-Order Derivative of the Simplified Gielis Equation EPE Calculation of the Ordinate For an Arbitrary Point on the Preston Curve in the Plane ETE Calculation of the Ordinate For an Arbitrary Point on the Troscianko Curve in the Plane GE Calculation of the Polar Radius of the Gielis Curve LeafSizeDist Leaf size distribution of _Shibataea chinensis_ MBriereE Modified Briere Equation MLRFE Modified Lobry-Rosso-Flandrois (LRF) Equation MPerformanceE Modified Performance Equation MbetaE Modified Beta Equation NRGE The Narushin-Romanov-Griffin Equation (NRGE) Neocinnamomum Leaf Boundary Data of Seven Species of _Neocinnamomum_ PE Calculation of the Abscissa, Ordinate and Distance From the Origin For an Arbitrary Point on the Preston Curve SCSE Sarabia-Castillo-Slottje Equation (SCSE) SHE Sitthiyot-Holasut Equation SarabiaE Sarabia Equation SurfaceAreaEPE Calculation of the Surface Area of An Egg Based on the Explicit Preston Equation SurfaceAreaETE Calculation of the Surface Area of An Egg Based on the Explicit Troscianko Equation SurfaceAreaNRGE Calculation of the Surface Area of An Egg Based on the Narushin-Romanov-Griffin Equation SurfaceAreaSGE Calculation of the Surface Area of An Egg Based on the Simplified Gielis Equation TE The Troscianko Equation (TE) TGE Calculation of the Polar Radius of the Twin Gielis Curve TSE The Todd-Smart Equation (TSE) VolumeEPE Calculation of the Volume of An Egg Based on the Explicit Preston Equation VolumeETE Calculation of the Volume of An Egg Based on the Explicit Troscianko Equation VolumeNRGE Calculation of the Volume of An Egg Based on the Narushin-Romanov-Griffin Equation VolumeSGE Calculation of the Volume of An Egg Based on the Simplified Gielis Equation adjdata Boundary Data Adjustment of A Polygon areaGE Area Calculation for the Gielis Curve Within [0, 2pi) areaovate Area Calculation for an Ovate Polygon bambooleaves Leaf Boundary Data of _Phyllostachys incarnata_ T. H. Wen (Poaceae: Bambusoideae) bilat Measure of the Extent of Bilateral Symmetry of A Polygon biogeom Biological Geometries curveEPE Drawing the Preston Curve Produced by the the Explicit Preston Equation curveETE Drawing the Troscianko Curve Produced by the Explicit Troscianko Equation curveGE Drawing the Gielis Curve curveNRGE Drawing the Egg Shape Predicted by the Narushin-Romanov-Griffin Equation curveovate Drawing the Ovate Leaf-Shape Curve eggs Egg Boundary Data of Nine Species of Birds fitEPE Data-Fitting Function for the Explicit Preston Equation fitETE Data-Fitting Function for the Explicit Troscianko Equation fitGE Data-Fitting Function for the Gielis Equation fitLorenz Data-Fitting Function for the Rotated and Right-Shifted Lorenz Curve fitNRGE Parameter Estimation for the Narushin-Romanov-Griffin Equation fitSuper Data-Fitting Function for the Superellipse Equation fitovate Data-Fitting Function for the Ovate Leaf-Shape Equation fitsigmoid Data-Fitting Function for the Sigmoid Growth Equation fracdim Calculation of Fractal Dimension of Lef Veins Based on the Box-Counting Method ginkgoseed Boundary Data of the Side Projections of _Ginkgo biloba_ Seeds kp Boundary Data of the Vertical Projections of _Koelreuteria paniculata_ Fruit lmPE Parameter Estimation for the Todd-Smart Equation lmTE Parameter Estimation for the Troscianko Equation shoots Height Growth Data of Bamboo Shoots sigmoid Sigmoid Growth Equation starfish Boundary Data of Eight Sea Stars veins Leaf Vein Data of _Michelia compressa_ whitespruce Planar Coordinates of _Picea glauca_ Tree Rings
We are deeply thankful to Cang Hui, Yang Li, Uwe Ligges, Valeriy G. Narushin, Ülo Niinemets, Karl J. Nikas, Honghua Ruan, David A. Ratkowsky, Julian Schrader, Rolf Turner, Lin Wang, and Victoria Wimmer for their valuable help during creating this package. This work was supported by the National Key Research and Development Program of China (Grant No. 2021YFD02200403) and Simon Stevin Institute for Geometry (Antwerpen, Belguim).
Peijian Shi [aut, cre], Johan Gielis [aut], Brady K. Quinn [aut]
Maintainer: Peijian Shi <[email protected]>
Gielis, J. (2003) A generic geometric transformation that unifies a wide range of natural
and abstract shapes. American Journal of Botany 90, 333338. doi:10.3732/ajb.90.3.333
Nelder, J.A., Mead, R. (1965). A simplex method for function minimization.
Computer Journal 7, 308313. doi:10.1093/comjnl/7.4.308
Shi, P., Fan, M., Ratkowsky, D.A., Huang, J., Wu, H., Chen, L., Fang, S.,
Zhang, C. (2017) Comparison of two ontogenetic growth equations for animals and plants.
Ecological Modelling 349, 110. doi:10.1016/j.ecolmodel.2017.01.012
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Ratkowsky, D.A., Gielis, J. (2020) The generalized Gielis geometric equation and its application. Symmetry 12, 645. doi:10.3390/sym12040645
Shi, P., Yu, K., Niklas, K.J., Schrader, J., Song, Y., Zhu, R., Li, Y., Wei, H., Ratkowsky, D.A. (2021) A general model for describing the ovate leaf shape. Symmetry, 13, 1524. doi:10.3390/sym13081524
curveEPE
is used to draw the Preston curve that is produced by the explicit Preston equation.
curveEPE(P, np = 5000, simpver = NULL, fig.opt = FALSE, deform.fun = NULL, Par = NULL, xlim = NULL, ylim = NULL, unit = NULL, main="")
curveEPE(P, np = 5000, simpver = NULL, fig.opt = FALSE, deform.fun = NULL, Par = NULL, xlim = NULL, ylim = NULL, unit = NULL, main="")
P |
the three location parameters and the parameters of the explicit Preston equation or one of its simplified versions. |
np |
the number of data points on the Preston curve. |
simpver |
an optional argument to use the simplfied version of the explicit Preston equation. |
fig.opt |
an optional argument to draw the Preston curve. |
deform.fun |
the deformation function used to describe the deviation from a theoretical Preston curve. |
Par |
the parameter(s) of the deformation function. |
xlim |
the range of the |
ylim |
the range of the |
unit |
the units of the |
main |
the main title of the figure. |
The first three elements of P
are location parameters. The first two are the planar coordinates of the transferred origin,
and the third is the angle between the major axis of the curve and the -axis. Here, the major axis is a straight line through
the two ends of an egg's profile (i.e., the mid-line of the egg's profile). The other arguments in
P
(except these first three location parameters), and simpver
should correspond to those of P
in EPE
.
deform.fun
should take the form as: deform.fun <- function(Par, z){...}
, where z
is
a two-dimensional matrix related to the and
values.
And the return value of
deform.fun
should be a list
with two variables x
and y
.
x |
the |
y |
the |
When the rotation angle is zero (i.e., the third element in P
is zero), np
data points
are distributed counterclockwise on the Preston curve from the rightmost end of the egg's profile to itself.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Preston, F.W. (1953) The shapes of birds' eggs. The Auk 70, 160182.
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Wang, L., Quinn, B.K., Gielis, J. (2023) A new program to estimate the parameters of Preston's equation, a general formula for describing the egg shape of birds. Symmetry 15, 231. doi:10.3390/sym15010231
Todd, P.H., Smart, I.H.M. (1984) The shape of birds' eggs. Journal of Theoretical Biology
106, 239243. doi:10.1016/0022-5193(84)90021-3
Para1 <- c(0, 0, 0, 10, 6, 0.325, -0.0415) curveEPE(P=Para1, simpver=1, fig.opt=TRUE) Para2 <- c(0, 0, pi, 10, 6, -0.325, -0.0415) curveEPE(P=Para2, simpver=1, fig.opt=TRUE) Para3 <- c(0, 0, 0, 10, 6, 0.325, -0.0415, 0.2) curveEPE(P=Para3, simpver=NULL, fig.opt=TRUE) Para4 <- c(0, 0, pi, 10, 6, -0.325, -0.0415, 0.2) curveEPE(P=Para4, simpver=NULL, fig.opt=TRUE) Para5 <- c(0, 0, pi/4, 10, 6, 0.325, -0.0415) curveEPE(P=Para5, simpver=1, fig.opt=TRUE, main="A rotated egg shape") # There is an example that introduces a deformation function in the egg-shape equation myfun <- function(Par, z){ x <- z[,1] y <- z[,2] k1 <- Par[1] k2 <- Par[2] y <- y - k1*(y+k2)^2 list(x=x, y=y) } deform.op <- curveEPE(P=Para1, np=5000, simpver=1, fig.opt=TRUE, deform.fun=myfun, Par=c(0.05, 8)) graphics.off()
Para1 <- c(0, 0, 0, 10, 6, 0.325, -0.0415) curveEPE(P=Para1, simpver=1, fig.opt=TRUE) Para2 <- c(0, 0, pi, 10, 6, -0.325, -0.0415) curveEPE(P=Para2, simpver=1, fig.opt=TRUE) Para3 <- c(0, 0, 0, 10, 6, 0.325, -0.0415, 0.2) curveEPE(P=Para3, simpver=NULL, fig.opt=TRUE) Para4 <- c(0, 0, pi, 10, 6, -0.325, -0.0415, 0.2) curveEPE(P=Para4, simpver=NULL, fig.opt=TRUE) Para5 <- c(0, 0, pi/4, 10, 6, 0.325, -0.0415) curveEPE(P=Para5, simpver=1, fig.opt=TRUE, main="A rotated egg shape") # There is an example that introduces a deformation function in the egg-shape equation myfun <- function(Par, z){ x <- z[,1] y <- z[,2] k1 <- Par[1] k2 <- Par[2] y <- y - k1*(y+k2)^2 list(x=x, y=y) } deform.op <- curveEPE(P=Para1, np=5000, simpver=1, fig.opt=TRUE, deform.fun=myfun, Par=c(0.05, 8)) graphics.off()
curveETE
is used to draw the Troscianko curve that is produced by the explicit Troscianko equation.
curveETE(P, np = 5000, fig.opt = FALSE, deform.fun = NULL, Par = NULL, xlim = NULL, ylim = NULL, unit = NULL, main="")
curveETE(P, np = 5000, fig.opt = FALSE, deform.fun = NULL, Par = NULL, xlim = NULL, ylim = NULL, unit = NULL, main="")
P |
the three location parameters and the parameters of the explicit Troscianko equation. |
np |
the number of data points on the Troscianko curve. |
fig.opt |
an optional argument to draw the Troscianko curve. |
deform.fun |
the deformation function used to describe the deviation from a theoretical Troscianko curve. |
Par |
the parameter(s) of the deformation function. |
xlim |
the range of the |
ylim |
the range of the |
unit |
the units of the |
main |
the main title of the figure. |
The first three elements of P
are location parameters. The first two are the planar coordinates of the transferred origin,
and the third is the angle between the major axis of the curve and the -axis. Here, the major axis is a straight line through
the two ends of an egg's profile (i.e., the mid-line of the egg's profile). The other arguments in
P
(except these first three location parameters) should correspond to those of P
in ETE
.
deform.fun
should take the form as: deform.fun <- function(Par, z){...}
, where z
is
a two-dimensional matrix related to the and
values.
And the return value of
deform.fun
should be a list
with two variables x
and y
.
x |
the |
y |
the |
When the rotation angle is zero (i.e., the third element in P
is zero), np
data points
are distributed counterclockwise on the Troscianko curve from the rightmost end of the egg's profile to itself.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Biggins, J.D., Montgomeries, R.M., Thompson, J.E., Birkhead, T.R. (2022) Preston’s universal formula for avian egg shape. Ornithology 139, ukac028. doi:10.1093/ornithology/ukac028
Biggins, J.D., Thompson, J.E., Birkhead, T.R. (2018) Accurately quantifying
the shape of birds' eggs. Ecology and Evolution 8, 97289738. doi:10.1002/ece3.4412
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Wang, L., Quinn, B.K., Gielis, J. (2023) A new program to estimate the parameters of Preston's equation, a general formula for describing the egg shape of birds. Symmetry 15, 231. doi:10.3390/sym15010231
Para1 <- c(0, 0, 0, 2.25, -0.377, -0.29, -0.16) curveETE(P=Para1, fig.opt=TRUE) # There is an example that introduces a deformation function in the egg-shape equation myfun <- function(Par, z){ x <- z[,1] y <- z[,2] k1 <- Par[1] k2 <- Par[2] y <- y - k1*(y+k2)^2 list(x=x, y=y) } deform.op <- curveETE(P=Para1, np=5000, fig.opt=TRUE, deform.fun=myfun, Par=c(0.05, 8)) graphics.off()
Para1 <- c(0, 0, 0, 2.25, -0.377, -0.29, -0.16) curveETE(P=Para1, fig.opt=TRUE) # There is an example that introduces a deformation function in the egg-shape equation myfun <- function(Par, z){ x <- z[,1] y <- z[,2] k1 <- Par[1] k2 <- Par[2] y <- y - k1*(y+k2)^2 list(x=x, y=y) } deform.op <- curveETE(P=Para1, np=5000, fig.opt=TRUE, deform.fun=myfun, Par=c(0.05, 8)) graphics.off()
curveGE
is used to draw the Gielis curve.
curveGE(expr, P, phi = seq(0, 2*pi, len = 2000), m = 1, simpver = NULL, nval = 1, fig.opt = FALSE, deform.fun = NULL, Par = NULL, xlim = NULL, ylim = NULL, unit = NULL, main="")
curveGE(expr, P, phi = seq(0, 2*pi, len = 2000), m = 1, simpver = NULL, nval = 1, fig.opt = FALSE, deform.fun = NULL, Par = NULL, xlim = NULL, ylim = NULL, unit = NULL, main="")
expr |
the original (or twin) Gielis equation or one of its simplified versions. |
P |
the three location parameters and the parameters of the original (or twin) Gielis equation or one of its simplified versions. |
phi |
the given polar angles at which we want to draw the Gielis curve. |
m |
the given |
simpver |
an optional argument to use the simplfied version of the original (or twin) Gielis equation. |
nval |
the specified value for |
fig.opt |
an optional argument to draw the Gielis curve. |
deform.fun |
the deformation function used to describe the deviation from a theoretical Gielis curve. |
Par |
the parameter(s) of the deformation function. |
xlim |
the range of the |
ylim |
the range of the |
unit |
the units of the |
main |
the main title of the figure. |
The first three elements of P
are location parameters. The first two are the planar coordinates of the transferred polar point,
and the third is the angle between the major axis of the curve and the -axis. The other arguments in
P
(except these first three location parameters), m
, simpver
, and nval
should correspond
to expr
(i.e., GE
or TGE
).
Please note the differences in the simplified version number and the
number of parameters between GE
and TGE
.
deform.fun
should take the form as: deform.fun <- function(Par, z){...}
, where z
is
a two-dimensional matrix related to the and
values.
And the return value of
deform.fun
should be a list
with two variables x
and y
.
x |
the |
y |
the |
r |
the polar radii of the Gielis curve corresponding to the given polar angles |
simpver
in GE
is different from that in TGE
.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Gielis, J. (2003) A generic geometric transformation that unifies a wide range of natural and abstract shapes. American Journal of Botany 90, 333-338. doi:10.3732/ajb.90.3.333
Li, Y., Quinn, B.K., Gielis, J., Li, Y., Shi, P. (2022) Evidence that supertriangles exist in nature from the vertical projections of Koelreuteria paniculata fruit. Symmetry 14, 23. doi:10.3390/sym14010023
Shi, P., Gielis, J., Niklas, K.J. (2022) Comparison of a universal (but complex) model for avian egg
shape with a simpler model. Annals of the New York Academy of Sciences 1514, 3442. doi:10.1111/nyas.14799
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Ratkowsky, D.A., Gielis, J. (2020) The generalized Gielis geometric equation and its application. Symmetry 12, 645. doi:10.3390/sym12040645
Shi, P., Xu, Q., Sandhu, H.S., Gielis, J., Ding, Y., Li, H., Dong, X. (2015) Comparison of dwarf bamboos (Indocalamus sp.) leaf parameters to determine relationship between spatial density of plants and total leaf area per plant. Ecology and Evolution 5, 4578-4589. doi:10.1002/ece3.1728
GE.par <- c(2, 1, 4, 6, 3) phi.vec <- seq(0, 2*pi, len=2000) r.theor <- GE(P=GE.par, phi=phi.vec, m=5) dev.new() plot( phi.vec, r.theor, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(phi)), ylab=expression(italic("r")), type="l", col=4 ) curve.par <- c(1, 1, pi/4, GE.par) GE.res <- curveGE(GE, P=curve.par, fig.opt=TRUE, deform.fun=NULL, Par=NULL, m=5) # GE.res$r GE.res <- curveGE( GE, P=c(0, 0, 0, 2, 4, 20), m=1, simpver=1, fig.opt=TRUE ) # GE.res$r GE.res <- curveGE( GE, P=c(1, 1, pi/4, 2, 1, 3), m=5, simpver=1, fig.opt=TRUE ) # GE.res$r GE.res <- curveGE( GE, P=c(1, 1, pi/4, 2, 1, 3), m=2, simpver=1, fig.opt=TRUE ) # GE.res$r GE.res <- curveGE( GE, P=c(1, 1, pi/4, 2, 0.05), m=1, simpver=2, fig.opt=TRUE ) # GE.res$r GE.res <- curveGE( GE, P=c(1, 1, pi/4, 2), m=4, simpver=3, nval=2, fig.opt=TRUE ) # GE.res$r GE.res <- curveGE( GE, P=c(1, 1, pi/4, 2, 0.6), m=4, simpver=8, nval=2, fig.opt=TRUE ) # GE.res$r graphics.off()
GE.par <- c(2, 1, 4, 6, 3) phi.vec <- seq(0, 2*pi, len=2000) r.theor <- GE(P=GE.par, phi=phi.vec, m=5) dev.new() plot( phi.vec, r.theor, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(phi)), ylab=expression(italic("r")), type="l", col=4 ) curve.par <- c(1, 1, pi/4, GE.par) GE.res <- curveGE(GE, P=curve.par, fig.opt=TRUE, deform.fun=NULL, Par=NULL, m=5) # GE.res$r GE.res <- curveGE( GE, P=c(0, 0, 0, 2, 4, 20), m=1, simpver=1, fig.opt=TRUE ) # GE.res$r GE.res <- curveGE( GE, P=c(1, 1, pi/4, 2, 1, 3), m=5, simpver=1, fig.opt=TRUE ) # GE.res$r GE.res <- curveGE( GE, P=c(1, 1, pi/4, 2, 1, 3), m=2, simpver=1, fig.opt=TRUE ) # GE.res$r GE.res <- curveGE( GE, P=c(1, 1, pi/4, 2, 0.05), m=1, simpver=2, fig.opt=TRUE ) # GE.res$r GE.res <- curveGE( GE, P=c(1, 1, pi/4, 2), m=4, simpver=3, nval=2, fig.opt=TRUE ) # GE.res$r GE.res <- curveGE( GE, P=c(1, 1, pi/4, 2, 0.6), m=4, simpver=8, nval=2, fig.opt=TRUE ) # GE.res$r graphics.off()
curveNRGE
is used to draw the egg shape predicted by the Narushin-Romanov-Griffin equation.
curveNRGE(P, np = 5000, fig.opt = FALSE, deform.fun = NULL, Par = NULL, xlim = NULL, ylim = NULL, unit = NULL, main = "")
curveNRGE(P, np = 5000, fig.opt = FALSE, deform.fun = NULL, Par = NULL, xlim = NULL, ylim = NULL, unit = NULL, main = "")
P |
the three location parameters and the four parameters of the Narushin-Romanov-Griffin equation (Narushin et al., 2021). |
np |
the number of data points on the Narushin-Romanov-Griffin curve. |
fig.opt |
an optional argument to draw the Narushin-Romanov-Griffin curve. |
deform.fun |
the deformation function used to describe the deviation from a theoretical Narushin-Romanov-Griffin curve. |
Par |
the parameter(s) of the deformation function. |
xlim |
the range of the |
ylim |
the range of the |
unit |
the units of the |
main |
the main title of the figure. |
The first three elements of P
are location parameters. The first two are the planar coordinates of the transferred origin,
and the third is the angle between the major axis of the curve and the -axis. The other arguments in
P
should be the same as those in NRGE
.
deform.fun
should take the form as: deform.fun <- function(Par, z){...}
, where z
is
a two-dimensional matrix related to the and
values.
And the return value of
deform.fun
should be a list
with two variables x
and y
.
x |
the |
y |
the |
When the rotation angle is zero (i.e., the third element in P
is zero), np
data points
are distributed counterclockwise on the Narushin-Romanov-Griffin curve from the rightmost end of the egg's profile to itself.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Narushin, V.G., Romanov, M.N., Griffin, D.K. (2021) Egg and math: introducing a universal formula for egg shape.
Annals of the New York Academy of Sciences 1505, 169177. doi:10.1111/nyas.14680
Shi, P., Gielis, J., Niklas, K.J. (2022) Comparison of a universal (but complex) model for avian egg
shape with a simpler model. Annals of the New York Academy of Sciences 1514, 3442. doi:10.1111/nyas.14799
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Tian, F., Wang, Y., Sandhu, H.S., Gielis, J., Shi, P. (2020) Comparison of seed morphology of two ginkgo cultivars.
Journal of Forestry Research 31, 751758. doi:10.1007/s11676-018-0770-y
PA <- c(1, 1, pi/4, 11.5, 7.8, 1.1, 5.6) resA <- curveNRGE(PA, np=5000, fig.opt=TRUE) resB <- curveNRGE(PA, np=5000, fig.opt=TRUE, xlim=c(-6, 6), ylim=c(-6, 6), main="A pear-shaped egg") cbind(resB$x, resB$y) graphics.off()
PA <- c(1, 1, pi/4, 11.5, 7.8, 1.1, 5.6) resA <- curveNRGE(PA, np=5000, fig.opt=TRUE) resB <- curveNRGE(PA, np=5000, fig.opt=TRUE, xlim=c(-6, 6), ylim=c(-6, 6), main="A pear-shaped egg") cbind(resB$x, resB$y) graphics.off()
curveovate
is used to draw the ovate leaf-shape curve.
curveovate(expr, P, x, fig.opt = FALSE, deform.fun = NULL, Par = NULL, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
curveovate(expr, P, x, fig.opt = FALSE, deform.fun = NULL, Par = NULL, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
expr |
the simplified version 1 of a performance equation. |
P |
the three location parameters and the parameters of the simplified version 1 of a performance equation. |
x |
the given |
fig.opt |
an optional argument to draw the ovate leaf-shape curve. |
deform.fun |
the deformation function used to describe the deviation from a theoretical ovate leaf-shape curve. |
Par |
the parameter(s) of the deformation function. |
xlim |
the range of the |
ylim |
the range of the |
unit |
the units of the |
main |
the main title of the figure. |
P
has two types of elements: three location parameters, and model parameters.
This means that expr
is limited to be the simplified version 1 (where )
in
MbetaE
, MBriereE
, MLRFE
, and MPerformanceE
.
The first three elements of P
are location parameters, among which the first two
are the planar coordinates of the transferred origin,
and the third is the angle between the major axis of the curve and the -axis.
deform.fun
should take the form as: deform.fun <- function(Par, z){...}
, where z
is
a two-dimensional matrix related to the and
values.
And the return value of
deform.fun
should be a list
with two variables x
and y
.
x |
the |
y |
the |
The number of elements in P
here has additional three location parameters than that
in MbetaE
, MBriereE
, MLRFE
, MPerformanceE
.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Jin, J., Quinn, B.K., Shi, P. (2022) The modified Brière equation and its applications. Plants 11, 1769. doi:10.3390/plants11131769
Huey, R.B., Stevenson, R.D. (1979) Integrating thermal physiology and ecology of ectotherms:
a discussion of approaches. American Zoologist 19, 357366. doi:10.1093/icb/19.1.357
Li, Y., Zheng, Y., Ratkowsky, D.A., Wei, H., Shi, P. (2022) Application of an ovate leaf shape model to evaluate leaf bilateral asymmetry and calculate lamina centroid location. Frontiers in Plant Science 12, 822907. doi:10.3389/fpls.2021.822907
Lian, M., Shi, P., Zhang, L., Yao, W., Gielis, J., Niklas, K.J. (2023) A generalized performance equation
and its application in measuring the Gini index of leaf size inequality.
Trees - Structure and Function 37, 15551565. doi:10.1007/s00468-023-02448-8
Shi, P., Fan, M., Ratkowsky, D.A., Huang, J., Wu, H., Chen, L., Fang, S.,
Zhang, C. (2017) Comparison of two ontogenetic growth equations for animals and plants.
Ecological Modelling 349, 110. doi:10.1016/j.ecolmodel.2017.01.012
Shi, P., Ge, F., Sun, Y., Chen, C. (2011) A simple model for describing
the effect of temperature on insect developmental rate. Journal of Asia-Pacific Entomology
14, 1520. doi:10.1016/j.aspen.2010.11.008
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Yu, K., Niklas, K.J., Schrader, J., Song, Y., Zhu, R., Li, Y., Wei, H., Ratkowsky, D.A. (2021) A general model for describing the ovate leaf shape. Symmetry 13, 1524. doi:10.3390/sym13081524
areaovate
, fitovate
, MbetaE
,
MBriereE
, MLRFE
, MPerformanceE
P1 <- c(1, 1, pi/4, 2, 3, 10, 4) RE1 <- curveovate(MLRFE, P=P1, x=seq(0, 10, by=0.1), fig.opt=TRUE) RE2 <- curveovate(MbetaE, P=P1, x=seq(0, 10, by=0.1), fig.opt=TRUE) dev.new() plot(RE1$x, RE1$y, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic(x)), ylab=expression(italic(y))) lines(RE2$x, RE2$y, col=4) P3 <- c(1, 1, pi/4, 2.4, 0.96, 0.64, 7.75, 1.76, 3.68) RE3 <- curveovate(MPerformanceE, P=P3, x=seq(0, 7.75, by=0.01), fig.opt=TRUE) dev.new() plot(RE3$x, RE3$y, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic(x)), ylab=expression(italic(y))) graphics.off()
P1 <- c(1, 1, pi/4, 2, 3, 10, 4) RE1 <- curveovate(MLRFE, P=P1, x=seq(0, 10, by=0.1), fig.opt=TRUE) RE2 <- curveovate(MbetaE, P=P1, x=seq(0, 10, by=0.1), fig.opt=TRUE) dev.new() plot(RE1$x, RE1$y, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic(x)), ylab=expression(italic(y))) lines(RE2$x, RE2$y, col=4) P3 <- c(1, 1, pi/4, 2.4, 0.96, 0.64, 7.75, 1.76, 3.68) RE3 <- curveovate(MPerformanceE, P=P3, x=seq(0, 7.75, by=0.01), fig.opt=TRUE) dev.new() plot(RE3$x, RE3$y, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic(x)), ylab=expression(italic(y))) graphics.off()
DEPE
is used to calculate the first-order derivative of the explicit Preston equation at a given x-value.
DEPE(P, x, simpver = NULL)
DEPE(P, x, simpver = NULL)
P |
the parameters of the explicit Preston equation or one of its simplified versions. |
x |
the x-value used in the explicit Preston equation. |
simpver |
an optional argument to use the simplified version of the explicit Preston equation. |
When simpver = NULL
, the first-order derivative of
the explicit Preston equation at a given x-value is selected:
where P
has five parameters: ,
,
,
, and
.
When
simpver = 1
, the first-order derivative of the simplified version 1 is selected:
where P
has four parameters: ,
,
, and
.
When
simpver = 2
, the first-order derivative of the simplified version 2 is selected:
where P
has three parameters: ,
, and
.
When
simpver = 3
, the first-order derivative of the simplified version 3 is selected:
where P
has three parameters: ,
, and
.
The argument P
in the DEPE
function has the same parameters, as those in the
EPE
function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
Par3 <- c(4.27, 2.90, 0.0868, 0.0224, -0.0287) xx1 <- seq(-4.27, 4.27, by=0.001) f1 <- DEPE(P=Par3, x=xx1, simpver=NULL) f2 <- -DEPE(P=Par3, x=xx1, simpver=NULL) dev.new() plot(xx1, f1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlim=c(-5, 5), ylim=c(-35, 35), xlab=expression(italic(x)), ylab=expression(paste(italic(f), "(", italic(x), ")", sep=""))) lines(xx1, f2, col=2) graphics.off()
Par3 <- c(4.27, 2.90, 0.0868, 0.0224, -0.0287) xx1 <- seq(-4.27, 4.27, by=0.001) f1 <- DEPE(P=Par3, x=xx1, simpver=NULL) f2 <- -DEPE(P=Par3, x=xx1, simpver=NULL) dev.new() plot(xx1, f1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlim=c(-5, 5), ylim=c(-35, 35), xlab=expression(italic(x)), ylab=expression(paste(italic(f), "(", italic(x), ")", sep=""))) lines(xx1, f2, col=2) graphics.off()
DETE
is used to calculate the first-order derivative of the explicit Troscianko equation at a given x-value.
DETE(P, x)
DETE(P, x)
P |
the parameters of the explicit Troscianko equation. |
x |
the x-value used in the explicit Troscianko equation. |
The first-order derivative of the explicit Troscianko equation at a given x-value is:
where P
has four parameters: ,
,
, and
.
The argument P
in the DETE
function has the same parameters, as those in the
ETE
function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
Par5 <- c(2.25, -0.38, -0.29, -0.16) xx2 <- seq(-2.25, 2.25, by=0.001) h1 <- DETE(P=Par5, x=xx2) h2 <- -DETE(P=Par5, x=xx2) ind <- which(is.na(h1) | is.na(h2)) xx2 <- xx2[-ind] h1 <- h1[-ind] h2 <- h2[-ind] dev.new() plot(xx2, h1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlim=c(-2.25, 2.25), ylim=c(-30, 30), xlab=expression(italic(x)), ylab=expression(paste(italic(h), "(", italic(x), ")", sep=""))) lines(xx2, h2, col=2) graphics.off()
Par5 <- c(2.25, -0.38, -0.29, -0.16) xx2 <- seq(-2.25, 2.25, by=0.001) h1 <- DETE(P=Par5, x=xx2) h2 <- -DETE(P=Par5, x=xx2) ind <- which(is.na(h1) | is.na(h2)) xx2 <- xx2[-ind] h1 <- h1[-ind] h2 <- h2[-ind] dev.new() plot(xx2, h1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlim=c(-2.25, 2.25), ylim=c(-30, 30), xlab=expression(italic(x)), ylab=expression(paste(italic(h), "(", italic(x), ")", sep=""))) lines(xx2, h2, col=2) graphics.off()
DNRGE
is used to calculate the first-order derivative of the Narushin-Romanov-Griffin equation at a given x-value.
DNRGE(P, x)
DNRGE(P, x)
P |
the parameters of the Narushin-Romanov-Griffin equation. |
x |
the x-value used in the Narushin-Romanov-Griffin equation. |
Let us define:
and then the first-order derivative of the Narushin-Romanov-Griffin equation at a given x-value is:
where P
has four parameters: ,
,
, and
.
The argument P
in the DNRGE
function has the same parameters, as those in the
NRGE
function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Narushin, V.G., Romanov, M.N., Griffin, D.K. (2021) Egg and math: introducing a universal formula for egg shape.
Annals of the New York Academy of Sciences 1505, 169177. doi:10.1111/nyas.14680
Narushin, V.G., Romanov, M.N., Mishra, B., Griffin, D.K. (2022) Mathematical progression of
avian egg shape with associated area and volume determinations.
Annals of the New York Academy of Sciences 1513, 6578. doi:10.1111/nyas.14771
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
fitNRGE
, NRGE
, SurfaceAreaNRGE
Par6 <- c(4.51, 3.18, 0.1227, 2.2284) xx3 <- seq(-4.51/2, 4.51/2, len=2000) J1 <- DNRGE(P=Par6, x=xx3) J2 <- -DNRGE(P=Par6, x=xx3) ind <- which(is.na(J1) | is.na(J2)) xx3 <- xx3[-ind] J1 <- J1[-ind] J2 <- J2[-ind] dev.new() plot(xx3, J1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlim=c(-4.51/2, 4.51/2), ylim=c(-20, 20), xlab=expression(italic(x)), ylab=expression(paste(italic(J), "(", italic(x), ")", sep=""))) lines(xx3, J2, col=2) graphics.off()
Par6 <- c(4.51, 3.18, 0.1227, 2.2284) xx3 <- seq(-4.51/2, 4.51/2, len=2000) J1 <- DNRGE(P=Par6, x=xx3) J2 <- -DNRGE(P=Par6, x=xx3) ind <- which(is.na(J1) | is.na(J2)) xx3 <- xx3[-ind] J1 <- J1[-ind] J2 <- J2[-ind] dev.new() plot(xx3, J1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlim=c(-4.51/2, 4.51/2), ylim=c(-20, 20), xlab=expression(italic(x)), ylab=expression(paste(italic(J), "(", italic(x), ")", sep=""))) lines(xx3, J2, col=2) graphics.off()
DSGE
is used to calculate the first-order derivative of the simplified Gielis equation at a given -value.
DSGE(P, phi)
DSGE(P, phi)
P |
the parameters of the simplified Gielis equation, including |
phi |
the |
The first-order derivative of the simplified Gielis equation with arguments simpver = 1
and m = 1
at a given -value is:
where P
has three parameters: ,
, and
.
The argument P
in the DSGE
function only has
the three parameters: ,
, and
.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Chen, Z. (2012) Volume and area of revolution under polar coordinate system.
Studies in College Mathematics 15(6), 911.
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
Par7 <- c(1.124, 14.86, 49.43) phi1 <- seq(0, pi, len=2000) g1 <- DSGE(P=Par7, phi=phi1) dev.new() plot(phi1, g1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(varphi)), ylab=expression(paste(italic(g), "(", italic(varphi), ")", sep=""))) graphics.off()
Par7 <- c(1.124, 14.86, 49.43) phi1 <- seq(0, pi, len=2000) g1 <- DSGE(P=Par7, phi=phi1) dev.new() plot(phi1, g1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(varphi)), ylab=expression(paste(italic(g), "(", italic(varphi), ")", sep=""))) graphics.off()
The data consist of the egg boundary data of nine species of birds.
data(eggs)
data(eggs)
In the data set, there are four columns of variables: Code
, LatinName
, x
, and y
.
Code
saves the codes of individual eggs;
LatinName
saves the Latin names of the nine species of birds;
x
saves the coordinates of the egg boundary in the Cartesian coordinate system (cm);
and
y
saves the coordinates of the egg boundary in the Cartesian coordinate system (cm).
In
Code
, codes 1-9 represent Strix uralensis, Dromaius novaehollandiae,
Turdus philomelos, Gallus gallus, Pandion haliaetus, Uria aalge,
Uria lomvia, Gallinago media, and Aptenodytes patagonicus, respectively.
Narushin, V.G., Romanov, M.N., Griffin, D.K. (2021) Egg and math: introducing a universal formula for egg shape.
Annals of the New York Academy of Sciences 1505, 169177. doi:10.1111/nyas.14680
Shi, P., Gielis, J., Niklas, K.J. (2022) Comparison of a universal (but complex) model for avian egg
shape with a simpler model. Annals of the New York Academy of Sciences 1514, 3442. doi:10.1111/nyas.14799
Tian, F., Wang, Y., Sandhu, H.S., Gielis, J., Shi, P. (2020) Comparison of seed morphology of two ginkgo cultivars.
Journal of Forestry Research 31, 751758. doi:10.1007/s11676-018-0770-y
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 8 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=1000, times=1.2, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( x1, y1, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic("x")), ylab=expression(italic("y")), pch=1, col=1 ) Res2 <- adjdata(x0, y0, ub.np=60, times=1, len.pro=1/2, index.sp=20) x2 <- Res2$x y2 <- Res2$y Res3 <- adjdata(x0, y0, ub.np=60, times=1, len.pro=1/2, index.sp=100) x3 <- Res3$x y3 <- Res3$y dev.new() plot( x2, y2, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")), pch=1, col=4 ) points( x3, y3, col=2) graphics.off()
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 8 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=1000, times=1.2, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( x1, y1, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic("x")), ylab=expression(italic("y")), pch=1, col=1 ) Res2 <- adjdata(x0, y0, ub.np=60, times=1, len.pro=1/2, index.sp=20) x2 <- Res2$x y2 <- Res2$y Res3 <- adjdata(x0, y0, ub.np=60, times=1, len.pro=1/2, index.sp=100) x3 <- Res3$x y3 <- Res3$y dev.new() plot( x2, y2, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")), pch=1, col=4 ) points( x3, y3, col=2) graphics.off()
EPE
is used to calculate the y-value for an arbitrary point on the Preston curve
that was generated by the explicit Preston equation or one of its simplified versions for a given x-value.
EPE(P, x, simpver = NULL)
EPE(P, x, simpver = NULL)
P |
the parameters of the explicit Preston equation or one of its simplified versions. |
x |
the x-value used in the explicit Preston equation. |
simpver |
an optional argument to use the simplified version of the explicit Preston equation. |
When simpver = NULL
, the explicit Preston equation is selected:
where P
has five parameters: ,
,
,
, and
.
When
simpver = 1
, the simplified version 1 is selected:
where P
has four parameters: ,
,
, and
.
When
simpver = 2
, the simplified version 2 is selected:
where P
has three parameters: ,
, and
.
When
simpver = 3
, the simplified version 3 is selected:
where P
has three parameters: ,
, and
.
The values predicted by the explicit Preston equation.
We only considered the upper part of the egg-shape curve in the above expressions because
the lower part is symmetrical to the upper part around the x-axis.
The mid-line of an egg's profile in EPE
is aligned to
the x-axis, while the mid-line of an egg's profile in PE
is aligned to the y-axis. The EPE
function has the same parameters,
P
, as those in the PE
function.
The explicit Preston equation is used for calculating an egg's volume and surface area,
when the parameters, which are saved in the P
vector,
are obtained using the fitEPE
function
or the lmPE
function based on the TSE
function.
In addition, the values in x
> (i.e., the first element in
P
)
are forced to be , and the values in
x
< will be forced to be
.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
Shi, P., Wang, L., Quinn, B.K., Gielis, J. (2023) A new program to estimate the parameters of Preston's equation, a general formula for describing the egg shape of birds. Symmetry 15, 231. doi:10.3390/sym15010231
curveEPE
, fitEPE
, PE
, SurfaceAreaEPE
, VolumeEPE
Par3 <- c(4.27, 2.90, 0.0868, 0.0224, -0.0287) xx1 <- seq(-4.27, 4.27, by=0.001) yy1 <- EPE(P=Par3, x=xx1, simpver=NULL) yy2 <- -EPE(P=Par3, x=xx1, simpver=NULL) dev.new() plot(xx1, yy1, asp=1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlim=c(-5, 5), ylim=c(-5, 5), xlab=expression(italic(x)), ylab=expression(italic(y))) lines(xx1, yy2, col=2) graphics.off()
Par3 <- c(4.27, 2.90, 0.0868, 0.0224, -0.0287) xx1 <- seq(-4.27, 4.27, by=0.001) yy1 <- EPE(P=Par3, x=xx1, simpver=NULL) yy2 <- -EPE(P=Par3, x=xx1, simpver=NULL) dev.new() plot(xx1, yy1, asp=1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlim=c(-5, 5), ylim=c(-5, 5), xlab=expression(italic(x)), ylab=expression(italic(y))) lines(xx1, yy2, col=2) graphics.off()
ETE
is used to calculate the y-value for an arbitrary point on the Troscianko curve
that was generated by the explicit Troscianko equation.
ETE(P, x)
ETE(P, x)
P |
the parameters of the explicit Troscianko equation, including |
x |
the x-value used in the explicit Troscianko equation. |
The explicit Troscianko equation is recommended as:
where and
represent the abscissa and ordinate of an arbitrary point on the explicit Troscianko curve;
,
,
, and
are parameters to be estimated.
The values predicted by the explicit Troscianko equation.
We only considered the upper part of the egg-shape curve in the above expressions because
the lower part is symmetrical to the upper part around the x-axis.
The mid-line of an egg's profile in ETE
is aligned to
the x-axis. The argument, P
, in the ETE
function has the same parameters,
,
, and
, as those in the
TE
function.
However, the former has an additional parameter than the latter, which represents half
the egg's length. The
lmTE
function is based on the TE
function, while
the fitETE
function is based on the ETE
function here.
In addition, the values in x
> (i.e., the first element in
P
)
are forced to be , and the values in
x
< will be forced to be
.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Biggins, J.D., Montgomeries, R.M., Thompson, J.E., Birkhead, T.R. (2022) Preston’s universal formula for avian egg shape. Ornithology 139, ukac028. doi:10.1093/ornithology/ukac028
Biggins, J.D., Thompson, J.E., Birkhead, T.R. (2018) Accurately quantifying
the shape of birds' eggs. Ecology and Evolution 8, 97289738. doi:10.1002/ece3.4412
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Wang, L., Quinn, B.K., Gielis, J. (2023) A new program to estimate the parameters of Preston's equation, a general formula for describing the egg shape of birds. Symmetry 15, 231. doi:10.3390/sym15010231
Troscianko, J. (2014). A simple tool for calculating egg shape, volume and surface area from digital images.
Ibis, 156, 874878. doi:10.1111/ibi.12177
curveETE
, fitETE
, SurfaceAreaETE
, VolumeETE
Par5 <- c(2.25, -0.38, -0.29, -0.16) xx2 <- seq(-2.25, 2.25, len=2000) yy3 <- ETE(P=Par5, x=xx2) yy4 <- -ETE(P=Par5, x=xx2) dev.new() plot(xx2, yy3, asp=1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlim=c(-3, 3), ylim=c(-3, 3), xlab=expression(italic(x)), ylab=expression(italic(y))) lines(xx2, yy4, col=2) graphics.off()
Par5 <- c(2.25, -0.38, -0.29, -0.16) xx2 <- seq(-2.25, 2.25, len=2000) yy3 <- ETE(P=Par5, x=xx2) yy4 <- -ETE(P=Par5, x=xx2) dev.new() plot(xx2, yy3, asp=1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlim=c(-3, 3), ylim=c(-3, 3), xlab=expression(italic(x)), ylab=expression(italic(y))) lines(xx2, yy4, col=2) graphics.off()
fitEPE
is used to estimate the parameters of the explicit Preston equation
or one of its simplified versions.
fitEPE(x, y, ini.val, simpver = NULL, control = list(), par.list = FALSE, stand.fig = TRUE, angle = NULL, fig.opt = FALSE, np = 2000, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
fitEPE(x, y, ini.val, simpver = NULL, control = list(), par.list = FALSE, stand.fig = TRUE, angle = NULL, fig.opt = FALSE, np = 2000, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
x |
the |
y |
the |
ini.val |
the list of initial values for the model parameters. |
simpver |
an optional argument to use the simplified version of the explicit Preston equation. |
control |
the list of control parameters for using the |
par.list |
the option of showing the list of parameters on the screen. |
stand.fig |
the option of drawing the observed and predicted profiles of an egg at the standard state
(i.e., the egg's centre is located at (0, 0), and the mid-line is aligned to the |
angle |
the angle between the mid-line and the |
fig.opt |
an optional argument of drawing the observed and predicted profiles of an egg
at arbitrary angle between the major axis and the |
np |
the number of data points on the predicted explicit Preston curve. |
xlim |
the range of the |
ylim |
the range of the |
unit |
the unit of the |
main |
the main title of the figure. |
The simpver
argument should correspond to EPE
. Here, the major axis
(i.e., the mid-line of an egg's profile) is the straight line trhough the two ends of the egg's length.
The Nelder-Mead algorithm (Nelder and Mead, 1965) is used to carry out the optimization of minimizing
the residual sum of squares (RSS) between the observed and predicted values.
The
optim
function in package stats was used to carry out the Nelder-Mead algorithm.
When angle = NULL
, the observed egg's profile will be shown at its initial angle in the scanned image;
when angle
is a numerical value (e.g., ) defined by the user, it indicates that the major axis
is rotated by the amount (
) counterclockwise from the
-axis.
par |
the estimates of the model parameters. |
scan.length |
the observed length of the egg's profile. |
scan.width |
the observed width of the egg's profile. |
scan.area |
the observed area of the egg's profile. |
scan.perimeter |
the observed perimeter of the egg's profile. |
r.sq |
the coefficient of determination between the observed and predicted |
RSS |
the residual sum of squares between the observed and predicted |
sample.size |
the number of data points used in the data fitting. |
x.stand.obs |
the observed |
y.stand.obs |
the observed |
y.stand.pred |
the predicted |
x.obs |
the observed |
y.obs |
the observed |
y.pred |
the predicted |
In the outputs, there are no x.stand.pred
and x.pred
, because y.stand.obs
and
y.stand.pred
share the same values (i.e.,
x.stand.obs
), and y.obs
and
y.pred
share the same values (i.e.,
x.obs
).
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Nelder, J.A., Mead, R. (1965) A simplex method for function minimization.
Computer Journal 7, 308313. doi:10.1093/comjnl/7.4.308
Preston, F.W. (1953) The shapes of birds' eggs. The Auk 70, 160182.
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Wang, L., Quinn, B.K., Gielis, J. (2023) A new program to estimate the parameters of Preston's equation, a general formula for describing the egg shape of birds. Symmetry 15, 231. doi:10.3390/sym15010231
Todd, P.H., Smart, I.H.M. (1984) The shape of birds' eggs. Journal of Theoretical Biology
106, 239243. doi:10.1016/0022-5193(84)90021-3
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 8 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=2000, times=1.2, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( x1, y1, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", col=4, xlab=expression(italic("x")), ylab=expression(italic("y")) ) simpver <- NULL res1 <- lmPE( x1, y1, simpver=simpver, dev.angle=seq(-0.05, 0.05, by=0.0001), unit="cm", fig.opt=FALSE ) x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- res1$theta a.ini <- res1$scan.length / 2 b.ini <- res1$scan.width / 2 c1.ini <- res1$par[2] / res1$par[1] c2.ini <- res1$par[3] / res1$par[1] c3.ini <- res1$par[4] / res1$par[1] ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, b.ini, c1.ini, c2.ini, c3.ini) res0 <- fitEPE( x=x1, y=y1, ini.val=ini.val, simpver=simpver, unit="cm", par.list=FALSE, stand.fig=FALSE, angle=NULL, fig.opt=FALSE, control=list(reltol=1e-30, maxit=50000), np=2000 ) n.loop <- 12 Show <- FALSE for(i in 1:n.loop){ ini.val <- res0$par if(i==n.loop) Show <- TRUE print(paste(i, "/", n.loop, sep="")) res0 <- fitEPE( x=x1, y=y1, ini.val=ini.val, simpver=simpver, unit="cm", par.list=FALSE, stand.fig=Show, angle=pi/4, fig.opt=Show, control=list(reltol=1e-30, maxit=50000), np=2000 ) } # The numerical values of the location and model parameters res0$par # The root-mean-square error (RMSE) between # the observed and predicted y values sqrt(res0$RSS/res0$sample.size) sqrt(sum((res0$y.stand.obs-res0$y.stand.pred)^2)/length(res0$y.stand.obs)) # To calculate the volume of the egg VolumeEPE(P=res0$par[4:8]) # To calculate the surface area of the egg SurfaceAreaEPE(P=res0$par[4:8]) graphics.off()
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 8 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=2000, times=1.2, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( x1, y1, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", col=4, xlab=expression(italic("x")), ylab=expression(italic("y")) ) simpver <- NULL res1 <- lmPE( x1, y1, simpver=simpver, dev.angle=seq(-0.05, 0.05, by=0.0001), unit="cm", fig.opt=FALSE ) x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- res1$theta a.ini <- res1$scan.length / 2 b.ini <- res1$scan.width / 2 c1.ini <- res1$par[2] / res1$par[1] c2.ini <- res1$par[3] / res1$par[1] c3.ini <- res1$par[4] / res1$par[1] ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, b.ini, c1.ini, c2.ini, c3.ini) res0 <- fitEPE( x=x1, y=y1, ini.val=ini.val, simpver=simpver, unit="cm", par.list=FALSE, stand.fig=FALSE, angle=NULL, fig.opt=FALSE, control=list(reltol=1e-30, maxit=50000), np=2000 ) n.loop <- 12 Show <- FALSE for(i in 1:n.loop){ ini.val <- res0$par if(i==n.loop) Show <- TRUE print(paste(i, "/", n.loop, sep="")) res0 <- fitEPE( x=x1, y=y1, ini.val=ini.val, simpver=simpver, unit="cm", par.list=FALSE, stand.fig=Show, angle=pi/4, fig.opt=Show, control=list(reltol=1e-30, maxit=50000), np=2000 ) } # The numerical values of the location and model parameters res0$par # The root-mean-square error (RMSE) between # the observed and predicted y values sqrt(res0$RSS/res0$sample.size) sqrt(sum((res0$y.stand.obs-res0$y.stand.pred)^2)/length(res0$y.stand.obs)) # To calculate the volume of the egg VolumeEPE(P=res0$par[4:8]) # To calculate the surface area of the egg SurfaceAreaEPE(P=res0$par[4:8]) graphics.off()
fitETE
is used to estimate the parameters of the explicit Troscianko equation.
fitETE(x, y, ini.val, control = list(), par.list = FALSE, stand.fig = TRUE, angle = NULL, fig.opt = FALSE, np = 2000, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
fitETE(x, y, ini.val, control = list(), par.list = FALSE, stand.fig = TRUE, angle = NULL, fig.opt = FALSE, np = 2000, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
x |
the |
y |
the |
ini.val |
the list of initial values for the model parameters. |
control |
the list of control parameters for using the |
par.list |
the option of showing the list of parameters on the screen. |
stand.fig |
the option of drawing the observed and predicted profiles of an egg at the standard state
(i.e., the egg's centre is located at (0, 0), and the mid-line is aligned to the |
angle |
the angle between the mid-line and the |
fig.opt |
an optional argument of drawing the observed and predicted profiles of an egg
at arbitrary angle between the major axis and the |
np |
the number of data points on the predicted explicit Troscianko curve. |
xlim |
the range of the |
ylim |
the range of the |
unit |
the unit of the |
main |
the main title of the figure. |
Here, the major axis (i.e., the mid-line of an egg's profile) is the straight line trhough the two ends of
the egg's length. The Nelder-Mead algorithm (Nelder and Mead, 1965) is used to carry out the optimization
of minimizing the residual sum of squares (RSS) between the observed and predicted values.
The
optim
function in package stats was used to carry out the Nelder-Mead algorithm.
When angle = NULL
, the observed egg's profile will be shown at its initial angle in the scanned image;
when angle
is a numerical value (e.g., ) defined by the user, it indicates that the major axis
is rotated by the amount (
) counterclockwise from the
-axis.
par |
the estimates of the model parameters. |
scan.length |
the observed length of the egg's profile. |
scan.width |
the observed width of the egg's profile. |
scan.area |
the observed area of the egg's profile. |
scan.perimeter |
the observed perimeter of the egg's profile. |
r.sq |
the coefficient of determination between the observed and predicted |
RSS |
the residual sum of squares between the observed and predicted |
sample.size |
the number of data points used in the data fitting. |
x.stand.obs |
the observed |
y.stand.obs |
the observed |
y.stand.pred |
the predicted |
x.obs |
the observed |
y.obs |
the observed |
y.pred |
the predicted |
In the outputs, there are no x.stand.pred
and x.pred
, because y.stand.obs
and
y.stand.pred
share the same values (i.e.,
x.stand.obs
), and y.obs
and
y.pred
share the same values (i.e.,
x.obs
).
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Biggins, J.D., Montgomeries, R.M., Thompson, J.E., Birkhead, T.R. (2022) Preston’s universal formula for avian egg shape. Ornithology 139, ukac028. doi:10.1093/ornithology/ukac028
Biggins, J.D., Thompson, J.E., Birkhead, T.R. (2018) Accurately quantifying
the shape of birds' eggs. Ecology and Evolution 8, 97289738. doi:10.1002/ece3.4412
Nelder, J.A., Mead, R. (1965) A simplex method for function minimization.
Computer Journal 7, 308313. doi:10.1093/comjnl/7.4.308
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Wang, L., Quinn, B.K., Gielis, J. (2023) A new program to estimate the parameters of Preston's equation, a general formula for describing the egg shape of birds. Symmetry 15, 231. doi:10.3390/sym15010231
Troscianko, J. (2014). A simple tool for calculating egg shape, volume and surface area from digital images.
Ibis, 156, 874878. doi:10.1111/ibi.12177
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 8 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=2000, times=1.2, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( x1, y1, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", col=4, xlab=expression(italic("x")), ylab=expression(italic("y")) ) res1 <- lmTE( x1, y1, unit="cm", fig.opt=FALSE ) if(FALSE){ P0 <- c(res1$scan.length/2, res1$par) xx <- seq(-res1$scan.length/2, res1$scan.length/2, len=2000) yy1 <- ETE(P0, xx) yy2 <- -ETE(P0, xx) dev.new() plot( xx, yy1, cex.lab=1.5, cex.axis=1.5, asp=1, col=2, ylim=c(-res1$scan.length/2, res1$scan.length/2), type="l", xlab=expression(x), ylab=expression(y) ) lines( xx, yy2, col=4 ) } x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- res1$theta a.ini <- res1$scan.length / 2 alpha0.ini <- res1$par[1] alpha1.ini <- res1$par[2] alpha2.ini <- res1$par[3] ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, alpha0.ini, alpha1.ini, alpha2.ini) res0 <- fitETE( x=x1, y=y1, ini.val=ini.val, unit="cm", par.list=FALSE, stand.fig=FALSE, angle=NULL, fig.opt=FALSE, control=list(reltol=1e-30, maxit=50000), np=2000 ) n.loop <- 12 Show <- FALSE for(i in 1:n.loop){ ini.val <- res0$par if(i==n.loop) Show <- TRUE print(paste(i, "/", n.loop, sep="")) res0 <- fitETE( x=x1, y=y1, ini.val=ini.val, unit="cm", par.list=FALSE, stand.fig=Show, angle=pi/4, fig.opt=Show, control=list(reltol=1e-30, maxit=50000), np=2000 ) } # The numerical values of the location and model parameters res0$par # The root-mean-square error (RMSE) between # the observed and predicted y values sqrt(res0$RSS/res0$sample.size) sqrt(sum((res0$y.stand.obs-res0$y.stand.pred)^2)/length(res0$y.stand.obs)) # To calculate the volume of the egg VolumeETE(P=res0$par[4:7]) # To calculate the surface area of the egg SurfaceAreaETE(P=res0$par[4:7]) graphics.off()
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 8 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=2000, times=1.2, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( x1, y1, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", col=4, xlab=expression(italic("x")), ylab=expression(italic("y")) ) res1 <- lmTE( x1, y1, unit="cm", fig.opt=FALSE ) if(FALSE){ P0 <- c(res1$scan.length/2, res1$par) xx <- seq(-res1$scan.length/2, res1$scan.length/2, len=2000) yy1 <- ETE(P0, xx) yy2 <- -ETE(P0, xx) dev.new() plot( xx, yy1, cex.lab=1.5, cex.axis=1.5, asp=1, col=2, ylim=c(-res1$scan.length/2, res1$scan.length/2), type="l", xlab=expression(x), ylab=expression(y) ) lines( xx, yy2, col=4 ) } x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- res1$theta a.ini <- res1$scan.length / 2 alpha0.ini <- res1$par[1] alpha1.ini <- res1$par[2] alpha2.ini <- res1$par[3] ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, alpha0.ini, alpha1.ini, alpha2.ini) res0 <- fitETE( x=x1, y=y1, ini.val=ini.val, unit="cm", par.list=FALSE, stand.fig=FALSE, angle=NULL, fig.opt=FALSE, control=list(reltol=1e-30, maxit=50000), np=2000 ) n.loop <- 12 Show <- FALSE for(i in 1:n.loop){ ini.val <- res0$par if(i==n.loop) Show <- TRUE print(paste(i, "/", n.loop, sep="")) res0 <- fitETE( x=x1, y=y1, ini.val=ini.val, unit="cm", par.list=FALSE, stand.fig=Show, angle=pi/4, fig.opt=Show, control=list(reltol=1e-30, maxit=50000), np=2000 ) } # The numerical values of the location and model parameters res0$par # The root-mean-square error (RMSE) between # the observed and predicted y values sqrt(res0$RSS/res0$sample.size) sqrt(sum((res0$y.stand.obs-res0$y.stand.pred)^2)/length(res0$y.stand.obs)) # To calculate the volume of the egg VolumeETE(P=res0$par[4:7]) # To calculate the surface area of the egg SurfaceAreaETE(P=res0$par[4:7]) graphics.off()
fitGE
is used to estimate the parameters of the original (or twin) Gielis equation
or one of its simplified versions.
fitGE(expr, x, y, ini.val, m = 1, simpver = NULL, nval = nval, control = list(), par.list = FALSE, stand.fig = TRUE, angle = NULL, fig.opt = FALSE, np = 2000, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
fitGE(expr, x, y, ini.val, m = 1, simpver = NULL, nval = nval, control = list(), par.list = FALSE, stand.fig = TRUE, angle = NULL, fig.opt = FALSE, np = 2000, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
expr |
the original (or twin) Gielis equation or one of its simplified versions. |
x |
the |
y |
the |
ini.val |
the list of initial values for the model parameters. |
m |
the given |
simpver |
an optional argument to use the simplified version of the original (or twin) Gielis equation. |
nval |
the specified value for |
control |
the list of control parameters for using the |
par.list |
the option of showing the list of parameters on the screen. |
stand.fig |
the option of drawing the observed and predicted polygons at the standard state
(i.e., the polar point is located at (0, 0), and the major axis overlaps with the |
angle |
the angle between the major axis and the |
fig.opt |
an optional argument of drawing the observed and predicted polygons at arbitrary angle
between the major axis and the |
np |
the number of data points on the predicted Gielis curve. |
xlim |
the range of the |
ylim |
the range of the |
unit |
the unit of the |
main |
the main title of the figure. |
The arguments of m
, simpver
, and nval
should correspond
to expr
(i.e., GE
or TGE
). Please note the differences in the simplified version number and the
number of parameters between GE
and TGE
.
The Nelder-Mead algorithm (Nelder and Mead, 1965) is used to carry out the optimization of minimizing
the residual sum of squares (RSS) between the observed and predicted radii.
The optim
function in package stats was used to carry out the Nelder-Mead algorithm.
When angle = NULL
, the observed polygon will be shown at its initial angle in the scanned image;
when angle
is a numerical value (e.g., ) defined by the user, it indicates that the major axis
is rotated by the amount (
) counterclockwise from the
-axis.
par |
the estimates of the model parameters. |
scan.length |
the observed length of the polygon. |
scan.width |
the observed width of the polygon. |
scan.area |
the observed area of the polygon. |
r.sq |
the coefficient of determination between the observed and predicted polar radii. |
RSS |
the residual sum of squares between the observed and predicted polar radii. |
sample.size |
the number of data points used in the data fitting. |
phi.stand.obs |
the polar angles at the standard state. |
phi.trans |
the transferred polar angles rotated as defined by the user. |
r.stand.obs |
the observed polar radii at the standard state. |
r.stand.pred |
the predicted polar radii at the standard state. |
x.stand.obs |
the observed |
x.stand.pred |
the predicted |
y.stand.obs |
the observed |
y.stand.pred |
the predicted |
r.obs |
the observed polar radii at the transferred polar angles as defined by the user. |
r.pred |
the predicted polar radii at the transferred polar angles as defined by the user. |
x.obs |
the observed |
x.pred |
the predicted |
y.obs |
the observed |
y.pred |
the predicted |
simpver
in GE
is different from that in TGE
.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Gielis, J. (2003) A generic geometric transformation that unifies a wide range of natural
and abstract shapes. American Journal of Botany 90, 333338. doi:10.3732/ajb.90.3.333
Li, Y., Quinn, B.K., Gielis, J., Li, Y., Shi, P. (2022) Evidence that supertriangles exist in nature from the vertical projections of Koelreuteria paniculata fruit. Symmetry 14, 23. doi:10.3390/sym14010023
Nelder, J.A., Mead, R. (1965) A simplex method for function minimization.
Computer Journal 7, 308313. doi:10.1093/comjnl/7.4.308
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Ratkowsky, D.A., Gielis, J. (2020) The generalized Gielis geometric equation and its application. Symmetry 12, 645. doi:10.3390/sym12040645
Shi, P., Xu, Q., Sandhu, H.S., Gielis, J., Ding, Y., Li, H., Dong, X. (2015) Comparison of dwarf bamboos (Indocalamus sp.) leaf parameters to determine relationship between spatial density of plants and total leaf area per plant. Ecology and Evolution 5, 4578-4589. doi:10.1002/ece3.1728
areaGE
, curveGE
, DSGE
, GE
,
SurfaceAreaSGE
, TGE
, VolumeSGE
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 1 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=200, times=1.2, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y Res2 <- adjdata(x0, y0, ub.np=40, times=1, len.pro=1/2, index.sp=20) x2 <- Res2$x y2 <- Res2$y Res3 <- adjdata(x0, y0, ub.np=100, times=1, len.pro=1/2, index.sp=100) x3 <- Res3$x y3 <- Res3$y dev.new() plot( x2, y2, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")), pch=1, col=4 ) points( x3, y3, col=2) x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- pi a.ini <- sqrt(2) * max( y0.ini-min(y1), x0.ini-min(x1) ) n1.ini <- c(5, 25) n2.ini <- c(15, 25) if(ind == 2){ n1.ini <- c(0.5, 1) n2.ini <- c(6, 12) } ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, n1.ini, n2.ini) Res4 <- fitGE( GE, x=x1, y=y1, ini.val=ini.val, m=1, simpver=1, nval=1, unit="cm", par.list=FALSE, fig.opt=TRUE, angle=NULL, control=list(reltol=1e-20, maxit=20000), np=2000 ) Res4$par sqrt(sum((Res4$y.stand.obs-Res4$y.stand.pred)^2)/Res4$sample.size) xx <- Res4$x.stand.obs yy <- Res4$y.stand.obs library(spatstat.geom) poly0 <- as.polygonal(owin(poly=list(x=xx, y=yy))) area(poly0) areaGE(GE, P = Res4$par[4:6], m=1, simpver=1) # The following code is used to # calculate the root-mean-square error (RMSE) in the y-coordinates ind1 <- which(yy >= 0) ind2 <- which(yy < 0) xx1 <- xx[ind1] # The upper part of the egg yy1 <- yy[ind1] xx2 <- xx[ind2] # The lower part of the egg yy2 <- yy[ind2] Para <- c(0, 0, 0, Res4$par[4:length(Res4$par)]) PartU <- curveGE(GE, P=Para, phi=seq(0, pi, len=100000), m=1, simpver=1, fig.opt=FALSE) xv1 <- PartU$x yv1 <- PartU$y PartL <- curveGE(GE, P=Para, phi=seq(pi, 2*pi, len=100000), m=1, simpver=1, fig.opt=FALSE) xv2 <- PartL$x yv2 <- PartL$y ind3 <- c() for(q in 1:length(xx1)){ ind.temp <- which.min(abs(xx1[q]-xv1)) ind3 <- c(ind3, ind.temp) } ind4 <- c() for(q in 1:length(xx2)){ ind.temp <- which.min(abs(xx2[q]-xv2)) ind4 <- c(ind4, ind.temp) } RSS <- sum((yy1-yv1[ind3])^2) + sum((yy2-yv2[ind4])^2) RMSE <- sqrt( RSS/length(yy) ) # To calculate the volume of the Gielis egg when simpver=1 & m=1 VolumeSGE(P=Res4$par[4:6]) # To calculate the surface area of the Gielis egg when simpver=1 & m=1 SurfaceAreaSGE(P=Res4$par[4:6]) graphics.off()
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 1 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=200, times=1.2, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y Res2 <- adjdata(x0, y0, ub.np=40, times=1, len.pro=1/2, index.sp=20) x2 <- Res2$x y2 <- Res2$y Res3 <- adjdata(x0, y0, ub.np=100, times=1, len.pro=1/2, index.sp=100) x3 <- Res3$x y3 <- Res3$y dev.new() plot( x2, y2, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")), pch=1, col=4 ) points( x3, y3, col=2) x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- pi a.ini <- sqrt(2) * max( y0.ini-min(y1), x0.ini-min(x1) ) n1.ini <- c(5, 25) n2.ini <- c(15, 25) if(ind == 2){ n1.ini <- c(0.5, 1) n2.ini <- c(6, 12) } ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, n1.ini, n2.ini) Res4 <- fitGE( GE, x=x1, y=y1, ini.val=ini.val, m=1, simpver=1, nval=1, unit="cm", par.list=FALSE, fig.opt=TRUE, angle=NULL, control=list(reltol=1e-20, maxit=20000), np=2000 ) Res4$par sqrt(sum((Res4$y.stand.obs-Res4$y.stand.pred)^2)/Res4$sample.size) xx <- Res4$x.stand.obs yy <- Res4$y.stand.obs library(spatstat.geom) poly0 <- as.polygonal(owin(poly=list(x=xx, y=yy))) area(poly0) areaGE(GE, P = Res4$par[4:6], m=1, simpver=1) # The following code is used to # calculate the root-mean-square error (RMSE) in the y-coordinates ind1 <- which(yy >= 0) ind2 <- which(yy < 0) xx1 <- xx[ind1] # The upper part of the egg yy1 <- yy[ind1] xx2 <- xx[ind2] # The lower part of the egg yy2 <- yy[ind2] Para <- c(0, 0, 0, Res4$par[4:length(Res4$par)]) PartU <- curveGE(GE, P=Para, phi=seq(0, pi, len=100000), m=1, simpver=1, fig.opt=FALSE) xv1 <- PartU$x yv1 <- PartU$y PartL <- curveGE(GE, P=Para, phi=seq(pi, 2*pi, len=100000), m=1, simpver=1, fig.opt=FALSE) xv2 <- PartL$x yv2 <- PartL$y ind3 <- c() for(q in 1:length(xx1)){ ind.temp <- which.min(abs(xx1[q]-xv1)) ind3 <- c(ind3, ind.temp) } ind4 <- c() for(q in 1:length(xx2)){ ind.temp <- which.min(abs(xx2[q]-xv2)) ind4 <- c(ind4, ind.temp) } RSS <- sum((yy1-yv1[ind3])^2) + sum((yy2-yv2[ind4])^2) RMSE <- sqrt( RSS/length(yy) ) # To calculate the volume of the Gielis egg when simpver=1 & m=1 VolumeSGE(P=Res4$par[4:6]) # To calculate the surface area of the Gielis egg when simpver=1 & m=1 SurfaceAreaSGE(P=Res4$par[4:6]) graphics.off()
fitLorenz
is used to estimate the parameters of the rotated and right-shifted Lorenz curve
using version 4 or 5 of MPerformanceE
, or the Lorenz equations including
SarabiaE
, SCSE
, and SHE
.
fitLorenz(expr, z, ini.val, simpver = 4, control = list(), par.list = FALSE, fig.opt = FALSE, np = 2000, xlab=NULL, ylab=NULL, main = NULL, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL, par.limit = TRUE)
fitLorenz(expr, z, ini.val, simpver = 4, control = list(), par.list = FALSE, fig.opt = FALSE, np = 2000, xlab=NULL, ylab=NULL, main = NULL, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL, par.limit = TRUE)
expr |
version 4 or 5 of |
z |
the observations of size distribution (i.e., the household income distribution, the leaf size distribution). |
ini.val |
the initial values of the model parameters. |
simpver |
an optional argument to use version 4 or 5 of |
control |
the list of control parameters for using the |
par.list |
the option of showing the list of parameters on the screen. |
fig.opt |
an optional argument to draw the original and rotated Lorenz curves. |
np |
the number of data points to draw the predicted original and rotated Lorenz curves. |
xlab |
the label of the |
ylab |
the label of the |
main |
the main title of the figure. |
subdivisions |
please see the arguments for the |
rel.tol |
please see the arguments for the |
abs.tol |
please see the arguments for the |
stop.on.error |
please see the arguments for the |
keep.xy |
please see the arguments for the |
aux |
please see the arguments for the |
par.limit |
an optional argument to limit the numerical ranges of model parameters of the three Lorenz
equations including |
Here, ini.val
only includes the initial values of the model parameters as a list.
The Nelder-Mead algorithm (Nelder and Mead, 1965) is used to carry out the optimization of minimizing the residual
sum of squares (RSS) between the observed and predicted values. The
optim
function in package stats was used to carry out the Nelder-Mead algorithm.
Here, versions 4 and 5 of MPerformanceE
and the Lorenz equations including
SarabiaE
, SCSE
, and SHE
can be used to fit the rotated
and right-shifted Lorenz curve.
When
simpver = 4
, the simplified version 4 of MPerformanceE
is selected:
There are five elements in P
, representing
the values of ,
,
,
, and
, respectively.
When
simpver = 5
, the simplified version 5 of MPerformanceE
is selected:
There are three elements in P
, representing
the values of ,
, and
, respectively.
For the Lorenz functions, the user can define any formulae that follow the below form:
Lorenz.fun <- function(P, x){...}
, where P
is the vector of parameter(s), x
is the preditor that
ranges between 0 and 1 representing the cumulative proportion of the number of individuals in a statistical
unit, and Lorenz.fun
is the name of a Lorenz function defined by the user, which also ranges between 0 and 1
representing the cumulative
proportion of the income or size in a statistical unit. This package provides three representative Lorenz
functions: SarabiaE
, SCSE
, and SHE
.
Here, the Gini coefficient (GC) is calculated as follows when
MPerformanceE
is selected:
where and
are the independent and dependent variables
in version 4 or 5 of
MPerformanceE
, respectively.
However, the Gini coefficient (GC) is calculated as follows when a Lorenz function, e.g.,
SCSE
, is selected:
where and
are the independent and dependent variables
in the Lorenz function, respectively.
For
SarabiaE
and SHE
, there are explicit formulae for GC
(Sarabia, 1997; Sitthiyot and Holasut, 2023).
x1 |
the cumulative proportion of the number of an entity of interest, i.e., the number of households of a city, the number of leaves of a plant. |
y1 |
the cumulative proportion of the size of an entity of interest. |
x |
the |
y |
the |
par |
the estimates of the model parameters. |
r.sq |
the coefficient of determination between the observed and predicted |
RSS |
the residual sum of squares between the observed and predicted |
sample.size |
the number of data points used in the data fitting. |
GC |
the calculated Gini coefficient. |
When MPerformanceE
is selected, the estimates of the model parameters
denote those in MPerformanceE
rather than
being obtained by directly fitting the y1
versus x1
data;
when a Lorenz function is selected, the estimates of the model parameters
denote those in the Lorenz function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Huey, R.B., Stevenson, R.D. (1979) Integrating thermal physiology and ecology of ectotherms:
a discussion of approaches. American Zoologist 19, 357366. doi:10.1093/icb/19.1.357
Lian, M., Shi, P., Zhang, L., Yao, W., Gielis, J., Niklas, K.J. (2023) A generalized performance equation
and its application in measuring the Gini index of leaf size inequality.
Trees Structure and Function 37, 1555
1565. doi:10.1007/s00468-023-02448-8
Lorenz, M.O. (1905) Methods of measuring the concentration of wealth.
Journal of the American Statistical Association 9(70), 209219. doi:10.2307/2276207
Nelder, J.A., Mead, R. (1965) A simplex method for function minimization.
Computer Journal 7, 308313. doi:10.1093/comjnl/7.4.308
Sarabia, J.-M. (1997) A hierarchy of Lorenz curves based on the generalized Tukey's lambda distribution.
Econometric Reviews 16, 305320. doi:10.1080/07474939708800389
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Sitthiyot, T., Holasut, K. (2023) A universal model for the Lorenz curve with novel applications for datasets containing zeros and/or exhibiting extreme inequality. Scientific Reports 13, 4729. doi:10.1038/s41598-023-31827-x
LeafSizeDist
, MPerformanceE
, SarabiaE
,
SCSE
, SHE
data(LeafSizeDist) CulmNumber <- c() for(i in 1:length(LeafSizeDist$Code)){ temp <- as.numeric( strsplit(LeafSizeDist$Code[i], "-", fixed=TRUE)[[1]][1] ) CulmNumber <- c(CulmNumber, temp) } uni.CN <- sort( unique(CulmNumber) ) ind <- CulmNumber==uni.CN[1] A0 <- LeafSizeDist$A[ind] ini.val1 <- list(0.5, 0.1, c(0.01, 0.1, 1, 5, 10), 1, 1) ini.val2 <- list(0.5, 0.1, c(0.01, 0.1, 1, 5, 10)) resu1 <- fitLorenz(MPerformanceE, z=A0, ini.val=ini.val1, simpver=4, fig.opt=TRUE) resu2 <- fitLorenz(MPerformanceE, z=A0, ini.val=ini.val2, simpver=5, fig.opt=TRUE) resu1$par resu2$par ini.val3 <- list(0.9, c(10, 50, 100, 500), 1, 0) resu3 <- fitLorenz( SarabiaE, z=A0, ini.val=ini.val3, par.limit=TRUE, fig.opt=TRUE, control=list(reltol=1e-20, maxit=10000) ) resu3$par resu3$GC graphics.off()
data(LeafSizeDist) CulmNumber <- c() for(i in 1:length(LeafSizeDist$Code)){ temp <- as.numeric( strsplit(LeafSizeDist$Code[i], "-", fixed=TRUE)[[1]][1] ) CulmNumber <- c(CulmNumber, temp) } uni.CN <- sort( unique(CulmNumber) ) ind <- CulmNumber==uni.CN[1] A0 <- LeafSizeDist$A[ind] ini.val1 <- list(0.5, 0.1, c(0.01, 0.1, 1, 5, 10), 1, 1) ini.val2 <- list(0.5, 0.1, c(0.01, 0.1, 1, 5, 10)) resu1 <- fitLorenz(MPerformanceE, z=A0, ini.val=ini.val1, simpver=4, fig.opt=TRUE) resu2 <- fitLorenz(MPerformanceE, z=A0, ini.val=ini.val2, simpver=5, fig.opt=TRUE) resu1$par resu2$par ini.val3 <- list(0.9, c(10, 50, 100, 500), 1, 0) resu3 <- fitLorenz( SarabiaE, z=A0, ini.val=ini.val3, par.limit=TRUE, fig.opt=TRUE, control=list(reltol=1e-20, maxit=10000) ) resu3$par resu3$GC graphics.off()
fitNRGE
is used to estimate the parameters of the Narushin-Romanov-Griffin equation.
fitNRGE(x, y, dev.angle = NULL, ini.C = c(-1, 0.1, 0.5, 1), strip.num = 2000, control = list(), fig.opt = TRUE, np = 2000, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
fitNRGE(x, y, dev.angle = NULL, ini.C = c(-1, 0.1, 0.5, 1), strip.num = 2000, control = list(), fig.opt = TRUE, np = 2000, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
x |
the |
y |
the |
dev.angle |
the angle of deviation for the axis associated with the maximum distance between two points on an egg's profile from the mid-line of the egg's profile. |
ini.C |
the initial value(s) of parameter |
strip.num |
the number of equidistant strips intersecting with the egg's boundary that are horizontally placed. See Shi et al. (2018, 2020) for details. |
control |
the list of control parameters for using the |
fig.opt |
an optional argument to draw the observed and predicted egg's boundaries. |
np |
the number of data points on the predicted Narushin-Romanov-Griffin curve. |
xlim |
the range of the |
ylim |
the range of the |
unit |
the units of the |
main |
the main title of the figure. |
The NRGE (see NRGE
) has a complex model structure with four parameters (i.e., ,
,
, and
).
Because three out of four parameters of NRGE have clear biological and geometric meanings (i.e.,
,
, and
),
their values could be estimated by means of numerical calculation. After obtaining the numerical values of the three parameters,
the Nelder-Mead algorithm (Nelder and Mead, 1965) was used to estimate
.
Because of the failure of the optimization method to estimate the major axis (i.e., the mid-line) and model parameters of NRGE,
it was difficult to define the egg length axis, although it is essential for calculating
,
, and
.
For this reason, two methods were used to obtain the major axis: the maximum distance method, and
the longest axis adjustment method. In the first method, the straight line through two points forming the maximum distance
on the egg's profile is defined as the major axis (i.e., the mid-line). In the second method, we assume that there is an angle
of deviation for the longest axis (i.e., the axis associated with the maximum distance between two points on an egg's profile)
from the mid-line of the egg's profile. That is to say,
the mid-line of an egg's profile is not the axis associated with the maximum distance between two points on the egg's profile.
When
angle = NULL
, the maximum distance method is used; when angle
is a numerical value or a numerical vector,
the longest axis adjustment method is used. Here, the numerical value of dev.angle
is not
the angle of deviation for the major axis of an egg's profile from the -axis, and instead it is the angle of deviation
for the longest axis (associated with the maximum distance between two points on the egg's profile) from the mid-line of the egg's profile.
Once the major axis is established, the distance of the major axis can be calculated as the estimate of
.
Using the maximum distance method,
equals the maximum distance. Using the longest axis adjustment method,
may be slightly smaller than the maximum distance. After rotating the major axis to make it overlap with the
-axis,
a large number of equidistant strips can be used (e.g., 2000) from the egg base to egg tip to intersect the egg's boundary.
This methodology makes it easy to obtain the maximum egg's breadth (i.e.,
) and
. The residual sum of squares (RSS)
between the observed and predicted
values can be minimized using an optimization method (Nelder and Mead, 1965) to estimate
.
Despite the complex structure of NRGE (see
NRGE
), the optimization method for estimating the remaining parameter
becomes feasible after the other three parameters have been numerically estimated. Please see Shi et al. (2022) for details.
theta |
the angle between the longest axis of an egg's profile (i.e.,
the axis associated with the maximum distance between two points on the egg's profile) and the |
epsilon |
the optimal angle of deviation for the longest axis (associated with the maximum distance between two points on an egg's profile)
from the mid-line of the egg's profile, when |
RSS.vector |
the vector of residual sum of squares corresponding to |
x.obs |
the observed |
y.obs |
the observed |
y.pred |
the predicted |
par |
the estimates of the four model parameters in the Narushin-Romanov-Griffin equation. |
scan.length |
the length of the egg's boundary. The default is the maximum distance between two points on the egg's boundary. |
scan.width |
the maximum width of the egg's boundary. |
scan.area |
the area of the egg's boundary. |
scan.perimeter |
the perimeter of the egg's boundary based on all data points on the egg's boundary. |
RSS |
the residual sum of squares between the observed and predicted |
sample.size |
the number of data points used in the numerical calculation. |
RMSE |
the root-mean-square errors between the observed and predicted |
theta
is the calculated angle between the longest axis (i.e., the axis associated with the maximum distance
between two points on an egg's profile) and the -axis, and
epsilon
is the calculated angle of deviation for the longest
axis from the mid-line of the egg's profile. This means that the angle between the mid-line and the -axis is equal to
theta
+ epsilon
. In the outputs, there is no x.pred
, because y.obs
and
y.pred
share the same values (i.e.,
x.obs
).
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Narushin, V.G., Romanov, M.N., Griffin, D.K. (2021) Egg and math: introducing a universal formula for egg shape.
Annals of the New York Academy of Sciences 1505, 169177. doi:10.1111/nyas.14680
Nelder, J.A., Mead, R. (1965). A simplex method for function minimization.
Computer Journal 7, 308313. doi:10.1093/comjnl/7.4.308
Shi, P., Gielis, J., Niklas, K.J. (2022) Comparison of a universal (but complex) model for avian egg
shape with a simpler model. Annals of the New York Academy of Sciences 1514, 3442. doi:10.1111/nyas.14799
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Niinemets, Ü., Hui, C., Niklas, K.J., Yu, X., Hölscher, D. (2020) Leaf bilateral symmetry and the scaling of the perimeter vs. the surface area in 15 vine species. Forests 11, 246. doi:10.3390/f11020246
Shi, P., Zheng, X., Ratkowsky, D.A., Li, Y., Wang, P., Cheng, L. (2018) A simple method for measuring the bilateral symmetry of leaves. Symmetry 10, 118. doi:10.3390/sym10040118
Tian, F., Wang, Y., Sandhu, H.S., Gielis, J., Shi, P. (2020) Comparison of seed morphology of two ginkgo cultivars.
Journal of Forestry Research 31, 751758. doi:10.1007/s11676-018-0770-y
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 8 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=3000, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) Res2 <- fitNRGE(x1, y1, dev.angle=NULL, ini.C=c(-1, -0.1, seq(0.1, 1, by=0.05)), strip.num=2000, fig.opt=TRUE) dev.new() plot(Res2$x.obs, Res2$y.obs, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")), type="l", col=4) lines( Res2$x.obs, Res2$y.pred, col=2) Res3 <- fitNRGE(x1, y1, dev.angle=seq(-0.05, 0.05, by=0.01), ini.C=c(-1, -0.1, seq(0.1, 1, by=0.05)), strip.num=2000, fig.opt=TRUE) zeta <- Res3$theta + Res3$epsilon x2 <- x1*cos(zeta) + y1*sin(zeta) y2 <- y1*cos(zeta) - x1*sin(zeta) plot( x2-min(x2), y2-min(y2), asp=1, col="grey70", cex=1, xlab=expression(italic("x")), ylab=expression(italic("y")) ) lines(Res3$x.obs-min(Res3$x.obs), Res3$y.obs-min(Res3$y.obs), col=4) lines(Res3$x.obs-min(Res3$x.obs), Res3$y.pred-min(Res3$y.obs), col=2) RMSE <- sqrt( Res3$RSS / Res3$sample.size ) RMSE # To calculate the volume of the egg VolumeNRGE(P=Res3$par) # To calculate the surface area of the egg SurfaceAreaNRGE(P=Res3$par) graphics.off()
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 8 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=3000, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) Res2 <- fitNRGE(x1, y1, dev.angle=NULL, ini.C=c(-1, -0.1, seq(0.1, 1, by=0.05)), strip.num=2000, fig.opt=TRUE) dev.new() plot(Res2$x.obs, Res2$y.obs, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")), type="l", col=4) lines( Res2$x.obs, Res2$y.pred, col=2) Res3 <- fitNRGE(x1, y1, dev.angle=seq(-0.05, 0.05, by=0.01), ini.C=c(-1, -0.1, seq(0.1, 1, by=0.05)), strip.num=2000, fig.opt=TRUE) zeta <- Res3$theta + Res3$epsilon x2 <- x1*cos(zeta) + y1*sin(zeta) y2 <- y1*cos(zeta) - x1*sin(zeta) plot( x2-min(x2), y2-min(y2), asp=1, col="grey70", cex=1, xlab=expression(italic("x")), ylab=expression(italic("y")) ) lines(Res3$x.obs-min(Res3$x.obs), Res3$y.obs-min(Res3$y.obs), col=4) lines(Res3$x.obs-min(Res3$x.obs), Res3$y.pred-min(Res3$y.obs), col=2) RMSE <- sqrt( Res3$RSS / Res3$sample.size ) RMSE # To calculate the volume of the egg VolumeNRGE(P=Res3$par) # To calculate the surface area of the egg SurfaceAreaNRGE(P=Res3$par) graphics.off()
fitovate
is used to estimate the parameters of a simplified performance equation.
fitovate(expr, x, y, ini.val, par.list = FALSE, stand.fig = TRUE, control = list(), angle = NULL, fig.opt = FALSE, index.xmax = 3, np = 2000, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
fitovate(expr, x, y, ini.val, par.list = FALSE, stand.fig = TRUE, control = list(), angle = NULL, fig.opt = FALSE, index.xmax = 3, np = 2000, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
expr |
the simplified version 1 of the performance equations. |
x |
the |
y |
the |
ini.val |
the initial values of the simplified version 1 of a performance equation. |
par.list |
an optional argument to show the list of parameters on the screen. |
stand.fig |
an optional argument to draw the observed and predicted polygons' boundaries at the standard state
(i.e., the origin is located at (0, 0), and the major axis overlaps with the |
control |
the list of control parameters for using the |
angle |
the angle between the major axis of the polygon and the |
fig.opt |
an optional argument to draw the observed and predicted polygons at an arbitrary angle
between the major axis and the |
index.xmax |
the specified index in parameters representing |
np |
the number of data points on the predicted ovate leaf-shape curve. |
xlim |
the range of the |
ylim |
the range of the |
unit |
the units of the |
main |
the main title of the figure. |
ini.val
is a list for two types of parameters: three location parameters,
and model parameters. This means that expr
is limited to
being the simplified version 1 (where )
in
MbetaE
, MBriereE
, MLRFE
, and MPerformanceE
.
The initial values for the first three parameters in ini.val
are location parameters,
among which the first two are the planar coordinates of the transferred origin,
and the third is the angle between the major axis of the polygon and the -axis.
The Nelder-Mead algorithm (Nelder and Mead, 1965) is used to carry out the optimization of minimizing
the residual sum of squares (RSS) between the observed and predicted
-axis.
The
optim
function in package stats was used to carry out the Nelder-Mead algorithm.
When angle = NULL
, the observed polygon will be shown at its initial angle in the scanned image;
when angle
is a numerical value (e.g., ) defined by the user, it indicates that the major axis
is rotated by the amount (
) counterclockwise from the
-axis.
par |
the estimates of the model parameters. |
r.sq |
the coefficient of determination between the observed and predicted |
RSS |
the residual sum of squares between the observed and predicted |
sample.size |
the number of data points on the polygon's boundary in the data fitting. |
scan.length |
the observed length of the polygon's boundary. |
scan.width |
the observed width of the polygon's boundary. |
scan.perimeter |
the observed perimeter of the polygon's boundary. |
scan.area |
the observed area of the polygon's boundary. |
pred.length |
the predicted length of the polygon's boundary. |
pred.width |
the predicted width of the polygon's boundary. |
pred.perimeter |
the predicted perimeter of the polygon's boundary. |
pred.area |
the predicted area of the polygon's boundary. |
x.stand.obs |
the observed |
x.stand.pred |
the predicted |
y.stand.obs |
the observed |
y.stand.pred |
the predicted |
x.obs |
the observed |
x.pred |
the predicted |
y.obs |
the observed |
y.pred |
the predicted |
There are two types of parameters (i.e., three location parameters and model parameters)
for the value of par
. The transferred
angle denotes the angle between the major axis and the -axis. For the argument
index.xmax
,
the default is 3, which corresponds to the order of the model parameter of in
MbetaE
, MBriereE
, and MLRFE
. However, in MPerformanceE
,
index.xmax
should be 4 rather than 3.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Jin, J., Quinn, B.K., Shi, P. (2022) The modified Brière equation and its applications. Plants 11, 1769. doi:10.3390/plants11131769
Huey, R.B., Stevenson, R.D. (1979) Integrating thermal physiology and ecology of ectotherms:
a discussion of approaches. American Zoologist 19, 357366. doi:10.1093/icb/19.1.357
Li, Y., Zheng, Y., Ratkowsky, D.A., Wei, H., Shi, P. (2022) Application of an ovate leaf shape model to evaluate leaf bilateral asymmetry and calculate lamina centroid location. Frontiers in Plant Science 12, 822907. doi:10.3389/fpls.2021.822907
Lian, M., Shi, P., Zhang, L., Yao, W., Gielis, J., Niklas, K.J. (2023) A generalized performance equation
and its application in measuring the Gini index of leaf size inequality.
Trees Structure and Function 37, 1555
1565. doi:10.1007/s00468-023-02448-8
Nelder, J.A., Mead, R. (1965) A simplex method for function minimization.
Computer Journal 7, 308313. doi:10.1093/comjnl/7.4.308
Shi, P., Fan, M., Ratkowsky, D.A., Huang, J., Wu, H., Chen, L., Fang, S.,
Zhang, C. (2017) Comparison of two ontogenetic growth equations for animals and plants.
Ecological Modelling 349, 110. doi:10.1016/j.ecolmodel.2017.01.012
Shi, P., Ge, F., Sun, Y., Chen, C. (2011) A simple model for describing
the effect of temperature on insect developmental rate. Journal of Asia-Pacific Entomology
14, 1520. doi:10.1016/j.aspen.2010.11.008
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Yu, K., Niklas, K.J., Schrader, J., Song, Y., Zhu, R., Li, Y., Wei, H., Ratkowsky, D.A. (2021) A general model for describing the ovate leaf shape. Symmetry 13, 1524. doi:10.3390/sym13081524
areaovate
, curveovate
, MbetaE
,
MBriereE
, MLRFE
, MPerformanceE
data(Neocinnamomum) uni.C <- sort( unique(Neocinnamomum$Code) ) ind <- 2 Data <- Neocinnamomum[Neocinnamomum$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=200, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) x0.ini <- min( x1 ) y0.ini <- min( y1 ) theta.ini <- pi/4 len.max <- max( max(y1)-min(y1), max(x1)-min(x1) ) *2/sqrt(2) a.ini <- c(0.1, 0.01, 0.001, 0.0001) m.ini <- c(0.1, 0.5, 1, 2) x2.ini <- len.max delta.ini <- c(0.5, 1) ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, m.ini, x2.ini, delta.ini) Res1 <- fitovate(MBriereE, x=x1, y=y1, ini.val=ini.val, par.list=FALSE, fig.opt=TRUE, angle=pi/6, control=list(reltol=1e-20, maxit=20000), np=2000, unit=NULL) Res1$RSS x0.ini <- min( x1 ) y0.ini <- min( y1 ) theta.ini <- pi/4 len.max <- max( max(y1)-min(y1), max(x1)-min(x1) ) *2/sqrt(2) yc.ini <- len.max/3 xc.ini <- 1/4*len.max x2.ini <- len.max delta.ini <- c(0.5, seq(1, 5, by=5)) ini.val <- list(x0.ini, y0.ini, theta.ini, yc.ini, xc.ini, x2.ini, delta.ini) Res2 <- fitovate( MbetaE, x=x1, y=y1, ini.val=ini.val, par.list=TRUE, fig.opt=TRUE, angle=pi/3, control=list(reltol=1e-20, maxit=20000), np=2000, unit=NULL ) Res2$RSS Res3 <- fitovate( MLRFE, x=x1, y=y1, ini.val=ini.val, unit=NULL, par.list=FALSE, fig.opt=TRUE, angle=NULL, control=list(reltol=1e-20, maxit=20000), np=2000) Res3$RSS x0.ini <- min( x1 ) y0.ini <- min( y1 ) theta.ini <- pi/4 len.max <- max( max(y1)-min(y1), max(x1)-min(x1) ) *2/sqrt(2) c.ini <- 1/5*len.max K1.ini <- c(0.1, 1, 5, 10) K2.ini <- 1 x2.ini <- len.max a.ini <- 1 b.ini <- 1 ini.val <- list(x0.ini, y0.ini, theta.ini, c.ini, K1.ini, K2.ini, x2.ini, a.ini, b.ini) Res4 <- fitovate( MPerformanceE, x=x1, y=y1, ini.val=ini.val, par.list=TRUE, fig.opt=TRUE, index.xmax=4, angle=pi/3, control=list(reltol=1e-20, maxit=20000), np=2000, unit=NULL ) Res4$RSS graphics.off()
data(Neocinnamomum) uni.C <- sort( unique(Neocinnamomum$Code) ) ind <- 2 Data <- Neocinnamomum[Neocinnamomum$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=200, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) x0.ini <- min( x1 ) y0.ini <- min( y1 ) theta.ini <- pi/4 len.max <- max( max(y1)-min(y1), max(x1)-min(x1) ) *2/sqrt(2) a.ini <- c(0.1, 0.01, 0.001, 0.0001) m.ini <- c(0.1, 0.5, 1, 2) x2.ini <- len.max delta.ini <- c(0.5, 1) ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, m.ini, x2.ini, delta.ini) Res1 <- fitovate(MBriereE, x=x1, y=y1, ini.val=ini.val, par.list=FALSE, fig.opt=TRUE, angle=pi/6, control=list(reltol=1e-20, maxit=20000), np=2000, unit=NULL) Res1$RSS x0.ini <- min( x1 ) y0.ini <- min( y1 ) theta.ini <- pi/4 len.max <- max( max(y1)-min(y1), max(x1)-min(x1) ) *2/sqrt(2) yc.ini <- len.max/3 xc.ini <- 1/4*len.max x2.ini <- len.max delta.ini <- c(0.5, seq(1, 5, by=5)) ini.val <- list(x0.ini, y0.ini, theta.ini, yc.ini, xc.ini, x2.ini, delta.ini) Res2 <- fitovate( MbetaE, x=x1, y=y1, ini.val=ini.val, par.list=TRUE, fig.opt=TRUE, angle=pi/3, control=list(reltol=1e-20, maxit=20000), np=2000, unit=NULL ) Res2$RSS Res3 <- fitovate( MLRFE, x=x1, y=y1, ini.val=ini.val, unit=NULL, par.list=FALSE, fig.opt=TRUE, angle=NULL, control=list(reltol=1e-20, maxit=20000), np=2000) Res3$RSS x0.ini <- min( x1 ) y0.ini <- min( y1 ) theta.ini <- pi/4 len.max <- max( max(y1)-min(y1), max(x1)-min(x1) ) *2/sqrt(2) c.ini <- 1/5*len.max K1.ini <- c(0.1, 1, 5, 10) K2.ini <- 1 x2.ini <- len.max a.ini <- 1 b.ini <- 1 ini.val <- list(x0.ini, y0.ini, theta.ini, c.ini, K1.ini, K2.ini, x2.ini, a.ini, b.ini) Res4 <- fitovate( MPerformanceE, x=x1, y=y1, ini.val=ini.val, par.list=TRUE, fig.opt=TRUE, index.xmax=4, angle=pi/3, control=list(reltol=1e-20, maxit=20000), np=2000, unit=NULL ) Res4$RSS graphics.off()
fitsigmoid
is used to estimate the parameters of a sigmoid growth equation based on the integral of
a performance equation or one of its simplified versions.
fitsigmoid(expr, x, y, ini.val, simpver = 1, control = list(), par.list = FALSE, fig.opt = FALSE, xlim = NULL, ylim = NULL, xlab = NULL, ylab = NULL, main = NULL, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
fitsigmoid(expr, x, y, ini.val, simpver = 1, control = list(), par.list = FALSE, fig.opt = FALSE, xlim = NULL, ylim = NULL, xlab = NULL, ylab = NULL, main = NULL, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
expr |
a performance equation or one of its simplified versions that is used to build a sigmoid growth equation. |
x |
the observed investigation times. |
y |
the observed |
ini.val |
the initial values of the model parameters. |
simpver |
an optional argument to use the simplified version of the performance equation. |
control |
the list of control parameters for using the |
par.list |
the option of showing the list of parameters on the screen. |
fig.opt |
an optional argument to draw the observations and the predicted sigmoid curve. |
xlim |
the range of the |
ylim |
the range of the |
xlab |
the label of the |
ylab |
the label of the |
main |
the main title of the figure. |
subdivisions |
please see the arguments for the |
rel.tol |
please see the arguments for the |
abs.tol |
please see the arguments for the |
stop.on.error |
please see the arguments for the |
keep.xy |
please see the arguments for the |
aux |
please see the arguments for the |
Here, ini.val
only includes the initial values of the model parameters as a list.
The Nelder-Mead algorithm (Nelder and Mead, 1965) is used to carry out the optimization of minimizing the residual
sum of squares (RSS) between the observed and predicted values. The
optim
function in package stats was used to carry out the Nelder-Mead algorithm.
The performance equations denote MbetaE
, MBriereE
,
MLRFE
, MPerformanceE
and their simplified versions.
The arguments of P
and simpver
should correspond
to expr
(i.e., MbetaE
or MBriereE
or MLRFE
or MPerformanceE
).
The sigmoid equation is the integral of a performance equation or one of its simplified versions.
par |
the estimates of the model parameters. |
r.sq |
the coefficient of determination between the observed and predicted |
RSS |
the residual sum of squares between the observed and predicted |
sample.size |
the number of data points used in the data fitting. |
x |
the observed |
y |
the observed |
y.pred |
the predicted |
Here, the user can define other performance equations, but new equations or
their simplified versions should include the lower and upper thresholds on
the -axis corresponding to
, whose indices should
be the same as those in
MbetaE
or MBriereE
or MLRFE
or MPerformanceE
.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Jin, J., Quinn, B.K., Shi, P. (2022) The modified Brière equation and its applications. Plants 11, 1769. doi:10.3390/plants11131769
Lian, M., Shi, P., Zhang, L., Yao, W., Gielis, J., Niklas, K.J. (2023) A generalized performance equation
and its application in measuring the Gini index of leaf size inequality.
Trees Structure and Function 37, 1555
1565. doi:10.1007/s00468-023-02448-8
Nelder, J.A., Mead, R. (1965) A simplex method for function minimization.
Computer Journal 7, 308313. doi:10.1093/comjnl/7.4.308
Shi, P., Fan, M., Ratkowsky, D.A., Huang, J., Wu, H., Chen, L., Fang, S.,
Zhang, C. (2017) Comparison of two ontogenetic growth equations for animals and plants.
Ecological Modelling 349, 110. doi:10.1016/j.ecolmodel.2017.01.012
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
areaovate
, MbetaE
, MBriereE
, MLRFE
,
MPerformanceE
, sigmoid
# The shrimp growth data(See the supplementary table in West et al., 2001) # West, G.B., Brown, J.H., Enquist, B.J. (2001) A general model for ontogenetic growth. # Nature 413, 628-631. t0 <- c(3, 60, 90, 120, 150, 180, 384) m0 <- c(0.001, 0.005, 0.018, 0.037, 0.06, 0.067, 0.07) dev.new() plot( t0, m0, cex.lab=1.5, cex.axis=1.5, col=4, xlab=expression(italic(x)), ylab=expression(italic(y)) ) xopt0 <- seq(100, 150, by=5) ini.val <- list(0.035, xopt0, 200, 1) resu1 <- fitsigmoid(MLRFE, x=t0, y=m0, ini.val=ini.val, simpver=1, fig.opt=TRUE, par.list=TRUE) xopt0 <- seq(100, 150, by=5) ini.val <- list(0.035, xopt0, 200, 1) resu1 <- fitsigmoid(MbetaE, x=t0, y=m0, ini.val=ini.val, simpver=1, fig.opt=TRUE) m.ini <- c(0.5, 1, 2, 3, 4, 5, 10, 20) ini.val <- list(1e-8, m.ini, 200, 1) resu2 <- fitsigmoid(MBriereE, x=t0, y=m0, ini.val=ini.val, simpver=1, fig.opt=TRUE, control=list(reltol=1e-20, maxit=20000, trace=FALSE), subdivisions=100L, rel.tol=.Machine$double.eps^0.25, abs.tol=.Machine$double.eps^0.25, stop.on.error=TRUE, keep.xy=FALSE, aux=NULL) delta0 <- c(0.5, 1, 2, 5, 10, 20) ini.val <- list(0.035, 150, -100, 200, delta0) resu3 <- fitsigmoid(MLRFE, x=t0, y=m0, ini.val=ini.val, simpver=NULL, fig.opt=TRUE, control=list(reltol=1e-20, maxit=20000), subdivisions = 100L, rel.tol=.Machine$double.eps^0.25, abs.tol=.Machine$double.eps^0.25, stop.on.error=TRUE, keep.xy=FALSE, aux=NULL) a.ini <- c(0.1, 1, 10, 100, 200) b.ini <- 200 ini.val <- list(0.001, 0.02, 0.15, 0, 200, a.ini, b.ini) resu4 <- fitsigmoid(MPerformanceE, x=t0, y=m0, ini.val=ini.val, simpver=NULL, fig.opt=TRUE, control=list(reltol=1e-20, maxit=20000, trace=FALSE), subdivisions=100L, rel.tol=.Machine$double.eps^0.25, abs.tol=.Machine$double.eps^0.25, stop.on.error=TRUE, keep.xy=FALSE, aux=NULL) resu5 <- fitsigmoid(MPerformanceE, x=t0, y=m0, ini.val=resu4$par, simpver=NULL, fig.opt=TRUE, control=list(reltol=1e-30, maxit=200000, trace=FALSE)) ini.val <- list(0.001, 0.01, c(0.1, 1, 10), 0, 200) resu6 <- fitsigmoid(MPerformanceE, x=t0, y=m0, ini.val=ini.val, simpver=2, fig.opt=TRUE, control=list(reltol=1e-20, maxit=20000, trace=FALSE), subdivisions=100L, rel.tol=.Machine$double.eps^0.25, abs.tol=.Machine$double.eps^0.25, stop.on.error=TRUE, keep.xy=FALSE, aux=NULL) graphics.off()
# The shrimp growth data(See the supplementary table in West et al., 2001) # West, G.B., Brown, J.H., Enquist, B.J. (2001) A general model for ontogenetic growth. # Nature 413, 628-631. t0 <- c(3, 60, 90, 120, 150, 180, 384) m0 <- c(0.001, 0.005, 0.018, 0.037, 0.06, 0.067, 0.07) dev.new() plot( t0, m0, cex.lab=1.5, cex.axis=1.5, col=4, xlab=expression(italic(x)), ylab=expression(italic(y)) ) xopt0 <- seq(100, 150, by=5) ini.val <- list(0.035, xopt0, 200, 1) resu1 <- fitsigmoid(MLRFE, x=t0, y=m0, ini.val=ini.val, simpver=1, fig.opt=TRUE, par.list=TRUE) xopt0 <- seq(100, 150, by=5) ini.val <- list(0.035, xopt0, 200, 1) resu1 <- fitsigmoid(MbetaE, x=t0, y=m0, ini.val=ini.val, simpver=1, fig.opt=TRUE) m.ini <- c(0.5, 1, 2, 3, 4, 5, 10, 20) ini.val <- list(1e-8, m.ini, 200, 1) resu2 <- fitsigmoid(MBriereE, x=t0, y=m0, ini.val=ini.val, simpver=1, fig.opt=TRUE, control=list(reltol=1e-20, maxit=20000, trace=FALSE), subdivisions=100L, rel.tol=.Machine$double.eps^0.25, abs.tol=.Machine$double.eps^0.25, stop.on.error=TRUE, keep.xy=FALSE, aux=NULL) delta0 <- c(0.5, 1, 2, 5, 10, 20) ini.val <- list(0.035, 150, -100, 200, delta0) resu3 <- fitsigmoid(MLRFE, x=t0, y=m0, ini.val=ini.val, simpver=NULL, fig.opt=TRUE, control=list(reltol=1e-20, maxit=20000), subdivisions = 100L, rel.tol=.Machine$double.eps^0.25, abs.tol=.Machine$double.eps^0.25, stop.on.error=TRUE, keep.xy=FALSE, aux=NULL) a.ini <- c(0.1, 1, 10, 100, 200) b.ini <- 200 ini.val <- list(0.001, 0.02, 0.15, 0, 200, a.ini, b.ini) resu4 <- fitsigmoid(MPerformanceE, x=t0, y=m0, ini.val=ini.val, simpver=NULL, fig.opt=TRUE, control=list(reltol=1e-20, maxit=20000, trace=FALSE), subdivisions=100L, rel.tol=.Machine$double.eps^0.25, abs.tol=.Machine$double.eps^0.25, stop.on.error=TRUE, keep.xy=FALSE, aux=NULL) resu5 <- fitsigmoid(MPerformanceE, x=t0, y=m0, ini.val=resu4$par, simpver=NULL, fig.opt=TRUE, control=list(reltol=1e-30, maxit=200000, trace=FALSE)) ini.val <- list(0.001, 0.01, c(0.1, 1, 10), 0, 200) resu6 <- fitsigmoid(MPerformanceE, x=t0, y=m0, ini.val=ini.val, simpver=2, fig.opt=TRUE, control=list(reltol=1e-20, maxit=20000, trace=FALSE), subdivisions=100L, rel.tol=.Machine$double.eps^0.25, abs.tol=.Machine$double.eps^0.25, stop.on.error=TRUE, keep.xy=FALSE, aux=NULL) graphics.off()
fitSuper
is used to estimate the parameters of the superellipse equation.
fitSuper(x, y, ini.val, control = list(), par.list = FALSE, stand.fig = TRUE, angle = NULL, fig.opt = FALSE, np = 2000, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
fitSuper(x, y, ini.val, control = list(), par.list = FALSE, stand.fig = TRUE, angle = NULL, fig.opt = FALSE, np = 2000, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
x |
the |
y |
the |
ini.val |
the list of initial values for the model parameters. |
control |
the list of control parameters for using the |
par.list |
the option of showing the list of parameters on the screen. |
stand.fig |
the option of drawing the observed and predicted polygons at the standard state
(i.e., the polar point is located at (0, 0), and the major axis overlaps with the |
angle |
the angle between the major axis and the |
fig.opt |
an optional argument of drawing the observed and predicted polygons at arbitrary angle
between the major axis and the |
np |
the number of data points on the predicted superellipse curve. |
xlim |
the range of the |
ylim |
the range of the |
unit |
the unit of the |
main |
the main title of the figure. |
The superellipse equation has its mathematical expression:
where represents the polar radius at the polar angle
.
represents the semi-major axis of the superellipse;
is the ratio of the semi-minor axis to the semi-major axis of the superellipse; and
determines
the curvature of the superellipse curve. This function is basically equal to the
fitGE
function with the arguments m = 4
and simpver = 9
.
However, this function can make the estimated value of the parameter be smaller than or equal to 1.
Apart from the above three model parameters, there are three additional location parameters, i.e., the
planar coordinates of the polar point (
and
), and the angle between the major axis
of the superellipse and the
-axis (
). The input order of
ini.val
is ,
,
,
,
, and
. The fitted parameters will be shown after running this function
in the same order. The Nelder-Mead algorithm (Nelder and Mead, 1965)
is used to carry out the optimization of minimizing the residual sum of squares (RSS) between the observed and predicted radii.
The
optim
function in package stats was used to carry out the Nelder-Mead algorithm.
When angle = NULL
, the observed polygon will be shown at its initial angle in the scanned image;
when angle
is a numerical value (e.g., ) defined by the user, it indicates that the major axis
is rotated by the amount (
) counterclockwise from the
-axis.
par |
the estimated values of the parameters in the superellipse equation. |
scan.length |
the observed length of the polygon. |
scan.width |
the observed width of the polygon. |
scan.area |
the observed area of the polygon. |
r.sq |
the coefficient of determination between the observed and predicted polar radii. |
RSS |
the residual sum of squares between the observed and predicted polar radii. |
sample.size |
the number of data points used in the data fitting. |
phi.stand.obs |
the polar angles at the standard state. |
phi.trans |
the transferred polar angles rotated as defined by the user. |
r.stand.obs |
the observed polar radii at the standard state. |
r.stand.pred |
the predicted polar radii at the standard state. |
x.stand.obs |
the observed |
x.stand.pred |
the predicted |
y.stand.obs |
the observed |
y.stand.pred |
the predicted |
r.obs |
the observed polar radii at the transferred polar angles as defined by the user. |
r.pred |
the predicted polar radii at the transferred polar angles as defined by the user. |
x.obs |
the observed |
x.pred |
the predicted |
y.obs |
the observed |
y.pred |
the predicted |
The output values of running this function can be used as those of running the fitGE
function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Huang, J., Hui, C., Grissino-Mayer, H.D., Tardif, J., Zhai, L., Wang, F., Li, B. (2015) Capturing spiral radial growth of conifers using the superellipse to model tree-ring geometric shape. Frontiers in Plant Science 6, 856. doi:10.3389/fpls.2015.00856
data(whitespruce) uni.C <- sort( unique(whitespruce$Code) ) Data <- whitespruce[whitespruce$Code==uni.C[12], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=200, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y plot( x1, y1, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", col="grey73", lwd=2, xlab=expression(italic("x")), ylab=expression(italic("y")) ) x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- c(0, pi/4, pi/2) a.ini <- mean(c(max(x1)-min(x1), max(y1)-min(y1)))/2 k.ini <- 1 n.ini <- c(1.5, 2, 2.5) ini.val <- list( x0.ini, y0.ini, theta.ini, a.ini, k.ini, n.ini ) Res2 <- fitSuper(x=x1, y=y1, ini.val=ini.val, unit="cm", par.list=FALSE, stand.fig=FALSE, angle=NULL, fig.opt=TRUE, control=list(reltol=1e-20, maxit=20000), np=2000) Res2$par Res2$r.sq graphics.off()
data(whitespruce) uni.C <- sort( unique(whitespruce$Code) ) Data <- whitespruce[whitespruce$Code==uni.C[12], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=200, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y plot( x1, y1, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", col="grey73", lwd=2, xlab=expression(italic("x")), ylab=expression(italic("y")) ) x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- c(0, pi/4, pi/2) a.ini <- mean(c(max(x1)-min(x1), max(y1)-min(y1)))/2 k.ini <- 1 n.ini <- c(1.5, 2, 2.5) ini.val <- list( x0.ini, y0.ini, theta.ini, a.ini, k.ini, n.ini ) Res2 <- fitSuper(x=x1, y=y1, ini.val=ini.val, unit="cm", par.list=FALSE, stand.fig=FALSE, angle=NULL, fig.opt=TRUE, control=list(reltol=1e-20, maxit=20000), np=2000) Res2$par Res2$r.sq graphics.off()
fracdim
is used to calculate the fractal dimension of leaf veins
based on the box-counting method.
fracdim(x, y, frac.fig = TRUE, denomi.range = seq(8, 30, by=1), ratiox = 0.02, ratioy = 0.08, main = NULL)
fracdim(x, y, frac.fig = TRUE, denomi.range = seq(8, 30, by=1), ratiox = 0.02, ratioy = 0.08, main = NULL)
x |
the |
y |
the |
frac.fig |
the option of drawing the results of the linear fitting. |
denomi.range |
the number of equidistant segments of the maximum range
between the range of the |
ratiox |
the the |
ratioy |
the the |
main |
the main title of the figure. |
The box-counting approach uses a group of boxes (squares for simplicity) with different
sizes () to divide the leaf vein image into different parts. Let
represent the number
of boxes that include at least one pixel of leaf vein.
The maximum of the range of the
coordinates and the range of the
coordinates
for leaf-vein pixels is defined as
. Let
represent the vector of
/
denomi.range
. Then, we used the following equation to calculate the fractal
dimension of leaf veins:
where is the theoretical value of the fractal dimension. We can use its estimate as the
numerical value of the fractal dimension for a leaf venation network.
a |
the estimate of the intercept. |
sd.a |
the standard deviation of the estimated intercept. |
lci.a |
the lower bound of the 95% confidence interval of the estimated intercept. |
uci.a |
the upper bound of the 95% confidence interval of the estimated intercept. |
b |
the estimate of the slope. |
sd.b |
the standard deviation of the estimated slope. |
lci.a |
the lower bound of the 95% confidence interval of the estimated slope. |
uci.a |
the upper bound of the 95% confidence interval of the estimated slope. |
r.sq |
the coefficient of determination. |
delta |
the vector of box sizes. |
N |
the number of boxes that include at least one pixel of leaf vein. |
Here, x
and y
cannot be adjusted by the adjdata
function
because the leaf veins are not the leaf's boundary data.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Yu, K., Niinemets, Ü., Gielis, J. (2021) Can leaf shape be represented by the ratio of leaf width to length? Evidence from nine species of Magnolia and Michelia (Magnoliaceae). Forests 12, 41. doi:10.3390/f12010041
Vico, P.G., Kyriacos, S., Heymans, O., Louryan, S., Cartilier, L. (1998)
Dynamic study of the extraembryonic vascular network of the
chick embryo by fractal analysis. Journal of Theoretical Biology 195, 525532.
doi:10.1006/jtbi.1998.0810
data(veins) dev.new() plot(veins$x, veins$y, cex=0.01, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y"))) fracdim(veins$x, veins$y) graphics.off()
data(veins) dev.new() plot(veins$x, veins$y, cex=0.01, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y"))) fracdim(veins$x, veins$y) graphics.off()
GE
is used to calculate polar radii of the original Gielis equation
or one of its simplified versions at given polar angles.
GE(P, phi, m = 1, simpver = NULL, nval = 1)
GE(P, phi, m = 1, simpver = NULL, nval = 1)
P |
the parameters of the original Gielis equation or one of its simplified versions. |
phi |
the polar angle(s). |
m |
the given |
simpver |
an optional argument to use the simplified version of the original Gielis equation. |
nval |
the specified value for |
When simpver = NULL
, the original Gielis equation is selected:
where represents the polar radius at the polar angle
;
determines the number of angles within
; and
,
,
,
, and
need to be provided in
P
.
When
simpver = 1
, the simplified version 1 is selected:
where ,
, and
need to be provided in
P
.
When
simpver = 2
, the simplified version 2 is selected:
where and
need to be provided in
P
, and should be specified in
nval
.
When
simpver = 3
, the simplified version 3 is selected:
where needs to be provided in
P
, and
should be specified in
nval
.
When
simpver = 4
, the simplified version 4 is selected:
where and
need to be provided in
P
.
When
simpver = 5
, the simplified version 5 is selected:
where ,
,
, and
need to be provided in
P
.
When
simpver = 6
, the simplified version 6 is selected:
where ,
,
, and
need to be provided in
P
.
When
simpver = 7
, the simplified version 7 is selected:
where ,
, and
need to be provided in
P
, and
should be specified in
nval
.
When
simpver = 8
, the simplified version 8 is selected:
where and
are parameters that need to be provided in
P
, and
should be specified in
nval
.
When
simpver = 9
, the simplified version 9 is selected:
where ,
, and
need to be provided in
P
.
The polar radii predicted by the original Gielis equation or one of its simplified versions.
simpver
here is different from that in the TGE
function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Gielis, J. (2003) A generic geometric transformation that unifies a wide range of natural
and abstract shapes. American Journal of Botany 90, 333338. doi:10.3732/ajb.90.3.333
Li, Y., Quinn, B.K., Gielis, J., Li, Y., Shi, P. (2022) Evidence that supertriangles exist in nature from the vertical projections of Koelreuteria paniculata fruit. Symmetry 14, 23. doi:10.3390/sym14010023
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Ratkowsky, D.A., Gielis, J. (2020) The generalized Gielis geometric equation and its application. Symmetry 12, 645. doi:10.3390/sym12040645
Shi, P., Xu, Q., Sandhu, H.S., Gielis, J., Ding, Y., Li, H., Dong, X. (2015) Comparison
of dwarf bamboos (Indocalamus sp.) leaf parameters to determine relationship
between spatial density of plants and total leaf area per plant. Ecology and Evolution 5,
45784589. doi:10.1002/ece3.1728
areaGE
, curveGE
, DSGE
, fitGE
,
SurfaceAreaSGE
, TGE
, VolumeSGE
GE.par <- c(2, 1, 4, 6, 3) varphi.vec <- seq(0, 2*pi, len=2000) r.theor <- GE(P=GE.par, phi=varphi.vec, m=5) dev.new() plot( varphi.vec, r.theor, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(varphi)), ylab=expression(italic("r")), type="l", col=4 ) graphics.off()
GE.par <- c(2, 1, 4, 6, 3) varphi.vec <- seq(0, 2*pi, len=2000) r.theor <- GE(P=GE.par, phi=varphi.vec, m=5) dev.new() plot( varphi.vec, r.theor, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(varphi)), ylab=expression(italic("r")), type="l", col=4 ) graphics.off()
The data consist of the boundary data of four side projections of G. biloba (Cultivar 'Fozhi') seeds sampled at Nanjing Forestry University campus on September 23, 2021.
data(ginkgoseed)
data(ginkgoseed)
In the data set, there are three columns of variables: Code
, x
, and y
.
Code
saves the codes of individual fruit;
x
saves the coordinates of the side projections of seeds in the Cartesian coordinate system (cm);
and
y
saves the coordinates of the side projections of seeds in the Cartesian coordinate system (cm).
Tian, F., Wang, Y., Sandhu, H.S., Gielis, J., Shi, P. (2020) Comparison of seed morphology of two ginkgo cultivars.
Journal of Forestry Research 31, 751758. doi:10.1007/s11676-018-0770-y
data(ginkgoseed) uni.C <- sort( unique(ginkgoseed$Code) ) ind <- 1 Data <- ginkgoseed[ginkgoseed$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=2000, len.pro=1/20) dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic("x")), ylab=expression(italic("y")) ) x1 <- Res1$x y1 <- Res1$y x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- pi/4 a.ini <- 1 n1.ini <- seq(0.6, 1, by=0.1) n2.ini <- 1 n3.ini <- 1 ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, n1.ini, n2.ini, n3.ini) Res2 <- fitGE( GE, x=x1, y=y1, ini.val=ini.val, m=2, simpver=5, nval=1, unit="cm", par.list=FALSE, fig.opt=TRUE, angle=NULL, control=list(reltol=1e-20, maxit=20000), np=2000 ) graphics.off()
data(ginkgoseed) uni.C <- sort( unique(ginkgoseed$Code) ) ind <- 1 Data <- ginkgoseed[ginkgoseed$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=2000, len.pro=1/20) dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic("x")), ylab=expression(italic("y")) ) x1 <- Res1$x y1 <- Res1$y x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- pi/4 a.ini <- 1 n1.ini <- seq(0.6, 1, by=0.1) n2.ini <- 1 n3.ini <- 1 ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, n1.ini, n2.ini, n3.ini) Res2 <- fitGE( GE, x=x1, y=y1, ini.val=ini.val, m=2, simpver=5, nval=1, unit="cm", par.list=FALSE, fig.opt=TRUE, angle=NULL, control=list(reltol=1e-20, maxit=20000), np=2000 ) graphics.off()
The data consist of the boundary data of four vertical projections of K. paniculata fruit sampled at Nanjing Forestry University campus in early October 2021.
data(kp)
data(kp)
In the data set, there are three columns of variables: Code
, x
, and y
.
Code
saves the codes of individual fruit;
x
saves the coordinates of the vertical projections of fruit in the Cartesian coordinate system (cm);
and
y
saves the coordinates of the vertical projections of fruit in the Cartesian coordinate system (cm).
Li, Y., Quinn, B.K., Gielis, J., Li, Y., Shi, P. (2022) Evidence that supertriangles exist in nature from the vertical projections of Koelreuteria paniculata fruit. Symmetry 14, 23. doi:10.3390/sym14010023
data(kp) uni.C <- sort( unique(kp$Code) ) ind <- 1 Data <- kp[kp$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=200, len.pro=1/20) dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic("x")), ylab=expression(italic("y")) ) x1 <- Res1$x y1 <- Res1$y x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- pi a.ini <- 0.9 n1.ini <- c(1, 4) n2.ini <- 5 n3.ini <- c(5, 10, 15) ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, n1.ini, n2.ini, n3.ini) Res2 <- fitGE( GE, x=x1, y=y1, ini.val=ini.val, m=3, simpver=5, nval=1, unit="cm", par.list=FALSE, fig.opt=TRUE, angle=NULL, control=list(reltol=1e-20, maxit=20000), np=2000 ) graphics.off()
data(kp) uni.C <- sort( unique(kp$Code) ) ind <- 1 Data <- kp[kp$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=200, len.pro=1/20) dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic("x")), ylab=expression(italic("y")) ) x1 <- Res1$x y1 <- Res1$y x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- pi a.ini <- 0.9 n1.ini <- c(1, 4) n2.ini <- 5 n3.ini <- c(5, 10, 15) ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, n1.ini, n2.ini, n3.ini) Res2 <- fitGE( GE, x=x1, y=y1, ini.val=ini.val, m=3, simpver=5, nval=1, unit="cm", par.list=FALSE, fig.opt=TRUE, angle=NULL, control=list(reltol=1e-20, maxit=20000), np=2000 ) graphics.off()
The data consist of the leaf size measures of 12 individuad culms of Shibataea chinensis.
data(LeafSizeDist)
data(LeafSizeDist)
In the data set, there are five columns of variables: Code
, L
,
W
, A
, and M
. Code
saves the bamboo number (the number before the hyphen)
and the leaf number (the number after the hyphen) on each bamboo. L
saves the length of each
leaf lamina (cm); W
saves the width of each leaf lamina (cm);
A
saves the area of each leaf lamina (); and
M
saves
the dry mass of each leaf lamina (g).
Lian, M., Shi, P., Zhang, L., Yao, W., Gielis, J., Niklas, K.J. (2023) A generalized performance equation
and its application in measuring the Gini index of leaf size inequality.
Trees Structure and Function 37, 1555
1565. doi:10.1007/s00468-023-02448-8
data(LeafSizeDist) CulmNumber <- c() for(i in 1:length(LeafSizeDist$Code)){ temp <- as.numeric( strsplit(LeafSizeDist$Code[i],"-", fixed=TRUE)[[1]][1] ) CulmNumber <- c(CulmNumber, temp) } uni.CN <- sort( unique(CulmNumber) ) ind <- CulmNumber==uni.CN[1] A0 <- LeafSizeDist$A[ind] A1 <- sort( A0 ) x1 <- 1:length(A1)/length(A1) y1 <- cumsum(A1)/sum(A1) x1 <- c(0, x1) y1 <- c(0, y1) M0 <- LeafSizeDist$M[ind] M1 <- sort( M0 ) x2 <- 1:length(M1)/length(M1) y2 <- cumsum(M1)/sum(M1) x2 <- c(0, x2) y2 <- c(0, y2) dev.new() plot(x1, y1, cex.lab=1.5, cex.axis=1.5, type="l", xlab="Cumulative proportion of the number of leaves", ylab="Cumulative proportion of leaf area") lines(c(0,1), c(0,1), col=4) dev.new() plot(x2, y2, cex.lab=1.5, cex.axis=1.5, type="l", xlab="Cumulative proportion of the number of leaves", ylab="Cumulative proportion of leaf dry mass") lines(c(0,1), c(0,1), col=4) graphics.off()
data(LeafSizeDist) CulmNumber <- c() for(i in 1:length(LeafSizeDist$Code)){ temp <- as.numeric( strsplit(LeafSizeDist$Code[i],"-", fixed=TRUE)[[1]][1] ) CulmNumber <- c(CulmNumber, temp) } uni.CN <- sort( unique(CulmNumber) ) ind <- CulmNumber==uni.CN[1] A0 <- LeafSizeDist$A[ind] A1 <- sort( A0 ) x1 <- 1:length(A1)/length(A1) y1 <- cumsum(A1)/sum(A1) x1 <- c(0, x1) y1 <- c(0, y1) M0 <- LeafSizeDist$M[ind] M1 <- sort( M0 ) x2 <- 1:length(M1)/length(M1) y2 <- cumsum(M1)/sum(M1) x2 <- c(0, x2) y2 <- c(0, y2) dev.new() plot(x1, y1, cex.lab=1.5, cex.axis=1.5, type="l", xlab="Cumulative proportion of the number of leaves", ylab="Cumulative proportion of leaf area") lines(c(0,1), c(0,1), col=4) dev.new() plot(x2, y2, cex.lab=1.5, cex.axis=1.5, type="l", xlab="Cumulative proportion of the number of leaves", ylab="Cumulative proportion of leaf dry mass") lines(c(0,1), c(0,1), col=4) graphics.off()
lmPE
is used to estimate the parameters of the Todd-Smart equation using the multiple linear regression.
lmPE(x, y, simpver = NULL, dev.angle = NULL, weights = NULL, fig.opt = TRUE, prog.opt = TRUE, xlim = NULL, ylim = NULL, unit = NULL, main = NULL, extr.method = "Shi")
lmPE(x, y, simpver = NULL, dev.angle = NULL, weights = NULL, fig.opt = TRUE, prog.opt = TRUE, xlim = NULL, ylim = NULL, unit = NULL, main = NULL, extr.method = "Shi")
x |
the |
y |
the |
simpver |
an optional argument to use the simplified version of the original Todd-Smart equation. |
dev.angle |
the angle of deviation for the axis associated with the maximum distance between two points on an egg's profile from the mid-line of the egg's profile. |
weights |
the weights for the multiple linear regression. |
fig.opt |
an optional argument to draw the observed and predicted egg's boundaries. |
prog.opt |
an optional argument to show the running progress for different values of |
xlim |
the range of the |
ylim |
the range of the |
unit |
the units of the |
main |
the main title of the figure. |
extr.method |
an optional argument to fit the planar coordinate data of an egg's profile extracted using different methods. |
There are two methods to obtain the major axis (i.e., the mid-line) of an egg's profile: the maximum distance method, and
the longest axis adjustment method. In the first method, the straight line through two points forming the maximum distance
on the egg's boundary is defined as the major axis. In the second method, we assume that there is an angle
of deviation for the longest axis (i.e., the axis associated with the maximum distance between two points on an egg's profile)
from the mid-line of the egg's profile. That is to say,
the mid-line of an egg's profile is not the axis associated with the maximum distance between two points on the egg's profile.
When dev.angle = NULL
, the maximum distance method is used; when dev.angle
is a numerical value or a numerical vector,
the longest axis adjustment method is used. Here, the numerical value of dev.angle
is not
the angle of deviation for the major axis of an egg's profile from the -axis, and instead it is the angle of deviation
for the longest axis (associated with the maximum distance between two points on the egg's profile) from the mid-line of the egg's profile.
It is better to take the option of
extr.method = "Shi"
for correctly fitting the planar coordinate data of an egg's profile extracted
using the protocols proposed by Shi et al. (2015, 2018) (and also see Su et al. (2019)),
while it is better to take the option of extr.method = "Biggins"
for correctly fitting the planar coordinate data
of an egg's profile extracted using the protocols proposed by Biggins et al. (2018). For the planar coordinate data extracted using
the protocols of Biggins et al. (2018), there are fewer data points on the two ends of the mid-line than other parts of the egg's profile,
which means that the range of the observed values might be smaller than the actual egg's length. A group of equidistant
values are set along the mid-line, and each
value corresponds to two
values that are respectively located at the upper
and lower sides of the egg's profile. Because of the difference in the curvature for differnt parts of the egg's profile,
the equidistant
values cannot render the extracted data points on the egg's profile to be regular. For the planar coordinate data
extracted using the protocols of Shi et al. (2015, 2018), the data points are more regularly distributed on the egg's profile (perimeter)
than those of Biggins et al. (2018), although the
values of the data points along the mid-line are not equidistant.
lm.tse |
the fitted results of the multiple linear regression. |
par |
the estimates of the four model parameters in the Todd-Smart equation. |
theta |
the angle between the longest axis of an egg's profile (i.e.,
the axis associated with the maximum distance between two points on the egg's profile) and the |
epsilon |
the optimal angle of deviation for the longest axis (associated with the maximum distance between two points on an egg's profile)
from the mid-line of the egg's profile, when |
RSS.vector |
the vector of residual sum of squares corresponding to |
x.obs |
the observed |
y.obs |
the observed |
y.pred |
the predicted |
x.stand.obs |
the observed |
y.stand.obs |
the observed |
y.stand.pred |
the predicted |
scan.length |
the length of the egg's boundary. The default is the maximum distance between two points on the egg's boundary. |
scan.width |
the maximum width of the egg's boundary. |
scan.area |
the area of the egg's boundary. |
scan.perimeter |
the perimeter of the egg's boundary based on all data points on the egg's boundary. |
RSS.scaled |
the residual sum of squares between the observed and predicted |
RSS |
the residual sum of squares between the observed and predicted |
sample.size |
the number of data points used in the numerical calculation. |
RMSE.scaled |
the root-mean-square errors between the observed and predicted |
RMSE |
the root-mean-square errors between the observed and predicted |
theta
is the calculated angle between the longest axis (i.e., the axis associated with the maximum distance
between two points on an egg's profile) and the -axis, and
epsilon
is the calculated angle of deviation for the longest
axis from the mid-line of the egg's profile. This means that the angle between the mid-line and the -axis is equal to
theta
+ epsilon
.
Here, RSS
, and RMSE
are for the observed and predicted coordinates of the egg's profile, not for those
when the egg's length is scaled to 2. There are two figures when
fig.opt = TRUE
: (i) the observed
and predicted egg's boundaries when the egg's length is scaled to 2, and (ii) the observed
and predicted egg's boundaries at their actual scales.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Biggins, J.D., Montgomeries, R.M., Thompson, J.E., Birkhead, T.R. (2022) Preston’s universal formula for avian egg shape. Ornithology 139, ukac028. doi:10.1093/ornithology/ukac028
Biggins, J.D., Thompson, J.E., Birkhead, T.R. (2018) Accurately quantifying
the shape of birds' eggs. Ecology and Evolution 8, 97289738. doi:10.1002/ece3.4412
Nelder, J.A., Mead, R. (1965). A simplex method for function minimization.
Computer Journal 7, 308313. doi:10.1093/comjnl/7.4.308
Preston, F.W. (1953) The shapes of birds' eggs. The Auk 70, 160182.
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Huang, J., Hui, C., Grissino-Mayer, H.D., Tardif, J., Zhai, L., Wang, F., Li, B. (2015) Capturing spiral radial growth of conifers using the superellipse to model tree-ring geometric shape. Frontiers in Plant Science 6, 856. doi:10.3389/fpls.2015.00856
Shi, P., Niinemets, Ü., Hui, C., Niklas, K.J., Yu, X., Hölscher, D. (2020) Leaf bilateral symmetry and the scaling of the perimeter vs. the surface area in 15 vine species. Forests 11, 246. doi:10.3390/f11020246
Shi, P., Ratkowsky, D.A., Li, Y., Zhang, L., Lin, S., Gielis, J. (2018) General leaf-area geometric formula exists for plants - Evidence from the simplified Gielis equation. Forests 9, 714. doi:10.3390/f9110714
Shi, P., Wang, L., Quinn, B.K., Gielis, J. (2023) A new program to estimate the parameters of Preston's equation, a general formula for describing the egg shape of birds. Symmetry 15, 231. doi:10.3390/sym15010231
Su, J., Niklas, K.J., Huang, W., Yu, X., Yang, Y., Shi, P. (2019) Lamina shape does not correlate with lamina surface area: An analysis based on the simplified Gielis equation. Global Ecology and Conservation 19, e00666. doi:10.1016/j.gecco.2019.e00666
Todd, P.H., Smart, I.H.M. (1984) The shape of birds' eggs. Journal of Theoretical Biology
106, 239243. doi:10.1016/0022-5193(84)90021-3
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 8 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=3000, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) Res2 <- lmPE(x1, y1, simpver=NULL, dev.angle=NULL, unit="cm") summary( Res2$lm.tse ) Res2$RMSE.scaled / 2 if(FALSE){ dev.new() xg1 <- seq(-1, 1, len=1000) yg1 <- TSE(P=Res2$par, x=xg1, simpver=NULL) xg2 <- seq(1, -1, len=1000) yg2 <- -TSE(P=Res2$par, x=xg2, simpver=NULL) plot(xg1, yg1, asp=1, type="l", col=2, ylim=c(-1,1), cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y))) lines(xg2, yg2, col=4) dev.new() plot(Res2$x.obs, Res2$y.obs, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y)), type="l") lines(Res2$x.obs, Res2$y.pred, col=2) dev.new() plot(Res2$x.stand.obs, Res2$y.stand.obs, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y)), type="l") lines(Res2$x.stand.obs, Res2$y.stand.pred, col=2) } Res3 <- lmPE(x1, y1, simpver=NULL, dev.angle=seq(-0.05, 0.05, by=0.0001), unit="cm") summary( Res3$lm.tse ) Res3$epsilon Res3$RMSE.scaled / 2 Res4 <- lmPE(x1, y1, simpver=1, dev.angle=NULL, unit="cm") summary( Res4$lm.tse ) graphics.off()
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 8 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=3000, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) Res2 <- lmPE(x1, y1, simpver=NULL, dev.angle=NULL, unit="cm") summary( Res2$lm.tse ) Res2$RMSE.scaled / 2 if(FALSE){ dev.new() xg1 <- seq(-1, 1, len=1000) yg1 <- TSE(P=Res2$par, x=xg1, simpver=NULL) xg2 <- seq(1, -1, len=1000) yg2 <- -TSE(P=Res2$par, x=xg2, simpver=NULL) plot(xg1, yg1, asp=1, type="l", col=2, ylim=c(-1,1), cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y))) lines(xg2, yg2, col=4) dev.new() plot(Res2$x.obs, Res2$y.obs, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y)), type="l") lines(Res2$x.obs, Res2$y.pred, col=2) dev.new() plot(Res2$x.stand.obs, Res2$y.stand.obs, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y)), type="l") lines(Res2$x.stand.obs, Res2$y.stand.pred, col=2) } Res3 <- lmPE(x1, y1, simpver=NULL, dev.angle=seq(-0.05, 0.05, by=0.0001), unit="cm") summary( Res3$lm.tse ) Res3$epsilon Res3$RMSE.scaled / 2 Res4 <- lmPE(x1, y1, simpver=1, dev.angle=NULL, unit="cm") summary( Res4$lm.tse ) graphics.off()
lmTE
is used to estimate the parameters of the Troscianko equation using the multiple linear regression,
and the estimated values of the parameters are only used as the initial values for using the fitETE
function
lmTE(x, y, dev.angle = NULL, weights = NULL, fig.opt = TRUE, prog.opt = TRUE, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
lmTE(x, y, dev.angle = NULL, weights = NULL, fig.opt = TRUE, prog.opt = TRUE, xlim = NULL, ylim = NULL, unit = NULL, main = NULL)
x |
the |
y |
the |
dev.angle |
the angle of deviation for the axis associated with the maximum distance between two points on an egg's profile from the mid-line of the egg's profile. |
weights |
the weights for the multiple linear regression. |
fig.opt |
an optional argument to draw the observed and predicted egg's boundaries. |
prog.opt |
an optional argument to show the running progress for different values of |
xlim |
the range of the |
ylim |
the range of the |
unit |
the units of the |
main |
the main title of the figure. |
The estimated values of the parameters using the lmTE
function tend to be NOT globally optimal, and the values
are only used as the initial values for using the fitETE
function.
There are two methods to obtain the major axis (i.e., the mid-line) of an egg's profile: the maximum distance method, and
the longest axis adjustment method. In the first method, the straight line through two points forming the maximum distance
on the egg's boundary is defined as the major axis. In the second method, we assume that there is an angle
of deviation for the longest axis (i.e., the axis associated with the maximum distance between two points on an egg's profile)
from the mid-line of the egg's profile. That is to say,
the mid-line of an egg's profile is not the axis associated with the maximum distance between two points on the egg's profile.
When dev.angle = NULL
, the maximum distance method is used; when dev.angle
is a numerical value or a numerical vector,
the longest axis adjustment method is used. Here, the numerical value of dev.angle
is not
the angle of deviation for the major axis of an egg's profile from the -axis, and instead it is the angle of deviation
for the longest axis (associated with the maximum distance between two points on the egg's profile) from the mid-line of the egg's profile.
The planar coordinate data of an egg's profile are extracted using the protocols proposed by Shi et al. (2015, 2018)
(and also see Su et al. (2019)). For the planar coordinate data extracted using the protocols of Shi et al. (2015, 2018),
the data points are more regularly distributed on the egg's profile (perimeter),
although the
values of the data points along the mid-line are not equidistant.
lm.te |
the fitted results of the multiple linear regression. |
par |
the estimates of the four model parameters in the Troscianko equation. |
theta |
the angle between the longest axis of an egg's profile (i.e.,
the axis associated with the maximum distance between two points on the egg's profile) and the |
epsilon |
the optimal angle of deviation for the longest axis (associated with the maximum distance between two points on an egg's profile)
from the mid-line of the egg's profile, when |
RSS.vector |
the vector of residual sum of squares corresponding to |
x.obs |
the observed |
y.obs |
the observed |
y.pred |
the predicted |
x.stand.obs |
the observed |
y.stand.obs |
the observed |
y.stand.pred |
the predicted |
scan.length |
the length of the egg's boundary. The default is the maximum distance between two points on the egg's boundary. |
scan.width |
the maximum width of the egg's boundary. |
scan.area |
the area of the egg's boundary. |
scan.perimeter |
the perimeter of the egg's boundary based on all data points on the egg's boundary. |
RSS.scaled |
the residual sum of squares between the observed and predicted |
RSS |
the residual sum of squares between the observed and predicted |
sample.size |
the number of data points used in the numerical calculation. |
RMSE.scaled |
the root-mean-square errors between the observed and predicted |
RMSE |
the root-mean-square errors between the observed and predicted |
theta
is the calculated angle between the longest axis (i.e., the axis associated with the maximum distance
between two points on an egg's profile) and the -axis, and
epsilon
is the calculated angle of deviation for the longest
axis from the mid-line of the egg's profile. This means that the angle between the mid-line and the -axis is equal to
theta
+ epsilon
.
Here, RSS
, and RMSE
are for the observed and predicted coordinates of the egg's profile, not for those
when the egg's length is scaled to 2. There are two figures when
fig.opt = TRUE
: (i) the observed
and predicted egg's boundaries when the egg's length is scaled to 2, and (ii) the observed
and predicted egg's boundaries at their actual scales.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Biggins, J.D., Montgomeries, R.M., Thompson, J.E., Birkhead, T.R. (2022) Preston’s universal formula for avian egg shape. Ornithology 139, ukac028. doi:10.1093/ornithology/ukac028
Biggins, J.D., Thompson, J.E., Birkhead, T.R. (2018) Accurately quantifying
the shape of birds' eggs. Ecology and Evolution 8, 97289738. doi:10.1002/ece3.4412
Nelder, J.A., Mead, R. (1965). A simplex method for function minimization.
Computer Journal 7, 308313. doi:10.1093/comjnl/7.4.308
Shi, P., Huang, J., Hui, C., Grissino-Mayer, H.D., Tardif, J., Zhai, L., Wang, F., Li, B. (2015) Capturing spiral radial growth of conifers using the superellipse to model tree-ring geometric shape. Frontiers in Plant Science 6, 856. doi:10.3389/fpls.2015.00856
Shi, P., Ratkowsky, D.A., Li, Y., Zhang, L., Lin, S., Gielis, J. (2018) General leaf-area geometric formula exists for plants - Evidence from the simplified Gielis equation. Forests 9, 714. doi:10.3390/f9110714
Shi, P., Wang, L., Quinn, B.K., Gielis, J. (2023) A new program to estimate the parameters of Preston's equation, a general formula for describing the egg shape of birds. Symmetry 15, 231. doi:10.3390/sym15010231
Su, J., Niklas, K.J., Huang, W., Yu, X., Yang, Y., Shi, P. (2019) Lamina shape does not correlate with lamina surface area: An analysis based on the simplified Gielis equation. Global Ecology and Conservation 19, e00666. doi:10.1016/j.gecco.2019.e00666
Troscianko, J. (2014). A simple tool for calculating egg shape, volume and surface area from digital images.
Ibis, 156, 874878. doi:10.1111/ibi.12177
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 8 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=3000, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) Res2 <- lmTE(x1, y1, dev.angle=NULL, unit="cm") summary( Res2$lm.te ) Res2$RMSE.scaled / 2 if(FALSE){ dev.new() xg1 <- seq(-1, 1, len=1000) yg1 <- TE(P=Res2$par, x=xg1) xg2 <- seq(1, -1, len=1000) yg2 <- -TE(P=Res2$par, x=xg2) plot(xg1, yg1, asp=1, type="l", col=2, ylim=c(-1,1), cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y))) lines(xg2, yg2, col=4) dev.new() plot(Res2$x.obs, Res2$y.obs, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y)), type="l") lines(Res2$x.obs, Res2$y.pred, col=2) dev.new() plot(Res2$x.stand.obs, Res2$y.stand.obs, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y)), type="l") lines(Res2$x.stand.obs, Res2$y.stand.pred, col=2) } Res3 <- lmTE(x1, y1, dev.angle=seq(-0.05, 0.05, by=0.0001), unit="cm") summary( Res3$lm.te ) Res3$epsilon Res3$RMSE.scaled / 2 graphics.off()
data(eggs) uni.C <- sort( unique(eggs$Code) ) ind <- 8 Data <- eggs[eggs$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=3000, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) Res2 <- lmTE(x1, y1, dev.angle=NULL, unit="cm") summary( Res2$lm.te ) Res2$RMSE.scaled / 2 if(FALSE){ dev.new() xg1 <- seq(-1, 1, len=1000) yg1 <- TE(P=Res2$par, x=xg1) xg2 <- seq(1, -1, len=1000) yg2 <- -TE(P=Res2$par, x=xg2) plot(xg1, yg1, asp=1, type="l", col=2, ylim=c(-1,1), cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y))) lines(xg2, yg2, col=4) dev.new() plot(Res2$x.obs, Res2$y.obs, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y)), type="l") lines(Res2$x.obs, Res2$y.pred, col=2) dev.new() plot(Res2$x.stand.obs, Res2$y.stand.obs, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y)), type="l") lines(Res2$x.stand.obs, Res2$y.stand.pred, col=2) } Res3 <- lmTE(x1, y1, dev.angle=seq(-0.05, 0.05, by=0.0001), unit="cm") summary( Res3$lm.te ) Res3$epsilon Res3$RMSE.scaled / 2 graphics.off()
MbetaE
is used to calculate values at given
values
using the modified beta equation or one of its simplified versions.
MbetaE(P, x, simpver = 1)
MbetaE(P, x, simpver = 1)
P |
the parameters of the modified beta equation or one of its simplified versions. |
x |
the given |
simpver |
an optional argument to use the simplified version of the modified beta equation. |
When simpver = NULL
, the modified beta equation is selected:
Here, and
represent the independent and dependent variables, respectively;
,
,
,
, and
are constants to be estimated;
represents the maximum
, and
is the
value associated with
the maximum
(i.e.,
);
and
and
represent the
lower and upper intersections between the curve and the
-axis.
is defined as 0
when
or
. There are five elements in
P
, representing
the values of ,
,
,
, and
, respectively.
When
simpver = 1
, the simplified version 1 is selected:
There are four elements in P
, representing
the values of ,
,
, and
, respectively.
When
simpver = 2
, the simplified version 2 is selected:
There are four elements in P
, representing
the values of ,
,
, and
, respectively.
When
simpver = 3
, the simplified version 3 is selected:
There are three elements in P
, representing
the values of ,
, and
, respectively.
The values predicted by the modified beta equation or one of its simplified versions.
We have added a parameter in the original beta equation (i.e.,
simpver = 2
) to increase the flexibility for data fitting.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Shi, P., Fan, M., Ratkowsky, D.A., Huang, J., Wu, H., Chen, L., Fang, S.,
Zhang, C. (2017) Comparison of two ontogenetic growth equations for animals and plants.
Ecological Modelling 349, 110. doi:10.1016/j.ecolmodel.2017.01.012
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
areaovate
, curveovate
, fitovate
, fitsigmoid
,
MBriereE
, MLRFE
, MPerformanceE
, sigmoid
x1 <- seq(-5, 15, len=2000) Par1 <- c(3, 3, 10, 2) y1 <- MbetaE(P=Par1, x=x1, simpver=1) dev.new() plot( x1, y1,cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic(x)), ylab=expression(italic(y)) ) graphics.off()
x1 <- seq(-5, 15, len=2000) Par1 <- c(3, 3, 10, 2) y1 <- MbetaE(P=Par1, x=x1, simpver=1) dev.new() plot( x1, y1,cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic(x)), ylab=expression(italic(y)) ) graphics.off()
MBriereE
is used to calculate values at given
values using
the modified Brière equation or one of its simplified versions.
MBriereE(P, x, simpver = 1)
MBriereE(P, x, simpver = 1)
P |
the parameters of the modified Brière equation or one of its simplified versions. |
x |
the given |
simpver |
an optional argument to use the simplified version of the modified Brière equation. |
When simpver = NULL
, the modified Brière equation is selected:
Here, and
represent the independent and dependent variables, respectively;
and
,
,
,
, and
are constants to be estimated,
where
and
represents the
lower and upper intersections between the curve and the
-axis.
is defined as 0
when
or
. There are five elements in
P
, representing
the values of ,
,
,
, and
, respectively.
When
simpver = 1
, the simplified version 1 is selected:
There are four elements in P
, representing
the values of ,
,
, and
, respectively.
When
simpver = 2
, the simplified version 2 is selected:
There are four elements in P
representing
the values of ,
,
, and
, respectively.
When
simpver = 3
, the simplified version 3 is selected:
There are three elements in P
representing
the values of ,
, and
, respectively.
The values predicted by the modified Brière equation or one of its simplified versions.
We have added a parameter in the original Brière equation (i.e.,
simpver = 2
) to increase the flexibility for data fitting.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Brière, J.-F., Pracros, P, Le Roux, A.-Y., Pierre, J.-S. (1999) A novel rate
model of temperature-dependent development for arthropods. Environmental
Entomology 28, 2229. doi:10.1093/ee/28.1.22
Cao, L., Shi, P., Li, L., Chen, G. (2019) A new flexible sigmoidal growth model. Symmetry 11, 204. doi:10.3390/sym11020204
Jin, J., Quinn, B.K., Shi, P. (2022) The modified Brière equation and its applications. Plants 11, 1769. doi:10.3390/plants11131769
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
areaovate
, curveovate
, fitovate
, fitsigmoid
,
MbetaE
, MLRFE
, MPerformanceE
, sigmoid
x2 <- seq(-5, 15, len=2000) Par2 <- c(0.01, 3, 0, 10, 1) y2 <- MBriereE(P=Par2, x=x2, simpver=NULL) dev.new() plot( x2, y2, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic(x)), ylab=expression(italic(y)) ) graphics.off()
x2 <- seq(-5, 15, len=2000) Par2 <- c(0.01, 3, 0, 10, 1) y2 <- MBriereE(P=Par2, x=x2, simpver=NULL) dev.new() plot( x2, y2, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic(x)), ylab=expression(italic(y)) ) graphics.off()
MLRFE
is used to calculate values at given
values
using the modified LRF equation or one of its simplified versions.
MLRFE(P, x, simpver = 1)
MLRFE(P, x, simpver = 1)
P |
the parameters of the modified LRF equation or one of its simplified versions. |
x |
the given |
simpver |
an optional argument to use the simplified version of the modified LRF equation. |
When simpver = NULL
, the modified LRF equation is selected:
Here, and
represent the independent and dependent variables, respectively;
,
,
,
, and
are constants to be estimated;
represents the maximum
, and
is the
value associated with
the maximum
(i.e.,
);
and
and
represents the
lower and upper intersections between the curve and the
-axis.
There are five elements in
P
, representing
the values of ,
,
,
, and
, respectively.
When
simpver = 1
, the simplified version 1 is selected:
There are four elements in P
, representing
the values of ,
,
, and
, respectively.
When
simpver = 2
, the simplified version 2 is selected:
There are four elements in P
, representing
the values of ,
,
, and
, respectively.
When
simpver = 3
, the simplified version 3 is selected:
There are three elements in P
, representing
the values of ,
, and
, respectively.
The values predicted by the modified LRF equation or one of its simplified versions.
We have added n parameter in the original LRF equation (i.e.,
simpver = 2
) to increase the flexibility for data fitting.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Shi, P., Fan, M., Ratkowsky, D.A., Huang, J., Wu, H., Chen, L., Fang, S.,
Zhang, C. (2017) Comparison of two ontogenetic growth equations for animals and plants.
Ecological Modelling 349, 110. doi:10.1016/j.ecolmodel.2017.01.012
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
areaovate
, curveovate
, fitovate
, fitsigmoid
,
MbetaE
, MBriereE
, MPerformanceE
, sigmoid
x3 <- seq(-5, 15, len=2000) Par3 <- c(3, 3, 10, 2) y3 <- MbetaE(P=Par3, x=x3, simpver=1) dev.new() plot( x3, y3, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic(x)), ylab=expression(italic(y)) ) graphics.off()
x3 <- seq(-5, 15, len=2000) Par3 <- c(3, 3, 10, 2) y3 <- MbetaE(P=Par3, x=x3, simpver=1) dev.new() plot( x3, y3, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic(x)), ylab=expression(italic(y)) ) graphics.off()
MPerformanceE
is used to calculate values at given
values using
the modified performance equation or one of its simplified versions.
MPerformanceE(P, x, simpver = 1)
MPerformanceE(P, x, simpver = 1)
P |
the parameters of the modified performance equation or one of its simplified versions. |
x |
the given |
simpver |
an optional argument to use the simplified version of the modified performance equation. |
When simpver = NULL
, the modified performance equation is selected:
Here, and
represent the independent and dependent variables, respectively;
and
,
,
,
,
,
,
and
are constants to be estimated,
where
and
represents the
lower and upper intersections between the curve and the
-axis.
is defined as 0
when
or
. There are seven elements in
P
, representing
the values of ,
,
,
,
,
, and
, respectively.
When
simpver = 1
, the simplified version 1 is selected:
There are six elements in P
, representing
the values of ,
,
,
,
, and
respectively.
When
simpver = 2
, the simplified version 2 is selected:
There are five elements in P
representing
the values of ,
,
,
, and
, respectively.
When
simpver = 3
, the simplified version 3 is selected:
There are four elements in P
representing
the values of ,
,
, and
, respectively.
When
simpver = 4
, the simplified version 4 is selected:
There are five elements in P
, representing
the values of ,
,
,
, and
, respectively.
When
simpver = 5
, the simplified version 5 is selected:
There are three elements in P
, representing
the values of ,
, and
, respectively.
The values predicted by the modified performance equation or one of its simplified versions.
We have added two parameters and
in the original performance
equation (i.e.,
simpver = 2
) to increase the flexibility for data fitting.
The cases of simpver = 4
and simpver = 5
are used to describe the rotated and right-shifted
Lorenz curve (see Lian et al. [2023] for details).
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Huey, R.B., Stevenson, R.D. (1979) Integrating thermal physiology and ecology of ectotherms:
a discussion of approaches. American Zoologist 19, 357366. doi:10.1093/icb/19.1.357
Lian, M., Shi, P., Zhang, L., Yao, W., Gielis, J., Niklas, K.J. (2023) A generalized performance equation
and its application in measuring the Gini index of leaf size inequality.
Trees Structure and Function 37, 1555
1565. doi:10.1007/s00468-023-02448-8
Shi, P., Ge, F., Sun, Y., Chen, C. (2011) A simple model for describing
the effect of temperature on insect developmental rate. Journal of Asia-Pacific Entomology
14, 1520. doi:10.1016/j.aspen.2010.11.008
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
areaovate
, curveovate
, fitLorenz
,
fitovate
, fitsigmoid
, MbetaE
,
MBriereE
, MLRFE
, sigmoid
x4 <- seq(0, 40, len=2000) Par4 <- c(0.117, 0.090, 0.255, 5, 35, 1, 1) y4 <- MPerformanceE(P=Par4, x=x4, simpver=NULL) dev.new() plot( x4, y4, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic(x)), ylab=expression(italic(y)) ) graphics.off()
x4 <- seq(0, 40, len=2000) Par4 <- c(0.117, 0.090, 0.255, 5, 35, 1, 1) y4 <- MPerformanceE(P=Par4, x=x4, simpver=NULL) dev.new() plot( x4, y4, cex.lab=1.5, cex.axis=1.5, type="l", xlab=expression(italic(x)), ylab=expression(italic(y)) ) graphics.off()
The data consist of the leaf boundary data of seven species of Neocinnamomum.
data(Neocinnamomum)
data(Neocinnamomum)
In the data set, there are four columns of variables: Code
, LatinName
, x
, and y
.
Code
saves the codes of individual leaves;
LatinName
saves the Latin names of the seven species of Neocinnamomum;
x
saves the coordinates of the leaf boundary in the Cartesian coordinate system (cm);
and
y
saves the coordinates of the leaf boundary in the Cartesian coordinate system (cm).
Shi, P., Fan, M., Ratkowsky, D.A., Huang, J., Wu, H., Chen, L., Fang, S.,
Zhang, C. (2017) Comparison of two ontogenetic growth equations for animals and plants.
Ecological Modelling 349, 110. doi:10.1016/j.ecolmodel.2017.01.012
Shi, P., Yu, K., Niklas, K.J., Schrader, J., Song, Y., Zhu, R., Li, Y., Wei, H., Ratkowsky, D.A. (2021) A general model for describing the ovate leaf shape. Symmetry, 13, 1524. doi:10.3390/sym13081524
data(Neocinnamomum) uni.C <- sort( unique(Neocinnamomum$Code) ) ind <- 2 Data <- Neocinnamomum[Neocinnamomum$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y length(x0) Res1 <- adjdata(x0, y0, ub.np=200, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y length(x1) dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) graphics.off()
data(Neocinnamomum) uni.C <- sort( unique(Neocinnamomum$Code) ) ind <- 2 Data <- Neocinnamomum[Neocinnamomum$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y length(x0) Res1 <- adjdata(x0, y0, ub.np=200, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y length(x1) dev.new() plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) graphics.off()
NRGE
is used to calculate values at given
values using
the Narushin-Romanov-Griffin equation (NRGE).
NRGE(P, x)
NRGE(P, x)
P |
the four parameters (i.e., |
x |
the given |
The Narushin-Romanov-Griffin equation (Narushin et al., 2021) has four parameters in total, among which three parameters have clear geometric meanings.
Here, is the Narushin-Romanov-Griffin equation, which is used to predict the
coordinates
at the given
coordinates;
represents the egg's length;
represents the egg's maximum breadth;
is a parameter
to be estimated, and it can be expressed as
, where
is a parameter to
be estimated;
represents the egg's breadth associated with
from the egg base (to the egg tip)
on the egg length axis (which can be regarded as the major axis of the egg shape).
The values predicted by the Narushin-Romanov-Griffin equation.
Here, parameter C
is a parameter to be estimated, which can be directly calculated numerically based on the egg-shape data.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Narushin, V.G., Romanov, M.N., Griffin, D.K. (2021) Egg and math: introducing a universal formula for egg shape.
Annals of the New York Academy of Sciences 1505, 169177. doi:10.1111/nyas.14680
Shi, P., Gielis, J., Niklas, K.J. (2022) Comparison of a universal (but complex) model for avian egg
shape with a simpler model. Annals of the New York Academy of Sciences 1514, 3442. doi:10.1111/nyas.14799
Tian, F., Wang, Y., Sandhu, H.S., Gielis, J., Shi, P. (2020) Comparison of seed morphology of two ginkgo cultivars.
Journal of Forestry Research 31, 751758. doi:10.1007/s11676-018-0770-y
curveNRGE
, fitNRGE
, SurfaceAreaNRGE
, VolumeNRGE
P0 <- c(11.5, 7.8, 1.1, 5.6) x <- seq(-11.5/2, 11.5/2, len=2000) y1 <- NRGE(P=P0, x=x) y2 <- -NRGE(P=P0, x=x) dev.new() plot(x, y1, cex.lab=1.5, cex.axis=1.5, type="l", col=4, ylim=c(-4, 4), asp=1, xlab=expression(italic(x)), ylab=expression(italic(y)) ) lines(x, y2, col=2) graphics.off()
P0 <- c(11.5, 7.8, 1.1, 5.6) x <- seq(-11.5/2, 11.5/2, len=2000) y1 <- NRGE(P=P0, x=x) y2 <- -NRGE(P=P0, x=x) dev.new() plot(x, y1, cex.lab=1.5, cex.axis=1.5, type="l", col=4, ylim=c(-4, 4), asp=1, xlab=expression(italic(x)), ylab=expression(italic(y)) ) lines(x, y2, col=2) graphics.off()
PE
is used to calculate the abscissa, ordinate and distance from the origin for an arbitrary point on
the Preston curve that was generated by the original Preston equation or one of its simplified versions at a given angle.
PE(P, zeta, simpver = NULL)
PE(P, zeta, simpver = NULL)
P |
the parameters of the original Preston equation or one of its simplified versions. |
zeta |
the angle(s) used in the Preston equation. |
simpver |
an optional argument to use the simplified version of the original Preston equation. |
When simpver = NULL
, the original Preston equation is selected:
where and
represent the abscissa and ordinate of an arbitrary point on the Preston curve
corresponding to an angle
;
represents the distance of the point from the origin;
,
,
,
, and
are parameters to be estimated.
When
simpver = 1
, the simplified version 1 is selected:
where and
represent the abscissa and ordinate of an arbitrary point on the Preston curve
corresponding to an angle
;
represents the distance of the point from the origin;
,
,
,
and
are parameters to be estimated.
When
simpver = 2
, the simplified version 2 is selected:
where and
represent the abscissa and ordinate of an arbitrary point on the Preston curve
corresponding to an angle
;
represents the distance of the point from the origin;
,
, and
are parameters to be estimated.
When
simpver = 3
, the simplified version 3 is selected:
where and
represent the abscissa and ordinate of an arbitrary point on the Preston curve
corresponding to an angle
;
represents the distance of the point from the origin;
,
, and
are parameters to be estimated.
x |
the abscissa(s) of the Preston curve corresponding to the given angle(s). |
y |
the ordinate(s) of the Preston curve corresponding to the given angle(s). |
r |
the distance(s) of the Preston curve corresponding to the given angle(s) from the origin. |
is NOT the polar angle corresponding to
, i.e.,
Let be the polar angle corresponding to
. We have:
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Biggins, J.D., Montgomeries, R.M., Thompson, J.E., Birkhead, T.R. (2022) Preston's universal formula for avian egg shape. Ornithology 139, ukac028. doi:10.1093/ornithology/ukac028
Biggins, J.D., Thompson, J.E., Birkhead, T.R. (2018) Accurately quantifying
the shape of birds' eggs. Ecology and Evolution 8, 97289738. doi:10.1002/ece3.4412
Preston, F.W. (1953) The shapes of birds' eggs. The Auk 70, 160182.
Shi, P., Wang, L., Quinn, B.K., Gielis, J. (2023) A new program to estimate the parameters of Preston's equation, a general formula for describing the egg shape of birds. Symmetry 15, 231. doi:10.3390/sym15010231
Todd, P.H., Smart, I.H.M. (1984) The shape of birds' eggs. Journal of Theoretical Biology
106, 239243. doi:10.1016/0022-5193(84)90021-3
zeta <- seq(0, 2*pi, len=2000) Par1 <- c(10, 6, 0.325, -0.0415) Res1 <- PE(P=Par1, zeta=zeta, simpver=1) Par2 <- c(10, 6, -0.325, -0.0415) Res2 <- PE(P=Par2, zeta=zeta, simpver=1) dev.new() plot(Res1$x, Res1$y, asp=1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y))) lines(Res2$x, Res2$y, col=2) dev.new() plot(Res1$r, Res2$r, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(paste(italic(r), ""[1], sep="")), ylab=expression(paste(italic(r), ""[2], sep=""))) abline(0, 1, col=4) graphics.off()
zeta <- seq(0, 2*pi, len=2000) Par1 <- c(10, 6, 0.325, -0.0415) Res1 <- PE(P=Par1, zeta=zeta, simpver=1) Par2 <- c(10, 6, -0.325, -0.0415) Res2 <- PE(P=Par2, zeta=zeta, simpver=1) dev.new() plot(Res1$x, Res1$y, asp=1, type="l", col=4, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y))) lines(Res2$x, Res2$y, col=2) dev.new() plot(Res1$r, Res2$r, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(paste(italic(r), ""[1], sep="")), ylab=expression(paste(italic(r), ""[2], sep=""))) abline(0, 1, col=4) graphics.off()
SarabiaE
is used to calculate values at given
values using
the Sarabia equation. The equation describes the
coordinates of the Lorenz curve.
SarabiaE(P, x)
SarabiaE(P, x)
P |
the parameters of the Sarabia equation. |
x |
the given |
Here, and
represent the independent and dependent variables, respectively;
and
,
,
and
are constants to be estimated, where
,
,
,
, and
.
There are four elements in
P
, representing
the values of ,
,
and
, respectively.
The values predicted by the Sarabia equation.
The numerical range of should range between 0 and 1.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Sarabia, J.-M. (1997) A hierarchy of Lorenz curves based on the generalized Tukey's lambda distribution.
Econometric Reviews 16, 305320. doi:10.1080/07474939708800389
Sitthiyot, T., Holasut, K. (2023) A universal model for the Lorenz curve with novel applications for datasets containing zeros and/or exhibiting extreme inequality. Scientific Reports 13, 4729. doi:10.1038/s41598-023-31827-x
fitLorenz
, MPerformanceE
, SCSE
, SHE
X1 <- seq(0, 1, len=2000) Pa1 <- c(0.295, 101.485, 0.705, 0.003762) Y1 <- SarabiaE(P=Pa1, x=X1) dev.new() plot( X1, Y1, cex.lab=1.5, cex.axis=1.5, type="l", asp=1, xaxs="i", yaxs="i", xlim=c(0, 1), ylim=c(0, 1), xlab="Cumulative proportion of the number of infructescences", ylab="Cumulative proportion of the infructescence length" ) graphics.off()
X1 <- seq(0, 1, len=2000) Pa1 <- c(0.295, 101.485, 0.705, 0.003762) Y1 <- SarabiaE(P=Pa1, x=X1) dev.new() plot( X1, Y1, cex.lab=1.5, cex.axis=1.5, type="l", asp=1, xaxs="i", yaxs="i", xlim=c(0, 1), ylim=c(0, 1), xlab="Cumulative proportion of the number of infructescences", ylab="Cumulative proportion of the infructescence length" ) graphics.off()
SCSE
is used to calculate values at given
values using
the Sarabia-Castillo-Slottje equation. The equation describes the
coordinates of the Lorenz curve.
SCSE(P, x)
SCSE(P, x)
P |
the parameters of the Sarabia-Castillo-Slottje equation. |
x |
the given |
Here, and
represent the independent and dependent variables, respectively;
and
,
and
are constants to be estimated, where
,
, and
.
There are three elements in
P
, representing
the values of ,
and
, respectively.
The values predicted by the Sarabia-Castillo-Slottje equation.
The numerical range of should range between 0 and 1.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Sarabia, J.-M., Castillo, E., Slottje, D.J. (1999) An ordered family of Lorenz curves.
Journal of Econometrics. 91, 4360. doi:10.1016/S0304-4076(98)00048-7
Sitthiyot, T., Holasut, K. (2023) A universal model for the Lorenz curve with novel applications for datasets containing zeros and/or exhibiting extreme inequality. Scientific Reports 13, 4729. doi:10.1038/s41598-023-31827-x
fitLorenz
, MPerformanceE
, SarabiaE
, SHE
X1 <- seq(0, 1, len=2000) Pa2 <- c(0, 0.790, 1.343) Y2 <- SCSE(P=Pa2, x=X1) dev.new() plot( X1, Y2, cex.lab=1.5, cex.axis=1.5, type="l", asp=1, xaxs="i", yaxs="i", xlim=c(0, 1), ylim=c(0, 1), xlab="Cumulative proportion of the number of infructescences", ylab="Cumulative proportion of the infructescence length" ) graphics.off()
X1 <- seq(0, 1, len=2000) Pa2 <- c(0, 0.790, 1.343) Y2 <- SCSE(P=Pa2, x=X1) dev.new() plot( X1, Y2, cex.lab=1.5, cex.axis=1.5, type="l", asp=1, xaxs="i", yaxs="i", xlim=c(0, 1), ylim=c(0, 1), xlab="Cumulative proportion of the number of infructescences", ylab="Cumulative proportion of the infructescence length" ) graphics.off()
SHE
is used to calculate values at given
values using
the Sitthiyot-Holasut equation. The equation describes the
coordinates of the Lorenz curve.
SHE(P, x)
SHE(P, x)
P |
the parameters of the Sitthiyot-Holasut equation. |
x |
the given |
Here, and
represent the independent and dependent variables, respectively;
and
,
,
and
are constants to be estimated, where
,
,
, and
.
There are four elements in
P
, representing
the values of ,
,
and
, respectively.
The values predicted by the Sitthiyot-Holasut equation.
The numerical range of should range between 0 and 1.
When
, the
value is assigned to be
.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Sitthiyot, T., Holasut, K. (2023) A universal model for the Lorenz curve with novel applications for datasets containing zeros and/or exhibiting extreme inequality. Scientific Reports 13, 4729. doi:10.1038/s41598-023-31827-x
fitLorenz
, MPerformanceE
, SarabiaE
, SCSE
X1 <- seq(0, 1, len=2000) Pa3 <- c(0, 1, 0.446, 1.739) Y3 <- SHE(P=Pa3, x=X1) dev.new() plot( X1, Y3, cex.lab=1.5, cex.axis=1.5, type="l", asp=1, xaxs="i", yaxs="i", xlim=c(0, 1), ylim=c(0, 1), xlab="Cumulative proportion of the number of infructescences", ylab="Cumulative proportion of the infructescence length" ) graphics.off()
X1 <- seq(0, 1, len=2000) Pa3 <- c(0, 1, 0.446, 1.739) Y3 <- SHE(P=Pa3, x=X1) dev.new() plot( X1, Y3, cex.lab=1.5, cex.axis=1.5, type="l", asp=1, xaxs="i", yaxs="i", xlim=c(0, 1), ylim=c(0, 1), xlab="Cumulative proportion of the number of infructescences", ylab="Cumulative proportion of the infructescence length" ) graphics.off()
The height data of four species of bamboo at Nanjing Forestry University campus in 2016.
data(shoots)
data(shoots)
In the data set, there are four columns of variables: Code
, LatinName
, x
, and y
.
Code
saves the number codes of different bamboo species;
LatinName
saves the Latin names of different bamboo species;
x
saves the investigation times (days from a specific starting time of growth,
and where every bamboo has a different starting time of growth);
and y
saves the measured aboveground height values (cm).
Code = 1
represents Phyllostachys iridescens, and the starting time (namely time = 0
) was defined as 12:00, 3rd April, 2016;
Code = 2
represents Phyllostachys mannii, and the starting time (namely time = 0
) was defined as 12:00, 4th April, 2016;
Code = 3
represents Pleioblastus maculatus, and the starting time (namely time = 0
) was defined as 12:00, 29th April, 2016;
Code = 4
represents Sinobambusa tootsik, and the starting time (namely time = 0
) was defined as 12:00, 18th April, 2016.
Shi, P., Fan, M., Ratkowsky, D.A., Huang, J., Wu, H., Chen, L., Fang, S.,
Zhang, C. (2017) Comparison of two ontogenetic growth equations for animals and plants.
Ecological Modelling 349, 110. doi:10.1016/j.ecolmodel.2017.01.012
data(shoots) attach(shoots) # Choose a species # 1: Phyllostachys iridescens; 2: Phyllostachys mannii; # 3: Pleioblastus maculatus; 4: Sinobambusa tootsik ind <- 4 x1 <- x[Code == ind] y1 <- y[Code == ind] dev.new() plot(x1, y1, cex=1.5, cex.lab=1.5, cex.axis=1.5, xlab="Time (days)", ylab="Height (cm)") delta0 <- c(0.5, 1, 2, 5, 10, 20) ini.val <- list(600, 25, 0, 40, delta0) resu1 <- fitsigmoid(MLRFE, x=x1, y=y1, ini.val=ini.val, simpver=NULL, fig.opt=TRUE, control=list(reltol=1e-20, maxit=20000), subdivisions = 100L, rel.tol=.Machine$double.eps^0.25, abs.tol=.Machine$double.eps^0.25, stop.on.error=TRUE, keep.xy=FALSE, aux=NULL) graphics.off()
data(shoots) attach(shoots) # Choose a species # 1: Phyllostachys iridescens; 2: Phyllostachys mannii; # 3: Pleioblastus maculatus; 4: Sinobambusa tootsik ind <- 4 x1 <- x[Code == ind] y1 <- y[Code == ind] dev.new() plot(x1, y1, cex=1.5, cex.lab=1.5, cex.axis=1.5, xlab="Time (days)", ylab="Height (cm)") delta0 <- c(0.5, 1, 2, 5, 10, 20) ini.val <- list(600, 25, 0, 40, delta0) resu1 <- fitsigmoid(MLRFE, x=x1, y=y1, ini.val=ini.val, simpver=NULL, fig.opt=TRUE, control=list(reltol=1e-20, maxit=20000), subdivisions = 100L, rel.tol=.Machine$double.eps^0.25, abs.tol=.Machine$double.eps^0.25, stop.on.error=TRUE, keep.xy=FALSE, aux=NULL) graphics.off()
sigmoid
is used to calculate the values (e.g., biomass, height, body length, and so on) at given investigation times.
sigmoid(expr, P, x, simpver = 1, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
sigmoid(expr, P, x, simpver = 1, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
expr |
a performance equation or one of its simplified versions. |
P |
the parameters of the performance equation or one of its simplified versions. |
x |
the given investigation times. |
simpver |
an optional argument to use the simplfied version of the performance equation. |
subdivisions |
please see the arguments for the |
rel.tol |
please see the arguments for the |
abs.tol |
please see the arguments for the |
stop.on.error |
please see the arguments for the |
keep.xy |
please see the arguments for the |
aux |
please see the arguments for the |
The performance equations denote MbetaE
, MBriereE
,
MLRFE
, and their simplified versions.
The arguments of P
and simpver
should correspond
to expr
(i.e., MbetaE
, MBriereE
,
MLRFE
, and MPerformanceE
).
The sigmoid curve is the integral of the performance equation or one of its simplified versions.
The values (i.e., biomass, height, body length, and so on) at given investigation times.
The growth euqation is actually an integral of the performance equation or one of its simplified versions.
Here, the user can define other performance equations, but new equations or
their simplified versions should include the lower and upper thresholds in
the -axis corresponding to
, whose indices of the parameters in
P
should
be the same as those in MbetaE
, MBriereE
,
MLRFE
, and MPerformanceE
).
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Jin, J., Quinn, B.K., Shi, P. (2022) The modified Brière equation and its applications. Plants 11, 1769. doi:10.3390/plants11131769
Lian, M., Shi, P., Zhang, L., Yao, W., Gielis, J., Niklas, K.J. (2023) A generalized performance equation
and its application in measuring the Gini index of leaf size inequality.
Trees Structure and Function 37, 1555
1565. doi:10.1007/s00468-023-02448-8
Shi, P., Fan, M., Ratkowsky, D.A., Huang, J., Wu, H., Chen, L., Fang, S.,
Zhang, C. (2017) Comparison of two ontogenetic growth equations for animals and plants.
Ecological Modelling 349, 110. doi:10.1016/j.ecolmodel.2017.01.012
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
fitsigmoid
, MbetaE
, MBriereE
,
MLRFE
, and MPerformanceE
Pa1 <- c(3, 3, 10, 1) xv1 <- seq(-5, 15, len=2000) yv1 <- sigmoid(MBriereE, P=Pa1, x=xv1, simpver=1) Pa2 <- c(3, 3, 2, 12, 1) yv2 <- sigmoid(MBriereE, P=Pa2, x=xv1, simpver=NULL) Pa3 <- c(200, 1.32, 0.85, 12, 3.7, 17) yv3 <- sigmoid(MPerformanceE, P=Pa3, x=xv1, simpver=1) dev.new() plot( xv1, yv2, cex.lab=1.5, cex.axis=1.5, type="l", col=4, xlab=expression(italic(x)), ylab=expression(italic(y)) ) lines( xv1, yv1, col=2 ) lines( xv1, yv3, col=3 ) graphics.off()
Pa1 <- c(3, 3, 10, 1) xv1 <- seq(-5, 15, len=2000) yv1 <- sigmoid(MBriereE, P=Pa1, x=xv1, simpver=1) Pa2 <- c(3, 3, 2, 12, 1) yv2 <- sigmoid(MBriereE, P=Pa2, x=xv1, simpver=NULL) Pa3 <- c(200, 1.32, 0.85, 12, 3.7, 17) yv3 <- sigmoid(MPerformanceE, P=Pa3, x=xv1, simpver=1) dev.new() plot( xv1, yv2, cex.lab=1.5, cex.axis=1.5, type="l", col=4, xlab=expression(italic(x)), ylab=expression(italic(y)) ) lines( xv1, yv1, col=2 ) lines( xv1, yv3, col=3 ) graphics.off()
The data consist of the boundary data of eight sea stars from five species.
data(starfish)
data(starfish)
In the data set, there are four columns of variables: Code
, LatinName
, x
, and y
.
Code
saves the codes of individual sea stars;
LatinName
saves the Latin names of the eight sea stars;
x
saves the coordinates of the eight sea stars in the Cartesian coordinate system (cm);
and
y
saves the coordinates of the eight sea stars in the Cartesian coordinate system (cm).
In
Code
, codes 1-9 represent Anthenoides tenuis, Culcita schmideliana sample 1,
Culcita schmideliana sample 2, Culcita schmideliana sample 3, Stellaster equestris, Tosia australis,
Tosia magnifica sample 1, and Tosia magnifica sample 2, respectively. See Table A1 published in Shi et al. (2020).
Shi, P., Ratkowsky, D.A., Gielis, J. (2020) The generalized Gielis geometric equation and its application. Symmetry 12, 645. doi:10.3390/sym12040645
data(starfish) uni.C <- sort( unique(starfish$Code) ) ind <- 2 Data <- starfish[starfish$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y dev.new() plot( x0, y0, asp=1, type="l", cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) length(x0) Res1 <- adjdata(x0, y0, ub.np=400, times=1.2, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( x1, y1, asp=1, type="l", cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) length(x1) graphics.off()
data(starfish) uni.C <- sort( unique(starfish$Code) ) ind <- 2 Data <- starfish[starfish$Code==uni.C[ind], ] x0 <- Data$x y0 <- Data$y dev.new() plot( x0, y0, asp=1, type="l", cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) length(x0) Res1 <- adjdata(x0, y0, ub.np=400, times=1.2, len.pro=1/20) x1 <- Res1$x y1 <- Res1$y dev.new() plot( x1, y1, asp=1, type="l", cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y")) ) length(x1) graphics.off()
SurfaceAreaEPE
is used to calculate the surface area of an egg that follows the explicit Preston equation.
SurfaceAreaEPE(P, simpver = NULL, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
SurfaceAreaEPE(P, simpver = NULL, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
P |
the parameters of the explicit Preston equation or one of its simplified versions. |
simpver |
an optional argument to use the simplified version of the explicit Preston equation. |
subdivisions |
please see the arguments for the |
rel.tol |
please see the arguments for the |
abs.tol |
please see the arguments for the |
stop.on.error |
please see the arguments for the |
keep.xy |
please see the arguments for the |
aux |
please see the arguments for the |
The formula of the surface area () of an egg based on the explicit Preston equation or one of its simplified versions is:
where denotes the explicit Preston equation (i.e.,
EPE
), and
denotes half the egg's length.
When
simpver = NULL
, P
has five parameters: ,
,
,
, and
;
when
simpver = 1
, P
has four parameters: ,
,
, and
;
when
simpver = 2
, P
has three parameters: ,
, and
;
when
simpver = 3
, P
has three parameters: ,
, and
.
The argument P
in the SurfaceAreaEPE
function has the same parameters, as those in the
EPE
function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Narushin, V.G., Romanov, M.N., Mishra, B., Griffin, D.K. (2022) Mathematical progression of
avian egg shape with associated area and volume determinations.
Annals of the New York Academy of Sciences 1513, 6578. doi:10.1111/nyas.14771
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
Par4 <- c(4.27, 2.90, 0.0868, 0.0224, -0.0287) SurfaceAreaEPE(P = Par4, simpver = NULL)
Par4 <- c(4.27, 2.90, 0.0868, 0.0224, -0.0287) SurfaceAreaEPE(P = Par4, simpver = NULL)
SurfaceAreaETE
is used to calculate the surface area of an egg that follows the explicit Troscianko equation.
SurfaceAreaETE(P, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
SurfaceAreaETE(P, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
P |
the parameters of the explicit Troscianko equation. |
subdivisions |
please see the arguments for the |
rel.tol |
please see the arguments for the |
abs.tol |
please see the arguments for the |
stop.on.error |
please see the arguments for the |
keep.xy |
please see the arguments for the |
aux |
please see the arguments for the |
The formula of the surface area () of an egg based on the explicit Troscianko equation is:
where denotes the explicit Troscianko equation (i.e.,
ETE
), and
denotes half the egg's length.
The argument P
in the SurfaceAreaETE
function has the same parameters, as those in the
ETE
function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Narushin, V.G., Romanov, M.N., Mishra, B., Griffin, D.K. (2022) Mathematical progression of
avian egg shape with associated area and volume determinations.
Annals of the New York Academy of Sciences 1513, 6578. doi:10.1111/nyas.14771
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
Par5 <- c(2.25, -0.38, -0.29, -0.16) SurfaceAreaETE(P = Par5)
Par5 <- c(2.25, -0.38, -0.29, -0.16) SurfaceAreaETE(P = Par5)
SurfaceAreaNRGE
is used to calculate the surface area of an egg that follows the Narushin-Romanov-Griffin equation.
SurfaceAreaNRGE(P, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
SurfaceAreaNRGE(P, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
P |
the parameters of the Narushin-Romanov-Griffin equation. |
subdivisions |
please see the arguments for the |
rel.tol |
please see the arguments for the |
abs.tol |
please see the arguments for the |
stop.on.error |
please see the arguments for the |
keep.xy |
please see the arguments for the |
aux |
please see the arguments for the |
The formula of the surface area () of an egg based on the Narushin-Romanov-Griffin equation is:
where denotes the Narushin-Romanov-Griffin equation (i.e.,
NRGE
), and
denotes the egg's length, which is the first element in the parameter vector,
P
.
The argument P
in the SurfaceAreaNRGE
function has the same parameters, as those in the
NRGE
function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Narushin, V.G., Romanov, M.N., Griffin, D.K. (2021) Egg and math: introducing a universal formula for egg shape.
Annals of the New York Academy of Sciences 1505, 169177. doi:10.1111/nyas.14680
Narushin, V.G., Romanov, M.N., Mishra, B., Griffin, D.K. (2022) Mathematical progression of
avian egg shape with associated area and volume determinations.
Annals of the New York Academy of Sciences 1513, 6578. doi:10.1111/nyas.14771
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
curveNRGE
, DNRGE
, fitNRGE
, NRGE
, VolumeNRGE
Par6 <- c(4.51, 3.18, 0.1227, 2.2284) SurfaceAreaNRGE(P = Par6)
Par6 <- c(4.51, 3.18, 0.1227, 2.2284) SurfaceAreaNRGE(P = Par6)
SurfaceAreaSGE
is used to calculate the surface area of an egg that follows the simplified Gielis equation.
SurfaceAreaSGE(P, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
SurfaceAreaSGE(P, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
P |
the parameters of the simplified Gielis equation, including |
subdivisions |
please see the arguments for the |
rel.tol |
please see the arguments for the |
abs.tol |
please see the arguments for the |
stop.on.error |
please see the arguments for the |
keep.xy |
please see the arguments for the |
aux |
please see the arguments for the |
The formula of the surface area () of an egg based on the simplified Gielis equation is:
where the polar raidus () is the function of the polar angle (
):
namely the simplified Gielis equation (i.e., GE
) with arguments
simpver = 1
and m = 1
.
The argument P
in the SurfaceAreaSGE
function only has
the three parameters: ,
, and
.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Chen, Z. (2012) Volume and area of revolution under polar coordinate system.
Studies in College Mathematics 15(6), 911.
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
Par7 <- c(1.124, 14.86, 49.43) SurfaceAreaSGE(P = Par7)
Par7 <- c(1.124, 14.86, 49.43) SurfaceAreaSGE(P = Par7)
TE
is used to calculate values at given
values using
the re-expression of Troscianko's egg-shape equation, which was proposed by Biggins et al. (2018, 2022).
TE(P, x)
TE(P, x)
P |
the parameters of the Troscianko equation, including |
x |
the given |
The Troscianko equation is recommended as (Biggins et al., 2022):
where and
represent the abscissa and ordinate of an arbitrary point on the Troscianko curve;
,
, and
are parameters to be estimated.
The values predicted by the Troscianko equation.
Here, and
in the Troscianko equation are actually equal to
and
, respectively, in the explicit Troscianko equation, where
represents half the egg length (See
ETE
for details).
This means that the egg length is scaled to be 2,
and the maximum egg width is correspondingly adjusted to keep the same scale.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Biggins, J.D., Montgomeries, R.M., Thompson, J.E., Birkhead, T.R. (2022) Preston's universal formula for avian egg shape. Ornithology 139, ukac028. doi:10.1093/ornithology/ukac028
Biggins, J.D., Thompson, J.E., Birkhead, T.R. (2018) Accurately quantifying
the shape of birds' eggs. Ecology and Evolution 8, 97289738. doi:10.1002/ece3.4412
Nelder, J.A., Mead, R. (1965). A simplex method for function minimization.
Computer Journal 7, 308313. doi:10.1093/comjnl/7.4.308
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Wang, L., Quinn, B.K., Gielis, J. (2023) A new program to estimate the parameters of Preston's equation, a general formula for describing the egg shape of birds. Symmetry 15, 231. doi:10.3390/sym15010231
Troscianko, J. (2014). A simple tool for calculating egg shape, volume and surface area from digital images.
Ibis, 156, 874878. doi:10.1111/ibi.12177
Par <- c(-0.377, -0.29, -0.16) xb1 <- seq(-1, 1, len=20000) yb1 <- TE(P=Par, x=xb1) xb2 <- seq(1, -1, len=20000) yb2 <- -TE(P=Par, x=xb2) dev.new() plot(xb1, yb1, asp=1, type="l", col=2, ylim=c(-1, 1), cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y))) lines(xb2, yb2, col=4) graphics.off()
Par <- c(-0.377, -0.29, -0.16) xb1 <- seq(-1, 1, len=20000) yb1 <- TE(P=Par, x=xb1) xb2 <- seq(1, -1, len=20000) yb2 <- -TE(P=Par, x=xb2) dev.new() plot(xb1, yb1, asp=1, type="l", col=2, ylim=c(-1, 1), cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y))) lines(xb2, yb2, col=4) graphics.off()
TGE
is used to calculate the polar radii of the twin Gielis equation or
one of its simplified versions at given polar angles.
TGE(P, phi, m = 1, simpver = NULL, nval = 1)
TGE(P, phi, m = 1, simpver = NULL, nval = 1)
P |
the parameters of the twin Gielis equation or one of its simplified versions. |
phi |
the polar angle(s). |
m |
the given |
simpver |
an optional argument to use the simplified version of the twin Gielis equation. |
nval |
the specified value for |
The general form of the twin Gielis equation can be represented as follows:
where represents the polar radius of the twin Gielis curve at the polar angle
, and
represents the elementary polar radius at the polar angle
. There is a hyperbolic
link function to link their log-transformations, i.e.,
The first three elements of P
are ,
, and
, and the remaining element(s) of
P
are the parameters of the elementary polar function, i.e., .
See Shi et al. (2020) for details.
When
simpver = NULL
, the original twin Gielis equation is selected:
where represents the elementary polar radius at the polar angle
;
determines the number of angles of the twin Gielis curve within
;
and
,
, and
are the fourth to the sixth elements in
P
. In total, there are
six elements in P
.
When
simpver = 1
, the simplified version 1 is selected:
where is the fourth element in
P
. There are
four elements in total in P
.
When
simpver = 2
, the simplified version 2 is selected:
where should be specified in
nval
, and P
only includes three elements, i.e.,
,
, and
.
When
simpver = 3
, the simplified version 3 is selected:
where and
are the fourth and fifth elements in
P
. There are
five elements in total in P
.
When
simpver = 4
, the simplified version 4 is selected:
where and
are the fourth and fifth elelments in
P
. There are
five elements in total in P
.
When
simpver = 5
, the simplified version 5 is selected:
where is the fourth elelment in
P
. There are
four elements in total in P
. should be specified in
nval
.
The polar radii predicted by the twin Gielis equation or one of its simplified versions.
simpver
here is different from that in the GE
function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Li, Y., Quinn, B.K., Gielis, J., Li, Y., Shi, P. (2022) Evidence that supertriangles exist in nature from the vertical projections of Koelreuteria paniculata fruit. Symmetry 14, 23. doi:10.3390/sym14010023
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Shi, P., Ratkowsky, D.A., Gielis, J. (2020) The generalized Gielis geometric equation and its application. Symmetry 12, 645. doi:10.3390/sym12040645
TGE.par <- c(2.88, 0.65, 1.16, 139) varphi.vec <- seq(0, 2*pi, len=2000) r2.theor <- TGE(P=TGE.par, phi=varphi.vec, simpver=1, m=5) dev.new() plot( varphi.vec, r2.theor, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(varphi)), ylab=expression(italic("r")), type="l", col=4 ) starfish4 <- curveGE(TGE, P=c(0, 0, 0, TGE.par), simpver=1, m=5, fig.opt=TRUE) graphics.off()
TGE.par <- c(2.88, 0.65, 1.16, 139) varphi.vec <- seq(0, 2*pi, len=2000) r2.theor <- TGE(P=TGE.par, phi=varphi.vec, simpver=1, m=5) dev.new() plot( varphi.vec, r2.theor, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(varphi)), ylab=expression(italic("r")), type="l", col=4 ) starfish4 <- curveGE(TGE, P=c(0, 0, 0, TGE.par), simpver=1, m=5, fig.opt=TRUE) graphics.off()
TSE
is used to calculate values at given
values using
the Todd and Smart's re-expression of Preston's universal egg shape.
TSE(P, x, simpver = NULL)
TSE(P, x, simpver = NULL)
P |
the parameters of the original Todd-Smart equation or one of its simplified versions. |
x |
the given |
simpver |
an optional argument to use the simplified version of the original Todd-Smart equation. |
When simpver = NULL
, the original Preston equation is selected:
where
Here, and
represent the abscissa and ordinate of an arbitrary point on the Todd-Smart curve;
,
,
, and
are parameters to be estimated.
When
simpver = 1
, the simplified version 1 is selected:
where and
represent the abscissa and ordinate of an arbitrary point on the Todd-Smart curve;
,
, and
are parameters to be estimated.
When
simpver = 2
, the simplified version 2 is selected:
where and
represent the abscissa and ordinate of an arbitrary point on the Todd-Smart curve;
, and
are parameters to be estimated.
When
simpver = 3
, the simplified version 3 is selected:
where and
represent the abscissa and ordinate of an arbitrary point on the Todd-Smart curve;
, and
are parameters to be estimated.
The values predicted by the Todd-Smart equation.
Here, and
in the Todd-Smart equation are actually equal to
and
, respectively, in the Preston equation (See
PE
for details).
Since represents half the egg length, this means that the egg length is fixed to be 2,
and the maximum egg width is correspondingly adjusted to keep the same scale.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Biggins, J.D., Montgomeries, R.M., Thompson, J.E., Birkhead, T.R. (2022) Preston's universal formula for avian egg shape. Ornithology 139, ukac028. doi:10.1093/ornithology/ukac028
Biggins, J.D., Thompson, J.E., Birkhead, T.R. (2018) Accurately quantifying
the shape of birds' eggs. Ecology and Evolution 8, 97289738. doi:10.1002/ece3.4412
Nelder, J.A., Mead, R. (1965). A simplex method for function minimization.
Computer Journal 7, 308313. doi:10.1093/comjnl/7.4.308
Preston, F.W. (1953) The shapes of birds' eggs. The Auk 70, 160182.
Shi, P., Gielis, J., Quinn, B.K., Niklas, K.J., Ratkowsky, D.A., Schrader, J., Ruan, H.,
Wang, L., Niinemets, Ü. (2022) 'biogeom': An R package for simulating and fitting natural
shapes. Annals of the New York Academy of Sciences 1516, 123134. doi:10.1111/nyas.14862
Todd, P.H., Smart, I.H.M. (1984) The shape of birds' eggs. Journal of Theoretical Biology
106, 239243. doi:10.1016/0022-5193(84)90021-3
Par <- c(0.695320398, -0.210538656, -0.070373518, 0.116839895) xb1 <- seq(-1, 1, len=20000) yb1 <- TSE(P=Par, x=xb1) xb2 <- seq(1, -1, len=20000) yb2 <- -TSE(P=Par, x=xb2) dev.new() plot(xb1, yb1, asp=1, type="l", col=2, ylim=c(-1, 1), cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y))) lines(xb2, yb2, col=4) graphics.off()
Par <- c(0.695320398, -0.210538656, -0.070373518, 0.116839895) xb1 <- seq(-1, 1, len=20000) yb1 <- TSE(P=Par, x=xb1) xb2 <- seq(1, -1, len=20000) yb2 <- -TSE(P=Par, x=xb2) dev.new() plot(xb1, yb1, asp=1, type="l", col=2, ylim=c(-1, 1), cex.lab=1.5, cex.axis=1.5, xlab=expression(italic(x)), ylab=expression(italic(y))) lines(xb2, yb2, col=4) graphics.off()
The data consist of the leaf vein data of a leaf of M. compressa sampled at Nanjing Forestry University campus in late July 2019.
data(veins)
data(veins)
In the data set, there are two columns of variables: x
and y
.
x
saves the coordinates of the leaf veins in the Cartesian coordinate system (cm);
y
saves the coordinates of the leaf veins in the Cartesian coordinate system (cm).
The data cannot be adjusted by the adjdata
function.
Shi, P., Yu, K., Niinemets, Ü., Gielis, J. (2021) Can leaf shape be represented by the ratio of leaf width to length? Evidence from nine species of Magnolia and Michelia (Magnoliaceae). Forests 12, 41. doi:10.3390/f12010041
data(veins) dev.new() plot(veins$x, veins$y, cex=0.01, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y"))) graphics.off()
data(veins) dev.new() plot(veins$x, veins$y, cex=0.01, asp=1, cex.lab=1.5, cex.axis=1.5, xlab=expression(italic("x")), ylab=expression(italic("y"))) graphics.off()
VolumeEPE
is used to calculate the volume of an egg that follows the explicit Preston equation.
VolumeEPE(P, simpver = NULL)
VolumeEPE(P, simpver = NULL)
P |
the parameters of the explicit Preston equation or one of its simplified versions. |
simpver |
an optional argument to use the simplified version of the explicit Preston equation. |
When simpver = NULL
, the volume formula () of the explicit Preston equation is selected:
where P
has five parameters: ,
,
,
, and
.
When
simpver = 1
, the volume formula of the simplified version 1 is selected:
where P
has four parameters: ,
,
, and
.
When
simpver = 2
, the volume formula of the simplified version 2 is selected:
where P
has three parameters: ,
, and
.
When
simpver = 3
, the volume formula of the simplified version 3 is selected:
where P
has three parameters: ,
, and
.
The argument P
in the VolumeEPE
function has the same parameters, as those in the
EPE
function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Narushin, V.G., Romanov, M.N., Mishra, B., Griffin, D.K. (2022) Mathematical progression of
avian egg shape with associated area and volume determinations.
Annals of the New York Academy of Sciences 1513, 6578. doi:10.1111/nyas.14771
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
Par3 <- c(4.27, 2.90, 0.0868, 0.0224, -0.0287) VolumeEPE(P=Par3, simpver=NULL) # Test the case when simpver = NULL a <- Par3[1] b <- Par3[2] c1 <- Par3[3] c2 <- Par3[4] c3 <- Par3[5] pi*4/315*a*b^2*(105+21*c1^2+42*c2+9*c2^2+18*c1*c3+5*c3^2) myfun <- function(x){ pi*EPE(P=Par3, x=x, simpver=NULL)^2 } integrate(myfun, -4.27, 4.27)$value
Par3 <- c(4.27, 2.90, 0.0868, 0.0224, -0.0287) VolumeEPE(P=Par3, simpver=NULL) # Test the case when simpver = NULL a <- Par3[1] b <- Par3[2] c1 <- Par3[3] c2 <- Par3[4] c3 <- Par3[5] pi*4/315*a*b^2*(105+21*c1^2+42*c2+9*c2^2+18*c1*c3+5*c3^2) myfun <- function(x){ pi*EPE(P=Par3, x=x, simpver=NULL)^2 } integrate(myfun, -4.27, 4.27)$value
VolumeETE
is used to calculate the volume of an egg that follows the explicit Troscianko equation.
VolumeETE(P, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
VolumeETE(P, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
P |
the parameters of the explicit Troscianko equation. |
subdivisions |
please see the arguments for the |
rel.tol |
please see the arguments for the |
abs.tol |
please see the arguments for the |
stop.on.error |
please see the arguments for the |
keep.xy |
please see the arguments for the |
aux |
please see the arguments for the |
The formula of the volume () of an egg based on the explicit Troscianko equation is:
where denotes the explicit Troscianko equation (i.e.,
ETE
), and
denotes half the egg's length.
The argument P
in the VolumeETE
function has the same parameters, as those in the
ETE
function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Narushin, V.G., Romanov, M.N., Mishra, B., Griffin, D.K. (2022) Mathematical progression of
avian egg shape with associated area and volume determinations.
Annals of the New York Academy of Sciences 1513, 6578. doi:10.1111/nyas.14771
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
Par5 <- c(2.25, -0.38, -0.29, -0.16) VolumeETE(P=Par5) myfun <- function(x){ pi*ETE(P=Par5, x=x)^2 } integrate(myfun, -2.25, 2.25)$value
Par5 <- c(2.25, -0.38, -0.29, -0.16) VolumeETE(P=Par5) myfun <- function(x){ pi*ETE(P=Par5, x=x)^2 } integrate(myfun, -2.25, 2.25)$value
VolumeNRGE
is used to calculate the volume of an egg that follows the Narushin-Romanov-Griffin equation.
VolumeNRGE(P, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
VolumeNRGE(P, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
P |
the parameters of the Narushin-Romanov-Griffin equation. |
subdivisions |
please see the arguments for the |
rel.tol |
please see the arguments for the |
abs.tol |
please see the arguments for the |
stop.on.error |
please see the arguments for the |
keep.xy |
please see the arguments for the |
aux |
please see the arguments for the |
The formula of the volume () of an egg based on the Narushin-Romanov-Griffin equation is:
where denotes the Narushin-Romanov-Griffin equation (i.e.,
NRGE
), and
denotes the egg's length, which is the first element in the parameter vector,
P
.
The argument P
in the VolumeNRGE
function has the same parameters, as those in the
NRGE
function.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Narushin, V.G., Romanov, M.N., Griffin, D.K. (2021) Egg and math: introducing a universal formula for egg shape.
Annals of the New York Academy of Sciences 1505, 169177. doi:10.1111/nyas.14680
Narushin, V.G., Romanov, M.N., Mishra, B., Griffin, D.K. (2022) Mathematical progression of
avian egg shape with associated area and volume determinations.
Annals of the New York Academy of Sciences 1513, 6578. doi:10.1111/nyas.14771
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
curveNRGE
, fitNRGE
, NRGE
, SurfaceAreaNRGE
Par6 <- c(4.51, 3.18, 0.1227, 2.2284) VolumeNRGE(P=Par6) myfun <- function(x){ pi*NRGE(P=Par6, x=x)^2 } integrate(myfun, -4.51/2, 4.51/2)$value
Par6 <- c(4.51, 3.18, 0.1227, 2.2284) VolumeNRGE(P=Par6) myfun <- function(x){ pi*NRGE(P=Par6, x=x)^2 } integrate(myfun, -4.51/2, 4.51/2)$value
VolumeSGE
is used to calculate the volume of an egg that follows the simplified Gielis equation.
VolumeSGE(P, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
VolumeSGE(P, subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
P |
the parameters of the simplified Gielis equation, including |
subdivisions |
please see the arguments for the |
rel.tol |
please see the arguments for the |
abs.tol |
please see the arguments for the |
stop.on.error |
please see the arguments for the |
keep.xy |
please see the arguments for the |
aux |
please see the arguments for the |
The formula of the volume () of an egg based on the simplified Gielis equation is:
where the polar raidus () is the function of the polar angle (
):
namely the simplified Gielis equation (i.e., GE
) with arguments
simpver = 1
and m = 1
.
The argument P
in the VolumeSGE
function only has
the three parameters: ,
, and
.
Peijian Shi [email protected], Johan Gielis [email protected], Brady K. Quinn [email protected].
Chen, Z. (2012) Volume and area of revolution under polar coordinate system.
Studies in College Mathematics 15(6), 911.
Shi, P., Chen, L., Quinn, B.K., Yu, K., Miao, Q., Guo, X., Lian, M., Gielis, J., Niklas, K.J. (2023)
A simple way to calculate the volume and surface area of avian eggs.
Annals of the New York Academy of Sciences 1524, 118131. doi:10.1111/nyas.15000
Par7 <- c(1.124, 14.86, 49.43) VolumeSGE(P = Par7)
Par7 <- c(1.124, 14.86, 49.43) VolumeSGE(P = Par7)
The data consist of the planar coordinates of Picea glauca tree rings.
data(whitespruce)
data(whitespruce)
In the data set, there are three columns of variables: Code
, x
, and y
.
Code
saves the age codes of tree rings from the 2nd year to the 44th year;
x
saves the coordinates of the tree rings in the Cartesian coordinate system (cm);
and
y
saves the coordinates of the tree rings in the Cartesian coordinate system (cm).
Shi, P., Huang, J., Hui, C., Grissino-Mayer, H.D., Tardif, J., Zhai, L., Wang, F., Li, B. (2015) Capturing spiral radial growth of conifers using the superellipse to model tree-ring geometric shape. Frontiers in Plant Science 6, 856. doi:10.3389/fpls.2015.00856
data(whitespruce) uni.C <- sort( unique(whitespruce$Code) ) Data <- whitespruce[whitespruce$Code==uni.C[10], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=2000, len.pro=1/20) plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", xlim=c(3, 13), ylim=c(3, 13), col="grey73", lwd=2, xlab=expression(italic(x)), ylab=expression(italic(y)) ) uni.C <- sort( unique(whitespruce$Code) ) for(i in 1:length(uni.C)){ Data <- whitespruce[whitespruce$Code==uni.C[i], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=200, len.pro=1/10) if(i == 1){ plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", xlim=c(3, 13), ylim=c(3, 13), col=1, lwd=1, xlab=expression(italic(x)), ylab=expression(italic(y)) ) } if(i > 1) lines(Res1$x, Res1$y, col=1, lwd=1) } uni.C <- sort( unique(whitespruce$Code) ) uni.C <- uni.C[1:12] Length <- c() results <- data.frame(Code=c(), x0=c(), y0=c(), theta=c(), a=c(), k=c(), n1=c(), r.sq=c(), RSS=c(), N=c()) for(i in 1:length(uni.C)){ Data <- whitespruce[whitespruce$Code==uni.C[i], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=200, len.pro=1/10) x1 <- Res1$x y1 <- Res1$y x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- c(0, pi/4, pi/2) a.ini <- 0.9 k.ini <- 1 n1.ini <- c(1.5, 2, 2.5) ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, k.ini, n1.ini) print(paste("Progress: ", i, "/", length(uni.C), sep="")) H <- NULL try( H <- fitGE(GE, x=x1, y=y1, ini.val=ini.val, m=4, simpver=9, unit="cm", par.list=FALSE, stand.fig=FALSE, angle=NULL, fig.opt=FALSE, control=list(reltol=1e-20, maxit=20000), np=2000), silent=TRUE ) if(is.null(H)){ RE <- data.frame(Code=uni.C[i], x0=NA, y0=NA, theta=NA, a=NA, k=NA, n1=NA, r.sq=NA, RSS=NA, N=NA) } if(!is.null(H)){ RE <- data.frame(Code=uni.C[i], x0=H$par[1], y0=H$par[2], theta=H$par[3], a=H$par[4], k=H$par[5], n1=H$par[6], r.sq=H$r.sq, RSS=H$RSS, N=H$sample.size) Length <- c(Length, max(max(H$y.stand.pred)[1]-min(H$y.stand.pred)[1], max(H$x.stand.pred)[1]-min(H$x.stand.pred)[1])[1]) if(i == 1){ plot(H$x.obs, H$y.obs, asp=1, xlim=c(7.4, 8.6), ylim=c(7.4, 8.6), cex.lab=1.5, cex.axis=1.5, type="l", lwd=2, col="grey70", xlab=expression(italic(x)), ylab=expression(italic(y))) lines(H$x.pred, H$y.pred, col=2) } if(i > 1){ lines(H$x.obs, H$y.obs, lwd=2, col="grey70") lines(H$x.pred, H$y.pred, col=2) } } results <- rbind(results, RE) } # To adjust the estimates of partial parameters to ensure k <= 1 results2 <- results Ind <- results$k > 1 results2$theta[Ind] <- results$theta[Ind] + pi/2 results2$a[Ind] <- results$a[Ind] * results$k[Ind]^(1/results$n1[Ind]) results2$k[Ind] <- 1/results$k[Ind] results2 Length/2
data(whitespruce) uni.C <- sort( unique(whitespruce$Code) ) Data <- whitespruce[whitespruce$Code==uni.C[10], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=2000, len.pro=1/20) plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", xlim=c(3, 13), ylim=c(3, 13), col="grey73", lwd=2, xlab=expression(italic(x)), ylab=expression(italic(y)) ) uni.C <- sort( unique(whitespruce$Code) ) for(i in 1:length(uni.C)){ Data <- whitespruce[whitespruce$Code==uni.C[i], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=200, len.pro=1/10) if(i == 1){ plot( Res1$x, Res1$y, asp=1, cex.lab=1.5, cex.axis=1.5, type="l", xlim=c(3, 13), ylim=c(3, 13), col=1, lwd=1, xlab=expression(italic(x)), ylab=expression(italic(y)) ) } if(i > 1) lines(Res1$x, Res1$y, col=1, lwd=1) } uni.C <- sort( unique(whitespruce$Code) ) uni.C <- uni.C[1:12] Length <- c() results <- data.frame(Code=c(), x0=c(), y0=c(), theta=c(), a=c(), k=c(), n1=c(), r.sq=c(), RSS=c(), N=c()) for(i in 1:length(uni.C)){ Data <- whitespruce[whitespruce$Code==uni.C[i], ] x0 <- Data$x y0 <- Data$y Res1 <- adjdata(x0, y0, ub.np=200, len.pro=1/10) x1 <- Res1$x y1 <- Res1$y x0.ini <- mean( x1 ) y0.ini <- mean( y1 ) theta.ini <- c(0, pi/4, pi/2) a.ini <- 0.9 k.ini <- 1 n1.ini <- c(1.5, 2, 2.5) ini.val <- list(x0.ini, y0.ini, theta.ini, a.ini, k.ini, n1.ini) print(paste("Progress: ", i, "/", length(uni.C), sep="")) H <- NULL try( H <- fitGE(GE, x=x1, y=y1, ini.val=ini.val, m=4, simpver=9, unit="cm", par.list=FALSE, stand.fig=FALSE, angle=NULL, fig.opt=FALSE, control=list(reltol=1e-20, maxit=20000), np=2000), silent=TRUE ) if(is.null(H)){ RE <- data.frame(Code=uni.C[i], x0=NA, y0=NA, theta=NA, a=NA, k=NA, n1=NA, r.sq=NA, RSS=NA, N=NA) } if(!is.null(H)){ RE <- data.frame(Code=uni.C[i], x0=H$par[1], y0=H$par[2], theta=H$par[3], a=H$par[4], k=H$par[5], n1=H$par[6], r.sq=H$r.sq, RSS=H$RSS, N=H$sample.size) Length <- c(Length, max(max(H$y.stand.pred)[1]-min(H$y.stand.pred)[1], max(H$x.stand.pred)[1]-min(H$x.stand.pred)[1])[1]) if(i == 1){ plot(H$x.obs, H$y.obs, asp=1, xlim=c(7.4, 8.6), ylim=c(7.4, 8.6), cex.lab=1.5, cex.axis=1.5, type="l", lwd=2, col="grey70", xlab=expression(italic(x)), ylab=expression(italic(y))) lines(H$x.pred, H$y.pred, col=2) } if(i > 1){ lines(H$x.obs, H$y.obs, lwd=2, col="grey70") lines(H$x.pred, H$y.pred, col=2) } } results <- rbind(results, RE) } # To adjust the estimates of partial parameters to ensure k <= 1 results2 <- results Ind <- results$k > 1 results2$theta[Ind] <- results$theta[Ind] + pi/2 results2$a[Ind] <- results$a[Ind] * results$k[Ind]^(1/results$n1[Ind]) results2$k[Ind] <- 1/results$k[Ind] results2 Length/2