R/lgm.R

In geostatsp: Geostatistical Modelling with Likelihood and Bayes

setGeneric('lgm', function(
	formula,data,grid,covariates, 
	buffer=0,
	shape=1, boxcox=1, nugget = 0, 
	expPred=FALSE, nuggetInPrediction=TRUE,
	reml=TRUE,mc.cores=1,
	aniso=FALSE,
	fixShape=TRUE,
	fixBoxcox=TRUE,
	fixNugget = FALSE,
	...) 
standardGeneric("lgm")
)

# sort out formula
# null formula
setMethod("lgm", 
	signature("missing", "ANY", "ANY", "ANY"), 
	function(formula, data, grid, 
		covariates, 
		buffer=0,
		shape=1, boxcox=1, nugget = 0, 
		expPred=FALSE, nuggetInPrediction=TRUE,
		reml=TRUE,mc.cores=1,
		aniso=FALSE,
		fixShape=TRUE,
		fixBoxcox=TRUE,
		fixNugget = FALSE,
		...) {

		formula =  1 

		callGeneric(formula, data, grid, 
			covariates, 
			buffer,
			shape, boxcox, nugget, 
			expPred, nuggetInPrediction,
			reml, mc.cores,
			aniso,
			fixShape,
			fixBoxcox,
			fixNugget,
			...)
	}
)


setMethod("lgm", 
	signature("numeric", "ANY", "ANY", "ANY"),
	function(formula, data, grid, 
		covariates,           
		buffer=0,
		shape=1, boxcox=1, nugget = 0, 
		expPred=FALSE, nuggetInPrediction=TRUE,
		reml=TRUE,mc.cores=1,
		aniso=FALSE,
		fixShape=TRUE,
		fixBoxcox=TRUE,
		fixNugget = FALSE,
		...) {

		formula = names(data)[formula]

		callGeneric()
	}
)

# change character to formula
setMethod("lgm", 
	signature("character", "ANY", "ANY", "ANY"),  
	function(formula, data, grid, 
		covariates, 
		buffer=0,
		shape=1, boxcox=1, nugget = 0, 
		expPred=FALSE, nuggetInPrediction=TRUE,
		reml=TRUE,mc.cores=1,
		aniso=FALSE,
		fixShape=TRUE,
		fixBoxcox=TRUE,
		fixNugget = FALSE,
		...) {

		if(length(names(covariates)))
			names(covariates) = gsub("[[:punct:]]|[[:space:]]","_", names(covariates))
		if(length(covariates) & !length(names(covariates))) 
			names(covariates) = paste("c", 1:length(covariates),sep="")			

		if(length(formula)==1)
			formula = unique(c(formula, names(covariates)))
		if(length(formula)==1)
			formula = c(formula, '1')

		formula = paste(formula[1] , "~",
			paste(formula[-1], collapse=" + ")
		)
		formula = as.formula(formula)

		callGeneric()
	}
)



# numeric cells, create raster from data bounding box

setMethod("lgm", 
	signature("formula", "SpatVector", "numeric", "ANY"),
	function(formula, data, grid, covariates, 
		buffer=0,
		shape=1, boxcox=1, nugget = 0, 
		expPred=FALSE, nuggetInPrediction=TRUE,
		reml=TRUE,mc.cores=1,
		aniso=FALSE,
		fixShape=TRUE,
		fixBoxcox=TRUE,
		fixNugget = FALSE,
		...) {

		grid = squareRaster(data, grid)

		callGeneric()
	}
)

# missing covariates, create empty list
setMethod("lgm", 
	signature("formula", "SpatVector", "SpatRaster", "missing"),
	function(formula, data, grid, covariates, 
		buffer=0,
		shape=1, boxcox=1, nugget = 0, 
		expPred=FALSE, nuggetInPrediction=TRUE,
		reml=TRUE,mc.cores=1,
		aniso=FALSE,
		fixShape=TRUE,
		fixBoxcox=TRUE,
		fixNugget = FALSE,
		...) {

		covariates = list()

		callGeneric(formula, data, grid, 
			covariates, 
			buffer,
			shape, boxcox, nugget, 
			expPred, nuggetInPrediction,
			reml, mc.cores,
			aniso,
			fixShape,
			fixBoxcox,
			fixNugget,
			...)
	}
)


setMethod("lgm",
	signature("formula", "SpatVector", "SpatRaster", "list"),
	function(formula, data, grid, covariates, 
		buffer=0,
		shape=1, boxcox=1, nugget = 0, 
		expPred=FALSE, nuggetInPrediction=TRUE,
		reml=TRUE,mc.cores=1,
		aniso=FALSE,
		fixShape=TRUE,
		fixBoxcox=TRUE,
		fixNugget = FALSE,
		...) {

		dataCov = #geostatsp:::
		gm.dataSpatial(
			formula, data, 
			grid, covariates, buffer)

		data=dataCov$data 
		grid=dataCov$grid 
		covariates=dataCov$covariates

		callGeneric()
	}
)

setMethod("lgm",
	signature("formula", "SpatVector", "SpatRaster", "SpatRaster"),
	function(formula, data, grid, covariates, 
		buffer=0,
		shape=1, boxcox=1, nugget = 0, 
		expPred=FALSE, nuggetInPrediction=TRUE,
		reml=TRUE,mc.cores=1,
		aniso=FALSE,
		fixShape=TRUE,
		fixBoxcox=TRUE,
		fixNugget = FALSE,
		...) {

		dataCov = gm.dataSpatial(
			formula, data, 
			grid, covariates, buffer)

		data=dataCov$data 
		grid=dataCov$grid 
		covariates=dataCov$covariates

		callGeneric()
	}
)

# the real work
setMethod("lgm", 
	signature("formula", "SpatVector", "SpatRaster","data.frame"), 
	function(formula, data, grid, covariates, 
		buffer=0,
		shape=1, boxcox=1, nugget = 0, 
		expPred=FALSE, nuggetInPrediction=TRUE,
		reml=TRUE,mc.cores=1,
		aniso=FALSE,
		fixShape=TRUE,
		fixBoxcox=TRUE,
		fixNugget = FALSE,
		...) {

		locations = grid

		dots <- list(...)  

		param=dots$param	
		if(!length(param)) {
			param=c()
		}

		paramToEstimate	= c(
			"variance", "range", "shape","nugget","boxcox"
		)[c(
			TRUE, TRUE, !fixShape, !fixNugget, !fixBoxcox
		)]
		
		range=NA
		Spar = c(shape=as.numeric(shape),
			nugget=as.numeric(nugget),
			range=NA,
			boxcox=as.numeric(boxcox))

		if(aniso) {
			Spar = c(Spar, anisoAngleRadians=NA,anisoRatio=NA)
			paramToEstimate = c(paramToEstimate,
				"anisoAngleDegrees","anisoRatio")		
		}

		Spar = Spar[!names(Spar) %in% names(param)]
		param = c(param, Spar)

# to do: make sure factors in rasters are set up correctly
	   # NA's for levels without data
	   # have most common level the baseline

# call likfit

		dots$param = param
		dots$formula=formula
		dots$data=data
		dots$paramToEstimate=paramToEstimate
		dots$reml = reml
		
		likRes = do.call(likfitLgm, dots)

			# call krige	
			krigeRes =  krigeLgm(
				formula=formula,data=data,
				grid=grid,
				covariates=covariates, param=likRes$param, 
				expPred=expPred,
				nuggetInPrediction=nuggetInPrediction,
				mc.cores=mc.cores
			)

		res = c(predict=krigeRes, likRes)

    # add confidence intervals for covariance parameters
		theInf=informationLgm(res)

		res$varBetaHat = list(beta=res$varBetaHat)
		names(res) = gsub("varBetaHat", "varParam", names(res))
		res$varParam$information = theInf$information


		res$summary = 	theInf$summary

		if(is.na(res$summary ['anisoAngleDegrees','ci0.05']) &
			!is.na(res$summary ['anisoAngleRadians','ci0.05']) ){
			ciCols = grep("^ci0\\.[[:digit:]]+$", colnames(res$summary ))
			res$summary ['anisoAngleDegrees',ciCols] =
			(360/(2*pi))*res$summary ['anisoAngleRadians',ciCols]
			res$summary ['anisoAngleDegrees','Estimated'] = 
			res$summary ['anisoAngleRadians','Estimated']
		}


		if(FALSE){
			for(Dvar in names(covariates)) {
				theLevels =levels(covariates[[Dvar]])[[1]]
				if(!is.null(nrow(theLevels))){
					for(D in 1:nrow(theLevels)) {
						rownames(res$summary) = gsub(
							paste("(factor)?(\\()?", Dvar, "(\\))?:?", 
								theLevels[D,1],"$",sep=""),
							paste(Dvar, ":",theLevels[D,2],sep=""), 
							rownames(res$summary))
					}
				}
			}
		}

		theOrder = c('sdNugget','sdSpatial', 'range', 'shape','anisoRatio',
			'anisoAngleRadians','anisoAngleDegrees', 'boxcox')  
		theOrder = na.omit(match(theOrder, rownames(res$summary)))
		notInOrder = (1:nrow(res$summary))[-theOrder]
		res$summary = res$summary[c(notInOrder,theOrder),]


	   # if range is very big, it's probably in metres, convert to km
		if(res$summary['range','estimate']>1000) {
			logicalCol = names(res$summary) == "Estimated"
			res$summary["range",!logicalCol] = 
			res$summary["range",!logicalCol] /1000
			rownames(res$summary) = gsub("^range$", "range/1000", 
				rownames(res$summary))
		}
		class(res) = c('lgm',class(res))    
		return(res)
	}

)
AIC.lgm = function(object, ..., k = 2) {

	AIC(logLik(object, ..., k))

}


logLik.lgm = function(object, ...){
	res = object$opt$logL
	res = res[grep('^logL', names(res))]

	srows = rownames(object$summary)
	df = sum(
		object$summary[
		grep("^anisoAngle", srows, invert=TRUE),
		'Estimated']
	) + any(
		object$summary[
		grep("^anisoAngle", srows, invert=FALSE),
		'Estimated']
	)
	attributes(res)$df = df
	attributes(res)$nobs = try(nrow(data.frame(object$data)), silent=TRUE)
	class(res)= 'logLik'
	res
}

summary.lgm = function(object, ...) {
	object$summary
}