diff --git a/.Rbuildignore b/.Rbuildignore index a82ff8a..92b7e9d 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -15,6 +15,7 @@ bigstudies .synctex.gz$ ^README.md$ ^devel/devel.R$ +^devel ^R/rhipe_fns.R$ ^cran-comments.md$ ^\.\.pdf$ diff --git a/NAMESPACE b/NAMESPACE index 5dc54a0..ac064f4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -66,6 +66,7 @@ importFrom(stats,model.frame) importFrom(stats,model.matrix) importFrom(stats,na.fail) importFrom(stats,optim) +importFrom(stats,rnorm) importFrom(stats,runif) importFrom(stats,sd) importFrom(stats,terms) diff --git a/R/AllGeneric.R b/R/AllGeneric.R index eb0447d..c7e609d 100644 --- a/R/AllGeneric.R +++ b/R/AllGeneric.R @@ -183,9 +183,11 @@ setGeneric("data.frame<-", function(x, value) standardGeneric("data.frame<-")) #' @param d data, that is, an object of class SpatialPointsDataFrame or SpatialPolygonsDataFrame. Provision of data implies that the domain is bounded (necessary with \code{real_line} and \code{plane} but not necessary with \code{sphere}) #' @param nonconvex_hull flag indicating whether \code{INLA} should be used to create a non-convex domain boundary #' @param convex convex parameter for the \code{INLA} function \code{inla.nonconvex.hull} used for smoothing an extended boundary when working on a finite domain (that is, when the object \code{d} is supplied), see details +#' @param xlims limits of the horizontal axis (overrides automatic selection). +#' @param ylims limits of the vertical axis (overrides automatic selection). #' @param ... currently unused #' @details This generic function is not called directly. Please refer to \code{auto_BAUs} for more details. -setGeneric("auto_BAU", function(manifold,type,cellsize,resl,d,nonconvex_hull,convex,...) standardGeneric("auto_BAU")) +setGeneric("auto_BAU", function(manifold,type,cellsize,resl,d,nonconvex_hull,convex,xlims,ylims,...) standardGeneric("auto_BAU")) #' @title Bin data into BAUs #' @description This is an internal function which bins data into BAUs or aggregates across BAUs if the data have a large footprint. If \code{est_error == TRUE}, the observation error is estimated using the variogram (see vignette for details). diff --git a/R/FRK.R b/R/FRK.R index b41468d..9715f15 100644 --- a/R/FRK.R +++ b/R/FRK.R @@ -1,21 +1,21 @@ -#' Fixed Rank Kriging -#' -#' Fixed Rank Kriging is a tool for spatial/spatio-temporal modelling and prediction with large datasets. The approach, discussed in Cressie and Johannesson (2008), decomposes the field, and hence the covariance function, using a fixed set of n basis functions, where n is typically much smaller than the number of data points (or polygons) m. The method naturally allows for non-stationary, anisotropic covariance functions and the use of observations with varying support (with known error variance). The projected field is a key building block of the Spatial Random Effects (SRE) model, on which this package is based. The package FRK provides helper functions to model, fit, and predict using an SRE with relative ease. Reference: Cressie, N. and Johannesson, G. (2008) . -#' @name FRK-package -#' @docType package -#' @useDynLib FRK, .registration=TRUE -#' @import methods -#' @import ggplot2 -#' @import Matrix -#' @import sp -#' @import spacetime -#' @import parallel -#' @import dplyr -#' @importFrom Hmisc round.POSIXt trunc.POSIXt ceil -#' @importFrom plyr ddply dlply rbind.fill -#' @importFrom digest digest -#' @importFrom Rcpp cppFunction -#' @importFrom grDevices chull -#' @importFrom stats .getXlevels coefficients dist kmeans lm median model.extract model.frame model.matrix na.fail optim runif sd terms var time -#' @importFrom utils data -NULL +#' Fixed Rank Kriging +#' +#' Fixed Rank Kriging is a tool for spatial/spatio-temporal modelling and prediction with large datasets. The approach, discussed in Cressie and Johannesson (2008), decomposes the field, and hence the covariance function, using a fixed set of n basis functions, where n is typically much smaller than the number of data points (or polygons) m. The method naturally allows for non-stationary, anisotropic covariance functions and the use of observations with varying support (with known error variance). The projected field is a key building block of the Spatial Random Effects (SRE) model, on which this package is based. The package FRK provides helper functions to model, fit, and predict using an SRE with relative ease. Reference: Cressie, N. and Johannesson, G. (2008) . +#' @name FRK-package +#' @docType package +#' @useDynLib FRK, .registration=TRUE +#' @import methods +#' @import ggplot2 +#' @import Matrix +#' @import sp +#' @import spacetime +#' @import parallel +#' @import dplyr +#' @importFrom Hmisc round.POSIXt trunc.POSIXt ceil +#' @importFrom plyr ddply dlply rbind.fill +#' @importFrom digest digest +#' @importFrom Rcpp cppFunction +#' @importFrom grDevices chull +#' @importFrom stats .getXlevels coefficients dist kmeans lm median model.extract model.frame model.matrix na.fail optim runif sd terms var time rnorm +#' @importFrom utils data +NULL diff --git a/R/FRK_wrapper.R b/R/FRK_wrapper.R index fb6ebda..be3c1a1 100644 --- a/R/FRK_wrapper.R +++ b/R/FRK_wrapper.R @@ -150,3 +150,39 @@ FRK <- function(f, # formula (compulsory) ## Return fitted SRE model S } + + +## The function below checks the arguments for the function FRK. The code is self-explanatory +## This is similar, but slightly different to, .check_args1() +.check_args_wrapper <- function(f,data,basis,BAUs,est_error) { + if(!is(f,"formula")) stop("f needs to be a formula.") + if(!is(data,"list")) + stop("Please supply a list of Spatial objects.") + if(!all(sapply(data,function(x) is(x,"Spatial") | is(x,"ST")))) + stop("All data list elements need to be of class Spatial or ST.") + if(!all(sapply(data,function(x) all.vars(f)[1] %in% names(x@data)))) + stop("All data list elements to have values for the dependent variable.") + if(!est_error & !all(sapply(data,function(x) "std" %in% names(x@data)))) + stop("If observational error is not going to be estimated, + please supply a field 'std' in the data objects.") + if(!(is.null(BAUs))) { + if(!(is(BAUs,"SpatialPolygonsDataFrame") | is(BAUs,"SpatialPixelsDataFrame") | is(BAUs,"STFDF"))) + stop("BAUs should be a SpatialPolygonsDataFrame or a STFDF object") + if(!all(sapply(data,function(x) identical(proj4string(x), proj4string(BAUs))))) + stop("Please ensure all data items and BAUs have the same coordinate reference system") + if(!(all(BAUs$fs >= 0))) + stop("fine-scale variation basis function needs to be nonnegative everywhere") + if(is(BAUs,"STFDF")) if(!(is(BAUs@sp,"SpatialPolygonsDataFrame") | is(BAUs@sp,"SpatialPixelsDataFrame"))) + stop("The spatial component of the BAUs should be a SpatialPolygonsDataFrame") + if(any(sapply(data,function(x) any(names(x@data) %in% names(BAUs@data))))) + stop("Please don't have overlapping variable names in data and BAUs. All covariates need to be in the BAUs.") + if(!all(all.vars(f)[-1] %in% c(names(BAUs@data),coordnames(BAUs)))) + stop("All covariates need to be in the SpatialPolygons BAU object.") + } + + if(!(is.null(basis))) { + if(!(is(basis,"Basis") | is(basis,"TensorP_Basis"))) + stop("basis needs to be of class Basis or TensorP_Basis (package FRK)") + } +} + diff --git a/R/SREutils.R b/R/SREutils.R index 38fde52..a5a7cd7 100644 --- a/R/SREutils.R +++ b/R/SREutils.R @@ -120,14 +120,12 @@ SRE <- function(f,data,basis,BAUs,est_error=TRUE,average_in_BAU = TRUE, ## Now estimate the measurement error using variogram methods this_data <- .est_obs_error(this_data,variogram.formula=f, - vgm_model = vgm_model) + vgm_model = vgm_model) ## Allocate the measurement standard deviation from variogram analysis data[[i]]$std <- this_data$std } - - ## The next step is to allocate all data (both point and polygon referenced) to BAUs. We can either average data points falling in ## the same BAU (average_in_BAU == TRUE) or not (average_in_BAU == FALSE) print("Binning data ...") data_proc <- map_data_to_BAUs(data[[i]], # data object @@ -139,6 +137,11 @@ SRE <- function(f,data,basis,BAUs,est_error=TRUE,average_in_BAU = TRUE, stop("NAs found when mapping data to BAUs. Do you have NAs in your data? If not, are you sure all your data are covered by BAUs?") + + + + ## The next step is to allocate all data (both point and polygon referenced) to BAUs. We can either average data points falling in + ## Extract information from the data using the .extract.from.formula internal function L <- .extract.from.formula(f,data=data_proc) X[[i]] <- as(L$X,"Matrix") # covariate information @@ -348,6 +351,7 @@ setMethod("summary",signature(object="SRE"), ## If zero fine-scale variation detected just make sure user knows. ## This can be symptomatic of poor fitting if(SRE_model@sigma2fshat == 0) + if(opts_FRK$get("verbose") > 0) message("sigma2fs is being estimated to zero. This might because of an incorrect binning procedure or because too much measurement error @@ -397,6 +401,7 @@ setMethod("summary",signature(object="SRE"), D <- Sm@sigma2fshat*Sm@Vfs + Sm@Ve # total variance of data if(isDiagonal(D)) { # if this is diagonal + D <- Diagonal(x = D@x) # cast to diagonal matrix cholD <- sqrt(D) # just compute sqrt cholDinvT <- solve(cholD) # and the inverse is just the reciprocal } else { @@ -406,7 +411,8 @@ setMethod("summary",signature(object="SRE"), ## Compute log-determinant. This is given by a formula in Section 2.2 S_Dinv_S <- crossprod(cholDinvT %*% S) - log_det_SigmaZ <- determinant(Kinv + S_Dinv_S,logarithm = TRUE)$modulus + + R <- chol(Kinv + S_Dinv_S) + log_det_SigmaZ <- logdet(R) + determinant(K,logarithm = TRUE)$modulus + logdet(cholD) # this computes the log-determinant of a matrix from its Cholesky factor @@ -415,8 +421,6 @@ setMethod("summary",signature(object="SRE"), # SigmaZ_inv <- Dinv - Dinv %*% S %*% solve(Kinv + S_Dinv_S) %*% t(S) %*% Dinv # SigmaZ_inv2 <- Dinv - tcrossprod(Dinv %*% S %*% solve(R)) - R <- chol(Kinv + S_Dinv_S) - ## Compute efficiently rDinv <- t(resid) %*% Dinv rDinv <- crossprod(cholDinvT %*% resid,cholDinvT) @@ -452,6 +456,7 @@ setMethod("summary",signature(object="SRE"), D <- sigma2fs*Sm@Vfs + Sm@Ve # total variance-covariance of Z if(isDiagonal(D)) { # if this is diagonal + D <- Diagonal(x=D@x) # cast to diagonal cholD <- sqrt(D) # then the Cholesky is the sqrt cholDinv <- solve(cholD) # the inverse Cholesky is the inverse-sqrt Dinv <- solve(D) # the inverse is just reciprocal of diagonal elements @@ -613,6 +618,7 @@ setMethod("summary",signature(object="SRE"), } D <- sigma2fs_new*Sm@Vfs + Sm@Ve # total data variance-covariance if(isDiagonal(D)) { # inverse of D (as above) + D <- Diagonal(x=D@x) # cast to Diagonal Dinv <- solve(D) } else { Dinv <- chol2inv(chol(D)) @@ -1014,6 +1020,7 @@ setMethod("summary",signature(object="SRE"), if(Sm@fs_model == "ind") { D <- sigma2fs*Sm@Vfs + Sm@Ve # compute total variance (data) if(isDiagonal(D)) { + D <- Diagonal(x=D@x) # cast to diagonal Dchol <- sqrt(D) # find the inverse of the variance-covariance matrix Dinv <- solve(D) # if Diagonal the Cholesky and inverse are just sqrt and reciprocal } else { @@ -1236,6 +1243,10 @@ setMethod("summary",signature(object="SRE"), ## to study variogram points close to the origin) cutoff <- sqrt(area * 100 / length(sp_pts_sub)) + ## Remove covariates since we are only interested at behaviour close to iring + ## and since the binning into BAUs hasn't (shouldn't have) occurred yet + variogram.formula <- .formula_no_covars(variogram.formula) + ## Extract data values from data object L <- .extract.from.formula(variogram.formula,data=sp_pts_sub) @@ -1256,51 +1267,64 @@ setMethod("summary",signature(object="SRE"), nugget = var(L$y)/2) ## Try to fit the model. - vgm.fit <- suppressWarnings(gstat::fit.variogram(v, model = vgm_model)) - - ## Check if the process of fitting generates a warning. If it did then OK == 0 - ## otherwise OK == 1 (this fits twice but since it's so quick it's not an issue) - OK <- tryCatch({vgm.fit <- gstat::fit.variogram(v, model = vgm_model); 1}, - warning=function(w) 0) - ## If the reporte psill is less or equal to zero, or fit.variogram reported a singularity, - ## or a Warning was thrown, then retry fitting using a linear model on just the first + ## Try fitting using a linear model on just the first two to ## four points (cf. Kang and Cressie) - if(vgm.fit$psill[1] <= 0 | attributes(vgm.fit)$singular | !OK) { - linfit <- lm(gamma~dist,data=v[1:4,]) # fit only first four points - vgm.fit$psill[1] <- coefficients(linfit)[1] # extract psill from intercept + OK <- FALSE # init. OK + num_points <- 2 # start with 2 points + while(!OK & num_points <= 4) { # while NOT OK + ## fit only to first 1:num_points points + linfit <- lm(gamma~dist,data=v[1:num_points,]) + + ## ans. will be returned in vgm.fit$psill later on + vgm.fit <- list(singular = 0) + vgm.fit$psill <- coefficients(linfit)[1] # extract psill from intercept + if(vgm.fit$psill > 0) OK <- TRUE # if variance > 0 OK = TRUE + num_points <- num_points + 1 # increment num_points } - ## If this still didn't work then try to fit an exponential model + + ## If we have estimated a negative variance, then try to fit a linear variogram + ## using more points if(vgm.fit$psill[1] <= 0) { - vgm_model <- gstat::vgm(psill = var(L$y)/2, - model = "Exp", - range = mean(v$dist), - nugget = var(L$y)/2) - OK <- tryCatch({vgm.fit = gstat::fit.variogram(v, model = vgm_model); OK <- 1},warning=function(w) 0) + vgm.fit <- suppressWarnings(gstat::fit.variogram(v, model = vgm_model)) - } - if(vgm.fit$psill[1] <= 0 | attributes(vgm.fit)$singular | !OK) { - ## Try with Gaussian, maybe process is very smooth or data has a large support - vgm_model <- gstat::vgm(var(L$y)/2, "Gau", mean(v$dist), var(L$y)/2) + ## Check if the process of fitting generates a warning. If it did then OK == 0 + ## otherwise OK == 1 (this fits twice but since it's so quick it's not an issue) + OK <- tryCatch({vgm.fit <- gstat::fit.variogram(v, model = vgm_model); 1}, + warning=function(w) 0) + ## If the reportef psill is less or equal to zero, or fit.variogram reported a singularity, + ## or a Warning was thrown, then retry fitting using an exponential model + if(vgm.fit$psill[1] <= 0 | attributes(vgm.fit)$singular | !OK) { + vgm_model <- gstat::vgm(psill = var(L$y)/2, + model = "Exp", + range = mean(v$dist), + nugget = var(L$y)/2) + OK <- tryCatch({vgm.fit = gstat::fit.variogram(v, model = vgm_model); OK <- 1},warning=function(w) 0) + vgm.fit <- suppressWarnings(gstat::fit.variogram(v, model = vgm_model)) + } - ## Try to fit the model. - vgm.fit <- suppressWarnings(gstat::fit.variogram(v, model = vgm_model)) + if(vgm.fit$psill[1] <= 0 | attributes(vgm.fit)$singular | !OK) { + ## Try with Gaussian, maybe process is very smooth or data has a large support + vgm_model <- gstat::vgm(var(L$y)/2, "Gau", mean(v$dist), var(L$y)/2) - ## Like above, we return OK = 0 if fit is still not good - OK <- tryCatch({vgm.fit = gstat::fit.variogram(v, model = vgm_model); OK <- 1},warning=function(w) 0) - } + ## Try to fit the model. + vgm.fit <- suppressWarnings(gstat::fit.variogram(v, model = vgm_model)) - ## If we still have problems, then just take the first point of the the empirical semivariogram and - ## throw a warning that this estimate is probably not very good - if(vgm.fit$psill[1] <= 0 | attributes(vgm.fit)$singular | !OK) { - vgm.fit$psill[1] <- v$gamma[1] - warning("Estimate of measurement error is probably inaccurate. + ## Like above, we return OK = 0 if fit is still not good + OK <- tryCatch({vgm.fit = gstat::fit.variogram(v, model = vgm_model); OK <- 1},warning=function(w) 0) + } + + ## If we still have problems, then just take the first point of the the empirical semivariogram and + ## throw a warning that this estimate is probably not very good + if(vgm.fit$psill[1] <= 0 | attributes(vgm.fit)$singular | !OK) { + vgm.fit$psill[1] <- v$gamma[1] + warning("Estimate of measurement error is probably inaccurate. Please consider setting it through the std variable in the data object if known.") + } } - print(paste0("sigma2e estimate = ",vgm.fit$psill[1])) ## Return the sqrt of the psill as the measurement error @@ -1309,44 +1333,6 @@ setMethod("summary",signature(object="SRE"), } -## The function below checks the arguments for the function FRK. The code is self-explanatory -## This is similar, but slightly different to, .check_args1() -.check_args_wrapper <- function(f,data,basis,BAUs,est_error) { - if(!is(f,"formula")) stop("f needs to be a formula.") - if(!is(data,"list")) - stop("Please supply a list of Spatial objects.") - if(!all(sapply(data,function(x) is(x,"Spatial") | is(x,"ST")))) - stop("All data list elements need to be of class Spatial or ST.") - if(!all(sapply(data,function(x) all.vars(f)[1] %in% names(x@data)))) - stop("All data list elements to have values for the dependent variable.") - if(!est_error & !all(sapply(data,function(x) "std" %in% names(x@data)))) - stop("If observational error is not going to be estimated, - please supply a field 'std' in the data objects.") - if(!(is.null(BAUs))) { - if(!(is(BAUs,"SpatialPolygonsDataFrame") | is(BAUs,"SpatialPixelsDataFrame") | is(BAUs,"STFDF"))) - stop("BAUs should be a SpatialPolygonsDataFrame or a STFDF object") - if(!all(sapply(data,function(x) identical(proj4string(x), proj4string(BAUs))))) - stop("Please ensure all data items and BAUs have the same coordinate reference system") - if(!(all(BAUs$fs >= 0))) - stop("fine-scale variation basis function needs to be nonnegative everywhere") - if(!("fs" %in% names(BAUs@data))) { - stop("BAUs should contain a field 'fs' containing a basis - function for fine-scale variation. ") - } - if(is(BAUs,"STFDF")) if(!(is(BAUs@sp,"SpatialPolygonsDataFrame") | is(BAUs@sp,"SpatialPixelsDataFrame"))) - stop("The spatial component of the BAUs should be a SpatialPolygonsDataFrame") - if(any(sapply(data,function(x) any(names(x@data) %in% names(BAUs@data))))) - stop("Please don't have overlapping variable names in data and BAUs. All covariates need to be in the BAUs.") - if(!all(all.vars(f)[-1] %in% c(names(BAUs@data),coordnames(BAUs)))) - stop("All covariates need to be in the SpatialPolygons BAU object.") - } - - if(!(is.null(basis))) { - if(!(is(basis,"Basis") | is(basis,"TensorP_Basis"))) - stop("basis needs to be of class Basis or TensorP_Basis (package FRK)") - } -} - ## Checks arguments for the SRE() function. Code is self-explanatory .check_args1 <- function(f,data,basis,BAUs,est_error) { diff --git a/R/geometryfns.R b/R/geometryfns.R index ef60fd6..9503a0d 100644 --- a/R/geometryfns.R +++ b/R/geometryfns.R @@ -286,6 +286,9 @@ setMethod("coordnames",signature(x="STIDF"),function(x) { #' @param nonconvex_hull flag indicating whether to use \code{INLA} to generate a non-convex hull. Otherwise a convex hull is used #' @param convex convex parameter used for smoothing an extended boundary when working on a bounded domain (that is, when the object \code{data} is supplied), see details. #' @param tunit temporal unit when requiring space-time BAUs. Can be either "secs", "mins", "hours" or "days". +#' @param xlims limits of the horizontal axis (overrides automatic selection). +#' @param ylims limits of the vertical axis (overrides automatic selection). +#' @param tunit temporal unit when requiring space-time BAUs. Can be either "secs", "mins", "hours" or "days". #' @param ... currently unused #' @details \code{auto_BAUs} constructs a set of Basic Areal Units (BAUs) used both for data pre-processing and for prediction. As such, the BAUs need to be of sufficienly fine resolution so that inferences are not affected due to binning. #' @@ -326,7 +329,7 @@ setMethod("coordnames",signature(x="STIDF"),function(x) { #' @export auto_BAUs <- function(manifold, type=NULL,cellsize = NULL, isea3h_res=NULL,data=NULL,nonconvex_hull=TRUE, - convex=-0.05,tunit=NULL,...) { + convex=-0.05,tunit=NULL,xlims=NULL,ylims=NULL,...) { ## Basic checks and setting of defaults if(!(is(data,"Spatial") | is(data,"ST") | is(data,"Date") | is.null(data))) @@ -400,6 +403,12 @@ auto_BAUs <- function(manifold, type=NULL,cellsize = NULL, stop("cellsize needs to be of length equal to dimension of manifold") } + ## Check xlims and ylims + if(!is.null(xlims)) + if(!(length(xlims) == 2) & !is.numeric(xlims)) stop("xlims need to be numeric and of length 2") + if(!is.null(ylims)) + if(!(length(ylims) == 2) & !is.numeric(ylims)) stop("ylims need to be numeric and of length 2") + ## Check and set tunit if we are in a space-time setting if(grepl("ST",class(manifold)) & is.null(data) & is.null(tunit)) stop("Need to specify tunit if data is not specified in ST case") @@ -407,20 +416,24 @@ auto_BAUs <- function(manifold, type=NULL,cellsize = NULL, tunit <- .choose_BAU_tunit_from_data(data) ## Call the internal function with checked arguments - auto_BAU(manifold=manifold,type=type,cellsize=cellsize,resl=resl, - d=data,nonconvex_hull=nonconvex_hull,convex=convex,tunit=tunit) + auto_BAU(manifold=manifold,type=type,cellsize=cellsize,resl=resl,d=data, + nonconvex_hull=nonconvex_hull,convex=convex,tunit=tunit,xlims=xlims,ylims=ylims) } ## Automatically generate BAUs on the real line setMethod("auto_BAU",signature(manifold="real_line"), - function(manifold,type="grid",cellsize = 1,resl=resl,d=NULL,...) { + function(manifold,type="grid",cellsize = 1,resl=resl,d=NULL,xlims=NULL,...) { if(is.null(d)) stop("Data must be supplied when generating BAUs on a plane") crs <- CRS(proj4string(d)) # CRS of data coords <- coordinates(d) # coordinates of data - xrange <- range(coords[,1]) # range of coordinates + + if(is.null(xlims)) # if x limits are not specified + xrange <- range(coords[,1]) # range of coordinates + else xrange <- xlims # else just allocate + drangex <- diff(xrange) # range of data ## Make a SpatialPoints object. Set y = 0 so all points are on @@ -511,10 +524,10 @@ auto_BAU_time <- function (manifold,type="grid",cellsize = 1,resl=resl,d=NULL,co return(tgrid) # Return time BAUs } -## Automatically generaying BAUs on the plane +## Automatically generating BAUs on the plane setMethod("auto_BAU",signature(manifold="plane"), function(manifold,type="grid",cellsize = c(1,1),resl=resl,d=NULL, - nonconvex_hull=TRUE,convex=-0.05,...) { + nonconvex_hull=TRUE,convex=-0.05,xlims=NULL,ylims=NULL,...) { ## To arrange BAUs in a nonconvex hull we need INLA to find the domain boundary if(nonconvex_hull) @@ -539,8 +552,14 @@ setMethod("auto_BAU",signature(manifold="plane"), function(i) coordinates(d@polygons[[i]]@Polygons[[1]]))) } coord_names <- coordnames(d) # extract coordinate names - xrange <- range(coords[,1]) # x-range of coordinates - yrange <- range(coords[,2]) # y-range of coordinates + + if(is.null(xlims)) # if xlims not specified + xrange <- range(coords[,1]) # find x-range of coordinates + else xrange <- xlims # else just allocate + + if(is.null(ylims)) # if ylims not specified + yrange <- range(coords[,2]) # y-range of coordinates + else yrange = ylims # else just allocate ## Increase convex until domain is contiguous and smooth ## (i.e., the distance betweeen successive points is small) @@ -581,13 +600,16 @@ setMethod("auto_BAU",signature(manifold="plane"), drangey <- diff(yrange) # range of y ## Create x and y grid with 20% buffer and selected cellsizes - xgrid <- seq(xrange[1] - drangex*0.2, - xrange[2] + drangex*0.2, + ## If the user has specified the limits do not do buffer + bufferx <- ifelse(is.null(xlims),0.2,0) # x buffer + buffery <- ifelse(is.null(ylims),0.2,0) # y buffer + + xgrid <- seq(xrange[1] - drangex*bufferx, + xrange[2] + drangex*bufferx, by=cellsize[1]) - ygrid <- seq(yrange[1] - drangey*0.2, - yrange[2] + drangey*0.2, + ygrid <- seq(yrange[1] - drangey*buffery, + yrange[2] + drangey*buffery, by=cellsize[2]) - ## Make a SpatialPoints grid xy <- SpatialPoints(expand.grid(x=xgrid,y=ygrid), proj4string = crs) @@ -630,8 +652,9 @@ setMethod("auto_BAU",signature(manifold="plane"), # Make sure to include all pixels that contain the data points idx3 <- over(d,xy) # cannot contain NAs by definition of how xy was constructed - ## Now take the union of all the indices - xy <- xy[union(union(idx1,idx2),idx3),] + ## Now take the union of all the indices but only if xlims and ylims were not specified + if(is.null(xlims) & is.null(ylims)) + xy <- xy[union(union(idx1,idx2),idx3),] ## Add UIDs row.names(xy) <- .UIDs(xy) @@ -650,7 +673,7 @@ setMethod("auto_BAU",signature(manifold="plane"), ## Constructing BAUs on the surface of the sphere setMethod("auto_BAU",signature(manifold="sphere"), - function(manifold,type="grid",cellsize = c(1,1),resl=2,d=NULL,...) { + function(manifold,type="grid",cellsize = c(1,1),resl=2,d=NULL,xlims=NULL,ylims=NULL,...) { ## For this function d (the data) may be NULL in which case the whole sphere is covered with BAUs if(is.null(d)) # set CRS if data not provided @@ -699,13 +722,20 @@ setMethod("auto_BAU",signature(manifold="sphere"), sphere_BAUs <- SpatialPolygonsDataFrame(isea3h_sp_pol,isea3h_df_info) } else if (type == "grid") { - ## If the user wants a grid - if(!is.null(d)) { - ## If there is data then find its extent - xrange <- range(coords$lon) - yrange <- range(coords$lat) + ## If the user has specified limits assign xmin,xmax,ymin and ymin clamped to + ## the spherical coordinate limits + if(!is.null(xlims) & !is.null(ylims)) { + xmin <- max(xlims[1],-180) + xmax <- min(xlims[2],180) + ymin <- max(ylims[1],-90) + ymax <- min(ylims[2],90) + + } else if(!is.null(d)) { + ## Else if there is data try to get limits from the data + xrange <- range(coords$lon) # find x-range of coordinates + yrange <- range(coords$lat) # y-range of coordinates ## And how long/wide it is in a lon/lat sense drangex <- diff(xrange) @@ -796,7 +826,7 @@ setMethod("auto_BAU",signature(manifold="sphere"), ## Constructing BAUs on the surface of the sphere x time setMethod("auto_BAU",signature(manifold = c("STmanifold")), function(manifold,type="grid",cellsize = c(1,1,1),resl=resl,d=NULL, - nonconvex_hull=TRUE,convex=-0.05,...) { + nonconvex_hull=TRUE,convex=-0.05,xlims=NULL,ylims=NULL,...) { ## In this function user can opt to just supply a Date object, in which case ## the whole surface of the sphere is covered and the temporal part of the BAUs @@ -834,7 +864,7 @@ setMethod("auto_BAU",signature(manifold = c("STmanifold")), ## Construct the spatial BAUs spatial_BAUs <- auto_BAU(manifold=spat_manifold,cellsize=cellsize_spat, resl=resl,type=type,d=space_part,nonconvex_hull=nonconvex_hull, - convex=convex,...) + convex=convex,xlims=xlims,ylims=ylims,...) ## Construct the temporal BAUs temporal_BAUs <- auto_BAU(manifold=real_line(), cellsize=cellsize_temp, @@ -927,8 +957,8 @@ SpatialPolygonsDataFrame_to_df <- function(sp_polys,vars = names(sp_polys)) { coords <- sp_polys@polygons[[i]]@Polygons[[1]]@coords # extract coordinates row.names(coords) <- NULL # set row names to NULL coords <- data.frame(coords) # convert to data frame - poldf <- cbind(coords,id=polynames[i]) # cbind the coordinates with the ID - + poldf <- cbind(coords,id=polynames[i], # cbind the coordinates with the ID + stringsAsFactors=FALSE) # ID not factor ## remove the rownames from the data frame rownames(poldf) <- NULL @@ -937,7 +967,7 @@ SpatialPolygonsDataFrame_to_df <- function(sp_polys,vars = names(sp_polys)) { poldf }) ## rbind the data frames for each polygon into one big data frame - df_polys <- data.frame(do.call("rbind",list_polys)) + df_polys <- bind_rows(list_polys) ## merge other information from the sp_polys with the data frame (merge by polygon ID) df_polys$id <- as.character(df_polys$id) @@ -1648,3 +1678,16 @@ process_isea3h <- function(isea3h,resl) { mean(x) } else { x[1] } } + +## Takes a formula including covariates, and returns a formula with only the intercept term +.formula_no_covars <- function(f) { + labs <- all.names(f) # extract all sub-strings in formula + dep_var <- all.vars(f)[1] # extract dep. variable + idx <- which(labs == dep_var) # first char is ~, second is the dep var + if(idx > 2) { # if we have a transformation of dep var + LHS <- paste0(paste0(labs[2:idx],collapse = "("), # concatenate all transfrmations + paste0(rep(")",idx-2,collapse=""), # and close the brackets + collapse="")) + } else { LHS <- dep_var } # otherwise it's just the dep var + newf <- formula(paste0(LHS,"~1")) # now return formula without covariates +} diff --git a/R/linalgfns.R b/R/linalgfns.R index 30d5503..b4c03d9 100644 --- a/R/linalgfns.R +++ b/R/linalgfns.R @@ -14,7 +14,7 @@ #' @examples #' require(Matrix) #' cholPermute(sparseMatrix(i=c(1,1,2,2),j=c(1,2,1,2),x=c(0.1,0.2,0.2,1))) -#' @references Davis T (2011). “SPARSEINV: a MATLAB toolbox for computing the sparse inverse subset using the Takahashi equations.” http://faculty.cse.tamu.edu/davis/suitesparse.html, Online: Last accessed 01 February 2016. +#' @references Davis T (2011). ``SPARSEINV: a MATLAB toolbox for computing the sparse inverse subset using the Takahashi equations.'' http://faculty.cse.tamu.edu/davis/suitesparse.html, Online: Last accessed 01 February 2016. #' Havard Rue and Leonhard Held (2005). Gaussian Markov Random Fields: Theory and Applications. Chapman & Hall/CRC Press, Boca Raton, FL. cholPermute <- function(Q,method="amd") { n <- nrow(Q) # matrix dimension @@ -138,8 +138,8 @@ cholsolveAQinvAT <- function(A,Lp,P) { #' X <- cholPermute(Q) #' S_partial = Takahashi_Davis(Q,cholQp = X$Qpermchol,P=X$P) #' @references Takahashi, K., Fagan, J., Chin, M.-S., 1973. Formation of a sparse bus impedance matrix and -#' its application to short circuit study. 8th PICA Conf. Proc.June 4–6, Minneapolis, Minn. -#' Davis T (2011). “SPARSEINV: a MATLAB toolbox for computing the sparse inverse subset using the Takahashi equations.” http://faculty.cse.tamu.edu/davis/suitesparse.html, Online: Last accessed 01 February 2016. +#' its application to short circuit study. 8th PICA Conf. Proc.June 4--6, Minneapolis, Minn. +#' Davis T (2011). ``SPARSEINV: a MATLAB toolbox for computing the sparse inverse subset using the Takahashi equations.'' http://faculty.cse.tamu.edu/davis/suitesparse.html, Online: Last accessed 01 February 2016. Takahashi_Davis <- function(Q,return_perm_chol = 0,cholQp = matrix(0,0,0),P=0) { n <- nrow(Q) # matrix dimension diff --git a/man/auto_BAUs.Rd b/man/auto_BAUs.Rd index 00178b7..8bd8565 100644 --- a/man/auto_BAUs.Rd +++ b/man/auto_BAUs.Rd @@ -5,7 +5,8 @@ \title{Automatic BAU generation} \usage{ auto_BAUs(manifold, type = NULL, cellsize = NULL, isea3h_res = NULL, - data = NULL, nonconvex_hull = TRUE, convex = -0.05, tunit = NULL, ...) + data = NULL, nonconvex_hull = TRUE, convex = -0.05, tunit = NULL, + xlims = NULL, ylims = NULL, ...) } \arguments{ \item{manifold}{object of class \code{manifold}} @@ -24,7 +25,13 @@ auto_BAUs(manifold, type = NULL, cellsize = NULL, isea3h_res = NULL, \item{tunit}{temporal unit when requiring space-time BAUs. Can be either "secs", "mins", "hours" or "days".} +\item{xlims}{limits of the horizontal axis (overrides automatic selection).} + +\item{ylims}{limits of the vertical axis (overrides automatic selection).} + \item{...}{currently unused} + +\item{tunit}{temporal unit when requiring space-time BAUs. Can be either "secs", "mins", "hours" or "days".} } \description{ This function calls the generic function \code{auto_BAU} (not exported) after a series of checks and is the easiest way to generate a set of Basic Areal Units (BAUs) on the manifold being used; see details. diff --git a/tests/testthat/test_BAUs.R b/tests/testthat/test_BAUs.R index 1158263..e76acbb 100644 --- a/tests/testthat/test_BAUs.R +++ b/tests/testthat/test_BAUs.R @@ -63,6 +63,21 @@ test_that("plane_BAUs",{ expect_equal(names(C),c("i_idx","j_idx")) expect_equal(length(C$i_idx),nrow(binned_data)) expect_equal(length(C$j_idx),nrow(binned_data)) + + ## Limited 2D grid + Grid2D_limited <- auto_BAUs(manifold = plane(), + type="grid", + cellsize = 0.5, + data=data, + nonconvex_hull = FALSE, + xlims=c(-2,2), + ylims=c(-2,2)) + expect_is(Grid2D_limited,"SpatialPixelsDataFrame") + expect_equal(names(Grid2D_limited),c("x","y")) + expect_equal(min(Grid2D_limited@data[,1]),-2) + expect_equal(max(Grid2D_limited@data[,1]),2) + expect_equal(min(Grid2D_limited@data[,2]),-2) + expect_equal(max(Grid2D_limited@data[,2]),2) }) @@ -85,6 +100,21 @@ test_that("sphere_BAUs",{ expect_equal(names(sphere_grid@data),c("lon","lat")) expect_true(grepl("+proj=longlat",proj4string(sphere_grid))) + sphere_grid_limited <- auto_BAUs(manifold=sphere(), + type="grid", + data=NULL, + cellsize=c(20,10), + xlims=c(-100,120), + ylims=c(-80,70)) + expect_is(sphere_grid_limited,"SpatialPolygonsDataFrame") + expect_equal(nrow(sphere_grid_limited@data),165) + expect_equal(names(sphere_grid_limited@data),c("lon","lat")) + expect_true(grepl("+proj=longlat",proj4string(sphere_grid_limited))) + expect_equal(min(sphere_grid_limited@data[,1]),-90) + expect_equal(max(sphere_grid_limited@data[,1]),110) + expect_equal(min(sphere_grid_limited@data[,2]),-75) + expect_equal(max(sphere_grid_limited@data[,2]),65) + }) diff --git a/tests/testthat/test_other.R b/tests/testthat/test_other.R index a954386..4267f9d 100644 --- a/tests/testthat/test_other.R +++ b/tests/testthat/test_other.R @@ -116,3 +116,17 @@ test_that("distR works", { D2 <- as.matrix(dist(x)) #expect_lt(sum(D1-D2), 1e-12) }) + +test_that("formula no covariates works", { + f1 <- y ~ x + 1 + f2 <- sin(y) ~ log(x) - 1 + f3 <- sin(cos(y)) ~ d + r + + f1B <- .formula_no_covars(f1) + f2B <- .formula_no_covars(f2) + f3B <- .formula_no_covars(f3) + + expect_equal(f1B,formula(y~1)) + expect_equal(f2B,formula(sin(y)~1)) + expect_equal(f3B,formula(sin(cos(y))~1)) +}) diff --git a/tests/testthat/test_wrapper.R b/tests/testthat/test_wrapper.R index 236ff37..c22572c 100644 --- a/tests/testthat/test_wrapper.R +++ b/tests/testthat/test_wrapper.R @@ -1,15 +1,15 @@ -context("BAUs") - -test_that("basic_wrapper",{ - library(sp) - Z <- data.frame(x = runif(1000), y= runif(1000)) - Z$z <- sin(8*Z$x) + cos(8*Z$y) + 0.5*rnorm(100) - coordinates(Z) = ~x+y - S <- FRK(f = z~1, - list(Z), - cellsize = c(0.1,0.1), - n_EM = 2, - nonconvex_hull = FALSE) - Pred <- SRE.predict(SRE_model = S, - obs_fs = TRUE) -}) +context("BAUs") + +test_that("basic_wrapper",{ + library(sp) + Z <- data.frame(x = runif(1000), y= runif(1000)) + Z$z <- sin(8*Z$x) + cos(8*Z$y) + 0.5*rnorm(100) + coordinates(Z) = ~x+y + S <- FRK(f = z~1, + list(Z), + cellsize = c(0.02,0.02), + n_EM = 2, + nonconvex_hull = FALSE) + Pred <- SRE.predict(SRE_model = S, + obs_fs = TRUE) +}) diff --git a/vignettes/FRK_intro.pdf b/vignettes/FRK_intro.pdf index 4a4c367..8f14290 100644 Binary files a/vignettes/FRK_intro.pdf and b/vignettes/FRK_intro.pdf differ