General Information

GENERAL INFORMATION AND INSTRUCTIONS

Program Structure

rootchem.f
This is the main program. It reads in and stores the chemical data, reads in or sets control variables, and invokes the simplex optimization. Programs rootchem and calcfx communicate both via the argument list of calcfx and via common blocks.
calcfx.f
This subroutine is passed values of parameters (ap, ao and xo), and uses thoses values plus the data to calculate the objective function, a measure of the strength of the relationship between soil-corrected estimates of plant composition (xp) and sample ash (as). This subroutine is called repeatedly by smplx as smplx tries to find the best parameter (ap, ao and xo) values. When called directly by rootchem, calcfx writes the output files.
smplx.f
This subroutine estimates the values of parameters (ap, ao and xo) which minimize the objective function (see below) defined by calcfx. It communicates with rootchem and calcfx entirely via argument lists.

Compiling and Running the Program
These programs were developed and tested under the UNIX operating system on a SUN workstation. Hopefully, the FORTRAN code doesn't have any system specific features. Under UNIX, the procedure for running and executing is:

compile the program:
`f77 rootchem.f calcfx.f smplx.f`
rename the object file if you wish:
`mv a.out root.rr`
put the data into file 'data.in':
`cp data.pasm data.in`
execute the program:
`root.rr < input.pasm`
rename the output files:
`mv data.out data.out.pasm mv parest parest.pasm mv smplx.out smplx.out.pasm`

'input' File
The program has two main modes of operation. In the first mode, it treats all three parameters (ap, ao and xo) as independent. In the second mode, it reduces to two independent parameters by assuming a fixed relationship between ao and xo. In either the three or two parameter mode, each parameter may be either fixed or estimated by fitting to the data. The 'input' file contains both the names of the variables and their values, and the program checks that it is reading the correct variable. The format is not fixed, and the delimiter is a space. There are two possible formats for the input file.

First format:

line # 1 = 1)'npar', 2) '3' (number of parameters)
line # 2 = 1) 'ap', 2) toggle variable determining whether parameter ap is to be estimated by optimization (1) or taken as a given (0), 3) value of ap
line # 3 = 1) 'ao', 2) toggle variable for ao, 3) value of ao
line # 4 = 1) 'xo', 2) toggle variable for xo, 3) value of xo
line # 5 = 1) 'iopt', 2) toggle variable determining whether (1) or not (0) the simplex routine will be invoked.
line # 6 = 1) 'iprd', 2) toggle variable determining whether (1) or not(0) output will be printed by the simplex routine

Second format:

line # 1 = 1)'npar', 2) '2' (number of parameters)
line # 2 = as above
line # 3 = as above
line # 4 = 1) 'xorg', 2) value of xorg which relates ao and xo according to:
xo = (1. - ao)xorg
line # 5 = as above
line # 6 = as above

'data.in' File:

line # 1 = 1)'ngrps', 2) number of groups within which to regress xp vs. as
line # 2 = 1) group identifier, up to 8 characters and beginning with alphabetic character, 2) number of replicates within the group
line # 3 to line # (2 + number of replicates in group) = 1) value of as, 2) corresponding value of xs
Group identifiers and data for replicates for additional groups follow directly behind the first group. Groups with only one observation may be included, but will be not be used for parameter estimation. The format is not fixed, and the delimiter is a space. The program does not count replicates, so the data file must be set up accurately. If the maximum numbers of groups (50) or replicates per group (60) are exceeded, the program will give an error message, and you will have to increase the size of program arrays. There are two data files for test cases, one with the P. smithii* root data used for Fig. 2 in the Ecology paper, and another set up to illustrate the operation of the program with fewer than three replicates per group.

Output Files

'smplx.out'
Output from subroutine smplx is written to this file when variable 'iprd' on the 'input' file equals 1. Use this file to examine the progress of the SIMPLEX optimization. There are two examples, the first for the P. smithii data and the second for the other data set.

'data.out'
This file contains 1) the original data (as and xs), fp (fraction of plant material), xp (soil-corrected) and xa (ash-free) estimates for each replicate, and 2) statistics for regressions within each group. There are two examples, the first generated using the P. smithii data and the second input format, and the second using the second data set with the first example input file.

'parest'
This file contains the initial and final parameter values, and the value of the objective function. There are two examples, the first for the P. smithii data and the second for the other data set.

Objective Function
The function minimized in our original work was the average R² (coefficient of determination) for the regression of xp vs. as within groups. This worked fine because all groups had about the same number of replicates, and there were no groups with as few as two replicates. However, R² for the regression in a group with only two replicates is 1.0 unless the slope is 0.0, in which case R² is zero. Thus R² is not useful as an objective function when there are only two replicates in a group. The above program uses a new objective function-- the sum over groups of the three-way product of 1) the absolute value of the slope , 2) R², and 3) the number of replicates minus one. In groups with two replicates, this new objective function will be be small when the slope is small. In groups with three or more replicates, the objective function will be small if ether the slope or the R² is small. Including the number of replicates minus one in the product serves to give greater weight to groups with more replicates. In addition, penalties are applied to the objective function to avoid regions in parameter space which yield unrealistic estimates of fp (fraction plant material in the sample) or xp. Refer to the listing of 'calcfx' for details of these penalties. The new objective function yields similar parameter estimates for the P. smithii data included here as did the original objective function used for the Ecology paper.

Initial Parameter Estimates and Local Minima
Most optimization algorithms, including the SIMPLEX routine utilized here, locate local minima, not the global minimum. The minimum to which the SIMPLEX routine converges may depend on the initial values assigned to parameters and on the peculiarities of the data set. Thus it is important to run the program several times with different initial values for the parameters to make sure you have found an acceptable local minimum. For the example data sets included here, there are several different combinations of values for parameters ap and ao that yield similar small values for the objective function. The value of xo is more critical than those of ap and ao. In my experience, these different sets of parameter values yield nearly the same values of xp. Thus the existence of local minima is more significant for estimating the composition of contaminating material than for the xp estimates. I recommend that the initial value of ap be set to less than the lowest value of as in the data set, and the initial value of ao be set to the highest value of as. The initial values of xo or xorg can be set to the best values available based on composition of whole soil.
Outliers and bad data points may have a large effect on soil-corrected estimates. Thus careful initial screening of the data is important.

Ash-free Estimates
The xa estimates are included on the 'data.out' file for comparison to xp estimates. To see the regressions of xa vs. as, fix ap on the 'input' file to 0, ao to 1, xo or xorg to 0, and iopt to 0.

back to home page