The likelihood methods of Kimura (1990) are used to fit any nonlinear equation, correct for
heterogeneity of variances, and estimate common or separate within-group variances depending on
user-specifications. A main assumption is errors are normally-distributed.
The results of the model fits can then be used in function compare.lrt.plus
to determine if parameterizations differ significantly from each other by using a likelihood ratio
and an F test.
Steps of the modeling process are as follows:
1) Specify the nonlinear model equation in the same
formula format as would be done in function nls. For example, the von Bertalanffy growth equation
is written as:
sl~Linf*(1-exp(-K*(age-t0)))
where sl is the variable name for length data, age is the variable name for age data,
and Linf, K and t0 are parameters to be estimated.
2) Specify the parameter formulae under params
. These formulae
are used to indicate that additional parameters based on some group variable should be estimated.
For example,
params=list(Linf~1,K~1,t0~1)
specifies single parameters are estimated for Linf, K and t0.
params=list(Linf~sex,K~1,t0~1)
specifies that separate Linfs are to be estimated for each sex and only single estimates
for K and t0.
params=list(Linf~sex,K~sex,t0~sex)
specifies that separate Linfs, Ks and t0s are to be estimated for each sex. Different group variables
for each parameter are not allowed.
3) Specify start values for all parameters. For example, if separate Linfs, Ks and t0s for a group
variable like sex (only two-levels: M and F), then 6 starting values must be given. When parameters are
based on a group variable, then the first estimate of a parameter will be the reference value (labeled as Intercept)
and the remaining parameters will be estimated as a deviation from that reference value. Reference values
are determined by alphanumeric order of levels within the group variable.
start=list(Linf=c(300,10),K=c(0.3,0.05),t0=c(0,-0.5))
is an example of the starting values for the 6 parameter model mentioned above. Warning messages are generated
if the number of start values is less than or greater than the number of parameters being estimated. Internally,
code will add (1/10th of first value) or truncate (last number(s) in list) start values in these cases.
However, the user should specify the appropriate number of values to ensure successful optimization.
4) Specify the within-group variance structure. If the within-group variance is assumed
the same among groups, then
within_grp_var=~1
which is the default specification. If within-group variances are suspected to differ among groups (e.g., sex), then
within_grp_var=~sex
Separate variances will be estimated for each level of the group variable. Whether or not better model fits can be
obtained by estimating separate group variances can be tested using the model comparison methods (see below).
When estimating thetas (correcting for heterogeneity), explore different starting values for the main parameters to
ensure global convergence.
5) Specify the correction for heterogenity (cfh
) argument(s) if needed. Initial curve fittings should be performed
and plots of residuals versus fitted values examined to determine if there is a change in residual variation with
increasing fitted values. If so, this indicates the presense of heterogeneity in variance which must be corrected to obtain
unbiased parameters estimates, standard errors, residual sum-of-squares, etc. Kimura (1990) uses the power function (same as the
varPower function in Pinheiro and Bates (2004)) and additional parameters known as theta are estimated.
If cfh
is NULL, then homogeneity of variance is assumed. If heterogeneity of variance needs to be accounted for, specify
cfh
as
cfh=list(form=~1,value=0,fixed=NULL)
form
is a formula and is 1 if a single theta is assumed equal among groups. If individual thetas are desired for groups (heterogenity is different
for each group), then a group variable is used (e.g.,form=~sex).
value
is the initial starting value(s) for theta(s). If more than 1 theta will
be estimated, provide the same number of starting values within c() as thetas.
fixed
is used to indicate whether the thetas will be estimated (default NULL) or assumed fixed to numeric values specified by the user.
cfh=list(form=~sex,value=0,fixed=c(0.5,0.9))
indicates that thetas for each sex (two-levels: M and F) will not be estimated, but fixed to values of 0.5 and 0.9
6) Run the model function. Parameter estimation is performed intially by using the optimization function nlminb.
The estimated parameters are then used as starting values and optimization is performed again by
using function optim to obtain the final parameter estimates and the Hessian matrix from
which standard errors are derived. Unlike estimation of thetas conducted in function gnls in package nlme, thetas herein
are estimated as parameters, standard errors are generated, and t-tests for significance are conducted. These extra parameters are counted in the
determination of residual and model degrees of freedom.
To convert a non-reference level estimate to the same scale as the reference level, the reference value and parameter estimate are added together.
For example, if estimates of Linf for two groups are 300 and 5, then adding them gives the Linf of 305 for the
non-reference group.
Model Comparisons
As in the growthlrt function based on Kimura (1980), growth curves are tested for differences by using likelihood ratio tests.
These tests assume homogeneity of variances among groups which is why heterogeneity must be corrected before proceeding. Unlike
the growthlrt function, growthlrt.plus does not automatically make the comparisons. The user must develop the model structures,
save each oject, and test for differences using the function compare.lrt.plus. Examples are provided below.