Description Usage Arguments Examples
All the arguments work as in likfit
, except the
additional argument dists.mat
, which takes a symmetric matrix of
distances between observation locations
1 2 3 4 5 6 | likfit(geodata, coords = geodata$coords, data = geodata$data,
trend = "cte", ini.cov.pars, fix.nugget = FALSE, nugget = 0,
fix.kappa = TRUE, kappa = 0.5, fix.lambda = TRUE, lambda = 1,
fix.psiA = TRUE, psiA = 0, fix.psiR = TRUE, psiR = 1, cov.model,
realisations, lik.method = "ML", components = TRUE, nospatial = TRUE,
limits = pars.limits(), dists.mat, print.pars = FALSE, messages, ...)
|
geodata |
a list containing elements |
coords |
an n x 2 matrix where each row has the 2-D
coordinates of the n data locations.
By default it takes the
component |
data |
a vector with n data values. By default it takes the
component |
trend |
specifies the mean part of the model. See documentation
of |
ini.cov.pars |
initial values for the covariance parameters:
sigma^2 (partial sill) and phi (range
parameter). Typically a vector with two components. However a
matrix can be used to provide several initial values. See
|
fix.nugget |
logical, indicating whether the parameter
tau^2 (nugget variance) should be regarded as fixed
( |
nugget |
value of the nugget parameter.
Regarded as a fixed value if |
fix.kappa |
logical, indicating whether the extra parameter
kappa should be regarded as fixed
( |
kappa |
value of the extra parameter kappa.
Regarded as a fixed value if |
fix.lambda |
logical, indicating whether the Box-Cox transformation parameter
lambda should be regarded as fixed
( |
lambda |
value of the Box-Cox transformation parameter
lambda.
Regarded as a fixed value if |
fix.psiA |
logical, indicating whether the anisotropy angle parameter
psi_R should be regarded as fixed
( |
psiA |
value (in radians) for the anisotropy angle parameter
psi_A.
Regarded as a fixed value if |
fix.psiR |
logical, indicating whether the anisotropy ratio parameter
psi_R should be regarded as fixed
( |
psiR |
value, always greater than 1, for the anisotropy ratio parameter
psi_R.
Regarded as a fixed value if |
cov.model |
a string specifying the model for the correlation
function. For further details see documentation for |
realisations |
optional. Logical or a vector indicating the number of replication
for each datum. For further information see |
lik.method |
(formely method.lik) options are |
components |
an n x 3 data-frame with fitted
values for the three model components: trend, spatial and residuals.
See the section |
nospatial |
logical. If |
limits |
values defining lower and upper limits for the model
parameters used in the numerical minimisation.
The auxiliary function |
dists.mat |
n x n symmetric matrix with cost-based distances between observations |
print.pars |
logical. If |
messages |
logical. Indicates whether status messages should be printed on the screen (or output device) while the function is running. |
... |
additional parameters to be passed to the minimisation
function. Typically arguments of the type |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 | ## geodata structure with transformed covariates
data(noise)
if (require(sp)) {
covarnames=sapply(1:3, function(x) paste("d2TV", x, sep=""))
obs.df <- data.frame(Leq=obs$Leq,
1/(1+(as.data.frame(obs)[covarnames]/20)^2))
obs.gd <- as.geodata(cbind(coordinates(obs), obs.df),
data.col="Leq",
covar.col=c('d2TV1','d2TV2','d2TV3'))
trend=~d2TV1*(d2TV2+d2TV3)
## compute euclidean and cost-based empirical variograms
vg.std <- variog(obs.gd, trend=trend)
vg.dmat <- variog(obs.gd, trend=trend, dists.mat=dd.distmat)
## fitting variogram models
vgmdl.std <- likfit(geodata = obs.gd, trend=trend,
ini = c(8,300), cov.model = "matern")
vgmdl.dmat <- likfit(geodata = obs.gd, trend=trend,
ini = c(8,300), cov.model = "matern",
dists.mat=dd.distmat)
## Fitted parameters
data.frame(
parameters=c("tausq","sigmasq","phi"),
Euclidean=c(round(vgmdl.std$tausq,2),round(vgmdl.std$sigmasq,2),round(vgmdl.std$phi,0)),
Cost_based=c(round(vgmdl.dmat$tausq,2),round(vgmdl.dmat$sigmasq,2),round(vgmdl.dmat$phi,0)))
## practical range
## defined as the value for which the correlation function
## decays to 5% of its value at 0
x=seq(0,800)
y=cov.spatial(x,cov.pars=vgmdl.std$cov.pars)
min(x[y<0.05*y[1]]) # 358
y=cov.spatial(x,cov.pars=vgmdl.dmat$cov.pars)
min(x[y<0.05*y[1]]) # 502
# Note that the cost-based analysis detects a
# longer-ranged correlation structure
## plotting and comparing empirical variograms
## (with classical and robust estimation)
## and fitted variogram models
op <- par(mfrow=c(2,2))
u = 13 # binning
for( est in c('classical','modulus') )
{
vg.std <- variog(obs.gd, trend=trend, estimator.type=est, uvec=u)
vg.dmat <- variog(obs.gd, trend=trend, dists.mat=dd.distmat, estimator.type=est, uvec=u)
plot(vg.std,pch=20, cex=1.2, max.dist=max(vg.dmat$u), ylim=c(0,20), col='gray')
lines(vg.std,pch=20, col='gray')
lines(vg.dmat,pch=20,cex=2,col='red')
legend('topleft',c('Euclidean','Cost'),lty=1,lwd=2,col=c('gray','red'))
title(paste('binning:',u,' estimator:',est))
plot(vg.std, pch=20, cex=1.2, max.dist=800, col='gray') # empirical standard
lines(vgmdl.std,lwd=2,col='gray') # fitted model standard
points(as.data.frame(vg.dmat[c(1,2)]),pch=20,cex=2,col='red') # empirical cost-based
lines(vgmdl.dmat,lwd=2,col='darkred') # fitted model cost-based
legend('topleft',c('Euclidean','Cost'),lty=1,lwd=2,col=c('gray','red'))
title(paste('binning:',u,' estimator:',est))
}
par(op)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.