likfitLgm: Likelihood Based Parameter Estimation for Gaussian Random...

In geostatsp: Geostatistical Modelling with Likelihood and Bayes

Description

Maximum likelihood (ML) or restricted maximum likelihood (REML) parameter estimation for (transformed) Gaussian random fields.

Usage

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)

Arguments

formula	A formula for the fixed effects portion of the model, specifying a response and covariates. Alternately, data can be a vector of observations and formula can be a model matrix.
data	An object of class SpatialPointsDataFrame, a vector of observations, or a data frame containing observations and covariates.
coordinates	A SpatialPoints object containing the locations of each observation, which defaults to data. Alternately, coordinates can be a symmetricMatrix-class or dist object reflecting the distance matrix of these coordinates (though this is only permitted if the model is isotropic).
param	A vector of model parameters, with named elements being amongst range, nugget, boxcox, shape, anisoAngleDegrees, anisoAngleRadians, anisoRatio, and possibly variance (see matern). When calling likfitLgm this vector is a combination of starting values for parameters to be estiamated and fixed values of parameters which will not be estimated. For loglikLgm, it is the covariance parameters for which the likelihood will be evaluated.
reml	Whether to use Restricted Likelihood rather than Likelihood, defaults to TRUE.
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 optim, defaults are used for parameters not specified.
upper	Upper bounds, as above.
parscale	Named vector of scaling of parameters passed as control=list(parscale=parscale) to optim.
minustwotimes	Return -2 times the log likelihood rather than the likelihood
moreParams	Vector of additional parameters, combined with param. Used for passing fixed parameters to loglikLgm from within optim.
verbose	if TRUE information is printed by optim.

Value

likfitLgm produces list with elements

parameters	Maximum Likelihood Estimates of model parameters
varBetaHat	Variance matrix of the estimated regression parameters
optim	results from optim
trend	Either formula for the fixed effects or names of the columns of the model matrix, depending on trend supplied.
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.

See Also

lgm

Examples

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))
}

Example output

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

geostatsp documentation built on Oct. 5, 2021, 9:10 a.m.