Description Usage Arguments Value See Also Examples
Maximum likelihood (ML) or restricted maximum likelihood (REML) parameter estimation for (transformed) Gaussian random fields.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
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.
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 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 | n=40
mydat = SpatialPointsDataFrame(
cbind(runif(n), seq(0,1,len=n)),
data=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)
mydat@data= cbind(mydat@data , RFsimulate(model=trueParam,x=mydat)@data)
# add fixed effects
mydat$Y = -3 + 0.5*mydat$cov1 + 0.2*mydat$cov2 +
mydat$sim + rnorm(length(mydat), 0, sd=sqrt(trueParam["nugget"]))
spplot(mydat, "sim", col.regions=rainbow(10), main="U")
spplot(mydat, "Y", col.regions=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]
if(interactive() | Sys.info()['user'] =='patrick') {
# 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"]))
par(mfrow=c(1,2))
oldmar = par("mar")
par(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")
par(mar=oldmar)
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(mfrow=c(1,2))
myraster = raster(nrows=30,ncols=30,xmn=0,xmx=1,ymn=0,ymx=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(mfrow=c(1,1))
}
|
Loading required package: Matrix
Loading required package: raster
Loading required package: sp
New output format of RFsimulate: S4 object of class 'RFsp';
for a bare, but faster array format use 'RFoptions(spConform=FALSE)'.
estimate stdErr ci0.005 ci0.995
range 0.2624976 NA NA NA
shape 2.0000000 NA NA NA
sdSpatial 1.7647809 NA NA NA
sdNugget 0.3887618 NA NA NA
anisoRatio 1.0000000 NA NA NA
anisoAngleRadians 0.0000000 NA NA NA
boxcox 1.0000000 NA NA NA
(Intercept) -3.0055478 0.31123617 -3.8072390 -2.2038566
cov1 0.5428786 0.07648083 0.3458770 0.7398802
cov2 0.1432875 0.12914746 -0.1893743 0.4759493
anisoAngleDegrees 0.0000000 NA NA NA
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.