likfitLgm {geostatsp} | R Documentation |
Maximum likelihood (ML) or restricted maximum likelihood (REML) parameter estimation for (transformed) Gaussian random fields.
likfitLgm(formula, data,
paramToEstimate = c("range","nugget"),
reml=TRUE,
coordinates=data,
param=NULL,
upper=NULL,lower=NULL, parscale=NULL,
verbose=FALSE)
loglikLgm(param,
data, formula, coordinates=data,
reml=TRUE,
minustwotimes=TRUE,
moreParams=NULL)
formula |
A formula for the fixed effects portion of the model, specifying a response and covariates.
Alternately, |
data |
An object of class |
coordinates |
A |
param |
A vector of model parameters, with named elements being amongst
|
reml |
Whether to use Restricted Likelihood rather than Likelihood, defaults to |
paramToEstimate |
Vector of names of model parameters to estimate, with parameters excluded from this list being fixed. The variance parameter and regression coefficients are always estimated even if not listed. |
lower |
Named vector of lower bounds for model parameters passed to |
upper |
Upper bounds, as above. |
parscale |
Named vector of scaling of parameters passed as |
minustwotimes |
Return -2 times the log likelihood rather than the likelihood |
moreParams |
Vector of additional parameters, combined with |
verbose |
if |
likfitLgm
produces list with elements
parameters |
Maximum Likelihood Estimates of model parameters |
varBetaHat |
Variance matrix of the estimated regression parameters |
optim |
results from |
trend |
Either formula for the fixed effects or names of the columns
of the model matrix, depending on |
summary |
a table of parameter estimates, standard errors, confidence intervals, p values, and a logical value indicating whether each parameter was estimated as opposed to fixed. |
resid |
residuals, being the observations minus the fixed effects, on the transformed scale. |
loglikLgm
returns a scalar value, either the log likelihood or -2 times the
log likelihood. Attributes of this result include the vector of
parameters (including the MLE's computed for the variance and coefficients),
and the variance matrix of the coefficient MLE's.
n=40
mydat = vect(
cbind(runif(n), seq(0,1,len=n)),
atts=data.frame(cov1 = rnorm(n), cov2 = rpois(n, 0.5))
)
# simulate a random field
trueParam = c(variance=2^2, range=0.35, shape=2, nugget=0.5^2)
set.seed(1)
oneSim = RFsimulate(model=trueParam,x=mydat)
values(mydat) = cbind(values(mydat) , values(oneSim))
# add fixed effects
mydat$Y = -3 + 0.5*mydat$cov1 + 0.2*mydat$cov2 +
mydat$sim + rnorm(length(mydat), 0, sd=sqrt(trueParam["nugget"]))
plot(mydat, "sim", col=rainbow(10), main="U")
plot(mydat, "Y", col=rainbow(10), main="Y")
myres = likfitLgm(
formula=Y ~ cov1 + cov2,
data=mydat,
param=c(range=0.1,nugget=0.1,shape=2),
paramToEstimate = c("range","nugget")
)
myres$summary[,1:4]
# plot variograms of data, true model, and estimated model
myv = variog(mydat, formula=Y ~ cov1 + cov2,option="bin", max.dist=0.5)
# myv will be NULL if geoR isn't installed
if(!is.null(myv)){
plot(myv, ylim=c(0, max(c(1.2*sum(trueParam[c("variance", "nugget")]),myv$v))),
main="variograms")
distseq = seq(0, 0.5, len=50)
lines(distseq,
sum(myres$param[c("variance", "nugget")]) - matern(distseq, param=myres$param),
col='blue', lwd=3)
lines(distseq,
sum(trueParam[c("variance", "nugget")]) - matern(distseq, param=trueParam),
col='red')
legend("bottomright", fill=c("black","red","blue"),
legend=c("data","true","MLE"))
}
# without a nugget
myresNoN = likfitLgm(
formula=Y ~ cov1 + cov2,
data=mydat,
param=c(range=0.1,nugget=0,shape=1),
paramToEstimate = c("range")
)
myresNoN$summary[,1:4]
# plot variograms of data, true model, and estimated model
myv = variog(mydat, formula=Y ~ cov1 + cov2,option="bin", max.dist=0.5)
if(!is.null(myv)){
plot(myv, ylim=c(0, max(c(1.2*sum(trueParam[c("variance", "nugget")]),myv$v))),
main="variograms")
distseq = seq(0, 0.5, len=50)
lines(distseq,
sum(myres$param[c("variance", "nugget")]) - matern(distseq, param=myres$param),
col='blue', lwd=3)
lines(distseq,
sum(trueParam[c("variance", "nugget")]) - matern(distseq, param=trueParam),
col='red')
lines(distseq,
sum(myresNoN$param[c("variance", "nugget")]) -
matern(distseq, param=myresNoN$param),
col='green', lty=2, lwd=3)
legend("bottomright", fill=c("black","red","blue","green"),
legend=c("data","true","MLE","no N"))
}
# calculate likelihood
temp=loglikLgm(param=myres$param,
data=mydat,
formula = Y ~ cov1 + cov2,
reml=FALSE, minustwotimes=FALSE)
# an anisotropic example
trueParamAniso = param=c(variance=2^2, range=0.2, shape=2,
nugget=0,anisoRatio=4,anisoAngleDegrees=10, nugget=0)
mydat$U = geostatsp::RFsimulate(trueParamAniso,mydat)$sim
mydat$Y = -3 + 0.5*mydat$cov1 + 0.2*mydat$cov2 +
mydat$U + rnorm(length(mydat), 0, sd=sqrt(trueParamAniso["nugget"]))
oldpar = par(no.readonly = TRUE)
par(mfrow=c(1,2), mar=rep(0.1, 4))
plot(mydat, col=as.character(cut(mydat$U, breaks=50, labels=heat.colors(50))),
pch=16, main="aniso")
plot(mydat, col=as.character(cut(mydat$Y, breaks=50, labels=heat.colors(50))),
pch=16,main="iso")
myres = likfitLgm(
formula=Y ~ cov1 + cov2,
data=mydat,
param=c(range=0.1,nugget=0,shape=2, anisoAngleDegrees=0, anisoRatio=2),
paramToEstimate = c("range","nugget","anisoRatio","anisoAngleDegrees")
)
myres$summary
par(oldpar)
par(mfrow=c(1,2))
myraster = rast(nrows=30,ncols=30,xmin=0,xmax=1,ymin=0,ymax=1)
covEst = matern(myraster, y=c(0.5, 0.5), par=myres$param)
covTrue = matern(myraster, y=c(0.5, 0.5), par=trueParamAniso)
plot(covEst, main="estimate")
plot(covTrue, main="true")
par(oldpar)