library("lgcp") INTERACTIVE <- FALSE # set to TRUE to use interactive tools in package load("v52i04.rda") x <- xyt$x y <- xyt$y t <- xyt$t win <-xyt$window data <- cbind(x,y,t) tlim <- c(0,100) tlim win xyt <- stppp(list(data=data,tlim=tlim,window=win)) xyt <- integerise(xyt) xyt { if (INTERACTIVE){ den <- lambdaEst(xyt,axes=TRUE) plot(den) } else{ den <- density.ppp(xyt,sigma=1) } } sar <- spatialAtRisk(den) { if (INTERACTIVE){ mut1 <- muEst(xyt) mut1 } } mut <- constantInTime(xyt) gin <- ginhomAverage(xyt,spatial.intensity=sar,temporal.intensity=mut) kin <- KinhomAverage(xyt,spatial.intensity=sar,temporal.intensity=mut) { if (INTERACTIVE){ sigmaphi1 <- spatialparsEst(gin,sigma.range=c(0,10),phi.range=c(0,10),spatial.covmodel="exponential") sigmaphi2 <- spatialparsEst(kin,sigma.range=c(0,10),phi.range=c(0,10),spatial.covmodel="exponential") theta <- thetaEst(xyt,spatial.intensity=sar,temporal.intensity=mut,sigma=1.6,phi=1.9) } } args(lgcpPredict) lgcppars(sigma=1.6,phi=1.9,theta=1.4) args(mcmcpars) args(andrieuthomsh) args(setoutput) args(dump2dir) gfun <- function(Y){ return(Y) } args(MonteCarloAverage) exceed <- exceedProbs(c(1.5,2,3)) tmpdr <- tempdir() lg <- lgcpPredict( xyt=xyt, T=50, laglength=5, model.parameters=lgcppars(sigma=1.6,phi=1.9,theta=1.4), cellwidth=2, spatial.intensity=sar, temporal.intensity=mut, mcmc.control=mcmcpars(mala.length=120000,burnin=20000, retain=100, adaptivescheme=andrieuthomsh(inith=1,alpha=0.5,C=1, targetacceptance=0.574)), output.control=setoutput(gridfunction=dump2dir(dirname=tmpdr,forceSave=TRUE), gridmeans=MonteCarloAverage("exceed"))) #save(list=ls(),file="new_results_for_paper.RData") # Wait about 2.77 hours ... # the lg file and associated NetCDF file are LARGE lg rastlg <- raster(intens(lg)) # a raster image of the intensity spdflg <- as.SpatialPixelsDataFrame(rr(lg)) { if (INTERACTIVE){ plot(lg,xlab="x coordinate",ylab="y coordinate") plot(lg,type="serr",xlab="x coordinate",ylab="y coordinate") plot(lg,type="intensity",xlab="x coordinate",ylab="y coordinate") } } { if (INTERACTIVE){ plot(hvals(lg)[20000:120000],type="l",xlab="Iteration",ylab="h") } } tr <- extract(lg,x=6,y=32,t=1,s=-1) { if (INTERACTIVE){ plot(tr,type="l",xlab="Iteration",ylab="Y") } } acor <- autocorr(lg,lags=c(1,5,10),inWindow=NULL) { if (INTERACTIVE){ plot(acor,zlim=c(-1,1),xlab="x-coords",ylab="y-coords") } } fcast <- lgcpForecast(lg,c(51,53,55,60),spatial.intensity=sar,temporal.intensity=function(x){return(100)}) subsamp <- extract(lg,x=c(4,10),y=c(32,35),t=1,s=-1) { if (INTERACTIVE){ win2 <- clickpoly() } else{ win2 <- owin(poly=list(x=c(440,480,480,440),y=c(80,80,120,120))) } } subsamp2 <- extract(lg,inWindow=win2,t=1) { if (INTERACTIVE){ plot(window(lg)) } } ex <- expectation(obj=lg,fun=exceed) qt <- quantile(lg,c(0.5,0.75,0.9),fun=exp) { if (INTERACTIVE){ plot(qt,xlab="X coords",ylab="y coords") plotExceed(lg, fun = "exceed") plotExceed(ex[[1]], fun = "exceed",lgcppredict=lg) } } set.seed(666) ncells <- 128 population <- exp(rgauss(ncells=ncells,model.parameters=lgcppars(sigma=2,phi=0.1))$grid[[1]]) ctsvar <- rgauss(ncells=ncells,model.parameters=lgcppars(sigma=1,phi=0.1))$grid[[1]] binvar <- matrix(as.numeric(rgauss(ncells=ncells,model.parameters=lgcppars(sigma=2,phi=0.05))$grid[[1]]>1),ncells,ncells) xvals <- seq(0,1,length.out=ncells) # centroids in x-direction yvals <- seq(0,1,length.out=ncells) # centroids in y-direction mut <- function(t){ return(exp(-7 + sin(2*pi*t/25)+cos(2*pi*t/25))) } param <- c(1,2) time <- seq(0:100) simd <- matrix(0,ncells,ncells) count <- c() for(t in 1:length(time)){ rate <- mut(time[t])*population*exp(ctsvar*param[1]+binvar*param[2]) dat <- matrix(rpois(ncells^2,rate),ncells,ncells) simd <- simd + dat count <- c(count,sum(dat)) } spatmod <- glm(as.vector(simd)~as.vector(ctsvar)+as.vector(binvar),family="poisson",offset=as.vector(log(population))) print(summary(spatmod)$coefficients) tempmod <- glm(count~sin(2*pi*time/25)+cos(2*pi*time/25),family="poisson") print(summary(tempmod)$coefficients) sar <- spatialAtRisk(list(X=xvals,Y=yvals,Zm=matrix(fitted(spatmod),ncells,ncells))) tar <- temporalAtRisk(fitted(tempmod),tlim=c(0,100),warn=FALSE) # note this a different example from the code in the paper roteffgain(xyt,cellwidth=24) W <- xyt$window tempfun <- function(t){return(100)} sim <- lgcpSim(owin=W, tlim=c(0,100), spatial.intensity=den, temporal.intensity=tempfun, cellwidth = 0.5, model.parameters=lgcppars(sigma=2,phi=5,theta=2))