R/profLlgm.R

In geostatsp: Geostatistical Modelling with Likelihood and Bayes

profLlgm = function(fit,mc.cores=1, ...) {
	
	dots = list(...)
	fit$parameters = fillParam(fit$parameters)
	varying = intersect(names(dots), names(fit$parameters))

	for(D in varying)
		dots[[D]] = sort(dots[[D]])
	
	
	nonLinearParams = c('boxcox','shape','nugget','variance',
		'anisoAngleRadians','anisoAngleDegrees','anisoRatio','range')
	
	reEstimate = rownames(fit$summary)[
	fit$summary[,"Estimated"]
	]
	reEstimate = gsub("sdNugget", "nugget", reEstimate)
	reEstimate = gsub("sdSpatial", "variance", reEstimate)
	reEstimate = gsub("range \\(km\\)", "range", reEstimate)
	reEstimate = gsub("/1000", "", reEstimate)
	reEstimate = intersect(reEstimate, nonLinearParams)
	reEstimate = reEstimate[!reEstimate %in% varying]
	
	
	baseParams = fit$parameters
	baseParams = baseParams[names(baseParams)%in%
	nonLinearParams]
	
	baseParams=baseParams[!names(baseParams)%in% varying]
	if(any(reEstimate=='variance')){
		baseParams['nugget'] = baseParams['nugget'] / baseParams['variance']
	}
	baseParams=baseParams[names(baseParams) != 'variance']
	
	
	if(length(grep("^anisoAngle", varying))){
		reEstimate = grep("^anisoAngle", reEstimate,value=TRUE,invert=TRUE)
		baseParams = baseParams[
		grep("^anisoAngle", names(baseParams),value=TRUE,invert=TRUE)
		]
	} 
	
	parValues = do.call(expand.grid, dots[varying])
	
	parList = apply(parValues,1,list)
	parList = lapply(parList, function(qq) c(unlist(qq), baseParams))
	
	if(any(grepl("^Spat", class(fit$data)))) {
		dataArg = as.data.frame(fit$data)
		coordinatesArg = terra::crds(fit$data)
	} else {
		dataArg = fit$data
		coordinatesArg = fit$data
	}
	formulaArg = fit$model$formula
	remlArg = fit$model$reml

	if(mc.cores > 1) {
		cl = parallel::makePSOCKcluster(mc.cores)
		on.exit(parallel::stopCluster(cl), add=TRUE)
		parallel::clusterEvalQ(cl, {
			loadNamespace("geostatsp")
			NULL
		})
		parallel::clusterExport(
			cl,
			varlist=c(
				"parList", "dataArg", "coordinatesArg",
				"formulaArg", "reEstimate", "remlArg"
			),
			envir=environment()
		)
		resL = parallel::parLapply(
			cl, X=parList,
			fun=function(paramHere) {
				get("likfitLgm", envir=asNamespace("geostatsp"))(
					param=paramHere,
					data=dataArg,
					formula=formulaArg,
					paramToEstimate=reEstimate,
					reml=remlArg,
					coordinates=coordinatesArg
				)
			}
		)
	} else {
		resL = mapply(
			likfitLgm, 
			param=parList,
			MoreArgs=list(
				data=dataArg, 
				formula=formulaArg,
				paramToEstimate=reEstimate,
				reml=remlArg,
				coordinates=coordinatesArg
			),
			SIMPLIFY=FALSE
		)
	}

	resL2 = lapply(resL, function(qq) qq$optim$logL)
	resL = simplify2array(resL2)

	if(length(varying)==1) {
		forNames =  c(
			dots[[varying]],
			fit$param[varying]
		)
		theorder = order(forNames)
		L = c(
			resL[grep("^m2logL",rownames(resL)),],
			fit$optim$logL[
			grep("^m2logL",
				names(fit$optim$logL))]
		)[theorder]
		dots[[varying]] = forNames = forNames[theorder]
		names(L) = paste(varying, forNames,sep="")
	} else {
		thedimnames=dots[varying]
		for(D in names(thedimnames)) {
			thedimnames[[D]] = paste(
				D, thedimnames[[D]],sep='_')
		}
		L = array(resL[grep("^m2logL",rownames(resL)),], 
			unlist(lapply(thedimnames, length)),
			dimnames=thedimnames)	
	} 
	

	Sprob = c(1, 0.999, 0.99, 0.95, 0.9, 0.8, 0.5, 0)
	Squant = qchisq(Sprob, df=length(varying))

  # get Scol from RColorBrewer
#	dput( rev(
#			RColorBrewer::brewer.pal(
#					length(Scontour)-1,
#					'Spectral')	))
	Scol=c("#3288BD", "#99D594", 
		"#E6F598", "#FFFFBF", "#FEE08B", "#FC8D59", 
		"#D53E4F")
	names(Scol) = as.character(Sprob[-length(Sprob)])
	
	res = list(
		logL=-L/2,
		legend = list(breaks=Sprob, col=Scol),
		prob = Sprob,
		col=Scol,
		MLE=fit$param[varying],
		basepars=baseParams
	)
	
	res$maxLogL = max(res$logL)
	res$breaks= res$maxLogL - Squant/2 	
	res$breaks[1] = min(c(res$breaks[2],min(res$logL)))-1

	res = c(dots[varying],res)
	
	if(length(varying)==1) {
		Skeep = seq(2, length(Sprob)-1)
		res$ci = cbind(
			prob=Sprob[Skeep],
			lower=NA, upper=NA)
		
		
		smaller = seq(1,which.max(res$logL))
		bigger = seq(which.max(res$logL), length(res$logL))
		
    # make the likelihood unimodal
		monotoneLik = res$logL
		
		if(length(smaller)>1) {
			for(D in 2:length(smaller)){
				monotoneLik[D] = max(monotoneLik[c(D-1,D)])
			}
			
			resCi = try(stats::spline(
				x=monotoneLik[smaller], 
				y=res[[1]][smaller],
				xout=res$breaks[Skeep],
				method= "hyman"), silent=TRUE)
			
			if(all(class(resCi)!='try-error'))
				res$ci[,'lower'] = resCi$y
		}

		if(sum(bigger)>1) {
			for(D in seq(length(smaller)+1, length(monotoneLik)-1)){
				monotoneLik[D] = min(monotoneLik[c(D-1,D)])
			}

			resCi = try(stats::spline(
				x=monotoneLik[bigger], 
				y=res[[1]][bigger], 
				xout=res$breaks[Skeep], method='hyman'),
			silent=TRUE)
			
			if(all(class(resCi)!='try-error'))
				res$ci[,'upper']= resCi$y     
		}
		
		res$ci = cbind(
			res$ci, 
			logL = res$breaks[match(res$ci[,'prob'], res$prob)]
		)
		res$ciLong = na.omit(
			reshape(as.data.frame(res$ci[,c('prob','upper','lower')]), 
				direction="long",
				varying=list(par=c('upper','lower')),
				v.names='par',
				times = c('upper','lower'),
				timevar=c('direction'),
				idvar='prob')	
		)
		res$ciLong = rbind(
			res$ciLong,
			data.frame(prob=0,direction='upper',
				par=res$MLE)[
			colnames(res$ciLong)
			])
		res$ciLong$quantile = (1-res$ciLong$prob)/2
		res$ciLong[res$ciLong$direction=='upper','quantile'] =
		1 - res$ciLong[res$ciLong$direction=='upper','quantile'] 
		
		
		res$ciLong = res$ciLong[order(res$ciLong$par),]
	}

	res
}