Function designed to be called by either eblupDyn or eblupRY
to produce EBLUP small area estimates of the dynamic or Rao-Yu time
series models through either ML or REML estimation of the variance
components. For completeness and as a reference to several optional
parameters, the function is documented here, but users are encouraged to
base applications on calls to the more convenient eblupDyn or eblupRY.
dynRYfit(y, X, M, TI, NV=1, vcov_e, maxiter=100,
iter.tol=.1e-5, ncolx=NULL, sig2_u=1, sig2_v=1,
rho=.8, rho_u =.4, delta=NULL, rho.fixed=NULL,
y.include=NULL, ids=NULL, contrast.matrix=NULL,
baby.steps=TRUE, dampening=NULL, iter.history=NULL,
sig2.min.factor=.0001, max.rho_u=.98, max.rho=NULL,
tol=.Machine$double.eps, y.rescale=NULL,
llike.only=FALSE, method=c("REML", "ML"),
model=c("dyn", "RY"))
In the univariate case, a vector of length M*nt with
the eblup estimates in the same sort order as y. In the multivariate
case, a matrix of M*nt rows and NV columns.
A list summarizing the fit of the model with the following:
model: form of the model: T - Dynamic or RaoYu; REML
or ML.
covergence: a logical value indicating whether the
convergence criterion was met.
A labelled vector with the estimated variance components, correlations, and number of iterations.
A labelled vector of coefficients for the fixed effects of the model or models.
A data frame with D rows and one or more columns of
numeric or character domain identifiers.
An ordered vector of the variance components (see above). It may be used as starting values for additional iterations.
MSE estimates for eblup.
The g1 term of the MSE estimate.
The g2 term of the MSE estimate.
The g3 term of the MSE estimate.
Estimates based on fixed effects only.
The variance-covariance matrix for the estimates in
coef.
Weights given to the direct estimate in forming eblup.
Weights given to the direct estimate, including effects through estimating the fixed effect coefficients.
Estimates requested by the specified contrasts.
MSE estimates for contrast.est.
The g1 term in the estimation of contrast.mse.
The g2 term in the estimation of contrast.mse.
The g3 term in the estimation of contrast.mse.
Contrast estimates based on the fixed effect model.
Variance estimates for the fixed effect model.
Weight wt1 given to the direct estimate in estimating the contrasts.
Weight wt2 in estimating the contrasts.
Information matrix for the components of delta.
Variance covariance matrix for coef.
Values of delta at each iteration.
Values of the log-likelihood (ML) or restricted log-likelihood (REML) at each iteration.
Number of cycles in the internal loop to determine
adj.factor within each iteration.
Values of inf.mat at each iteration.
Vector to be multiplied by the inverse information matrix to determine the change in the parameters.
List of parameters eligible for change at each iteration. Parameters with estimated changes out of bounds will not be eligible.
Adjustment to the vector change in the parameters at each iteration.
A 4-row matrix of warnings at each iteration, where warning 1 is set to 1 for the iteration if the algorithm has not found an increase in the restricted log likelihood or log likelihood, warning 2 is set to 1 if the maximum number of iterations is reached, warning 3 is set to 1 if the estimated variance-covariance matrix becomes singular, and warning 4 is set to 1 if the coefficients of the fixed effects cannot be estimated.
For a univariate model, the dependent variable sorted in ascending order by time within domain. For a multivariate model, the dependent variables sorted in ascending order by time within variable within domain.
A matrix of independent variables with the same number of
rows as the length of y.
The total number of domains, equivalent to D in
eblupDyn and eblupRY.
If TI is of length 1, it should specify the number of equally
spaced time instances, nt; for example, the number of years in an
unbroken sequence of years.
If some of the time instances in the time interval are to be omitted from the model,
the first element of TI should be 1, the last element should be the ending
instance, and the elements in between should indicate the other instances to include.
For example, over a period of 6 years, TI=c(1,2, 4:6) omits the 3rd year
from the model, giving nt=5.
The number of dependent variables.
For the univariate model, the sampling covariance matrix for
the direct estimates of the M*nt elements of the dependent
variable, where nt is based on TI. The covariance matrix should be
in the form of a square matrix with M*nt rows and columns.
Non-zero covariances between domains are not allowed, so the matrix must have
block diagonal form with M blocks, each of which is a square matrix
with nt rows and columns. Note that within domain, non-zero
covariances are allowed over time.
For the multivariate model, the square covariance matrix for the
M*NV*nt elements of the dependent variables. The matrix should be
in the form of a square matrix with M*NV*nt rows and columns. Time
should vary within variable, which should vary within domain.
Non-zero covariances between domains are not allowed, but non-zero
covariances may be present across time and between variables within each domain.
The maximum number of iterations allowed for the Fisher-scoring algorithm, with a default value of 100.
The convergence tolerance limit for the Fisher-scoring algorithm, with a default value of .000001.
For a univariate model, the number of columns of X. For a multivariate model, a vector of length NV must be specified giving the number of columns of X used for each dependent variable.
An initial starting value or values for the variance of the random increments.
An initial starting value or values for a domain level random effect. In the Rao-Yu model, the random effect is constant over time, whereas in the dynamic model it is an initial effect subject to dampening over time.
The correlation across time. This correlation is assumed to be the same for the dependent variables in the multivariate model.
For NV > 1 only, the (NV*(NV-1))/2 starting values for
the correlations between the random effects of the different dependent
variables. If a single value is given, it will be used for the
(NV*(NV-1))/2 components. The sort order corresponds to a lower triangle
of the covariance matrix.
The random effect components in the preferred internal order.
Specification of delta will override any specification of
sig2_u, sig2_v, rho, or rho_u. In the
univariate case, delta should contain sig2_u, sig2_v, and
rho. In the bivariate case, delta should be of length 6 with
sig2_u and sig2_v each of length 2, then rho, and
rho_u. Similarly, for 3 dependent variables, the length of
delta is 10, with the last 3 elements rho_u in lower
triangular form.
If TRUE, the value of rho imbedded in delta,
if specified, or else given by rho, will remain fixed during the
iterations. Among other features, this allows the likelihood function
for trial values of rho to be computed at the maximum over the
other random effect parameters.
If specified, vector of length M to indicate
which domains to include in the estimation, with 1 signalling inclusion
and 0 exclusion. Estimates for the excluded domains will be based on
the fixed effects model only.
A data frame with M rows giving ids for each of the
domains. The data frame is copied into the returned object.
A matrix of coefficients of contrasts. The matrix
must have nt*NV rows, but it can contain an arbitrary number
of columns. Within each domain, the coefficients are applied
to the nt*NV EBLUP estimates.
Unless specified as FALSE, the first five iterations of the
Fisher scoring algorithm are dampened by factors of
c(.0625, .125, .25, .5, .75). These heuristically derived
factors appear to lessen drastic overshooting of the true maximum
in the initial iterations.
A factor used to dampen the changes in the random
effect parameters. Unlike baby.steps, its effect persists
during all of the iterations until convergence. Note that the
"factory setting" of this parameter is 1 for the dynamic model but
.9 for the Rao-Yu model.
If TRUE, key values are saved during each iteration and included as additional items, described below, in the returned list:
delta.hist,
llikelihood.hist,
adj.hist,
s.hist,
ix.hist,
adj.factor.hist, and
warning.hist.
If iter.history is not specified, these items will be returned only
if the calculations are begun but not successfully completed.
A factor to multiply the minimum direct variance to use as a minimum value for any of the variance components. The iterations will be constrained not to go below the resulting bounds.
A maximum allowed value for the estimated rho_u. The
default value is .98.
A maximum allowed value for rho. By default,
eblupRY sets this value to .98.
A tolerance value used by matrix routines to prevent numerical instability. The value may be set to a lower value to encourage covergence, but appropriate caution should be applied.
In the univariate case, a scaler multiplier for all of
the y values. If the y values are either too small or
too large, the information matrix may become unstable. Setting this
value to 1 has no effect; setting it to 10 or 100 rescales very
small y values to a more appropriate range. Similarly,
positive values less than 1 may be used to rescale large y
values. The effect of rescaling is removed before normal return
from the function, within the limits of normal precision.
In the multivariate case, y.rescale may be a vector of
length NV to rescale each component separately.
Compute the log-likelihood (ML) or restricted
log-likelihood (REML) without further iteration, typically from
values specified by delta.
Use restricted maximum likelihood ("REML") or
maximum likelihood ("ML").
Dynamic ("dyn") or Rao-Yu ("RY").
Robert E. Fay, Mamadou Diallo
Many of arguments can be used to control the iterations if the defaults lead to convergence difficulties.
llike.only in combination with delta permits a
point-by-point investigation of the likelihood surface.
The primary functions eblupDyn and eblupRY determine X,
NV, colx, and model, but the remaining parameters can
be passed to dynRYfit through eblupDyn or eblupRY.