In villandre/ISMRAanalyses: Analyses for IS-MRA study
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
We start by simulating data.
if ("package:ISMRAanalyses" %in% search()) {
detach("package:ISMRAanalyses", unload = TRUE)
}
library(ISMRAanalyses)
library(MRAinla)
library(sp)
setwd("/home/luc/INLAMRAfiles/INLAMRApaper1/simulations")
# setwd("/home/luc/ISMRAfiles/")
The simulated data are based on satellite imagery data obtained in Maharashtra. We first define a rectangular zone in the longitude-latitude projection. Since MODIS data are in the sinusoidal projection, we project the zone we just defined in that coordinate system. Warnings can be safely ignored: they result from a change in the conventions in referring to projections, cf. details in sp::proj4string.
MaharashtraPolygonEdges <- rbind(c(18.8, 73.2), c(18.8, 73.4), c(18.6, 73.4), c(18.6, 73.2), c(18.8, 73.2))
MaharashtraPolygonEdges <- MaharashtraPolygonEdges[ , 2:1]
MaharashtraPolygon <- sp::SpatialPolygons(Srl = list(sp::Polygons(list(sp::Polygon(coords = MaharashtraPolygonEdges)), ID = "Mumbai")))
raster::crs(MaharashtraPolygon) <- sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
MaharashtraPolygonOtherCRS <- sp::spTransform(MaharashtraPolygon, CRSobj = sp::CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"))
We then import land cover and elevation values from the MODIS and ASTER data files we obtained from EarthData.
folderForEarthDataFiles <- "/store/luc/rawDataFiles/"
landCoverFiles <- list.files(folderForEarthDataFiles, pattern = "MCD*", full.names = TRUE)
landCoverRasterSinusoidal <- produceLandCover(landCoverFiles, regionPolygon = MaharashtraPolygonOtherCRS)
elevationFiles <- list.files(path = folderForEarthDataFiles, pattern = "*dem.tif", full.names = TRUE)
elevationRasterList <- lapply(elevationFiles, raster::raster)
From the imported data, we produce the covariate matrix used in the simulations.
dayOffset <- 121 # Days are listed numerically, from 1 to 366 (leap year). Day 121 corresponds to April 30.
dayRange <- 27:29 # We extract data for May 27 until May 29.
collectionDatesPOSIX <- as.POSIXct(paste("2012-05-", dayRange, sep = ""))
covariateData <- prepareCovariateDataForISMRA(elevationRasterList, landCoverRasterSinusoidal, collectionDatesPOSIX)
saveRDS(covariateData, file = "data/covariatesForSimulationStudy.rds", compress = TRUE)
We are now ready to simulate $100$ vectors from a spatiotemporal Gaussian process.
timepoints <- as.numeric(time(covariateData@time))/(3600*24) # Time is in days
# Parameter values based on real data analysis (cf. table 1 in paper)
# Spatial parameters
spSmoothness <- 1.5
spRange <- 5.66
spScale <- 4.14
# Time parameters
timeSmoothness <- 0.5
timeRange <- 3.60
timeScale <- 1
# Uncorrelated error
errorSD <- 0.5
covFunAlt <- function(spacetimeObject, timepointsVec) {
spatialDistances <- fields::RdistEarth(spacetimeObject@sp@coords, miles=FALSE)
timeDistances <- fields::rdist(timepointsVec, timepointsVec)
MRAinla::maternCov(d = spatialDistances, rho = spRange, smoothness = spSmoothness, scale = spScale) * MRAinla::maternCov(d = timeDistances, scale = timeScale, rho = timeRange, smoothness = timeSmoothness)
}
covarStruct <- covFunAlt(covariateData, timepoints) + errorSD^2 * diag(nrow(covariateData@data))
choleskyDecomp <- t(chol(covarStruct))
set.seed(10)
numReplicates <- 100
fieldValuesWithError <- replicate(n = numReplicates, expr = as.numeric(choleskyDecomp %*% rnorm(nrow(covarStruct))))
We then add fixed effects to the field values, yielding the data we are going to model using the methods we selected for the simulation study. We start by inputting parameter values.
landCoverEffects <- c(
Water = 0,
EvergreenNeedleleaf = 1.09,
EvergreenBroadleaf = 0.979,
DeciduousBroadleaf = 1.181,
MixedForest = 1.087,
ClosedShrublands = 2.097,
OpenShrublands = 1.448,
WoodySavannas = 1.329,
Savannas = 1.393,
Grasslands = 1.506,
PermanentWetlands = 0.705,
Croplands = 1.588,
Urban = 1.623,
CroplandNaturalMosaics = 1.430,
NonVegetated = 0.667)
# Water is the reference (landCover = 0)
referenceLandCover <- "Water"
referenceLandCoverIndex <- match(referenceLandCover, names(landCoverEffects))
timeEffects <- c(
May27 = -0.033,
May28 = -6.319,
May29 = 1.155) + 0.033 # We're moving the reference day from May 25 to May 27, hence + 0.033
elevationEffect <- -0.001
interceptValue <- 40 # This is arbitrary. It should not affect the analyses.
referenceTimepoint <- "May27"
referenceTimepointIndex <- match(referenceTimepoint, names(timeEffects))
fixedEffectsVec <- c(Intercept = interceptValue, landCoverEffects[-referenceLandCoverIndex], elevation = elevationEffect, timeEffects[-referenceTimepointIndex])
We now create the covariate matrix and add covariate effects to the previously-simulated responses.
landCoverFactor <- factor(covariateData@data$landCover, levels = 0:14, labels = names(landCoverEffects))
timepointsFactor <- factor(timepoints, levels = unique(timepoints), labels = names(timeEffects))
covariateMatrix <- cbind(model.matrix(object = ~ landCover, data = data.frame(landCover = landCoverFactor)), elevation = covariateData@data$elevation, model.matrix(object = ~ time, data = data.frame(time = timepointsFactor))[ , -1]) # We remove the intercept, which has already been created for land cover. Ocean and May 27 are reference categories.
colnames(covariateMatrix) <- replace(colnames(covariateMatrix), which(colnames(covariateMatrix) == "(Intercept)"), "Intercept")
colnames(covariateMatrix) <- gsub(pattern = "landCover", replacement = "", x = colnames(covariateMatrix))
colnames(covariateMatrix) <- gsub(pattern = "time", replacement = "", x = colnames(covariateMatrix))
# Getting the simulated values.
observedValues <- fieldValuesWithError + as.vector(covariateMatrix %*% fixedEffectsVec[colnames(covariateMatrix)])
# Fixing the covariate matrix
landCoverTypes <- names(landCoverEffects)[-referenceLandCoverIndex]
presentLandCoverTypes <- landCoverTypes[sapply(landCoverTypes, function(coverType) !all(covariateMatrix[ , coverType] == 0))]
covariateMatrixNoAbsentLandCovers <- cbind(covariateMatrix[ , -match(landCoverTypes, colnames(covariateMatrix))], covariateMatrix[ , presentLandCoverTypes])
simulatedDatasetsList <- list(
responses = observedValues,
covariates = covariateMatrixNoAbsentLandCovers[ , -match("Intercept", colnames(covariateMatrixNoAbsentLandCovers))],
coords = covariateData@sp@coords,
time = timepointsFactor,
timeNumeric = timepoints)
colnames(simulatedDatasetsList$coords) <- c("longitude", "latitude")
saveRDS(simulatedDatasetsList, file = "data/dataForSimulations.rds", compress = TRUE)
As a sanity check, we plot one of the simulated datasets.
numRasterRows <- numRasterCols <- 45 # Approximate, hard to recover because data were collected in the sinusoidal projection, which does not translate to a perfectly-spaced grid of coordinates in the longitude/latitude projection.
observationsToPlot <- which(simulatedDatasetsList$time == "May28")
jpeg("outputFiles/dataExample.jpeg")
fields::quilt.plot(
x = simulatedDatasetsList$coords[ , "x"][observationsToPlot],
y = simulatedDatasetsList$coords[ , "y"][observationsToPlot],
z = simulatedDatasetsList$responses[observationsToPlot, 1],
nx = numRasterRows, ny = numRasterCols)
dev.off()
We now begin the simulations per se.
simulatedDatasetsList <- readRDS("data/dataForSimulations.rds")
simulatedDatasetsList$coords <- as.data.frame(simulatedDatasetsList$coords)
longitudeQuantiles <- quantile(simulatedDatasetsList$coords$longitude, probs = c(0.25, 0.75))
latitudeQuantiles <- quantile(simulatedDatasetsList$coords$latitude, probs = c(0.25, 0.75))
obsIndicesForTraining <-
!((simulatedDatasetsList$time == "May28") &
(simulatedDatasetsList$coords$longitude <= longitudeQuantiles[[2]]) &
(simulatedDatasetsList$coords$longitude >= longitudeQuantiles[[1]]) &
(simulatedDatasetsList$coords$latitude <= latitudeQuantiles[[2]]) &
(simulatedDatasetsList$coords$latitude >= latitudeQuantiles[[1]]))
realValues <- simulatedDatasetsList$responses[!obsIndicesForTraining, 4]
fixedEffsContribs <- as.vector(simulatedDatasetsList$covariates[!obsIndicesForTraining, ] %*% fixedEffectsVec[colnames(simulatedDatasetsList$covariates)]) + interceptValue
# For Vecchia
# dataCovarianceMatrix <- matrix(0, nrow(simulatedDatasetsList$responses), nrow(simulatedDatasetsList$responses))
# allCoordsWithTime <- cbind(simulatedDatasetsList$coords, time = simulatedDatasetsList$timeNumeric)
#
# for (i in 1:nrow(dataCovarianceMatrix)) {
# dataCovarianceMatrix[i:nrow(dataCovarianceMatrix), i] <- customCovFct(allCoordsWithTime[i, ], allCoordsWithTime[i:nrow(dataCovarianceMatrix), ])
# }
# dataCovarianceMatrix[upper.tri(dataCovarianceMatrix)] <- dataCovarianceMatrix[lower.tri(dataCovarianceMatrix)]
fittedModelsByDataset <- lapply(
X = 1,
FUN = simulationFun,
responseMatrix = simulatedDatasetsList$responses,
covariateMatrix = simulatedDatasetsList$covariates,
coordinatesMatrix = as.matrix(simulatedDatasetsList$coords),
timeVecNumeric = simulatedDatasetsList$timeNumeric,
obsIndicesForTraining = obsIndicesForTraining,
funToFitSPDE = fitSPDE,
funToFitVecchia = fitVecchia,
funToFitISMRA = fitISMRA,
numThreads = 4,
controlForSPDE = list(
mesh.2d.offset = c(2,2),
mesh.2d.max.edge = 2,
sigma0 = 1,
kappa0 = sqrt(8 * 1.5)), # We assume that spatial smoothness is 1.5 and baseline range rho0 is 1, i.e. log(rho0) = 0.
controlForISMRA = list(control = list(
Mlon = 3,
Mlat = 3,
Mtime = 0,
numKnotsRes0 = 8,
numIterOptim = 20,
tipKnotsThinningRate = 1,
numValuesForIS = 50,
returnQmat = TRUE
)))
# saveDirectory = "outputFiles/")
We now obtain summaries for prediction accuracy across all simulations for both IS-MRA and SPDE.
simulationFiles <- list.files(path = "outputFilesSmallerM/", pattern = "ISMRAsimulationResults", full.names=T)
simulationResults <- lapply(simulationFiles, readRDS)
realFEvalues <- fixedEffectsVec[colnames(simulatedDatasetsList$covariates)]
realHyperparsLogScale <- c(space.rho = log(spRange), time.rho = log(timeRange), scale = log(spScale))
funToProduceGraphs <- function(dirName, dpi = 400, height = 8, width = 8) {
predictionSummariesByMethod <- analysePredResults(folderForSimResults = dirName, patternForFilename = "Dataset[[:digit:]]+\\.rds$", simulatedDataList = simulatedDatasetsList, obsIndicesForTraining = obsIndicesForTraining, shiftISMRApostPredSDs = 0)
FEandHyperGraphsAndCovProbs <- analyseParaEsts(folderForSimResults = dirName, patternForFilename = "Dataset[[:digit:]]+.rds$", simulatedDataList = simulatedDatasetsList, obsIndicesForTraining = obsIndicesForTraining, realFEs = realFEvalues, realHyperparsLogScale = realHyperparsLogScale)
# ggplot2::ggsave(plot = predictionSummariesByMethod$diffBoxPlots, filename = paste(dirName, "diffBoxPlots.jpeg", sep = "/"), device = "jpeg", width = width, height = height, units = "in", dpi = dpi)
# ggplot2::ggsave(plot = predictionSummariesByMethod$summaryBoxPlots, filename = paste(dirName, "boxPlotsByMethod.jpeg", sep = "/"), device = "jpeg", width = width, height = height, units = "in", dpi = dpi)
# lapply(names(FEandHyperGraphsAndCovProbs$FEgraphs), function(FEname) {
# ggplot2::ggsave(plot = FEandHyperGraphsAndCovProbs$FEgraphs[[FEname]], filename = paste(dirName, "/FEgraph_", FEname,".jpeg", sep = ""), device = "jpeg", width = width, height = height, units = "in", dpi = dpi)
# })
# lapply(names(FEandHyperGraphsAndCovProbs$hyperparGraphs), function(hyperName) {
# hyperNameForFilename <- stringr::str_replace(hyperName, pattern = "\\.", replacement = "") # An extra dot in the filename might cause problems on certain systems
# ggplot2::ggsave(plot = FEandHyperGraphsAndCovProbs$hyperparGraphs[[hyperName]], filename = paste(dirName, "/hyperparGraph_", hyperNameForFilename, ".jpeg", sep = ""), device = "jpeg", width = width, height = height, units = "in", dpi = dpi)
# })
list(FEandHyperSummary = FEandHyperGraphsAndCovProbs, predSummary = predictionSummariesByMethod)
}
predAndFittedParsSummaries <- lapply(c("outputFiles/", "outputFilesSmallerM/"), funToProduceGraphs)
absDiffMats <- list(
K6 = predAndFittedParsSummaries[[1]]$FEandHyperSummary$parsAbsDiffSummary,
K4 = predAndFittedParsSummaries[[2]]$FEandHyperSummary$parsAbsDiffSummary)
fixColnames <- function(matName) {
colnames(absDiffMats[[matName]]) <- stringr::str_replace(colnames(absDiffMats[[matName]]), "ISMRA", paste("ISMRA", matName, sep = "_"))
absDiffMats[[matName]]
}
absDiffMats <- lapply(names(absDiffMats), fixColnames)
tableForPaper <- cbind(absDiffMats[[1]], absDiffMats[[2]][ , grep(pattern = "ISMRA", x = colnames(absDiffMats[[2]]))])
tableForPaper <- tableForPaper[ , c(sort(grep("Mean", colnames(tableForPaper), value = T)), sort(grep("Min", colnames(tableForPaper), value = T)), sort(grep("Max", colnames(tableForPaper), value = T)))]
xtable::xtable(tableForPaper, digits = 4)
summaryPlotFrames <- list(
K6 = subset(predAndFittedParsSummaries[[1]]$predSummary$summaryBoxPlotFrame, methodName == "IS-MRA"),
K4 = subset(predAndFittedParsSummaries[[2]]$predSummary$summaryBoxPlotFrame, methodName == "IS-MRA"),
SPDE = subset(predAndFittedParsSummaries[[2]]$predSummary$summaryBoxPlotFrame, methodName == "INLA-SPDE")) # SPDE is the same fit in outputFiles and outputFilesSmallerM.
namesToIterOver <- c("K6", "K4")
summaryPlotFrames[namesToIterOver] <- lapply(namesToIterOver, function(nameToInput) {
K <- stringr::str_sub(nameToInput, start = -1)
summaryPlotFrames[[nameToInput]]$methodName <- replace(summaryPlotFrames[[nameToInput]]$methodName, which(summaryPlotFrames[[nameToInput]]$methodName == "IS-MRA"), paste("IS-MRA (K=", K, ")", sep = ""))
summaryPlotFrames[[nameToInput]]
})
boxPlotObject <- ggplot2::ggplot(do.call("rbind", summaryPlotFrames), ggplot2::aes(x = predStatName, y = Value, colour = factor(methodName))) + ggplot2::geom_boxplot(outlier.colour = "red", outlier.shape = 1) + ggplot2::theme_bw() + ggplot2::xlab("Prediction statistic") + ggplot2::scale_colour_hue(name = "Method") + ggplot2::theme(legend.position = c(0.825, 0.875), text = ggplot2::element_text(size = 25)) + ggplot2::scale_x_discrete(limits = c("MSPE", "MedSPE"), labels = c("MSPE", "Med-SPE")) + ggplot2::ylab("Prediction error")
# + ggplot2::scale_y_continuous(trans = "log")
ggplot2::ggsave(plot = boxPlotObject, filename = "outputFiles/summaryBoxPlotsK6_K4_SPDE.jpeg", device = "jpeg", width = 8, height = 8, dpi = 400)
FEandHyperPlotFrames <- list(K6 = predAndFittedParsSummaries[[1]]$FEandHyperSummary$parPlotFrames, K4 = predAndFittedParsSummaries[[2]]$FEandHyperSummary$parPlotFrames)
funToUpdateFrames <- function(depthName) {
K <- stringr::str_sub(depthName, start = -1)
mergedFrames <- as.data.frame(do.call("rbind", lapply(FEandHyperPlotFrames[[depthName]], function(x) x[ , c("dataIndex", "paraName", "Method", "CredInt2.5", "CredInt97.5")])))
mergedFrames$Method <- with(mergedFrames, replace(Method, which(Method == "ISMRA"), paste("IS-MRA (K=", K, ")", sep = "")))
hyperparIndices <- which(mergedFrames$paraName %in% c("space.rho", "time.rho", "scale"))
mergedFrames$Method <- replace(mergedFrames$Method, which(mergedFrames$Method == "SPDE"), "INLA-SPDE")
mergedFrames[hyperparIndices, c("CredInt2.5", "CredInt97.5")] <- exp(mergedFrames[hyperparIndices, c("CredInt2.5", "CredInt97.5")])
mergedFrames
}
updatedFrames <- lapply(names(FEandHyperPlotFrames), funToUpdateFrames)
frameForSegmentPlots <- do.call("rbind", updatedFrames)
keepIndices <- which(!duplicated(paste(frameForSegmentPlots$Method, frameForSegmentPlots$paraName, frameForSegmentPlots$dataIndex)))
frameForSegmentPlots <- frameForSegmentPlots[keepIndices, ]
numSimsForGraphs <- 50
equivForHyper <- list(scale = bquote(sigma), space.rho = bquote(rho), time.rho = bquote(phi))
equivForFixed <- c(elevation = "Elevation", May28 = "May 28", May29 = "May 29", EvergreenBroadleaf = "Evergreen broadleaf", MixedForest = "Mixed forest", ClosedShrublands = "Closed shrublands", Savannas = "Savannas", Grasslands = "Grasslands", PermanentWetlands = "Permanent wetlands", Croplands = "Croplands", Urban = "Urban", CroplandNaturalMosaics = "Croplands-natural mosaics", NonVegetated = "Non-vegetated")
realParValues <- c(realFEvalues, exp(realHyperparsLogScale))
funToGetParCredIntPlots <- function(parameterName, plotFrame, equivForFixed, realParValues, equivForHyper, numSimsForGraphs) {
# We won't plot credibility intervals for dataset 18, which produced a very unusually bad estimate, which throws off the graph scaling.
yAxisName <- NULL
if (parameterName %in% names(equivForHyper)) {
yAxisName <- equivForHyper[[parameterName]]
} else {
yAxisName <- bquote(beta ~ ": " ~ .(equivForFixed[[parameterName]]))
}
plotFrame <- subset(plotFrame, subset = (dataIndex != 18) & (dataIndex <= numSimsForGraphs + 1) & (paraName == parameterName)) # Result for sim. dataset 18 is odd
dataIndexShift <- c("IS-MRA (K=4)" = 0, "IS-MRA (K=6)" = max(plotFrame$dataIndex), "INLA-SPDE" = max(plotFrame$dataIndex) * 2)
plotFrame$dataIndex <- rank(plotFrame$dataIndex + dataIndexShift[plotFrame$Method])
segmentPlot <- ggplot2::ggplot(plotFrame, ggplot2::aes(x = dataIndex, group = Method, colour = Method)) + ggplot2::theme_bw() + ggplot2::theme(legend.position = c(0.825, 0.875), text = ggplot2::element_text(size = 23), axis.text.x= ggplot2::element_blank()) + ggplot2::scale_colour_manual(values = c("goldenrod", "blue", "firebrick1")) + ggplot2::geom_vline(xintercept = numSimsForGraphs + 0.5, linetype="dashed", color = "black") + ggplot2::geom_vline(xintercept = numSimsForGraphs*2 + 0.5, linetype="dashed", color = "black") + ggplot2::geom_hline(yintercept = realParValues[[parameterName]], linetype="dashed", color = "red") + ggplot2::xlab("") + ggplot2::ylab(yAxisName) +
#ggplot2::geom_errorbar(ggplot2::aes(ymin = CredInt2.5, ymax = CredInt97.5), width = .2, position = ggplot2::position_dodge(.95))
ggplot2::geom_errorbar(ggplot2::aes(ymin = CredInt2.5, ymax = CredInt97.5), width = .2)
if (parameterName %in% c("space.rho", "time.rho", "scale")) {
segmentPlot <- segmentPlot + ggplot2::coord_trans(y = "log")
if (parameterName == "space.rho") {
breaks <- 1:15
} else if (parameterName == "time.rho") {
breaks <- 1:6
} else {
breaks <- 1:14
}
segmentPlot <- segmentPlot + ggplot2::scale_y_continuous(breaks = breaks)
}
ggplot2::ggsave(filename = paste("outputFiles/CredIntPlot_", parameterName, ".jpeg", sep = ""), plot = segmentPlot, device = "jpeg", width = 8, height = 8, dpi = 400)
NULL
}
lapply(unique(frameForSegmentPlots$paraName), funToGetParCredIntPlots, plotFrame = frameForSegmentPlots, equivForFixed = equivForFixed, realParValues = realParValues, equivForHyper = equivForHyper, numSimsForGraphs = numSimsForGraphs)
formattedTables <- lapply(predAndFittedParsSummaries, function(x) {
methodNames <- c("SPDE", "ISMRA")
printList <- lapply(methodNames, function(methodName) {
stringToPrint <- xtable::xtable(x$predSummary$summaryStats[[methodName]], digits = 4)
stringToPrint
})
names(printList) <- methodNames
printList
})
coverageProbs <- cbind(predAndFittedParsSummaries[[1]]$FEandHyperSummary$coverageProbsMatrix, predAndFittedParsSummaries[[2]]$FEandHyperSummary$coverageProbsMatrix[, "ISMRA"])
colnames(coverageProbs) <- c("SPDE", "IS-MRA_K6", "IS-MRA_K4")
xtable::xtable(coverageProbs)
hyperparNames <- c(spatialRange = "space.rho", temporalRange = "time.rho", scale = "scale")
lapply(seq_along(hyperparNames), function(hyperparIndex) {
filename <- paste("outputFiles/", names(hyperparNames)[[hyperparIndex]], "CredIntsAcrossSims.jpeg", sep = "")
ggplot2::ggsave(filename = filename, predAndFittedParsSummaries[[1]]$FEandHyperSummary$hyperparGraphs[[hyperparNames[[hyperparIndex]]]], device = "jpeg", width = 8, height = 8, dpi = 400)
})
In order to reduce the memory footprint, the INLAMRA function does not return the full conditional precision matrix by default. However, in order to show its structure, we can make it return it.
ISMRAsampleResultWithQmatrix <- MRAinla::INLAMRA(
responseVec = simulatedDatasetsList$responses[obsIndicesForTraining, 1],
spatialCoordMat = simulatedDatasetsList$coords[obsIndicesForTraining, ],
predSpatialCoordMat = simulatedDatasetsList$coords[!obsIndicesForTraining, ],
covariateFrame = as.data.frame(simulatedDatasetsList$covariates[obsIndicesForTraining, ]),
predCovariateFrame = as.data.frame(simulatedDatasetsList$covariates[!obsIndicesForTraining, ]),
timePOSIXorNumericVec = simulatedDatasetsList$timeNumeric[obsIndicesForTraining],
predTimePOSIXorNumericVec = simulatedDatasetsList$timeNumeric[!obsIndicesForTraining],
spatialRangeList = list(log(4), c(log(4), 1)),
spatialSmoothnessList = list(log(1.5)),
timeRangeList = list(log(2), c(log(2), 1)),
timeSmoothnessList = list(log(0.5)),
scaleList = list(log(3), c(log(3), 1)),
errorSDlist = list(log(0.5)),
fixedEffSDlist = list(log(2)),
control = list(
Mlon = 2,
Mlat = 2,
Mtime = 0,
numKnotsRes0 = 20,
numIterOptim = 20,
tipKnotsThinningRate = 1,
numValuesForIS = 50,
returnQmat = TRUE,
numOpenMPthreads = 4))
library(Matrix)
jpeg("outputFiles/QmatrixExample.jpg", width = 1200, height = 1200)
image(ISMRAsampleResultWithQmatrix$Qmat)
dev.off()
In the following outputs, we have the posterior predictive distribution for INLA based on fitted values, which should also reflect the contribution from the uncorrelated error term.
controlForSPDE <- list(mesh.2d.max.n = c(100, 100),
mesh.2d.offset = c(5,5),
useFittedValues = TRUE)
updatedSimResults <- recomputePredictionsForSimOutputs(
folderForSimResults = "outputFiles/",
patternForFilename = "Dataset.+\\.rds$",
coordinatesMatrixTraining = as.matrix(simulatedDatasetsList$coords[obsIndicesForTraining, ]),
coordinatesMatrixTest = as.matrix(simulatedDatasetsList$coords[!obsIndicesForTraining, ]),
timeVecNumericTraining = as.numeric(simulatedDatasetsList$time[obsIndicesForTraining]),
timeVecNumericTest = as.numeric(simulatedDatasetsList$time[!obsIndicesForTraining]),
covariateMatrixTraining = as.matrix(simulatedDatasetsList$covariates[obsIndicesForTraining, ]),
covariateMatrixTest = as.matrix(simulatedDatasetsList$covariates[!obsIndicesForTraining, ]),
responseMatrixTraining = simulatedDatasetsList$responses[obsIndicesForTraining, ],
saveResult = FALSE,
controlForSPDE = controlForSPDE)
We produce a plotted example of the training data.
numRasterRows <- 45
numRasterCols <- 44 # Approximate, hard to recover because data were collected in the sinusoidal projection, which does not translate to a perfectly-spaced grid of coordinates in the longitude/latitude projection.
boundaries <- c(min(simulatedDatasetsList$responses[ , 1]) - 1e-300, max(simulatedDatasetsList$responses[ , 1]) + 1e-300)
colourBreaks <- round(seq(from = boundaries[[1]], to = boundaries[[2]], length.out = 8), digits = 2)
jpeg("outputFiles/sampleSimulatedDataset.jpeg", width = 1000, height = 1000, units = "px")
layout(mat = matrix(c(1, 3, 2, 4), 2, 2))
for (dayToPlot in c("May27", "May28", "May29")) {
observationsToPlot <- intersect(which(simulatedDatasetsList$time == dayToPlot), which(obsIndicesForTraining))
rasterToPlot <- raster::raster(nrows=numRasterRows, ncols=numRasterCols, xmn=min(simulatedDatasetsList$coords$longitude) - 0.0001, xmx=max(simulatedDatasetsList$coords$longitude) + 0.0001, ymn=min(simulatedDatasetsList$coords$latitude) - 0.0001, ymx=max(simulatedDatasetsList$coords$latitude) + 0.0001)
rasterToPlot <- raster::rasterize(x = sp::SpatialPoints(coords = simulatedDatasetsList$coords[observationsToPlot, ]), y = rasterToPlot, field = simulatedDatasetsList$responses[observationsToPlot, 1])
ecol <- mapmisc::colourScale(x = raster::values(rasterToPlot), breaks = 9, col = "Spectral")
ecol$col <- rev(ecol$col)
raster::plot(rasterToPlot, col = ecol$col, breaks = colourBreaks, legend = TRUE)
}
dev.off()
# copySPDEresults(originFolder = "outputFiles/", destinationFolder = "outputFilesSmallerM/", patternForFilename = "Dataset[[:digit:]]+.rds$")
copySPDEresults(originFolder = "outputFiles/", destinationFolder = "outputFilesSmallerM/", patternForFilename = "Dataset[[:digit:]]+\\.rds$")
villandre/ISMRAanalyses documentation built on Dec. 23, 2021, 3:11 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
We start by simulating data.
if ("package:ISMRAanalyses" %in% search()) { detach("package:ISMRAanalyses", unload = TRUE) } library(ISMRAanalyses) library(MRAinla) library(sp) setwd("/home/luc/INLAMRAfiles/INLAMRApaper1/simulations") # setwd("/home/luc/ISMRAfiles/")
The simulated data are based on satellite imagery data obtained in Maharashtra. We first define a rectangular zone in the longitude-latitude projection. Since MODIS data are in the sinusoidal projection, we project the zone we just defined in that coordinate system. Warnings can be safely ignored: they result from a change in the conventions in referring to projections, cf. details in sp::proj4string.
MaharashtraPolygonEdges <- rbind(c(18.8, 73.2), c(18.8, 73.4), c(18.6, 73.4), c(18.6, 73.2), c(18.8, 73.2)) MaharashtraPolygonEdges <- MaharashtraPolygonEdges[ , 2:1] MaharashtraPolygon <- sp::SpatialPolygons(Srl = list(sp::Polygons(list(sp::Polygon(coords = MaharashtraPolygonEdges)), ID = "Mumbai"))) raster::crs(MaharashtraPolygon) <- sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") MaharashtraPolygonOtherCRS <- sp::spTransform(MaharashtraPolygon, CRSobj = sp::CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"))
We then import land cover and elevation values from the MODIS and ASTER data files we obtained from EarthData.
folderForEarthDataFiles <- "/store/luc/rawDataFiles/" landCoverFiles <- list.files(folderForEarthDataFiles, pattern = "MCD*", full.names = TRUE) landCoverRasterSinusoidal <- produceLandCover(landCoverFiles, regionPolygon = MaharashtraPolygonOtherCRS) elevationFiles <- list.files(path = folderForEarthDataFiles, pattern = "*dem.tif", full.names = TRUE) elevationRasterList <- lapply(elevationFiles, raster::raster)
From the imported data, we produce the covariate matrix used in the simulations.
dayOffset <- 121 # Days are listed numerically, from 1 to 366 (leap year). Day 121 corresponds to April 30. dayRange <- 27:29 # We extract data for May 27 until May 29. collectionDatesPOSIX <- as.POSIXct(paste("2012-05-", dayRange, sep = "")) covariateData <- prepareCovariateDataForISMRA(elevationRasterList, landCoverRasterSinusoidal, collectionDatesPOSIX) saveRDS(covariateData, file = "data/covariatesForSimulationStudy.rds", compress = TRUE)
We are now ready to simulate $100$ vectors from a spatiotemporal Gaussian process.
timepoints <- as.numeric(time(covariateData@time))/(3600*24) # Time is in days # Parameter values based on real data analysis (cf. table 1 in paper) # Spatial parameters spSmoothness <- 1.5 spRange <- 5.66 spScale <- 4.14 # Time parameters timeSmoothness <- 0.5 timeRange <- 3.60 timeScale <- 1 # Uncorrelated error errorSD <- 0.5 covFunAlt <- function(spacetimeObject, timepointsVec) { spatialDistances <- fields::RdistEarth(spacetimeObject@sp@coords, miles=FALSE) timeDistances <- fields::rdist(timepointsVec, timepointsVec) MRAinla::maternCov(d = spatialDistances, rho = spRange, smoothness = spSmoothness, scale = spScale) * MRAinla::maternCov(d = timeDistances, scale = timeScale, rho = timeRange, smoothness = timeSmoothness) } covarStruct <- covFunAlt(covariateData, timepoints) + errorSD^2 * diag(nrow(covariateData@data)) choleskyDecomp <- t(chol(covarStruct)) set.seed(10) numReplicates <- 100 fieldValuesWithError <- replicate(n = numReplicates, expr = as.numeric(choleskyDecomp %*% rnorm(nrow(covarStruct))))
We then add fixed effects to the field values, yielding the data we are going to model using the methods we selected for the simulation study. We start by inputting parameter values.
landCoverEffects <- c( Water = 0, EvergreenNeedleleaf = 1.09, EvergreenBroadleaf = 0.979, DeciduousBroadleaf = 1.181, MixedForest = 1.087, ClosedShrublands = 2.097, OpenShrublands = 1.448, WoodySavannas = 1.329, Savannas = 1.393, Grasslands = 1.506, PermanentWetlands = 0.705, Croplands = 1.588, Urban = 1.623, CroplandNaturalMosaics = 1.430, NonVegetated = 0.667) # Water is the reference (landCover = 0) referenceLandCover <- "Water" referenceLandCoverIndex <- match(referenceLandCover, names(landCoverEffects)) timeEffects <- c( May27 = -0.033, May28 = -6.319, May29 = 1.155) + 0.033 # We're moving the reference day from May 25 to May 27, hence + 0.033 elevationEffect <- -0.001 interceptValue <- 40 # This is arbitrary. It should not affect the analyses. referenceTimepoint <- "May27" referenceTimepointIndex <- match(referenceTimepoint, names(timeEffects)) fixedEffectsVec <- c(Intercept = interceptValue, landCoverEffects[-referenceLandCoverIndex], elevation = elevationEffect, timeEffects[-referenceTimepointIndex])
We now create the covariate matrix and add covariate effects to the previously-simulated responses.
landCoverFactor <- factor(covariateData@data$landCover, levels = 0:14, labels = names(landCoverEffects)) timepointsFactor <- factor(timepoints, levels = unique(timepoints), labels = names(timeEffects)) covariateMatrix <- cbind(model.matrix(object = ~ landCover, data = data.frame(landCover = landCoverFactor)), elevation = covariateData@data$elevation, model.matrix(object = ~ time, data = data.frame(time = timepointsFactor))[ , -1]) # We remove the intercept, which has already been created for land cover. Ocean and May 27 are reference categories. colnames(covariateMatrix) <- replace(colnames(covariateMatrix), which(colnames(covariateMatrix) == "(Intercept)"), "Intercept") colnames(covariateMatrix) <- gsub(pattern = "landCover", replacement = "", x = colnames(covariateMatrix)) colnames(covariateMatrix) <- gsub(pattern = "time", replacement = "", x = colnames(covariateMatrix)) # Getting the simulated values. observedValues <- fieldValuesWithError + as.vector(covariateMatrix %*% fixedEffectsVec[colnames(covariateMatrix)]) # Fixing the covariate matrix landCoverTypes <- names(landCoverEffects)[-referenceLandCoverIndex] presentLandCoverTypes <- landCoverTypes[sapply(landCoverTypes, function(coverType) !all(covariateMatrix[ , coverType] == 0))] covariateMatrixNoAbsentLandCovers <- cbind(covariateMatrix[ , -match(landCoverTypes, colnames(covariateMatrix))], covariateMatrix[ , presentLandCoverTypes]) simulatedDatasetsList <- list( responses = observedValues, covariates = covariateMatrixNoAbsentLandCovers[ , -match("Intercept", colnames(covariateMatrixNoAbsentLandCovers))], coords = covariateData@sp@coords, time = timepointsFactor, timeNumeric = timepoints) colnames(simulatedDatasetsList$coords) <- c("longitude", "latitude") saveRDS(simulatedDatasetsList, file = "data/dataForSimulations.rds", compress = TRUE)
As a sanity check, we plot one of the simulated datasets.
numRasterRows <- numRasterCols <- 45 # Approximate, hard to recover because data were collected in the sinusoidal projection, which does not translate to a perfectly-spaced grid of coordinates in the longitude/latitude projection. observationsToPlot <- which(simulatedDatasetsList$time == "May28") jpeg("outputFiles/dataExample.jpeg") fields::quilt.plot( x = simulatedDatasetsList$coords[ , "x"][observationsToPlot], y = simulatedDatasetsList$coords[ , "y"][observationsToPlot], z = simulatedDatasetsList$responses[observationsToPlot, 1], nx = numRasterRows, ny = numRasterCols) dev.off()
We now begin the simulations per se.
simulatedDatasetsList <- readRDS("data/dataForSimulations.rds") simulatedDatasetsList$coords <- as.data.frame(simulatedDatasetsList$coords) longitudeQuantiles <- quantile(simulatedDatasetsList$coords$longitude, probs = c(0.25, 0.75)) latitudeQuantiles <- quantile(simulatedDatasetsList$coords$latitude, probs = c(0.25, 0.75)) obsIndicesForTraining <- !((simulatedDatasetsList$time == "May28") & (simulatedDatasetsList$coords$longitude <= longitudeQuantiles[[2]]) & (simulatedDatasetsList$coords$longitude >= longitudeQuantiles[[1]]) & (simulatedDatasetsList$coords$latitude <= latitudeQuantiles[[2]]) & (simulatedDatasetsList$coords$latitude >= latitudeQuantiles[[1]])) realValues <- simulatedDatasetsList$responses[!obsIndicesForTraining, 4] fixedEffsContribs <- as.vector(simulatedDatasetsList$covariates[!obsIndicesForTraining, ] %*% fixedEffectsVec[colnames(simulatedDatasetsList$covariates)]) + interceptValue # For Vecchia # dataCovarianceMatrix <- matrix(0, nrow(simulatedDatasetsList$responses), nrow(simulatedDatasetsList$responses)) # allCoordsWithTime <- cbind(simulatedDatasetsList$coords, time = simulatedDatasetsList$timeNumeric) # # for (i in 1:nrow(dataCovarianceMatrix)) { # dataCovarianceMatrix[i:nrow(dataCovarianceMatrix), i] <- customCovFct(allCoordsWithTime[i, ], allCoordsWithTime[i:nrow(dataCovarianceMatrix), ]) # } # dataCovarianceMatrix[upper.tri(dataCovarianceMatrix)] <- dataCovarianceMatrix[lower.tri(dataCovarianceMatrix)] fittedModelsByDataset <- lapply( X = 1, FUN = simulationFun, responseMatrix = simulatedDatasetsList$responses, covariateMatrix = simulatedDatasetsList$covariates, coordinatesMatrix = as.matrix(simulatedDatasetsList$coords), timeVecNumeric = simulatedDatasetsList$timeNumeric, obsIndicesForTraining = obsIndicesForTraining, funToFitSPDE = fitSPDE, funToFitVecchia = fitVecchia, funToFitISMRA = fitISMRA, numThreads = 4, controlForSPDE = list( mesh.2d.offset = c(2,2), mesh.2d.max.edge = 2, sigma0 = 1, kappa0 = sqrt(8 * 1.5)), # We assume that spatial smoothness is 1.5 and baseline range rho0 is 1, i.e. log(rho0) = 0. controlForISMRA = list(control = list( Mlon = 3, Mlat = 3, Mtime = 0, numKnotsRes0 = 8, numIterOptim = 20, tipKnotsThinningRate = 1, numValuesForIS = 50, returnQmat = TRUE ))) # saveDirectory = "outputFiles/")
We now obtain summaries for prediction accuracy across all simulations for both IS-MRA and SPDE.
simulationFiles <- list.files(path = "outputFilesSmallerM/", pattern = "ISMRAsimulationResults", full.names=T) simulationResults <- lapply(simulationFiles, readRDS) realFEvalues <- fixedEffectsVec[colnames(simulatedDatasetsList$covariates)] realHyperparsLogScale <- c(space.rho = log(spRange), time.rho = log(timeRange), scale = log(spScale)) funToProduceGraphs <- function(dirName, dpi = 400, height = 8, width = 8) { predictionSummariesByMethod <- analysePredResults(folderForSimResults = dirName, patternForFilename = "Dataset[[:digit:]]+\\.rds$", simulatedDataList = simulatedDatasetsList, obsIndicesForTraining = obsIndicesForTraining, shiftISMRApostPredSDs = 0) FEandHyperGraphsAndCovProbs <- analyseParaEsts(folderForSimResults = dirName, patternForFilename = "Dataset[[:digit:]]+.rds$", simulatedDataList = simulatedDatasetsList, obsIndicesForTraining = obsIndicesForTraining, realFEs = realFEvalues, realHyperparsLogScale = realHyperparsLogScale) # ggplot2::ggsave(plot = predictionSummariesByMethod$diffBoxPlots, filename = paste(dirName, "diffBoxPlots.jpeg", sep = "/"), device = "jpeg", width = width, height = height, units = "in", dpi = dpi) # ggplot2::ggsave(plot = predictionSummariesByMethod$summaryBoxPlots, filename = paste(dirName, "boxPlotsByMethod.jpeg", sep = "/"), device = "jpeg", width = width, height = height, units = "in", dpi = dpi) # lapply(names(FEandHyperGraphsAndCovProbs$FEgraphs), function(FEname) { # ggplot2::ggsave(plot = FEandHyperGraphsAndCovProbs$FEgraphs[[FEname]], filename = paste(dirName, "/FEgraph_", FEname,".jpeg", sep = ""), device = "jpeg", width = width, height = height, units = "in", dpi = dpi) # }) # lapply(names(FEandHyperGraphsAndCovProbs$hyperparGraphs), function(hyperName) { # hyperNameForFilename <- stringr::str_replace(hyperName, pattern = "\\.", replacement = "") # An extra dot in the filename might cause problems on certain systems # ggplot2::ggsave(plot = FEandHyperGraphsAndCovProbs$hyperparGraphs[[hyperName]], filename = paste(dirName, "/hyperparGraph_", hyperNameForFilename, ".jpeg", sep = ""), device = "jpeg", width = width, height = height, units = "in", dpi = dpi) # }) list(FEandHyperSummary = FEandHyperGraphsAndCovProbs, predSummary = predictionSummariesByMethod) } predAndFittedParsSummaries <- lapply(c("outputFiles/", "outputFilesSmallerM/"), funToProduceGraphs) absDiffMats <- list( K6 = predAndFittedParsSummaries[[1]]$FEandHyperSummary$parsAbsDiffSummary, K4 = predAndFittedParsSummaries[[2]]$FEandHyperSummary$parsAbsDiffSummary) fixColnames <- function(matName) { colnames(absDiffMats[[matName]]) <- stringr::str_replace(colnames(absDiffMats[[matName]]), "ISMRA", paste("ISMRA", matName, sep = "_")) absDiffMats[[matName]] } absDiffMats <- lapply(names(absDiffMats), fixColnames) tableForPaper <- cbind(absDiffMats[[1]], absDiffMats[[2]][ , grep(pattern = "ISMRA", x = colnames(absDiffMats[[2]]))]) tableForPaper <- tableForPaper[ , c(sort(grep("Mean", colnames(tableForPaper), value = T)), sort(grep("Min", colnames(tableForPaper), value = T)), sort(grep("Max", colnames(tableForPaper), value = T)))] xtable::xtable(tableForPaper, digits = 4) summaryPlotFrames <- list( K6 = subset(predAndFittedParsSummaries[[1]]$predSummary$summaryBoxPlotFrame, methodName == "IS-MRA"), K4 = subset(predAndFittedParsSummaries[[2]]$predSummary$summaryBoxPlotFrame, methodName == "IS-MRA"), SPDE = subset(predAndFittedParsSummaries[[2]]$predSummary$summaryBoxPlotFrame, methodName == "INLA-SPDE")) # SPDE is the same fit in outputFiles and outputFilesSmallerM. namesToIterOver <- c("K6", "K4") summaryPlotFrames[namesToIterOver] <- lapply(namesToIterOver, function(nameToInput) { K <- stringr::str_sub(nameToInput, start = -1) summaryPlotFrames[[nameToInput]]$methodName <- replace(summaryPlotFrames[[nameToInput]]$methodName, which(summaryPlotFrames[[nameToInput]]$methodName == "IS-MRA"), paste("IS-MRA (K=", K, ")", sep = "")) summaryPlotFrames[[nameToInput]] }) boxPlotObject <- ggplot2::ggplot(do.call("rbind", summaryPlotFrames), ggplot2::aes(x = predStatName, y = Value, colour = factor(methodName))) + ggplot2::geom_boxplot(outlier.colour = "red", outlier.shape = 1) + ggplot2::theme_bw() + ggplot2::xlab("Prediction statistic") + ggplot2::scale_colour_hue(name = "Method") + ggplot2::theme(legend.position = c(0.825, 0.875), text = ggplot2::element_text(size = 25)) + ggplot2::scale_x_discrete(limits = c("MSPE", "MedSPE"), labels = c("MSPE", "Med-SPE")) + ggplot2::ylab("Prediction error") # + ggplot2::scale_y_continuous(trans = "log") ggplot2::ggsave(plot = boxPlotObject, filename = "outputFiles/summaryBoxPlotsK6_K4_SPDE.jpeg", device = "jpeg", width = 8, height = 8, dpi = 400) FEandHyperPlotFrames <- list(K6 = predAndFittedParsSummaries[[1]]$FEandHyperSummary$parPlotFrames, K4 = predAndFittedParsSummaries[[2]]$FEandHyperSummary$parPlotFrames) funToUpdateFrames <- function(depthName) { K <- stringr::str_sub(depthName, start = -1) mergedFrames <- as.data.frame(do.call("rbind", lapply(FEandHyperPlotFrames[[depthName]], function(x) x[ , c("dataIndex", "paraName", "Method", "CredInt2.5", "CredInt97.5")]))) mergedFrames$Method <- with(mergedFrames, replace(Method, which(Method == "ISMRA"), paste("IS-MRA (K=", K, ")", sep = ""))) hyperparIndices <- which(mergedFrames$paraName %in% c("space.rho", "time.rho", "scale")) mergedFrames$Method <- replace(mergedFrames$Method, which(mergedFrames$Method == "SPDE"), "INLA-SPDE") mergedFrames[hyperparIndices, c("CredInt2.5", "CredInt97.5")] <- exp(mergedFrames[hyperparIndices, c("CredInt2.5", "CredInt97.5")]) mergedFrames } updatedFrames <- lapply(names(FEandHyperPlotFrames), funToUpdateFrames) frameForSegmentPlots <- do.call("rbind", updatedFrames) keepIndices <- which(!duplicated(paste(frameForSegmentPlots$Method, frameForSegmentPlots$paraName, frameForSegmentPlots$dataIndex))) frameForSegmentPlots <- frameForSegmentPlots[keepIndices, ] numSimsForGraphs <- 50 equivForHyper <- list(scale = bquote(sigma), space.rho = bquote(rho), time.rho = bquote(phi)) equivForFixed <- c(elevation = "Elevation", May28 = "May 28", May29 = "May 29", EvergreenBroadleaf = "Evergreen broadleaf", MixedForest = "Mixed forest", ClosedShrublands = "Closed shrublands", Savannas = "Savannas", Grasslands = "Grasslands", PermanentWetlands = "Permanent wetlands", Croplands = "Croplands", Urban = "Urban", CroplandNaturalMosaics = "Croplands-natural mosaics", NonVegetated = "Non-vegetated") realParValues <- c(realFEvalues, exp(realHyperparsLogScale)) funToGetParCredIntPlots <- function(parameterName, plotFrame, equivForFixed, realParValues, equivForHyper, numSimsForGraphs) { # We won't plot credibility intervals for dataset 18, which produced a very unusually bad estimate, which throws off the graph scaling. yAxisName <- NULL if (parameterName %in% names(equivForHyper)) { yAxisName <- equivForHyper[[parameterName]] } else { yAxisName <- bquote(beta ~ ": " ~ .(equivForFixed[[parameterName]])) } plotFrame <- subset(plotFrame, subset = (dataIndex != 18) & (dataIndex <= numSimsForGraphs + 1) & (paraName == parameterName)) # Result for sim. dataset 18 is odd dataIndexShift <- c("IS-MRA (K=4)" = 0, "IS-MRA (K=6)" = max(plotFrame$dataIndex), "INLA-SPDE" = max(plotFrame$dataIndex) * 2) plotFrame$dataIndex <- rank(plotFrame$dataIndex + dataIndexShift[plotFrame$Method]) segmentPlot <- ggplot2::ggplot(plotFrame, ggplot2::aes(x = dataIndex, group = Method, colour = Method)) + ggplot2::theme_bw() + ggplot2::theme(legend.position = c(0.825, 0.875), text = ggplot2::element_text(size = 23), axis.text.x= ggplot2::element_blank()) + ggplot2::scale_colour_manual(values = c("goldenrod", "blue", "firebrick1")) + ggplot2::geom_vline(xintercept = numSimsForGraphs + 0.5, linetype="dashed", color = "black") + ggplot2::geom_vline(xintercept = numSimsForGraphs*2 + 0.5, linetype="dashed", color = "black") + ggplot2::geom_hline(yintercept = realParValues[[parameterName]], linetype="dashed", color = "red") + ggplot2::xlab("") + ggplot2::ylab(yAxisName) + #ggplot2::geom_errorbar(ggplot2::aes(ymin = CredInt2.5, ymax = CredInt97.5), width = .2, position = ggplot2::position_dodge(.95)) ggplot2::geom_errorbar(ggplot2::aes(ymin = CredInt2.5, ymax = CredInt97.5), width = .2) if (parameterName %in% c("space.rho", "time.rho", "scale")) { segmentPlot <- segmentPlot + ggplot2::coord_trans(y = "log") if (parameterName == "space.rho") { breaks <- 1:15 } else if (parameterName == "time.rho") { breaks <- 1:6 } else { breaks <- 1:14 } segmentPlot <- segmentPlot + ggplot2::scale_y_continuous(breaks = breaks) } ggplot2::ggsave(filename = paste("outputFiles/CredIntPlot_", parameterName, ".jpeg", sep = ""), plot = segmentPlot, device = "jpeg", width = 8, height = 8, dpi = 400) NULL } lapply(unique(frameForSegmentPlots$paraName), funToGetParCredIntPlots, plotFrame = frameForSegmentPlots, equivForFixed = equivForFixed, realParValues = realParValues, equivForHyper = equivForHyper, numSimsForGraphs = numSimsForGraphs) formattedTables <- lapply(predAndFittedParsSummaries, function(x) { methodNames <- c("SPDE", "ISMRA") printList <- lapply(methodNames, function(methodName) { stringToPrint <- xtable::xtable(x$predSummary$summaryStats[[methodName]], digits = 4) stringToPrint }) names(printList) <- methodNames printList }) coverageProbs <- cbind(predAndFittedParsSummaries[[1]]$FEandHyperSummary$coverageProbsMatrix, predAndFittedParsSummaries[[2]]$FEandHyperSummary$coverageProbsMatrix[, "ISMRA"]) colnames(coverageProbs) <- c("SPDE", "IS-MRA_K6", "IS-MRA_K4") xtable::xtable(coverageProbs) hyperparNames <- c(spatialRange = "space.rho", temporalRange = "time.rho", scale = "scale") lapply(seq_along(hyperparNames), function(hyperparIndex) { filename <- paste("outputFiles/", names(hyperparNames)[[hyperparIndex]], "CredIntsAcrossSims.jpeg", sep = "") ggplot2::ggsave(filename = filename, predAndFittedParsSummaries[[1]]$FEandHyperSummary$hyperparGraphs[[hyperparNames[[hyperparIndex]]]], device = "jpeg", width = 8, height = 8, dpi = 400) })
In order to reduce the memory footprint, the INLAMRA function does not return the full conditional precision matrix by default. However, in order to show its structure, we can make it return it.
ISMRAsampleResultWithQmatrix <- MRAinla::INLAMRA( responseVec = simulatedDatasetsList$responses[obsIndicesForTraining, 1], spatialCoordMat = simulatedDatasetsList$coords[obsIndicesForTraining, ], predSpatialCoordMat = simulatedDatasetsList$coords[!obsIndicesForTraining, ], covariateFrame = as.data.frame(simulatedDatasetsList$covariates[obsIndicesForTraining, ]), predCovariateFrame = as.data.frame(simulatedDatasetsList$covariates[!obsIndicesForTraining, ]), timePOSIXorNumericVec = simulatedDatasetsList$timeNumeric[obsIndicesForTraining], predTimePOSIXorNumericVec = simulatedDatasetsList$timeNumeric[!obsIndicesForTraining], spatialRangeList = list(log(4), c(log(4), 1)), spatialSmoothnessList = list(log(1.5)), timeRangeList = list(log(2), c(log(2), 1)), timeSmoothnessList = list(log(0.5)), scaleList = list(log(3), c(log(3), 1)), errorSDlist = list(log(0.5)), fixedEffSDlist = list(log(2)), control = list( Mlon = 2, Mlat = 2, Mtime = 0, numKnotsRes0 = 20, numIterOptim = 20, tipKnotsThinningRate = 1, numValuesForIS = 50, returnQmat = TRUE, numOpenMPthreads = 4)) library(Matrix) jpeg("outputFiles/QmatrixExample.jpg", width = 1200, height = 1200) image(ISMRAsampleResultWithQmatrix$Qmat) dev.off()
In the following outputs, we have the posterior predictive distribution for INLA based on fitted values, which should also reflect the contribution from the uncorrelated error term.
controlForSPDE <- list(mesh.2d.max.n = c(100, 100), mesh.2d.offset = c(5,5), useFittedValues = TRUE) updatedSimResults <- recomputePredictionsForSimOutputs( folderForSimResults = "outputFiles/", patternForFilename = "Dataset.+\\.rds$", coordinatesMatrixTraining = as.matrix(simulatedDatasetsList$coords[obsIndicesForTraining, ]), coordinatesMatrixTest = as.matrix(simulatedDatasetsList$coords[!obsIndicesForTraining, ]), timeVecNumericTraining = as.numeric(simulatedDatasetsList$time[obsIndicesForTraining]), timeVecNumericTest = as.numeric(simulatedDatasetsList$time[!obsIndicesForTraining]), covariateMatrixTraining = as.matrix(simulatedDatasetsList$covariates[obsIndicesForTraining, ]), covariateMatrixTest = as.matrix(simulatedDatasetsList$covariates[!obsIndicesForTraining, ]), responseMatrixTraining = simulatedDatasetsList$responses[obsIndicesForTraining, ], saveResult = FALSE, controlForSPDE = controlForSPDE)
We produce a plotted example of the training data.
numRasterRows <- 45 numRasterCols <- 44 # Approximate, hard to recover because data were collected in the sinusoidal projection, which does not translate to a perfectly-spaced grid of coordinates in the longitude/latitude projection. boundaries <- c(min(simulatedDatasetsList$responses[ , 1]) - 1e-300, max(simulatedDatasetsList$responses[ , 1]) + 1e-300) colourBreaks <- round(seq(from = boundaries[[1]], to = boundaries[[2]], length.out = 8), digits = 2) jpeg("outputFiles/sampleSimulatedDataset.jpeg", width = 1000, height = 1000, units = "px") layout(mat = matrix(c(1, 3, 2, 4), 2, 2)) for (dayToPlot in c("May27", "May28", "May29")) { observationsToPlot <- intersect(which(simulatedDatasetsList$time == dayToPlot), which(obsIndicesForTraining)) rasterToPlot <- raster::raster(nrows=numRasterRows, ncols=numRasterCols, xmn=min(simulatedDatasetsList$coords$longitude) - 0.0001, xmx=max(simulatedDatasetsList$coords$longitude) + 0.0001, ymn=min(simulatedDatasetsList$coords$latitude) - 0.0001, ymx=max(simulatedDatasetsList$coords$latitude) + 0.0001) rasterToPlot <- raster::rasterize(x = sp::SpatialPoints(coords = simulatedDatasetsList$coords[observationsToPlot, ]), y = rasterToPlot, field = simulatedDatasetsList$responses[observationsToPlot, 1]) ecol <- mapmisc::colourScale(x = raster::values(rasterToPlot), breaks = 9, col = "Spectral") ecol$col <- rev(ecol$col) raster::plot(rasterToPlot, col = ecol$col, breaks = colourBreaks, legend = TRUE) } dev.off()
# copySPDEresults(originFolder = "outputFiles/", destinationFolder = "outputFilesSmallerM/", patternForFilename = "Dataset[[:digit:]]+.rds$") copySPDEresults(originFolder = "outputFiles/", destinationFolder = "outputFilesSmallerM/", patternForFilename = "Dataset[[:digit:]]+\\.rds$")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.