R Under development (unstable) (2024-02-19 r85946 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> havePackages = c(
+     'INLA' = requireNamespace('INLA', quietly=TRUE)
+ )
> 
> if(requireNamespace("INLA", quietly=TRUE) ) {
+   INLA::inla.setOption(num.threads=2)
+   # not all versions of INLA support blas.num.threads
+   try(INLA::inla.setOption(blas.num.threads=2), silent=TRUE)
+ }
> 
> print(havePackages)
 INLA 
FALSE 
> 
> # as in example
> require('geostatsp')
Loading required package: geostatsp
Loading required package: Matrix
Loading required package: terra
terra 1.7.71
> myPoints = vect(cbind(rbeta(100,2,2), rbeta(100,3,4)))
> 
> mycov = rast(matrix(rbinom(100, 1, 0.5), 10, 10), extent=ext(0, 1, 0, 1))
> names(mycov)="x1"
> 
> if(all(havePackages)) {
+   
+ res = lgcp(
+     formula=~factor(x1),
+     data=myPoints, 
+     grid=20, 
+     covariates=mycov,
+ 		priorCI=list(sd=c(u=0.1, alpha = 0.01), range=c(0.4, 0.41))
+ )
+ 
+ if(length(res$parameters)) {
+ knitr::kable(res$parameters$summary, digits=3)
+ 
+ plot(res$raster[["predict.exp"]])
+ plot(myPoints,add=TRUE,col="#0000FF30",cex=0.5)
+ }
+ 	
+ # intercept only
+ 
+ res = lgcp(
+     data=myPoints, 
+     grid=20, 
+     covariates=mycov,
+ 		formula=~1,
+ 		priorCI=list(sd=c(0.9, 1.1), range=c(0.4, 0.41))
+ )
+ 
+ if(length(res$parameters)) {
+ 
+ knitr::kable(res$parameters$summary, digits=3)
+ 
+ plot(res$raster[["predict.exp"]])
+ plot(myPoints,add=TRUE,col="#0000FF30",cex=0.5)
+ }
+ }
> 
> 
> # some missing values
> mycov2 = deepcopy(mycov)
> temp = values(mycov2)
> temp[1:4] = NA
> values(mycov2) = temp
> 
> if(all(havePackages)) {
+ 
+ res = lgcp(data=myPoints, grid=20, covariates=mycov2,
+ 		formula=~factor(x1),
+ 		priorCI=list(sd=c(0.9, 1.1), range=c(0.4, 0.41))
+ )
+ 
+ if(length(res$parameters)) {
+ 
+ plot(res$raster[["predict.exp"]])
+ plot(myPoints,add=TRUE,col="#0000FF30",cex=0.5)
+ }
+ }
> 
> 
> data('murder')
> murder = unwrap(murder)
> data('torontoPop')
> torontoPdens = unwrap(torontoPdens)
> torontoIncome = unwrap(torontoIncome)
> torontoBorder = unwrap(torontoBorder)
> 
> myCov = list(
+     pop=torontoPdens,
+     inc = log(torontoIncome)
+ )
> 
> formula = ~ inc + offset(pop, log=TRUE)
> 
> if(all(havePackages)) {
+ 
+ resL=lgcp(formula, data=murder, 
+     grid=squareRaster(murder, 30),
+     covariates=myCov,
+     border=torontoBorder)
+ 
+ resO = lgcp( ~ inc + pop, 
+     data=murder, 
+     grid=squareRaster(murder, 30),
+     covariates=list(inc=myCov$inc, pop=log(myCov$pop)),
+     border=torontoBorder)
+ 
+ if(length(resL$parameters)) {
+ 
+ 	rbind(resL$param$summary, resO$param$summary)
+ }
+ 
+ 
+ }	
> 
> 
> # check spdfToBrock
> 
> if(requireNamespace('diseasemapping', quietly=TRUE)){
+ 	require('diseasemapping')
+ 	
+ 	data('kentucky')
+ 	kentucky = unwrap(kentucky)
+ 	
+ 	popList = list(
+ 			'2002' = kentucky[,c('M.50', 'M.55', 'M.60')],
+ 			'2005' = kentucky[,c('F.50', 'F.55', 'F.60')]
+ 			)
+ 	for(D in names(popList))
+ 		names(popList[[D]]) = paste('expected_', seq(as.numeric(D)-1, len=3), sep='')
+ 	
+ 	popBrick = spdfToBrick(
+ 			x=popList,
+ 			template=squareRaster(kentucky, 10),
+     	logSumExpected=FALSE
+ 			)
+ 	sum(popList[['2002']]$expected_2001, na.rm=TRUE)		
+ 	sum(values(popBrick[['expected_2001']]), na.rm=TRUE)*prod(res(popBrick))
+ 	
+ 	popBrick2 = spdfToBrick(
+ 			x=popList,
+ 			template=squareRaster(kentucky, 10),
+     	logSumExpected=TRUE
+ 	)
+ 	
+ 	sum(unlist(lapply(popList, function(qq) apply(values(qq), 2, sum, na.rm=TRUE))))
+ 	sum(exp(values(popBrick2)), na.rm=TRUE)*prod(res(popBrick2))
+ 			
+ }
> 	
> 
> 
> 
> proc.time()
   user  system elapsed 
   4.54    0.20    4.75