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))
	
	
	resL = parallel::mcmapply(
		likfitLgm, 
		param=parList,
		MoreArgs=list(
			data=as.data.frame(fit$data), 
			formula=fit$model$formula,
			paramToEstimate=reEstimate,
			reml=fit$model$reml,
			coordinates=terra::crds(fit$data)
		),
		mc.cores=mc.cores,
		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
}