corrHLfit
performs the joint estimation of correlation parameters, fixed effect and dispersion parameters.
corrHLfit(formula, data, init.corrHLfit = list(), init.HLfit = list(), ranFix = list(), lower = list(), upper = list(), trace = list(file = NULL, append = TRUE), objective = "p_bv", resid.model = ~1, resid.formula, control.dist = list(), control.corrHLfit = list(), processed = NULL, family = gaussian(), ...)
list(rho=1,nu=1,lambda=1,phi=1)
where rho
and nu
are parameters of the Matérn family, and
lambda
and phi
are dispersion parameters (see Details in spaMM
for the meaning of these parameters).
All are optional, but giving values for
a dispersion parameter changes the ways it is estimated (see Details).
rho
may be a vector (see make_scaled_dist
) and, in that case, it is possible that some or all of its elements are NA
, for which corrHLfit
substitute automatically determined values.
HLfit
argument.
init.corrHLfit
, but specifying fixed values of the parameters not estimated. init.corrHLfit
, used as lower values in calls to optim
. See Details for default values.
lower
, but upper values.trace=list(file=,append=F)
, some trace information is written in the file 'filename'.
This file is written over by each new call of corrHLfit
unless append=T
.
Information is written for each HLCor
call. It includes the APHLs
, and all given fixed random effect
parameters in the HLCor call. Information about the variables is printed at the end of the file, but may be slightly inaccurate
(some amount of guesswork is expected from users venturing into trace
). The arguments of all HLCor calls are saved in RData files
(if this is felt inappropriate, then using spaMM.options(TRACE.UNLINK=TRUE)
will keep only the latest HLCor call).
optim
.
Either "p_bv"
for restricted likelihood or "p_v"
for marginal likelihood.
HLfit
arguments.control.dist
in HLCor
control.corrHLfit$optim
should be a list of arguments for optim, e.g. control.corrHLfit = list(optim=list(control=list(ndeps=rep(2e-3,3))))
HLCor
, HLfit
or designL.from.Corr
, for example the distMatrix
argument
of HLCor
. Arguments that do not fit within these functions are detected and a warning is issued.
HLCor
call, with additional attributes. The HLCor
call is evaluated at the estimated correlation parameter values. These values are included in the return object as its $corrPars
member. The attributes added by corrHLfit
include the original call of the function (which can be retrived by getCall
(corrHLfit
.
spaMM.options
. By default corrHLfit
will estimate correlation parameters by maximizing the objective
value returned by HLCor
calls wherein the dispersion parameters are estimated jointly with fixed effects for given correlation parameters. If dispersion parameters are specified in init.corrHLfit
, they will also be estimated by maximizing the objective
value, and HLCor calls
will not estimate them jointly with fixed effects. This means that in general the fixed effect estimates may vary depending on init.corrHLfit
when any form of REML correction is applied.
Correctly using corrHLfit
for likelihood ratio tests
of fixed effects may then by tricky. It is safe to perform full ML fits of all parameters (using objective="p_v",HLmethod="ML"
) for such tests (see Examples). The higher level function fixedLRT
is a safe interface for likelihood ratio tests using some form of REML estimation in corrHLfit
.
attr(
and ...$upper
gives the lower and upper bounds for optimization of correlation parameters. These are the default values if the user did not provide explicit values. For the adjacency model, the default values are the inverse of the maximum and minimum eigenvalues of the adjMatrix
. For the Matérn model, the default values are not so easily summarized: they are intended to cover the range of values for which there is statistical information to distinguish among them.
Loaloa
.
See fixedLRT
for likelihood ratio tests.
# Example with an adjacency matrix (autoregressive model):
# see 'adjacency' documentation page
#### Examples with Matérn correlations
## A likelihood ratio test based on the ML fits of a full and of a null model.
if (spaMM.getOption("example_maxtime")>18) {
data(blackcap)
fullfit <- corrHLfit(migStatus ~ means+ Matern(1|latitude+longitude),data=blackcap,
HLmethod="ML")
summary(fullfit)
nullfit <- corrHLfit(migStatus ~ 1 + Matern(1|latitude+longitude),data=blackcap,
HLmethod="ML")
summary(nullfit)
## p-value:
1-pchisq(2*(logLik(fullfit)-logLik(nullfit)),df=1)
}
## see data set Loaloa for additional examples
Run the code above in your browser using DataLab