%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{LGCP with PC priors}

\documentclass[12pt]{article}
\usepackage[margin=1in]{geometry}
\usepackage{caption}
\usepackage{subcaption}
\providecommand{\subfloat}[2][need a sub-caption]{\subcaptionbox{#1}{#2}}
\title{LGCP with PC priors}
\author{Patrick Brown}

\begin{document}

\maketitle

<<knitr, include=FALSE>>=
require('knitr')
opts_chunk$set(out.width='0.48\\textwidth', fig.align='default', fig.height=3, fig.width=6)
@    

<<packages>>=
library('mapmisc')
library("geostatsp")
data('murder')
data('torontoPop')
@



<<Covariates, tidy=TRUE>>=
murderT = spTransform(murder, omerc(murder, angle=-20))
covList = list(
    pop=torontoPdens,
    inc = log(torontoIncome)
)

borderT = spTransform(torontoBorder, projection(murderT))
borderC = crop(borderT, extent(-12700, 7000, -7500, 3100))

formula = ~ inc + offset(pop, log=TRUE)
@


\section*{LGCP with gamma priors on precision}

<<lgcpGamma, tidy=TRUE>>=
resG=lgcp(formula, data=murderT, 
    grid=squareRaster(borderC, 30),
    covariates=covList,
    border=borderC, buffer=2000,
	priorCI = list(sd = c(0.05, 0.5), range = c(1, 10)*1000)
)
@





\section*{LGCP with penalised complexity prior}

$pr(sd > 1) = 0.05$ and $pr(phi < 0.2) = 0.95$

<<lgcpPc, tidy=TRUE>>=
resP=lgcp(formula, data=murderT, 
    grid=squareRaster(borderC, 30),
    covariates=covList,
    border=borderC, buffer=2000,
	priorCI = list(sd = c(u=0.5, alpha=0.05), range = c(1, 10)*1000)
)
@




\section*{Parameters}

<<summary>>=
if(!is.null(resG$parameters))
	knitr::kable(resG$parameters$summary[,c(1,3,5)], digits=3)
@

<<summaryPc>>=
if(!is.null(resP$parameters))
	knitr::kable(resP$parameters$summary[,c(1,3,5)], digits=3)
@


<<priorPost, fig.cap="gamma prior", fig.height=4, fig.width=3, fig.subcap=c("sd", "range"), echo=FALSE>>=

if(!is.null(resG$parameters)) {
	
plot(resG$parameters$sd$posterior, type='l', 
		xlim=c(0,1), col='red')
lines(resG$parameters$sd$prior, col='blue')
lines(resP$parameters$sd$posterior, col='red', lty=2)
lines(resP$parameters$sd$prior, col='blue', lty=2)


plot(resG$parameters$range$posterior, type='l', 
		xlim=c(0,10000), col='red')
lines(resG$parameters$range$prior, col='blue')
lines(resP$parameters$range$posterior, col='red', lty=2)
lines(resP$parameters$range$prior, col='blue', lty=2)
legend('topright', lty=c(1,1,1,2), 
		col=c('red','blue', 'black','black'), 
		legend=c('posterior','prior', 'Gamma','PC'))

} else {
	plot(1:10)
	plot(1:10)
}
@



\section*{Maps}

<<maps, fig.cap='Random effects and fitted values', fig.subcap=c('gamma, fitted','pc fitted','gamma random','pc random'), echo=FALSE>>=

if(require('mapmisc', quietly=TRUE) & !is.null(resG$raster)) {
	
thecex=1.2	
	
colFit = colourScale(resG$raster[['predict.exp']],
		breaks=6, dec=8, style='equal',
		transform='log')


map.new(resG$raster, TRUE)
plot(resG$raster[['predict.exp']], 
		col=colFit$col,breaks=colFit$breaks,
		legend=FALSE, add=TRUE)
points(murderT, col='#00FF0050',cex=0.2)
legendBreaks('bottomright', colFit, cex=thecex)


map.new(resG$raster, TRUE)
plot(resP$raster[['predict.exp']], 
		col=colFit$col,breaks=colFit$breaks,
		legend=FALSE, add=TRUE)
points(murderT, col='#00FFFF50',cex=0.2)
legendBreaks('bottomright', colFit, cex=thecex)


colR = colourScale(resG$raster[['random.mean']],
		breaks=12, dec=0, style='equal')

map.new(resG$raster, TRUE)
plot(resG$raster[['random.mean']], 
		col=colR$col,breaks=colR$breaks,
		legend=FALSE, add=TRUE)
legendBreaks('bottomright', colR, cex=thecex)

map.new(resG$raster, TRUE)
plot(resP$raster[['random.mean']], 
		col=colR$col,breaks=colR$breaks,
		legend=FALSE, add=TRUE)
legendBreaks('bottomright', colR, cex=thecex)


} else {
	plot(1:10)
	plot(1:10)
	plot(1:10)
	plot(1:10)
	
}
@


\end{document}
  