# 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

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14``` ```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.

`lgm`

## Examples

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

### Example output

```Loading required package: Matrix
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