## package library("lgcp") ############ ## aegiss ## ############ # load data load("aegiss.RData") load("popshape_aegiss.RData") # minimum contrast estimation minc <- minimum.contrast.spatiotemporal(xyt, model = "exponential", method = "g", transform = log, spatial.dens = density.ppp(xyt), temporal.intens = muEst(xyt)) minc # select cell width chooseCellwidth(xyt, cwinit = 1400) CELLWIDTH <- 1400 # select extension EXT <- 2 # perform polygon overlay operations and compute computational grid polyolay <- getpolyol(data = xyt, regionalcovariates = popshape, cellwidth = CELLWIDTH, ext = EXT) # the statistical model for the main effets, $\beta$ FORM <- X ~ pop + IMD + dotw FORM.spatial <- X ~ pop + IMD TEMPORAL.FORMULA <- t ~ dotw - 1 # set the interpolation type for each variable popshape@data <- guessinterp(popshape@data) popshape@data <- assigninterp(df = popshape@data, vars = c("pop", "males", "females"), value = "ArealWeightedSum") class(popshape@data$pop) # perform interpolation of the covariates onto the Zmat <- getZmat(formula = FORM.spatial, data = xyt, regionalcovariates = popshape, cellwidth = CELLWIDTH, ext = EXT, overl = polyolay) Zmat[, "pop"] <- log(Zmat[, "pop"]) Zmat[, "pop"][is.infinite(Zmat[, "pop"])] <- min(Zmat[, "pop"][!is.infinite(Zmat[, "pop"])]) # plot the spatial interpolated covariates plot(Zmat, ask = FALSE) # construct dummy temporal data frame days <- c("Mon", "Tue", "Wed", "Thur", "Fri", "Sat", "Sun") tvec <- xyt$tlim[1]:xyt$tlim[2] da <- rep(days, length.out = length(tvec)) tdata <- data.frame(t = tvec, dotw = da) # choose last time point and number of proceeding time-points to include T <- 595 LAGLENGTH <- 13 # bolt on the temporal covariates ZmatList <- addTemporalCovariates(temporal.formula = TEMPORAL.FORMULA, T = T, laglength = LAGLENGTH, tdata = tdata, Zmat = Zmat) # define the priors lgprior <- PriorSpec(LogGaussianPrior(mean = log(c(1, 2000, 1)), variance = diag(0.2, 3))) gprior <- PriorSpec(GaussianPrior(mean = rep(0, 9), variance = diag(1e+06, 9))) priors <- lgcpPrior(etaprior = lgprior, betaprior = gprior) # set initial values for the algorithm INITS <- lgcpInits(etainit = log(c(sqrt(2.8), 1286, 0.5)), betainit = NULL) # choose the covariance function CF <- CovFunction(exponentialCovFct) # run the MCMC algorithm, note that as this takes a long time, we present a # shorter run here in the article ,we used mala.length=1000000, burnin=100000, # thin=900 DIRNAME <- getwd() lg <- lgcpPredictSpatioTemporalPlusPars(formula = FORM, xyt = xyt, T = T, laglength = LAGLENGTH, ZmatList = ZmatList, model.priors = priors, model.inits = INITS, spatial.covmodel = CF, cellwidth = CELLWIDTH, mcmc.control = mcmcpars(mala.length = 1000, burnin = 100, retain = 9, adaptivescheme = andrieuthomsh(inith = 1, alpha = 0.5, C = 1, targetacceptance = 0.574)), output.control = setoutput(gridfunction = dump2dir(dirname = file.path(DIRNAME, "aegiss"), forceSave = TRUE)), ext = EXT) save(list = ls(), file = file.path(DIRNAME, "aegiss", "aegissout.RData")) # plot the log target plot(ltar(lg), type = "s", xlab = "Iteration/900", ylab = "log target", ask = FALSE) # compute and plot autocorrelations in the latent field lagch <- c(1, 5, 15) Sacf <- autocorr(lg, lagch, inWindow = NULL) library("fields") for (i in 1:3) { image.plot(xvals(lg), yvals(lg), Sacf[, , i], zlim = c(-1, 1), axes = FALSE, xlab = "", ylab = "", asp = 1, sub = paste("Lag:", lagch[i]), ask = FALSE) plot(xyt$window, add = TRUE, ask = FALSE) scalebar(5000, label = "5 km") } # produce traceplots of beta and eta traceplots(lg, ask = FALSE) # produce autocorrelation plots of beta and eta parautocorr(lg, ask = FALSE) # a summary table of beta and eta parsum <- parsummary(lg) # the above converted to LaTeX library("miscFuncs") latextable(parsum, rownames = rownames(parsum), colnames = c("Parameter", colnames(parsum)), digits = 4) # a text summary of the model parameters textsummary(lg, digits = 4) # a plot of the prior and posterior densities priorpost(lg, ask = FALSE) # the posterior covariance function postcov(lg, ask = FALSE) # exceedance and lower-tail exceedance probabilities ep <- exceedProbs(c(1.5, 3)) ex <- lgcp:::expectation.lgcpPredict(lg, ep) plotExceed(ex[[1]], "ep", lg, zlim = c(0, 1), asp = 1, axes = FALSE, xlab = "", ylab = "", sub = "", ask = FALSE) scalebar(5000, label = "5 km") ######### ## BTB ## ######### # load data load("BTBppp.RData") load("farmspdf.RData") # find highly correlated variables library("caret") d <- farmspdf@data[, 30:40] d <- d[!is.na(d[, 1]), ] findCorrelation(cor(d)) # minimum contrast estimation W <- pppdata$window simpW <- simplify.owin(W, 1000) ttx <- pppdata$x[pppdata$marks == "x1" | pppdata$marks == "x4" | pppdata$marks == "x5" | pppdata$marks == "x7"] tty <- pppdata$y[pppdata$marks == "x1" | pppdata$marks == "x4" | pppdata$marks == "x5" | pppdata$marks == "x7"] ttm <- as.character(pppdata$marks[pppdata$marks == "x1" | pppdata$marks == "x4" | pppdata$marks == "x5" | pppdata$marks == "x7"]) tempppp <- ppp(x = ttx, y = tty, window = simpW, marks = as.factor(ttm)) denls <- lapply(as.character(levels(tempppp$marks)), function(x) { density.ppp(tempppp[tempppp$marks == x]) }) mn <- minimum.contrast(tempppp, model = "exponential", method = "g", intens = denls, transform = log) # select cell width chooseCellwidth(pppdata, cwinit = 1800) CELLWIDTH <- 1800 # select extension EXT <- 2 # perform polygon overlay operations and compute computational grid polyolay <- getpolyol(data = pppdata, pixelcovariates = farmspdf, cellwidth = CELLWIDTH, ext = EXT) # the statistical model for the main effets, $\beta$ fl <- formulaList(list(x1 ~ K207 + K208 + K209 + K210, x4 ~ K207 + K208 + K209 + K210, x5 ~ K207 + K208 + K209 + K210, x7 ~ K207 + K208 + K209 + K210)) # specify formula for purposes of interpolation FORM <- X ~ K207 + K208 + K209 + K210 # interpolation has already been set in farmspdf.RData # perform interpolation of the covariates onto the Zmat <- getZmat(formula = FORM, data = pppdata, pixelcovariates = farmspdf, cellwidth = CELLWIDTH, ext = EXT, overl = polyolay) Zmat[, "K207"] <- log(Zmat[, "K207"]) Zmat[, "K207"][is.infinite(Zmat[, "K207"])] <- min(Zmat[, "K207"][!is.infinite(Zmat[, "K207"])]) Zmat[, "K208"] <- log(Zmat[, "K208"]) Zmat[, "K208"][is.infinite(Zmat[, "K208"])] <- min(Zmat[, "K208"][!is.infinite(Zmat[, "K208"])]) Zmat[, "K209"] <- log(Zmat[, "K209"]) Zmat[, "K209"][is.infinite(Zmat[, "K209"])] <- min(Zmat[, "K209"][!is.infinite(Zmat[, "K209"])]) Zmat[, "K210"] <- log(Zmat[, "K210"]) Zmat[, "K210"][is.infinite(Zmat[, "K210"])] <- min(Zmat[, "K210"][!is.infinite(Zmat[, "K210"])]) # plot the interpolated covariates plot(Zmat, ask = FALSE) # define the priors pr.mn <- log(c(1, 1500)) pr.vr <- c(0.2, 0.05) priors <- list() for (i in 1:4) { priors[[i]] <- lgcpPrior(etaprior = PriorSpec(LogGaussianPrior(mean = pr.mn, variance = diag(pr.vr))), betaprior = PriorSpec(GaussianPrior(mean = rep(0, 5), variance = diag(10^6, 5)))) } priors[[5]] <- lgcpPrior(etaprior = PriorSpec(LogGaussianPrior(mean = pr.mn, variance = diag(pr.vr))), betaprior = NULL) # set initial values were not set for this example # choose the covariance function cfs <- list() for (i in 1:4) { cfs[[i]] <- CovFunction(exponentialCovFct) } cfs[[5]] <- CovFunction(RandomFieldsCovFct(model = "matern", additionalparameters = 1)) # run the MCMC algorithm, note that as this takes a long time, we present a # shorter run here in the article ,we used mala.length=1000000, burnin=100000, # thin=900 BASEDR <- getwd() lg <- lgcpPredictMultitypeSpatialPlusPars(formulaList = fl, sd = pppdata, Zmat = Zmat, model.priorsList = priors, spatial.covmodelList = cfs, cellwidth = CELLWIDTH, mcmc.control = mcmcpars(mala.length = 1000, burnin = 100, retain = 9, adaptivescheme = andrieuthomsh(inith = 1, alpha = 0.5, C = 1, targetacceptance = 0.574)), output.control = setoutput(gridfunction = dump2dir(dirname = file.path(BASEDR, "BTB"), forceSave = TRUE)), ext = EXT) save(list = ls(), file = file.path(BASEDR, "BTB", "BTB.RData")) # plot the log target plot(ltar(lg), type = "s", xlab = "Iteration/900", ylab = "log target") # compute and plot autocorrelations in the latent field for (i in 1:4) { Yi <- autocorrMultitype(lg, c(1, 5, 15), i, inWindow = NULL) plot(Yi, zlim = c(-1, 1), xlab = "", ylab = "", ask = FALSE) } # produce traceplots of beta and eta traceplots(lg, ask = FALSE) # produce autocorrelation plots of beta and eta parautocorr(lg, ask = FALSE) # a summary table of beta and eta parsum <- parsummary(lg) # the above converted to LaTeX library("miscFuncs") latextable(parsum, rownames = rownames(parsum), colnames = c("Parameter", colnames(parsum)), digits = 4) # a plot of the prior and posterior densities priorpost(lg, ask = FALSE) # the posterior covariance function postcov(lg, ask = FALSE) # conditional probabilities cp <- condProbs(lg) # segregation probabilities sr <- segProbs(lg, domprob = 0.8) ########### ## liver ## ########### ### aggregated # load data load("liver_spdf.RData") load("popshape_liver.RData") spdf # select cell width chooseCellwidth(spdf, cwinit = 300) CELLWIDTH <- 300 # select extension EXT <- 3 # perform polygon overlay operations and compute computational grid polyolay <- getpolyol(data = spdf, regionalcovariates = popshape, cellwidth = CELLWIDTH, ext = EXT) # the statistical model for the main effets, $\beta$ FORM <- X ~ pop + propmale + Income + Employment + Education + Barriers + Crime + Environment # set the interpolation type for each variable popshape@data <- guessinterp(popshape@data) popshape@data <- assigninterp(df = popshape@data, vars = c("pop", "males", "females"), value = "ArealWeightedSum") class(popshape@data$pop) # perform interpolation of the covariates onto the Zmat <- getZmat(formula = FORM, data = spdf, regionalcovariates = popshape, cellwidth = CELLWIDTH, ext = EXT, overl = polyolay) Zmat[, "pop"] <- log(Zmat[, "pop"]) Zmat[, "pop"][is.infinite(Zmat[, "pop"])] <- min(Zmat[, "pop"][!is.infinite(Zmat[, "pop"])]) # plot the interpolated covariates plot(Zmat, ask = FALSE) # define the priors priors <- lgcpPrior(etaprior = PriorSpec(LogGaussianPrior(mean = log(c(1, 500)), variance = diag(0.15, 2))), betaprior = PriorSpec(GaussianPrior(mean = rep(0, 9), variance = diag(10^6, 9)))) # no initial values for this algorithm # choose the covariance function cf <- CovFunction(exponentialCovFct) # run the MCMC algorithm, note that as this takes a long time, we present a # shorter run here in the article ,we used mala.length=3100000, burnin=100000, # thin=3000 BASEDR <- getwd() lg <- lgcpPredictAggregateSpatialPlusPars(formula = FORM, spdf = spdf, Zmat = Zmat, overlayInZmat = FALSE, model.priors = priors, spatial.covmodel = cf, cellwidth = CELLWIDTH, mcmc.control = mcmcpars(mala.length = 1000, burnin = 100, retain = 9, adaptivescheme = andrieuthomsh(inith = 1, alpha = 0.5, C = 1, targetacceptance = 0.574)), output.control = setoutput(gridfunction = dump2dir(dirname = paste(BASEDR, "/liveraggregated/", sep = ""), forceSave = TRUE)), ext = EXT) save(list = ls(), file = file.path(BASEDR, "liveraggregated", "liveraggregated.RData")) # plot the log target plot(ltar(lg), type = "s", xlab = "Iteration/3000", ylab = "log target", ask = FALSE) # compute and plot autocorrelations in the latent field lagch <- c(1, 5, 15) Sacf <- autocorr(lg, lagch, inWindow = NULL) library("fields") for (i in 1:3) { image.plot(xvals(lg), yvals(lg), Sacf[, , i], zlim = c(-1, 1), axes = FALSE, xlab = "", ylab = "", asp = 1, sub = paste("Lag:", lagch[i]), ask = FALSE) plot(spdf, add = TRUE, ask = FALSE) scalebar(5000, label = "5 km") } # produce traceplots of beta and eta traceplots(lg, ask = FALSE) # produce autocorrelation plots of beta and eta parautocorr(lg, ask = FALSE) # a summary table of beta and eta parsum <- parsummary(lg) # the above converted to LaTeX library("miscFuncs") latextable(parsum, rownames = rownames(parsum), colnames = c("Parameter", colnames(parsum)), digits = 4) # a text summary of the model parameters textsummary(lg, digits = 4) # a plot of the prior and posterior densities priorpost(lg, ask = FALSE) # the posterior covariance function postcov(lg, ask = FALSE) # exceedance and lower-tail exceedance probabilities ep <- exceedProbs(c(1.5, 2, 5, 10)) sp <- exceedProbs(c(2/3, 1/2, 1/5, 1/10), direction = "lower") ex <- lgcp:::expectation.lgcpPredict(lg, ep) su <- lgcp:::expectation.lgcpPredict(lg, sp) plotExceed(ex[[1]], "ep", lg, zlim = c(0, 1), asp = 1, axes = FALSE, xlab = "", ylab = "", sub = "", ask = FALSE) scalebar(5000, label = "5 km") plotExceed(su[[1]], "sp", lg, zlim = c(0, 1), asp = 1, axes = FALSE, xlab = "", ylab = "", sub = "", ask = FALSE) scalebar(5000, label = "5 km") ### spatial point process # load data load("sd_liver.RData") load("popshape_liver.RData") # minimum contrast estimation minimum.contrast(sd, model = "exponential", method = "g", intens = density(sd), transform = log) # select cell width chooseCellwidth(sd, cwinit = 300) CELLWIDTH <- 300 # select extension EXT <- 2 # perform polygon overlay operations and compute computational grid polyolay <- getpolyol(data = sd, regionalcovariates = popshape, cellwidth = CELLWIDTH, ext = EXT) # the statistical model for the main effets, $\beta$ FORM <- X ~ pop + propmale + Income + Employment + Education + Barriers + Crime + Environment # set the interpolation type for each variable popshape@data <- guessinterp(popshape@data) popshape@data <- assigninterp(df = popshape@data, vars = c("pop", "males", "females"), value = "ArealWeightedSum") class(popshape@data$pop) # perform interpolation of the covariates onto the Zmat <- getZmat(formula = FORM, data = sd, regionalcovariates = popshape, cellwidth = CELLWIDTH, ext = EXT, overl = polyolay) Zmat[, "pop"] <- log(Zmat[, "pop"]) Zmat[, "pop"][is.infinite(Zmat[, "pop"])] <- min(Zmat[, "pop"][!is.infinite(Zmat[, "pop"])]) # plot the interpolated covariates plot(Zmat, ask = FALSE) # define the priors priors <- lgcpPrior(etaprior = PriorSpec(LogGaussianPrior(mean = log(c(1, 500)), variance = diag(0.15, 2))), betaprior = PriorSpec(GaussianPrior(mean = rep(0, 9), variance = diag(10^6, 9)))) # set initial values for the algorithm INITS <- lgcpInits(etainit = log(c(sqrt(1.5), 275)), betainit = NULL) # choose the covariance function cf <- CovFunction(exponentialCovFct) # run the MCMC algorithm, note that this takes a long time, we present a shorter # run here in the article ,we used mala.length=1000000, burnin=100000, thin=900 BASEDR <- getwd() lg <- lgcpPredictSpatialPlusPars(formula = FORM, sd = sd, Zmat = Zmat, model.priors = priors, model.inits = INITS, spatial.covmodel = cf, cellwidth = CELLWIDTH, poisson.offset = NULL, mcmc.control = mcmcpars(mala.length = 1000, burnin = 100, retain = 9, adaptivescheme = andrieuthomsh(inith = 1, alpha = 0.5, C = 1, targetacceptance = 0.574)), output.control = setoutput(gridfunction = dump2dir(dirname = paste(BASEDR, "/liver/", sep = ""), forceSave = TRUE)), ext = EXT) save(list = ls(), file = file.path(BASEDR, "liver", "liver.RData")) # plot the log target plot(ltar(lg), type = "s", xlab = "Iteration/900", ylab = "log target", ask = FALSE) # compute and plot autocorrelations in the latent field lagch <- c(1, 5, 15) Sacf <- autocorr(lg, lagch, inWindow = NULL) library("fields") for (i in 1:3) { image.plot(xvals(lg), yvals(lg), Sacf[, , i], zlim = c(-1, 1), axes = FALSE, xlab = "", ylab = "", asp = 1, sub = paste("Lag:", lagch[i]), ask = FALSE) plot(sd$window, add = TRUE, ask = FALSE) scalebar(5000, label = "5 km") } # produce traceplots of beta and eta traceplots(lg, ask = FALSE) # produce autocorrelation plots of beta and eta parautocorr(lg, ask = FALSE) # a summary table of beta and eta parsum <- parsummary(lg) # the above converted to LaTeX library("miscFuncs") latextable(parsum, rownames = rownames(parsum), colnames = c("Parameter", colnames(parsum)), digits = 4) # a text summary of the model parameters textsummary(lg, digits = 4) # a plot of the prior and posterior densities priorpost(lg, ask = FALSE) # the posterior covariance function postcov(lg, ask = FALSE) # exceedance and lower-tail exceedance probabilities ep <- exceedProbs(c(1.5, 2, 5, 10)) sp <- exceedProbs(c(2/3, 1/2, 1/5, 1/10), direction = "lower") ex <- lgcp:::expectation.lgcpPredict(lg, ep) su <- lgcp:::expectation.lgcpPredict(lg, sp) plotExceed(ex[[1]], "ep", lg, zlim = c(0, 1), asp = 1, axes = FALSE, xlab = "", ylab = "", sub = "", ask = FALSE) scalebar(5000, label = "5 km") plotExceed(su[[1]], "sp", lg, zlim = c(0, 1), asp = 1, axes = FALSE, xlab = "", ylab = "", sub = "", ask = FALSE) scalebar(5000, label = "5 km") ############### ## Figure 19 ## ############### library("OpenStreetMap") ncmap <- openmap(upperLeft = c(lat = 50.94028, lon = -5.760212), lowerRight = c(lat = 49.94531, lon = -4.145465), type = "stamen-toner") sdmc <- osppp2merc(pppdata) cprob <- as.array(condProbs(lg)) for (i in 1:4) { fna <- paste("conditionalprobS_", i, sep = "") ra <- raster(lgcpgrid(cprob[, , i], xvals(lg), yvals(lg)), crs = CRS("+init=epsg:27700")) xdiff <- diff(lg$mcens[1:2]) ydiff <- diff(lg$ncens[1:2]) bb <- SpatialPoints(matrix(c(lg$mcens[1] - xdiff/2, lg$mcens[64] + xdiff/2, lg$ncens[64] + ydiff/2, lg$ncens[1] - ydiff/2), 2, 2), proj4string = CRS("+init=epsg:27700")) bb <- coordinates(spTransform(bb, CRS("+init=epsg:3857"))) probri <- raster(nrows = 64, ncols = 64, xmn = bb[1, 1], xmx = bb[2, 1], ymn = bb[2, 2], ymx = bb[1, 2], crs = CRS("+init=epsg:3857")) tp <- projectRaster(ra, probri) plot(sdmc$window, main = "", lwd = 0.5) plot(ncmap, add = TRUE) cpal <- rev(heat.colors(64, alpha = 0.5)) if (i < 4) { plot(tp, add = TRUE, legend = FALSE, col = cpal, zlim = c(0, 1)) } else { plot(tp, add = TRUE, col = cpal, zlim = c(0, 1), legend.mar = 3.1) } } ############### ## Figure 20 ## ############### thresh <- c(0.6, 0.7, 0.8, 0.9) nthresh <- length(thresh) nits <- nrow(lg$betarec) domprob <- 0.8 library("OpenStreetMap") ncmap <- openmap(upperLeft = c(lat = 50.94028, lon = -5.760212), lowerRight = c(lat = 49.94531, lon = -4.145465), type = "stamen-toner") sdmc <- osppp2merc(pppdata) segrec <- as.array(segProbs(lg, domprob = domprob)) segexceed <- array(NA, dim = c(lg$M, lg$N, 4, nthresh)) for (k in 1:4) { for (l in 1:nthresh) { segexceed[, , k, l] <- segrec[, , k] > thresh[l] } } cols <- c(transblue(0.5), transgreen(0.5), transblack(0.5), transred(0.5)) for (l in 1:nthresh) { for (k in 1:4) { ra <- raster(lgcpgrid(matrix(as.numeric(segexceed[, , k, l]) * k, lg$M, lg$N), xvals(lg), yvals(lg)), crs = CRS("+init=epsg:27700")) xdiff <- diff(lg$mcens[1:2]) ydiff <- diff(lg$ncens[1:2]) bb <- SpatialPoints(matrix(c(lg$mcens[1] - xdiff/2, lg$mcens[64] + xdiff/2, lg$ncens[64] + ydiff/2, lg$ncens[1] - ydiff/2), 2, 2), proj4string = CRS("+init=epsg:27700")) bb <- coordinates(spTransform(bb, CRS("+init=epsg:3857"))) probri <- raster(nrows = 64, ncols = 64, xmn = bb[1, 1], xmx = bb[2, 1], ymn = bb[2, 2], ymx = bb[1, 2], crs = CRS("+init=epsg:3857")) tp <- projectRaster(ra, probri) assign(paste("lay", k, sep = ""), tp) } bri <- brick(lay1, lay2, lay3, lay4) plot(sdmc$window, main = "", lwd = 0.5) plot(ncmap, add = TRUE) for (k in 1:4) { laytemp <- get(paste("lay", k, sep = "")) laytemp[laytemp == 0] <- NA plot(laytemp, add = TRUE, legend = FALSE, col = cols[k]) } if (l == nthresh) { legend(-625721.6, 6593320, pch = rep(15, 4), col = cols, legend = paste("Genotype", c(9, 12, 15, 20)), bg = "white") } }