diff --git a/.Rbuildignore b/.Rbuildignore index 3ef5878..3123049 100755 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -5,4 +5,5 @@ ^\.Rproj\.user$ ^doc$ ^Meta$ -^\.github/ +^\.lintr +^\.github diff --git a/.github/workflows/lint.yaml b/.github/workflows/lint.yaml index df2a9fb..a37aaf1 100644 --- a/.github/workflows/lint.yaml +++ b/.github/workflows/lint.yaml @@ -2,7 +2,7 @@ # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: - branches: [main, master] + branches: [main, master, build-workflow] pull_request: branches: [main, master] @@ -12,6 +12,7 @@ jobs: lint: runs-on: ubuntu-latest env: + R_LINTR_LINTER_FILE: .lintr GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - uses: actions/checkout@v3 @@ -26,8 +27,7 @@ jobs: needs: lint - name: Lint - run: lintr::lint_package(linters = lintr::linters_with_defaults(object_name_linter = NULL, - commented_code_linter = NULL, cyclocomp_linter = NULL)) + run: lintr::lint_package() shell: Rscript {0} env: LINTR_ERROR_ON_LINT: true diff --git a/.github/workflows/release.yaml b/.github/workflows/release.yaml new file mode 100644 index 0000000..f7bddc0 --- /dev/null +++ b/.github/workflows/release.yaml @@ -0,0 +1,109 @@ +# This workflow uses actions that are not certified by GitHub. +# They are provided by a third-party and are governed by +# separate terms of service, privacy policy, and support +# documentation. +# +# See https://github.com/r-lib/actions/tree/master/examples#readme for +# additional example workflows available for the R community. + +name: Release New version + +on: + push: + tags: + - '*' + branches: ["build-workflow"] + workflow_dispatch: + +jobs: + create_release: + name: Create release + runs-on: ubuntu-latest + outputs: + upload_url: ${{ steps.create_release.outputs.upload_url }} + steps: + - name: Create release + id: create_release + uses: ncipollo/release-action@v1 + with: + allowUpdates: true + build_upload_artefacts: + runs-on: ${{ matrix.config.os }} + name: ${{ matrix.config.os }} (${{ matrix.config.r }}) + strategy: + fail-fast: false + matrix: + config: + - {os: windows-latest, r: 'release'} + - {os: ubuntu-latest, r: 'release'} + - {os: macos-latest, r: 'release'} + # - {os: windows-latest, r: '4.2'} + # - {os: ubuntu-latest, r: '4.2'} + # - {os: macos-latest, r: '4.2'} + env: + R_KEEP_PKG_SOURCE: yes + steps: + # see this for details: https://msmith.de/2020/03/12/r-cmd-check-github-actions.html + - name: Configure git + run: git config --global core.autocrlf false + - uses: actions/checkout@v3 + - uses: r-lib/actions/setup-pandoc@v2 + - uses: r-lib/actions/setup-tinytex@v2 + - uses: r-lib/actions/setup-r@v2 + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + cache-version: 2 + extra-packages: | + any::ggplot2 + any::rcmdcheck + any::roxygen2 + needs: | + check + roxygen2 + - name: Read VERSION file + if: runner.os != 'macOs' + id: getversion + shell: bash + run: | + echo "VERSION=$(cat DESCRIPTION | grep -Po '(?<=Version\:\s).*')" >> $GITHUB_OUTPUT + - name: Read VERSION file (macOS) + if: runner.os == 'macOs' + id: getversion_mac + run: | + echo "VERSION=$(sed -n 's/Version:[[:space:]]*//p' DESCRIPTION | tr -d '[:space:]')" >> $GITHUB_OUTPUT + - name: Build package (Windows) + if: runner.os == 'Windows' + shell: cmd + run: R CMD build --no-build-vignettes . + - name: Build package + if: runner.os == 'Linux' || runner.os == 'macOs' + run: R CMD build --no-build-vignettes . + - name: Test Install (Windows) + if: runner.os == 'Windows' + shell: cmd + run: R CMD INSTALL --build PrestoGP_${{ steps.getversion.outputs.VERSION }}.tar.gz + - name: Test Install (Linux) + if: runner.os == 'Linux' + run: R CMD INSTALL --build PrestoGP_${{ steps.getversion.outputs.VERSION }}.tar.gz + - name: Test Install (macOs) + if: runner.os == 'macOs' + run: R CMD INSTALL --build PrestoGP_${{ steps.getversion_mac.outputs.VERSION }}.tar.gz + - uses: svenstaro/upload-release-action@v2 + if: runner.os == 'macOs' + with: + tag: ${{ github.ref }} + file: PrestoGP_${{ steps.getversion_mac.outputs.VERSION }}.tgz + asset_name: "PrestoGP_${{ steps.getversion_mac.outputs.VERSION }}-x86_64-macOs-R.${{ matrix.config.r }}.tgz" + - uses: svenstaro/upload-release-action@v2 + if: runner.os == 'Linux' + with: + tag: ${{ github.ref }} + file: "PrestoGP_${{ steps.getversion.outputs.VERSION }}_R_x86_64-pc-linux-gnu.tar.gz" + asset_name: PrestoGP_${{ steps.getversion.outputs.VERSION }}-x86_64-linux-R.${{ matrix.config.r }}.zip + - uses: svenstaro/upload-release-action@v2 + if: runner.os == 'Windows' + with: + tag: ${{ github.ref }} + file: PrestoGP_${{ steps.getversion.outputs.VERSION }}.zip + asset_name: PrestoGP_${{ steps.getversion.outputs.VERSION }}-windows-R.${{ matrix.config.r }}.zip + \ No newline at end of file diff --git a/.github/workflows/sanitizers.yaml b/.github/workflows/sanitizers.yaml index 786c154..3c964ff 100644 --- a/.github/workflows/sanitizers.yaml +++ b/.github/workflows/sanitizers.yaml @@ -84,12 +84,12 @@ jobs: run: rcmdcheck::rcmdcheck(build_args = "--no-build-vignettes", args = c("--no-codoc", "--no-examples", "--no-manual", "--ignore-vignettes"), error_on = "warning", check_dir = "check") shell: Rscript {0} - - name: Upload check results - if: failure() - uses: actions/upload-artifact@v3 - with: - name: ${{ runner.os }}-asan-rrelease-results - path: check + # - name: Upload check results + # if: failure() + # uses: actions/upload-artifact@v3 + # with: + # name: ${{ runner.os }}-asan-rrelease-results + # path: check sanitizer-check-usan: runs-on: ubuntu-latest env: @@ -162,9 +162,9 @@ jobs: run: rcmdcheck::rcmdcheck(build_args = "--no-build-vignettes", args = c("--no-codoc", "--no-examples", "--no-manual", "--ignore-vignettes"), error_on = "warning", check_dir = "check") shell: Rscript {0} - - name: Upload check results - if: failure() - uses: actions/upload-artifact@v3 - with: - name: ${{ runner.os }}-ubsan-rrelease-results - path: check \ No newline at end of file + # - name: Upload check results + # if: failure() + # uses: actions/upload-artifact@v3 + # with: + # name: ${{ runner.os }}-ubsan-rrelease-results + # path: check \ No newline at end of file diff --git a/.lintr b/.lintr new file mode 100644 index 0000000..6c1cc72 --- /dev/null +++ b/.lintr @@ -0,0 +1,13 @@ +linters: linters_with_defaults( + line_length_linter(160L), + commented_code_linter = NULL, + object_name_linter = NULL, + cyclocomp_linter = NULL, + object_length_linter = NULL, + indentation_linter( + indent = 2L, + hanging_indent_style = "never", + assignment_as_infix = FALSE + ) + ) +encoding: "UTF-8" diff --git a/DESCRIPTION b/DESCRIPTION index 9705ed4..6b3e42e 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: PrestoGP Type: Package Title: Penalized Regression for Spatio-Temporal Outcomes via Gaussian Processes -Version: 0.2.0.9021 +Version: 0.2.0.9022 Authors@R: c( person(given = "Eric", family = "Bair", @@ -45,8 +45,10 @@ Imports: rlang, mvtnorm, spam, - psych - doParallel + psych, + doParallel, + covr, + mvtnorm License: GPL-3 Encoding: UTF-8 VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index b523133..961e489 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,15 +2,11 @@ export(":=") export(.data) -export(FullSpatialModel) export(FullModel) export(Kr_pred) export(MultivariateVecchiaModel) export(PrestoGPModel) export(ST_Krig_Param_Avg) -export(SpatialModel) -export(SpatiotemporalFullModel) -export(SpatiotemporalModel) export(VecchiaModel) export(as_label) export(as_name) diff --git a/R/Log_Likelihood.R b/R/Log_Likelihood.R index d663715..9b66d60 100755 --- a/R/Log_Likelihood.R +++ b/R/Log_Likelihood.R @@ -12,8 +12,7 @@ #' #' @examples #' @noRd -negloglik_vecchia_ST <- function(logparms, res, vecchia.approx, param.seq, - scaling, nscale) { +negloglik_vecchia_ST <- function(logparms, res, vecchia.approx, param.seq, scaling, nscale) { parms <- unlog.params(logparms, param.seq, 1) locs.scaled <- vecchia.approx$locsord for (j in 1:nscale) { @@ -21,9 +20,13 @@ negloglik_vecchia_ST <- function(logparms, res, vecchia.approx, param.seq, parms[param.seq[2, 1] + j - 1] } vecchia.approx$locsord <- locs.scaled - -vecchia_likelihood(res, vecchia.approx, c(parms[1], 1, - parms[param.seq[3, 1]]), - parms[param.seq[4, 1]]) + -vecchia_likelihood( + res, vecchia.approx, c( + parms[1], 1, + parms[param.seq[3, 1]] + ), + parms[param.seq[4, 1]] + ) } #' negloglik_vecchia @@ -42,8 +45,10 @@ negloglik_vecchia_ST <- function(logparms, res, vecchia.approx, param.seq, #' @noRd negloglik_vecchia <- function(logparms, res, vecchia.approx, param.seq) { parms <- unlog.params(logparms, param.seq, 1) - -vecchia_likelihood(res, vecchia.approx, c(parms[1], parms[2], parms[3]), - parms[4]) + -vecchia_likelihood( + res, vecchia.approx, c(parms[1], parms[2], parms[3]), + parms[4] + ) } #' negloglik_full_ST @@ -95,8 +100,10 @@ negloglik.full <- function(logparams, d, y, param.seq) { params <- unlog.params(logparams, param.seq, 1) # d <- fields::rdist(locs) N <- nrow(d) - cov.mat <- params[1] * fields::Matern(d, range = params[2], - smoothness = params[3]) + + cov.mat <- params[1] * fields::Matern(d, + range = params[2], + smoothness = params[3] + ) + params[4] * diag(N) return(-1 * mvtnorm::dmvnorm(y, rep(0, N), cov.mat, log = TRUE)) } @@ -133,8 +140,7 @@ mvnegloglik <- function(logparams, vecchia.approx, y, param.seq, P) { ############################################################################## ### Flexible Spatiotemporal Multivariate Matern Negative Loglikelihood Function ########### -mvnegloglik_ST <- function(logparams, vecchia.approx, y, param.seq, P, scaling, - nscale) { +mvnegloglik_ST <- function(logparams, vecchia.approx, y, param.seq, P, scaling, nscale) { # Input- # logparams: A numeric vector of length (4*P)+(4*choose(P,2)). # To construct these parameters we unlist a list of the 7 covariance @@ -227,8 +233,7 @@ mvnegloglik.full <- function(logparams, locs, y, param.seq) { } ############################################################################## -create.cov.upper.flex <- function(P, marg.var, marg.range, marg.smooth, - nugget, R.corr) { +create.cov.upper.flex <- function(P, marg.var, marg.range, marg.smooth, nugget, R.corr) { # Create the symmetrical marginal+cross-covariance flexible matern from the # given parameters. Output is a list of the 4 Matern parameters as matrices sig2.mat <- diag(marg.var, P, P) @@ -237,20 +242,17 @@ create.cov.upper.flex <- function(P, marg.var, marg.range, marg.smooth, nugget.mat <- diag(nugget, P, P) if (P > 1) { combs <- gtools::combinations(P, 2) - for (iter in 1:nrow(combs)) { + for (iter in seq_len(nrow(combs))) { i <- combs[iter, 1] j <- combs[iter, 2] smoothness.mat[i, j] <- (marg.smooth[i] + marg.smooth[j]) / 2 - range.mat[i, j] <- 1 / sqrt(((1 / marg.range[i])^2 + - (1 / marg.range[j])^2) / 2) + range.mat[i, j] <- 1 / sqrt(((1 / marg.range[i])^2 + (1 / marg.range[j])^2) / 2) s1 <- sqrt(marg.var[i] * marg.var[j]) - s2 <- ((1 / marg.range[i])^marg.smooth[i] * - (1 / marg.range[j])^marg.smooth[j]) / + s2 <- ((1 / marg.range[i])^marg.smooth[i] * (1 / marg.range[j])^marg.smooth[j]) / ((1 / range.mat[i, j])^(2 * smoothness.mat[i, j])) - s3 <- gamma(smoothness.mat[i, j]) / (sqrt(gamma(marg.smooth[i])) * - sqrt(gamma(marg.smooth[j]))) + s3 <- gamma(smoothness.mat[i, j]) / (sqrt(gamma(marg.smooth[i])) * sqrt(gamma(marg.smooth[j]))) s4 <- R.corr[iter] sig2.mat[i, j] <- s1 * s2 * s3 * s4 } @@ -284,7 +286,7 @@ cat.covariances <- function(locs.list, sig2, range, smoothness, nugget) { l <- length(locs.list) combs <- gtools::combinations(l, 2, repeats.allowed = TRUE) - for (iter in 1:nrow(combs)) { + for (iter in seq_len(nrow(combs))) { i <- combs[iter, 1] j <- combs[iter, 2] # d <- fields::rdist.earth(locs.list[[i]],locs.list[[j]],miles = FALSE) @@ -292,12 +294,16 @@ cat.covariances <- function(locs.list, sig2, range, smoothness, nugget) { # Calculate the covariance matrix - if/then based on its location in the super-matrix N <- nrow(d) if (i == j) { # To accomodate varying size outcomes- the nugget is not included on cross-covariances - cov.mat.ij <- sig2[i, j] * geoR::matern(d, phi = range[i, j], kappa = - smoothness[i, j]) + + cov.mat.ij <- sig2[i, j] * geoR::matern(d, + phi = range[i, j], kappa = + smoothness[i, j] + ) + nugget[i, j] * diag(N) } else { - cov.mat.ij <- sig2[i, j] * geoR::matern(d, phi = range[i, j], kappa = - smoothness[i, j]) + cov.mat.ij <- sig2[i, j] * geoR::matern(d, + phi = range[i, j], kappa = + smoothness[i, j] + ) } @@ -328,8 +334,7 @@ cat.covariances <- function(locs.list, sig2, range, smoothness, nugget) { ############################################################################## ### Create the likelihood initial values ######### -create.initial.values.flex <- function(marg.var, marg.range, marg.smooth, - nugget, R.corr, P) { +create.initial.values.flex <- function(marg.var, marg.range, marg.smooth, nugget, R.corr, P) { # Log-transform the covariance parameters and arrange in the proper order # for the likelihood function logparams.init <- c( diff --git a/R/PrestoGP_CreateU_Multivariate.R b/R/PrestoGP_CreateU_Multivariate.R index bb099f7..721bb3a 100755 --- a/R/PrestoGP_CreateU_Multivariate.R +++ b/R/PrestoGP_CreateU_Multivariate.R @@ -31,14 +31,14 @@ max_min_ordering <- function(locs, dist_func) { # find the point closest to the mean of all points dists <- dist_func(center, locs) first <- which.min(dists) - unsolved <- 1:nrow(locs) + unsolved <- seq_len(nrow(locs)) unsolved <- unsolved[-first] order <- c(first) while (length(order) < nrow(locs)) { max_min <- 0 max_min_i <- unsolved[1] - in_order <- locs[order[1:length(order)], ] + in_order <- locs[order[seq_along(order)], ] dim(in_order) <- c(length(order), ncol(locs)) for (i in unsolved) { loc_i <- locs[i, ] @@ -66,8 +66,7 @@ max_min_ordering <- function(locs, dist_func) { #' @param dist_func Any distance function with a signature of dist(query_location, locations_matrix) #' #' @return A vector containing the indices of the neighbors -knn_indices <- function(ordered_locs, query, n_neighbors, - dist_func, dist_func_code) { +knn_indices <- function(ordered_locs, query, n_neighbors, dist_func, dist_func_code) { if (dist_func_code == "custom") { dists <- dist_func(query, ordered_locs) dists_order <- order(dists) @@ -90,9 +89,9 @@ knn_indices <- function(ordered_locs, query, n_neighbors, #' @param n_neighbors The number of neighbors to find (K) for each location #' @param dist_func Any distance function with a signature of dist(query_location, locations_matrix) #' -#' @return A list containing two matrices, each with one row per location: an indices matrix with the indices of nearest neighbors for each location, and a distance matrix with the associated distances -sparseNN <- function(ordered_locs, n_neighbors, - dist_func, dist_func_code, ordered_locs_pred = NULL) { +#' @return A list containing two matrices, each with one row per location: +#' an indices matrix with the indices of nearest neighbors for each location, and a distance matrix with the associated distances +sparseNN <- function(ordered_locs, n_neighbors, dist_func, dist_func_code, ordered_locs_pred = NULL) { ee <- min(apply(ordered_locs, 2, stats::sd)) n <- nrow(ordered_locs) ordered_locs <- ordered_locs + matrix( @@ -100,25 +99,33 @@ sparseNN <- function(ordered_locs, n_neighbors, stats::rnorm(n * ncol(ordered_locs)), n, ncol(ordered_locs) ) - indices_matrix <- matrix(data = NA, nrow = nrow(ordered_locs), - ncol = n_neighbors) - distances_matrix <- matrix(data = NA, nrow = nrow(ordered_locs), - ncol = n_neighbors) + indices_matrix <- matrix( + data = NA, nrow = nrow(ordered_locs), + ncol = n_neighbors + ) + distances_matrix <- matrix( + data = NA, nrow = nrow(ordered_locs), + ncol = n_neighbors + ) for (row in 1:n_neighbors) { # for the locations from 1 to n_neighbors, use the entire locs list to find the neighbors - nn <- knn_indices(ordered_locs[1: - (n_neighbors + 1), , drop = FALSE][-row, , - drop = FALSE], - ordered_locs[row, , drop = FALSE], n_neighbors, - dist_func, dist_func_code) + nn <- knn_indices( + ordered_locs[1:(n_neighbors + 1), , drop = FALSE][-row, , + drop = FALSE + ], + ordered_locs[row, , drop = FALSE], n_neighbors, + dist_func, dist_func_code + ) indices_matrix[row, 1:n_neighbors] <- nn$indices[1:n_neighbors] distances_matrix[row, 1:n_neighbors] <- nn$distances[1:n_neighbors] } for (row in (n_neighbors + 1):nrow(ordered_locs)) { # get the m nearest neighbors from the locs before this one in the max-min order - nn <- knn_indices(ordered_locs[1:(row - 1), , drop = FALSE], - ordered_locs[row, , drop = FALSE], n_neighbors, - dist_func, dist_func_code) + nn <- knn_indices( + ordered_locs[1:(row - 1), , drop = FALSE], + ordered_locs[row, , drop = FALSE], n_neighbors, + dist_func, dist_func_code + ) indices_matrix[row, 1:n_neighbors] <- nn$indices[1:n_neighbors] distances_matrix[row, 1:n_neighbors] <- nn$distances[1:n_neighbors] } @@ -131,7 +138,7 @@ sparseNN <- function(ordered_locs, n_neighbors, data = NA, nrow = nrow(ordered_locs_pred), ncol = n_neighbors ) - for (row in 1:nrow(ordered_locs_pred)) { + for (row in seq_len(nrow(ordered_locs_pred))) { nn <- knn_indices( ordered_locs, ordered_locs_pred[row, , drop = FALSE], n_neighbors, @@ -170,7 +177,7 @@ calc.q <- function(nn.obj, firstind.pred) { for (j in 2:m) { cur.k <- cur.q[j] cur.qy <- intersect(q.y[[cur.k]], cur.q) - if (length(cur.qy) > length(best.qy) & cur.k < firstind.pred) { + if (length(cur.qy) > length(best.qy) && cur.k < firstind.pred) { best.k <- cur.k best.qy <- cur.qy } @@ -188,9 +195,9 @@ calc.q <- function(nn.obj, firstind.pred) { #' @export vecchia_Mspecify <- function(locs.list, m, locs.list.pred = NULL, - dist.func = NULL, - ordering.pred = c("obspred", "general"), - pred.cond = c("independent", "general")) { + dist.func = NULL, + ordering.pred = c("obspred", "general"), + pred.cond = c("independent", "general")) { ordering.pred <- match.arg(ordering.pred) pred.cond <- match.arg(pred.cond) @@ -245,7 +252,7 @@ vecchia_Mspecify <- function(locs.list, m, locs.list.pred = NULL, loc.order <- max_min_ordering(locs.all, dist.func) loc.order <- c(unique(loc.order), setdiff(1:n, loc.order)) } else { - if (is.null(locs.list.pred) | ordering.pred == "general") { + if (is.null(locs.list.pred) || ordering.pred == "general") { loc.order <- GPvecchia::order_maxmin_exact(locs.all) # I am not sure why the next two lines are here. I added them because # similar code exists in the GPvecchia package. But I don't know why @@ -274,7 +281,7 @@ vecchia_Mspecify <- function(locs.list, m, locs.list.pred = NULL, # is non-deterministic, so there may be some slight differences # between the output of this function and the output of createU # in the GPvecchia package. - if (is.null(locs.list.pred) | pred.cond == "general") { + if (is.null(locs.list.pred) || pred.cond == "general") { nn.mat <- sparseNN(olocs, m, dist.func, dist.func.code) } else { nn.mat <- sparseNN( @@ -282,7 +289,7 @@ vecchia_Mspecify <- function(locs.list, m, locs.list.pred = NULL, olocs[-(1:n), , drop = FALSE] ) } - last.obs <- max((1:length(obs))[obs]) + last.obs <- max(which(obs)) q.list <- calc.q(nn.mat$indices, last.obs + 1) return(list( @@ -348,8 +355,7 @@ createUMultivariate <- function(vec.approx, params, cov_func = NULL) { ajj <- 1 / rangep[ondx[2]] aij <- sqrt((aii^2 + ajj^2) / 2) K1 <- rho.mat[ondx[1], ondx[2]] * sqrt(sig2[ondx[1]]) * sqrt(sig2[ondx[2]]) * - aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * - gamma(vjj))) * + aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * cov_func(dist_func(olocs[1, , drop = FALSE], olocs[2, , drop = FALSE], ), smoothness = vij, alpha = aij ) @@ -404,8 +410,7 @@ createUMultivariate <- function(vec.approx, params, cov_func = NULL) { ajj <- 1 / rangep[ondx[2]] aij <- sqrt((aii^2 + ajj^2) / 2) K1 <- rho.mat[ondx[1], ondx[2]] * sqrt(sig2[ondx[1]]) * sqrt(sig2[ondx[2]]) * - aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * - gamma(vjj))) * + aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * cov_func(dist_func(olocs[1, , drop = FALSE], olocs[2, , drop = FALSE], ), smoothness = vij, alpha = aij ) @@ -416,6 +421,7 @@ createUMultivariate <- function(vec.approx, params, cov_func = NULL) { U1[7, ] <- c(1, 3, -1 * bi * ri^(-1 / 2)) # U[3,3] <- ri^(-1/2) # U[1,3] <- -1*bi*ri^(-1/2) + i <- NULL # lintr requirement U2 <- foreach(i = 3:n, .combine = rbind) %dopar% { # U[2*i,2*i] <- nugget[ondx[i]]^(-1/2) # U[2*i-1,2*i] <- -1*U[2*i,2*i] @@ -443,8 +449,7 @@ createUMultivariate <- function(vec.approx, params, cov_func = NULL) { # positive definite. See equation (9) in Apanasovich (2011). K1[j] <- rho.mat[ondx[i], ondx[cur.q[j]]] * sqrt(sig2[ondx[i]]) * sqrt(sig2[ondx[cur.q[j]]]) * - aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * - sqrt(gamma(vii) * gamma(vjj))) * + aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * cov_func( dist_func( olocs[i, , drop = FALSE], diff --git a/R/PrestoGP_Model.R b/R/PrestoGP_Model.R index 5b5a26a..3c3fc28 100755 --- a/R/PrestoGP_Model.R +++ b/R/PrestoGP_Model.R @@ -68,8 +68,22 @@ setMethod("initialize", "PrestoGPModel", function(.Object, ...) { }) setGeneric("show_theta", function(object, Y_names) standardGeneric("show_theta")) -setGeneric("prestogp_fit", function(model, Y, X, locs, scaling = NULL, apanasovich = FALSE, covparams = NULL, beta.hat = NULL, tol = 0.999999, max_iters = 100, verbose = FALSE, optim.method = "Nelder-Mead", optim.control = list(trace = 0, reltol = 1e-3, maxit = 5000), parallel = FALSE, foldid = NULL) standardGeneric("prestogp_fit")) -setGeneric("prestogp_predict", function(model, X = "matrix", locs = "matrix", m = "numeric", ordering.pred = c("obspred", "general"), pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) standardGeneric("prestogp_predict")) +setGeneric( + "prestogp_fit", + function(model, Y, X, locs, scaling = NULL, apanasovich = FALSE, + covparams = NULL, beta.hat = NULL, tol = 0.999999, max_iters = 100, verbose = FALSE, + optim.method = "Nelder-Mead", optim.control = list(trace = 0, reltol = 1e-3, maxit = 5000), + parallel = FALSE, foldid = NULL) { + standardGeneric("prestogp_fit") + } +) +setGeneric( + "prestogp_predict", + function(model, X = "matrix", locs = "matrix", m = "numeric", ordering.pred = c("obspred", "general"), + pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) { + standardGeneric("prestogp_predict") + } +) setGeneric("calc_covparams", function(model, locs, Y, covparams) standardGeneric("calc_covparams")) setGeneric("specify", function(model, ...) standardGeneric("specify")) setGeneric("compute_residuals", function(model, Y, Y.hat) standardGeneric("compute_residuals")) @@ -96,7 +110,7 @@ setMethod( cat("Covariance Parameters:\n") Y_names <- colnames(object@Y_train) if (is.null(Y_names)) { - Y_names <- unlist(lapply(1:ncol(object@Y_train), function(x) { + Y_names <- unlist(lapply(seq_len(ncol(object@Y_train)), function(x) { paste("Outcome", x) })) } @@ -111,21 +125,21 @@ setMethod( # TODO compare to zero within a tolerance # nnz_betas <- lapply(object@beta, 2, function(x){which(x != 0.0)}) nnz_betas <- list() - for (col in 1:ncol(object@Y_train)) { + for (col in seq_len(ncol(object@Y_train))) { nnz_betas <- append(nnz_betas, list(which(object@beta[, col] != 0.0))) } X_names <- colnames(object@X_train) if (is.null(X_names)) { - X_names <- unlist(lapply(1:ncol(object@X_train), function(x) { + X_names <- unlist(lapply(seq_len(ncol(object@X_train)), function(x) { paste("Ind. Variable", x) })) } X_names <- append("Intercept", X_names) - for (i in 1:ncol(object@Y_train)) { + for (i in seq_len(ncol(object@Y_train))) { cat(Y_names[i], " Parameters:\n") beta_summary <- data.frame(matrix(ncol = 4, nrow = 0, dimnames = list(NULL, c("Parameter", "Estimate", "Standard Error", "Walds P-value")))) # for(nnz in nnz_betas[i]){ - for (j in 1:length(nnz_betas[[i]])) { + for (j in seq_along(nnz_betas[[i]])) { nnz <- nnz_betas[[i]][[j]] walds <- wald.test(covm * mse[i, i], object@beta[, i], Terms = nnz) std_err <- sqrt(diag(covm) * mse[i, i]) @@ -150,13 +164,13 @@ setMethod( function(object, Y_names) { theta_name_arr <- theta_names(object) theta_summary <- data.frame(matrix(ncol = ncol(object@Y_train) + 1, nrow = length(theta_name_arr), dimnames = list(NULL, c("Parameter", Y_names)))) - for (i in 1:length(theta_name_arr)) { + for (i in seq_along(theta_name_arr)) { theta_row <- object@covparams[((i - 1) * ncol(object@Y_train) + 1):(i * ncol(object@Y_train))] - for (j in 1:ncol(object@Y_train)) { + for (j in seq_len(ncol(object@Y_train))) { theta_summary[i, j + 1] <- theta_row[j] } } - for (j in 1:length(theta_name_arr)) { + for (j in seq_along(theta_name_arr)) { theta_summary[j, 1] <- theta_name_arr[j] } print(theta_summary, row.names = FALSE) @@ -201,19 +215,19 @@ setMethod( setMethod( "prestogp_fit", "PrestoGPModel", function(model, Y, X, locs, scaling = NULL, apanasovich = NULL, - covparams = NULL, beta.hat = NULL, tol = 0.999999, - max_iters = 100, verbose = FALSE, optim.method = "Nelder-Mead", - optim.control = list(trace = 0, reltol = 1e-3, maxit = 5000), - parallel = FALSE, foldid = NULL) { + covparams = NULL, beta.hat = NULL, tol = 0.999999, + max_iters = 100, verbose = FALSE, optim.method = "Nelder-Mead", + optim.control = list(trace = 0, reltol = 1e-3, maxit = 5000), + parallel = FALSE, foldid = NULL) { model <- check_input(model, Y, X, locs) if (!is.null(beta.hat)) { - if (!is.vector(beta.hat) | !is.numeric(beta.hat)) { - stop("beta.hat parameter must be a numeric vector") - } - if (length(beta.hat) != (ncol(model@X_train) + 1)) { - stop("Length of beta.hat must match the number of predictors") - } - beta.hat <- as.matrix(beta.hat) + if (!is.vector(beta.hat) | !is.numeric(beta.hat)) { + stop("beta.hat parameter must be a numeric vector") + } + if (length(beta.hat) != (ncol(model@X_train) + 1)) { + stop("Length of beta.hat must match the number of predictors") + } + beta.hat <- as.matrix(beta.hat) } if (!is.numeric(tol)) { stop("tol must be numeric") @@ -221,20 +235,20 @@ setMethod( if (length(tol) != 1) { stop("tol must be a scalar") } - if (tol<=0 | tol>1) { + if (tol <= 0 | tol > 1) { stop("tol must satisfy 0 1) { - model@linear_model <- cv.glmnet(as.matrix(model@X_tilde), as.matrix(model@y_tilde), family = "mgaussian", alpha = model@alpha, parallel = parallel, foldid = foldid) + model@linear_model <- cv.glmnet( + as.matrix(model@X_tilde), + as.matrix(model@y_tilde), + family = "mgaussian", + alpha = model@alpha, + parallel = parallel, + foldid = foldid + ) } else { model@linear_model <- cv.glmnet(as.matrix(model@X_tilde), as.matrix(model@y_tilde), alpha = model@alpha, parallel = parallel, foldid = foldid) } @@ -360,8 +378,8 @@ sparseToDenseBeta <- function(linear_model) { } beta_construct <- matrix(data = 0, nrow = coefs[[1]]@Dim[1], ncol = length(coefs)) # coefs[[1]]@Dim[1]+2s because dgCMatrix is 0 offset, and we want to include intercept - for (i in 1:length(coefs)) { - for (j in 1:length(coefs[[i]]@i)) { + for (i in seq_along(coefs)) { + for (j in seq_along(coefs[[i]]@i)) { k <- coefs[[i]]@i[j] # beta_construct[k+1,i] <- coefs[[i]]@x[j] beta_construct[k + 1, i] <- coefs[[i]]@x[j] @@ -418,64 +436,65 @@ setMethod("calc_covparams", "PrestoGPModel", function(model, locs, Y, covparams) } pseq <- create.param.sequence(P, model@nscale) if (is.null(covparams)) { - col.vars <- rep(NA, P) - D.sample.bar <- rep(NA, model@nscale * P) - for (i in 1:P) { - col.vars[i] <- var(Y[[i]]) - N <- length(Y[[i]]) - # TODO find a better way to compute initial spatial range - for (j in 1:model@nscale) { - d.sample <- sample(1:N, max(2, ceiling(N / 50)), replace = FALSE) - D.sample <- rdist(locs[[i]][d.sample, model@scaling == j]) - D.sample.bar[(i - 1) * model@nscale + j] <- mean(D.sample) / 4 - } + col.vars <- rep(NA, P) + D.sample.bar <- rep(NA, model@nscale * P) + for (i in 1:P) { + col.vars[i] <- var(Y[[i]]) + N <- length(Y[[i]]) + # TODO find a better way to compute initial spatial range + for (j in 1:model@nscale) { + d.sample <- sample(1:N, max(2, ceiling(N / 50)), replace = FALSE) + D.sample <- rdist(locs[[i]][d.sample, model@scaling == j]) + D.sample.bar[(i - 1) * model@nscale + j] <- mean(D.sample) / 4 } - model@logparams <- create.initial.values.flex( - c(0.9 * col.vars), # marginal variance - D.sample.bar, # range - rep(0.5, P), # smoothness - c(.1 * col.vars), # nuggets - rep(0, choose(P, 2)), - P - ) + } + model@logparams <- create.initial.values.flex( + c(0.9 * col.vars), # marginal variance + D.sample.bar, # range + rep(0.5, P), # smoothness + c(.1 * col.vars), # nuggets + rep(0, choose(P, 2)), + P + ) } else { - if (P==1) { - if (length(covparams) != pseq[4,2]) { - stop("Incorrect number of parameters in covparams") - } - } else { - if (length(covparams) != pseq[5,2]) { - stop("Incorrect number of parameters in covparams") - } - } - init.var <- covparams[pseq[1,1]:pseq[1,2]] - init.range <- covparams[pseq[2,1]:pseq[2,2]] - init.smooth <- covparams[pseq[3,1]:pseq[3,2]] - init.nugget <- covparams[pseq[4,1]:pseq[4,2]] - if (P>1) { - init.corr <- covparams[pseq[5,1]:pseq[5,2]] + if (P == 1) { + if (length(covparams) != pseq[4, 2]) { + stop("Incorrect number of parameters in covparams") } - else { - init.corr <- 0 - } - if (sum(init.var<=0)>0) { - stop("Initial variance estimates must be positive") - } - if (sum(init.range<=0)>0) { - stop("Initial range estimates must be positive") - } - if (sum(init.nugget<=0)>0) { - stop("Initial nugget estimates must be positive") - } - if (sum(init.smooth<=0)>0 | sum(init.smooth>=2.5)>0) { - stop("Initial smoothness estimates must be between 0 and 2.5") - } - if (sum(init.corr < -1)>0 | sum(init.corr > 1)>0) { - stop("Initial correlation estimates must be between -1 and 1") + } else { + if (length(covparams) != pseq[5, 2]) { + stop("Incorrect number of parameters in covparams") } - model@logparams <- create.initial.values.flex(init.var, init.range, - init.smooth, init.nugget, - init.corr, P) + } + init.var <- covparams[pseq[1, 1]:pseq[1, 2]] + init.range <- covparams[pseq[2, 1]:pseq[2, 2]] + init.smooth <- covparams[pseq[3, 1]:pseq[3, 2]] + init.nugget <- covparams[pseq[4, 1]:pseq[4, 2]] + if (P > 1) { + init.corr <- covparams[pseq[5, 1]:pseq[5, 2]] + } else { + init.corr <- 0 + } + if (sum(init.var <= 0) > 0) { + stop("Initial variance estimates must be positive") + } + if (sum(init.range <= 0) > 0) { + stop("Initial range estimates must be positive") + } + if (sum(init.nugget <= 0) > 0) { + stop("Initial nugget estimates must be positive") + } + if (sum(init.smooth <= 0) > 0 | sum(init.smooth >= 2.5) > 0) { + stop("Initial smoothness estimates must be between 0 and 2.5") + } + if (sum(init.corr < -1) > 0 | sum(init.corr > 1) > 0) { + stop("Initial correlation estimates must be between -1 and 1") + } + model@logparams <- create.initial.values.flex( + init.var, init.range, + init.smooth, init.nugget, + init.corr, P + ) } model@param_sequence <- pseq model <- transform_covariance_parameters(model) @@ -495,12 +514,11 @@ setMethod("scale_locs", "PrestoGPModel", function(model, locs) { return(locs) } else { locs.out <- locs - for (i in 1:length(locs)) { + for (i in seq_along(locs)) { for (j in 1:model@nscale) { locs.out[[i]][, model@scaling == j] <- locs[[i]][, model@scaling == j] / - model@covparams[model@param_sequence[2, 1] + - model@nscale * (i - 1) + j - 1] + model@covparams[model@param_sequence[2, 1] + model@nscale * (i - 1) + j - 1] } } return(locs.out) @@ -513,25 +531,20 @@ setMethod("transform_covariance_parameters", "PrestoGPModel", function(model) { model@covparams <- c( exp(model@logparams[1:model@param_sequence[2, 2]]), gtools::inv.logit( - model@logparams[model@param_sequence[3, 1]: - model@param_sequence[3, 2]], + model@logparams[model@param_sequence[3, 1]:model@param_sequence[3, 2]], 0, 2.5 ), - exp(model@logparams[model@param_sequence[4, 1]: - model@param_sequence[4, 2]]), - tanh(model@logparams[model@param_sequence[5, 1]: - model@param_sequence[5, 2]]) + exp(model@logparams[model@param_sequence[4, 1]:model@param_sequence[4, 2]]), + tanh(model@logparams[model@param_sequence[5, 1]:model@param_sequence[5, 2]]) ) } else { model@covparams <- c( exp(model@logparams[1:model@param_sequence[2, 2]]), gtools::inv.logit( - model@logparams[model@param_sequence[3, 1]: - model@param_sequence[3, 2]], + model@logparams[model@param_sequence[3, 1]:model@param_sequence[3, 2]], 0, 2.5 ), - exp(model@logparams[model@param_sequence[4, 1]: - model@param_sequence[4, 2]]), 1 + exp(model@logparams[model@param_sequence[4, 1]:model@param_sequence[4, 2]]), 1 ) } invisible(model) diff --git a/R/PrestoGP_Multivariate_Vecchia.R b/R/PrestoGP_Multivariate_Vecchia.R index 08dab5b..642c125 100755 --- a/R/PrestoGP_Multivariate_Vecchia.R +++ b/R/PrestoGP_Multivariate_Vecchia.R @@ -45,75 +45,81 @@ setMethod("initialize", "MultivariateVecchiaModel", function(.Object, n_neighbor #' prediction <- prestogp_predict(model, X.test, locs.test) #' Vec.mean <- prediction[[1]] #' Vec.sds <- prediction[[2]] -setMethod("prestogp_predict", "MultivariateVecchiaModel", function(model, X, locs, m = NULL, ordering.pred = c("obspred", "general"), pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) { - # validate parameters - ordering.pred <- match.arg(ordering.pred) - pred.cond <- match.arg(pred.cond) - return.values <- match.arg(return.values) - X <- check_input_pred(model, X, locs) - if (is.null(m)) { # m defaults to the value used for training - m <- model@n_neighbors - } - if (m < model@min_m) { - stop(paste("m must be at least ", model@min_m, sep = "")) - } - if (m >= nrow(model@X_train)) { - warning("Conditioning set size m chosen to be >=n. Changing to m=n-1") - m <- nrow(model@X_train) - 1 - } - - # Vecchia prediction at new locations - Vecchia.Pred <- predict(model@linear_model, newx = X, s = model@linear_model$lambda[model@lambda_1se_idx]) - # Vecchia trend prediction at observed data - Vecchia.hat <- predict(model@linear_model, newx = model@X_train, s = model@linear_model$lambda[model@lambda_1se_idx]) +setMethod("prestogp_predict", "MultivariateVecchiaModel", + function(model, X, locs, m = NULL, ordering.pred = c("obspred", "general"), pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) { + # validate parameters + ordering.pred <- match.arg(ordering.pred) + pred.cond <- match.arg(pred.cond) + return.values <- match.arg(return.values) + X <- check_input_pred(model, X, locs) + if (is.null(m)) { # m defaults to the value used for training + m <- model@n_neighbors + } + if (m < model@min_m) { + stop(paste("m must be at least ", model@min_m, sep = "")) + } + if (m >= nrow(model@X_train)) { + warning("Conditioning set size m chosen to be >=n. Changing to m=n-1") + m <- nrow(model@X_train) - 1 + } - # Test set prediction - res <- model@Y_train - Vecchia.hat + # Vecchia prediction at new locations + Vecchia.Pred <- predict(model@linear_model, newx = X, s = model@linear_model$lambda[model@lambda_1se_idx]) + # Vecchia trend prediction at observed data + Vecchia.hat <- predict(model@linear_model, newx = model@X_train, s = model@linear_model$lambda[model@lambda_1se_idx]) - locs.train.scaled <- scale_locs(model, model@locs_train) - locs.scaled <- scale_locs(model, locs) - vec.approx.test <- vecchia_Mspecify(locs.train.scaled, m, - locs.list.pred = locs.scaled, - ordering.pred = ordering.pred, - pred.cond = pred.cond - ) + # Test set prediction + res <- model@Y_train - Vecchia.hat - ## carry out prediction - if (!model@apanasovich) { - params <- model@covparams - param.seq <- model@param_sequence - pred <- vecchia_Mprediction(res, vec.approx.test, - c( - params[1:param.seq[1, 2]], - rep(1, param.seq[2, 2] - param.seq[2, 1] + 1), - params[param.seq[3, 1]: - param.seq[5, 2]] - ), - return.values = return.values + locs.train.scaled <- scale_locs(model, model@locs_train) + locs.scaled <- scale_locs(model, locs) + vec.approx.test <- vecchia_Mspecify(locs.train.scaled, m, + locs.list.pred = locs.scaled, + ordering.pred = ordering.pred, + pred.cond = pred.cond ) - } else { - pred <- vecchia_Mprediction(res, vec.approx.test, model@covparams, - return.values = return.values - ) - } - # prediction function can return both mean and sds - # returns a list with elements mu.pred,mu.obs,var.pred,var.obs,V.ord - Vec.mean <- pred$mu.pred + Vecchia.Pred # residual + mean trend - if (return.values == "mean") { - return.list <- list(means = Vec.mean) - } else { - warning("Variance estimates do not include model fitting variance and are anticonservative. Use with caution.") - vec.var <- pred$var.pred - for (i in 1:length(locs)) { - vec.sds[ndx.out == i] <- sqrt(vec.var[ndx.out == i] + - model@covparams[model@param_sequence[4, i]]) + ## carry out prediction + if (!model@apanasovich) { + params <- model@covparams + param.seq <- model@param_sequence + pred <- vecchia_Mprediction(res, vec.approx.test, + c( + params[1:param.seq[1, 2]], + rep(1, param.seq[2, 2] - param.seq[2, 1] + 1), + params[param.seq[3, 1]:param.seq[5, 2]] + ), + return.values = return.values + ) + } else { + pred <- vecchia_Mprediction(res, vec.approx.test, model@covparams, + return.values = return.values + ) } - return.list <- list(means = Vec.mean, sds = vec.sds) - } - return(return.list) -}) + # prediction function can return both mean and sds + # returns a list with elements mu.pred,mu.obs,var.pred,var.obs,V.ord + Vec.mean <- pred$mu.pred + Vecchia.Pred # residual + mean trend + if (return.values == "mean") { + return.list <- list(means = Vec.mean) + } else { + warning("Variance estimates do not include model fitting variance and are anticonservative. Use with caution.") + Vec.sds <- pred$var.pred + ndx.out <- NULL + for (i in seq_along(locs)) { + ndx.out <- c(ndx.out, rep(i, nrow(locs[[i]]))) + } + for (i in seq_along(locs)) { + Vec.sds[ndx.out == i] <- sqrt(Vec.sds[ndx.out == i] + + model@covparams[model@param_sequence[4, + i]]) + } + return.list <- list(means = Vec.mean, sds = Vec.sds) + } + + return(return.list) + } +) setMethod("check_input", "MultivariateVecchiaModel", function(model, Y, X, locs) { if (!is.list(locs)) { @@ -131,7 +137,7 @@ setMethod("check_input", "MultivariateVecchiaModel", function(model, Y, X, locs) if (length(locs) != length(X)) { stop("locs and X must have the same length") } - for (i in 1:length(locs)) { + for (i in seq_along(locs)) { if (!is.matrix(locs[[i]])) { stop("Each locs must be a matrix") } @@ -182,14 +188,12 @@ setMethod("check_input_pred", "MultivariateVecchiaModel", function(model, X, loc if (length(locs) != length(model@locs_train)) { stop("Training and test set locs must have the same length") } - for (i in 1:length(locs)) { + for (i in seq_along(locs)) { if (!is.matrix(locs[[i]])) { stop("Each locs must be a matrix") } - if (i > 1) { - if (ncol(locs[[i]]) != ncol(model@locs_train[[1]])) { - stop("All locs must have the same number of columns as locs_train") - } + if (ncol(locs[[i]]) != ncol(model@locs_train[[1]])) { + stop("All locs must have the same number of columns as locs_train") } if (!is.matrix(X[[i]])) { stop("Each X must be a matrix") @@ -199,13 +203,12 @@ setMethod("check_input_pred", "MultivariateVecchiaModel", function(model, X, loc } } if (length(X) == 1) { - X <- X[[1]] - } - else { - X <- psych::superMatrix(X) + X <- X[[1]] + } else { + X <- psych::superMatrix(X) } if (ncol(X) != ncol(model@X_train)) { - stop("X and X_train must have the same number of predictors") + stop("X and X_train must have the same number of predictors") } return(X) }) @@ -216,12 +219,11 @@ setMethod("specify", "MultivariateVecchiaModel", function(model) { model@vecchia_approx <- vecchia_Mspecify(locs.scaled, model@n_neighbors) if (!model@apanasovich) { olocs.scaled <- model@vecchia_approx$locsord - for (i in 1:length(locs)) { + for (i in seq_along(locs)) { for (j in 1:model@nscale) { olocs.scaled[model@vecchia_approx$ondx == i, model@scaling == j] <- olocs.scaled[model@vecchia_approx$ondx == i, model@scaling == j] * - model@covparams[model@param_sequence[2, 1] + - model@nscale * (i - 1) + j - 1] + model@covparams[model@param_sequence[2, 1] + model@nscale * (i - 1) + j - 1] } } model@vecchia_approx$locsord <- olocs.scaled @@ -292,12 +294,8 @@ setMethod("transform_data", "MultivariateVecchiaModel", function(model, Y, X) { vecchia.approx = vecchia.approx, c( params[1:param.seq[1, 2]], - rep( - 1, - param.seq[2, 2] - param.seq[2, 1] + 1 - ), - params[param.seq[3, 1]: - param.seq[5, 2]] + rep(1, param.seq[2, 2] - param.seq[2, 1] + 1), + params[param.seq[3, 1]:param.seq[5, 2]] ) ) } else { diff --git a/R/PrestoGP_Util_Functions.R b/R/PrestoGP_Util_Functions.R old mode 100644 new mode 100755 index 0a7db69..b77c14d --- a/R/PrestoGP_Util_Functions.R +++ b/R/PrestoGP_Util_Functions.R @@ -3,7 +3,7 @@ ################################################################################ SCAD_Penalty_Loglike <- function(beta.in, lambda) { penalty <- matrix(NA, nrow = nrow(beta.in), ncol = ncol(beta.in)) - for (j in 1:nrow(beta.in)) { + for (j in seq_len(nrow(beta.in))) { beta.j <- beta.in[j, ] idx1 <- abs(beta.j) <= lambda idx2 <- abs(beta.j) > lambda & abs(beta.j) <= 3.7 * lambda @@ -54,7 +54,7 @@ Kr_pred <- function(new_coords, obs_coords, Y_obs, cov.pars, NN) { df_new <- new_coords_scaled df_obs <- obs_coords_scaled - for (i in 1:nrow(new_coords)) { + for (i in seq_len(nrow(new_coords))) { # print(i) locs.test.i <- rep.row(df_new[i, ], nrow(df_obs)) dist.op <- fields::rdist.vec(locs.test.i, df_obs) @@ -139,16 +139,16 @@ transform_iid <- function(data, vecchia.approx, covparms, nuggets) { # compute transformed data in parts part1.ord <- Matrix::crossprod(U.z, data[U.obj$ord.z, ]) # C.hat^-1 - temp1 <- U.y %*% part1.ord # - revord <- nrow(temp1):1 + temp1 <- U.y %*% part1.ord + revord <- rev(seq_len(nrow(temp1))) temp2 <- spam::solve(V.ord, temp1[revord, ]) part2.rev <- spam::solve(Matrix::t(V.ord), temp2) part2.ord <- spam::crossprod(U.y, part2.rev[revord, ]) transform.ord <- part1.ord - part2.ord # return to original ordering - orig.order <- order(U.obj$ord) - transformed.data <- transform.ord[orig.order, ] + # orig.order <- order(U.obj$ord) + # transformed.data <- transform.ord[orig.order, ] return(transform.ord) } @@ -168,15 +168,15 @@ transform_miid <- function(data, vecchia.approx, params) { # compute transformed data in parts part1.ord <- Matrix::crossprod(U.z, data[U.obj$ord.z, ]) # C.hat^-1 temp1 <- U.y %*% part1.ord # - revord <- nrow(temp1):1 + revord <- rev(seq_len(nrow(temp1))) temp2 <- spam::solve(V.ord, temp1[revord, ]) part2.rev <- spam::solve(Matrix::t(V.ord), temp2) part2.ord <- spam::crossprod(U.y, part2.rev[revord, ]) transform.ord <- part1.ord - part2.ord # return to original ordering - orig.order <- order(U.obj$ord) - transformed.data <- transform.ord[orig.order, ] + # orig.order <- order(U.obj$ord) + # transformed.data <- transform.ord[orig.order, ] return(transform.ord) } @@ -218,11 +218,17 @@ U2V <- function(U.obj) { ########################################################### ## Reverse order of matrix rows,cols -revMat <- function(mat) mat[nrow(mat):1, ncol(mat):1, drop = FALSE] +revMat <- function(mat) { + if (nrow(mat) == 0 || ncol(mat) == 0) { + return(mat) + } + row_seq <- rev(seq_len(nrow(mat))) + col_seq <- rev(seq_len(ncol(mat))) + mat[row_seq, col_seq, drop = FALSE] +} #' @export -vecchia_Mprediction <- function(z, vecchia.approx, covparms, var.exact = NULL, - return.values = "mean") { +vecchia_Mprediction <- function(z, vecchia.approx, covparms, var.exact = NULL, return.values = "mean") { GPvecchia:::removeNAs() U.obj <- createUMultivariate(vecchia.approx, covparms) V.ord <- U2V(U.obj) @@ -233,11 +239,11 @@ vecchia_Mprediction <- function(z, vecchia.approx, covparms, var.exact = NULL, mu.pred = vecchia.mean$mu.pred, mu.obs = vecchia.mean$mu.obs, var.pred = NULL, var.obs = NULL, V.ord = NULL, U.obj = NULL ) - if (return.values == "meanmat" | return.values == "all") { + if (return.values == "meanmat" || return.values == "all") { return.list$V.ord <- V.ord return.list$U.obj <- U.obj } - if (return.values == "meanvar" | return.values == "all") { + if (return.values == "meanvar" || return.values == "all") { if (is.null(var.exact)) { var.exact <- (sum(!vecchia.approx$obs) < 4 * 10000) } diff --git a/R/PrestoGP_Vecchia.R b/R/PrestoGP_Vecchia.R index ebd28ec..34ad020 100755 --- a/R/PrestoGP_Vecchia.R +++ b/R/PrestoGP_Vecchia.R @@ -42,57 +42,70 @@ setMethod("initialize", "VecchiaModel", function(.Object, n_neighbors = 25, ...) #' prediction <- prestogp_predict(model, X.test, locs.test) #' Vec.mean <- prediction[[1]] #' Vec.sds <- prediction[[2]] -setMethod("prestogp_predict", "VecchiaModel", function(model, X, locs, m = NULL, ordering.pred = c("obspred", "general"), pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) { - # validate parameters - ordering.pred <- match.arg(ordering.pred) - pred.cond <- match.arg(pred.cond) - return.values <- match.arg(return.values) - model <- check_input_pred(model, X, locs) - if (is.null(m)) { # m defaults to the value used for training - m <- model@n_neighbors - } - if (m < model@min_m) { - stop(paste("m must be at least ", model@min_m, sep = "")) - } - if (m >= nrow(model@X_train)) { - warning("Conditioning set size m chosen to be >=n. Changing to m=n-1") - m <- nrow(model@X_train) - 1 - } +setMethod("prestogp_predict", "VecchiaModel", + function(model, X, locs, m = NULL, ordering.pred = c("obspred", "general"), pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) { + # validate parameters + ordering.pred <- match.arg(ordering.pred) + pred.cond <- match.arg(pred.cond) + return.values <- match.arg(return.values) + model <- check_input_pred(model, X, locs) + if (is.null(m)) { # m defaults to the value used for training + m <- model@n_neighbors + } + if (m < model@min_m) { + stop(paste("m must be at least ", model@min_m, sep = "")) + } + if (m >= nrow(model@X_train)) { + warning("Conditioning set size m chosen to be >=n. Changing to m=n-1") + m <- nrow(model@X_train) - 1 + } - # Vecchia prediction at new locations - # Vecchia.Pred <- predict(model@Vecchia_SCAD_fit[[1]], X = X, which = model@lambda_1se_idx[[1]]) - Vecchia.Pred <- predict(model@linear_model, newx = X, s = model@linear_model$lambda[model@lambda_1se_idx]) - # Vecchia trend prediction at observed data - # Vecchia.hat <- predict(model@Vecchia_SCAD_fit[[1]], X = model@X_train, which = model@lambda_1se_idx[[1]]) - Vecchia.hat <- predict(model@linear_model, newx = model@X_train, s = model@linear_model$lambda[model@lambda_1se_idx]) + # Vecchia prediction at new locations + # Vecchia.Pred <- predict(model@Vecchia_SCAD_fit[[1]], X = X, which = model@lambda_1se_idx[[1]]) + Vecchia.Pred <- predict(model@linear_model, newx = X, s = model@linear_model$lambda[model@lambda_1se_idx]) + # Vecchia trend prediction at observed data + # Vecchia.hat <- predict(model@Vecchia_SCAD_fit[[1]], X = model@X_train, which = model@lambda_1se_idx[[1]]) + Vecchia.hat <- predict(model@linear_model, newx = model@X_train, s = model@linear_model$lambda[model@lambda_1se_idx]) - # Test set prediction - res <- model@Y_train - Vecchia.hat + # Test set prediction + res <- model@Y_train - Vecchia.hat - locs.train.scaled <- scale_locs(model, model@locs_train)[[1]] - locs.scaled <- scale_locs(model, list(locs))[[1]] - vec.approx.test <- vecchia_specify(locs.train.scaled, m, locs.pred = locs.scaled, ordering.pred=ordering.pred, pred.cond=pred.cond) + locs.train.scaled <- scale_locs(model, model@locs_train)[[1]] + locs.scaled <- scale_locs(model, list(locs))[[1]] + vec.approx.test <- vecchia_specify(locs.train.scaled, m, locs.pred = locs.scaled, ordering.pred = ordering.pred, pred.cond = pred.cond) ## carry out prediction if (!model@apanasovich) { - pred <- vecchia_prediction(res, vec.approx.test, c(model@covparams[1], 1, model@covparams[3]), model@covparams[4], return.values=return.values) + pred <- vecchia_prediction( + res, + vec.approx.test, + c(model@covparams[1], 1, model@covparams[3]), + model@covparams[4], + return.values = return.values + ) } else { - pred <- vecchia_prediction(res, vec.approx.test, c(model@covparams[1], model@covparams[2], model@covparams[3]), model@covparams[4], return.values=return.values) + pred <- vecchia_prediction( + res, + vec.approx.test, + c(model@covparams[1], model@covparams[2], model@covparams[3]), + model@covparams[4], return.values = return.values + ) } - # prediction function can return both mean and sds - # returns a list with elements mu.pred,mu.obs,var.pred,var.obs,V.ord - Vec.mean <- pred$mu.pred + Vecchia.Pred # residual + mean trend - if (return.values == "mean") { + # prediction function can return both mean and sds + # returns a list with elements mu.pred,mu.obs,var.pred,var.obs,V.ord + Vec.mean <- pred$mu.pred + Vecchia.Pred # residual + mean trend + if (return.values == "mean") { return.list <- list(means = Vec.mean) - } else { + } else { warning("Variance estimates do not include model fitting variance and are anticonservative. Use with caution.") - Vec.sds <- sqrt(pred$var.pred + model@covparams[4]) - return.list <- list(means = Vec.mean, sds = vec.sds) - } + Vec.sds <- sqrt(pred$var.pred + model@covparams[4]) # nolint + return.list <- list(means = Vec.mean, sds = Vec.sds) # nolint + } - return(return.list) -}) + return(return.list) + } +) setMethod("check_input", "VecchiaModel", function(model, Y, X, locs) { if (!is.matrix(locs)) { @@ -130,13 +143,13 @@ setMethod("check_input_pred", "VecchiaModel", function(model, X, locs) { stop("X must be a matrix") } if (ncol(locs) != ncol(model@locs_train[[1]])) { - stop("locs must have the same number of columns as locs_train") + stop("locs must have the same number of columns as locs_train") } if (nrow(X) != nrow(locs)) { - stop("X must have the same number of rows as locs") + stop("X must have the same number of rows as locs") } if (ncol(X) != ncol(model@X_train)) { - stop("X and X_train must have the same number of predictors") + stop("X and X_train must have the same number of predictors") } invisible(model) }) diff --git a/R/RcppExports.R b/R/RcppExports.R old mode 100644 new mode 100755 diff --git a/inst/doc/test.R b/inst/doc/test.R index 4c0587d..5992b7d 100644 --- a/inst/doc/test.R +++ b/inst/doc/test.R @@ -1,3 +1,2 @@ ## ----------------------------------------------------------------------------- cat("test", "\n") - diff --git a/man/compute_error-PrestoGPModel-method.Rd b/man/compute_error-PrestoGPModel-method.Rd old mode 100644 new mode 100755 index ca1c37a..ba2f959 --- a/man/compute_error-PrestoGPModel-method.Rd +++ b/man/compute_error-PrestoGPModel-method.Rd @@ -13,6 +13,5 @@ The total error used in the main optimization loop } \description{ -Compute the error (log likelihood using the GP log likelihood -and penalty from the beta coefficients) +Compute the error (log likelihood using the GP log likelihood and penalty from the beta coefficients) } diff --git a/man/estimate_betas-PrestoGPModel-method.Rd b/man/estimate_betas-PrestoGPModel-method.Rd old mode 100644 new mode 100755 index 1ed4099..2bffacf --- a/man/estimate_betas-PrestoGPModel-method.Rd +++ b/man/estimate_betas-PrestoGPModel-method.Rd @@ -4,7 +4,7 @@ \alias{estimate_betas,PrestoGPModel-method} \title{estimate_betas} \usage{ -\S4method{estimate_betas}{PrestoGPModel}(model, parallel) +\S4method{estimate_betas}{PrestoGPModel}(model, parallel, foldid) } \arguments{ \item{model}{the model to estimate coeffients for} diff --git a/man/prestogp_fit-PrestoGPModel-method.Rd b/man/prestogp_fit-PrestoGPModel-method.Rd index 728f2c5..96f695e 100755 --- a/man/prestogp_fit-PrestoGPModel-method.Rd +++ b/man/prestogp_fit-PrestoGPModel-method.Rd @@ -17,8 +17,9 @@ max_iters = 100, verbose = FALSE, optim.method = "Nelder-Mead", - optim.control = list(trace = 0, reltol = 1e-04, maxit = 5000), - parallel = FALSE + optim.control = list(trace = 0, reltol = 0.001, maxit = 5000), + parallel = FALSE, + foldid = NULL ) } \arguments{ @@ -34,21 +35,17 @@ \item{beta.hat}{The initial beta parameters to use (optional).} -\item{tol}{The model is considered converged when error isn't less than -tol*previous_error (optional).} +\item{tol}{The model is considered converged when error isn't less than tol*previous_error (optional).} -\item{verbose}{If TRUE, additional information about model fit -will be printed.} +\item{verbose}{If TRUE, additional information about model fit will be printed.} -\item{m}{The number of neighboring datapoints to condition on in the -likelihood function (optional).} +\item{m}{The number of neighboring datapoints to condition on in the likelihood function (optional).} } \value{ An object containing model parameters for spatiotemporal prediction. } \description{ -This method fits any PrestoGP model given a matrix of locations, a matrix of -#' independent variables, and a matrix of dependent variables. +This method fits any PrestoGP model given a matrix of locations, a matrix of independent variables, and a matrix of dependent variables. } \examples{ diff --git a/man/prestogp_predict-VecchiaModel-method.Rd b/man/prestogp_predict-VecchiaModel-method.Rd old mode 100644 new mode 100755 index cef3331..52941a8 --- a/man/prestogp_predict-VecchiaModel-method.Rd +++ b/man/prestogp_predict-VecchiaModel-method.Rd @@ -4,7 +4,15 @@ \alias{prestogp_predict,VecchiaModel-method} \title{Make predictions using a previously fit spatiotemporal model.} \usage{ -\S4method{prestogp_predict}{VecchiaModel}(model, X, locs, m = NULL) +\S4method{prestogp_predict}{VecchiaModel}( + model, + X = "matrix", + locs = "matrix", + m = NULL, + ordering.pred = c("obspred", "general"), + pred.cond = c("independent", "general"), + return.values = c("mean", "meanvar") +) } \arguments{ \item{model}{A model object returned by prestogp_fit.} diff --git a/man/show_theta-PrestoGPModel-method.Rd b/man/show_theta-PrestoGPModel-method.Rd old mode 100644 new mode 100755 index 8a9e697..0004602 --- a/man/show_theta-PrestoGPModel-method.Rd +++ b/man/show_theta-PrestoGPModel-method.Rd @@ -9,8 +9,7 @@ \arguments{ \item{object}{the PrestoGP model object} -\item{Y_names}{the names of the different outcome variables -(may just be numbers if not provided in training input)} +\item{Y_names}{the names of the different outcome variables (may just be numbers if not provided in training input)} } \description{ Print the covariance parameters in a table diff --git a/man/sparseNN.Rd b/man/sparseNN.Rd old mode 100644 new mode 100755 index 13f7c63..777e84f --- a/man/sparseNN.Rd +++ b/man/sparseNN.Rd @@ -20,7 +20,8 @@ sparseNN( \item{dist_func}{Any distance function with a signature of dist(query_location, locations_matrix)} } \value{ -A list containing two matrices, each with one row per location: an indices matrix with the indices of nearest neighbors for each location, and a distance matrix with the associated distances +A list containing two matrices, each with one row per location: +an indices matrix with the indices of nearest neighbors for each location, and a distance matrix with the associated distances } \description{ Find the index of K nearest neighbors within a set of locations for a given query and distance function diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp old mode 100644 new mode 100755 diff --git a/src/createUMC.cpp b/src/createUMC.cpp index 24b3a7a..213581b 100644 --- a/src/createUMC.cpp +++ b/src/createUMC.cpp @@ -4,59 +4,47 @@ // [[Rcpp::plugins(cpp17)]] #include -double matern1(const double &x, const double &smooth, const double &alpha) -{ - if (smooth == 0.5) - { +double matern1(const double &x, const double &smooth, const double &alpha) { + if (smooth == 0.5) { return exp(-x * alpha); - } - else if (smooth == 1.5) - { + } else if (smooth == 1.5) { double normcon = (pow(alpha, 3) / sqrt(M_PI)) * (tgamma(2.0) / tgamma(1.5)); - return normcon * 0.5 * M_PI * pow(alpha, -3) * exp(-x * alpha) * (1 + x * alpha); - } - else if (smooth == 2.5) - { + return normcon * 0.5 * M_PI * pow(alpha, -3) * exp(-x * alpha) * + (1 + x * alpha); + } else if (smooth == 2.5) { double normcon = (pow(alpha, 5) / sqrt(M_PI)) * (tgamma(3.0) / tgamma(2.5)); double scaled_x = x * alpha; - return normcon * 0.125 * M_PI * pow(alpha, -5) * exp(-scaled_x) * (3 + 3 * scaled_x + scaled_x * scaled_x); - } - else if (x == 0.0) - { + return normcon * 0.125 * M_PI * pow(alpha, -5) * exp(-scaled_x) * + (3 + 3 * scaled_x + scaled_x * scaled_x); + } else if (x == 0.0) { return 1.0; - } - else - { + } else { double normcon = pow(2.0, 1.0 - smooth) / tgamma(smooth); double scaled_x = x * alpha; return normcon * pow(scaled_x, smooth) * R::bessel_k(scaled_x, smooth, 1); } } -// computes Euclidean distance between 2 vectors ([slightly slower] helper function for rdist in C++) -double dist1(const arma::rowvec &a, const arma::rowvec &b) -{ +// computes Euclidean distance between 2 vectors ([slightly slower] helper +// function for rdist in C++) +double dist1(const arma::rowvec &a, const arma::rowvec &b) { double temp = 0.0; // loop over each element and square the distance, add to the placeholder - for (arma::uword i = 0; i < a.n_elem; i++) - { + for (arma::uword i = 0; i < a.n_elem; i++) { temp += (a(i) - b(i)) * (a(i) - b(i)); } return sqrt(temp); } // [[Rcpp::export]] -arma::vec na_omit_c(arma::vec x) -{ +arma::vec na_omit_c(arma::vec x) { // make placeholder vector r arma::vec r(x.size()); // make counter for number of non-NAs int k = 0; - for (unsigned j = 0; j < x.size(); j++) - { + for (unsigned j = 0; j < x.size(); j++) { // if not NA, add 1 to counter and set r(j) = x(j) - if (x(j) == x(j)) - { + if (x(j) == x(j)) { r(j) = x(j); k++; } @@ -70,11 +58,12 @@ arma::vec na_omit_c(arma::vec x) // [[Rcpp::export]] arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, - const arma::mat &curqys, const arma::mat &curqzs, - const arma::mat &vijs, const arma::mat &aijs, - const arma::mat &full_const, const arma::vec &nugget, - const arma::vec &sig2, const arma::vec &U_beginning) -{ + const arma::mat &curqys, + const arma::mat &curqzs, const arma::mat &vijs, + const arma::mat &aijs, + const arma::mat &full_const, + const arma::vec &nugget, const arma::vec &sig2, + const arma::vec &U_beginning) { int n = arma::as_scalar(ondx.n_elem); int m = arma::as_scalar(curqys.n_rows); int n_inds = 2 * n * (m + 3); @@ -85,16 +74,14 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, // vals.rows(0, 6) = U_beginning; // feeder.cols(0, 6) = {{1, 1, 2, 1, 3, 3, 4}, // {1, 2, 2, 3, 3, 4, 4}}; - ndx.cols(0, 6) = {{0, 0, 1, 0, 2, 2, 3}, - {0, 1, 1, 2, 2, 3, 3}}; + ndx.cols(0, 6) = {{0, 0, 1, 0, 2, 2, 3}, {0, 1, 1, 2, 2, 3, 3}}; vals(arma::span(0, 6)) = U_beginning; // feeder.cols(0,6) = {{0, 0, 1, 0, 2, 2, 3}, // {0, 1, 1, 2, 2, 3, 3}, {1,1,1,1,1,1,1}}; // feeder(2, arma::span(0,6)) = U_beginning.t(); int ind = 7; // int ind = 0; - for (arma::uword i = 3; i < (ondx.n_elem + 1); i++) - { + for (arma::uword i = 3; i < (ondx.n_elem + 1); i++) { double temppi = ondx(i - 1); arma::vec cy = na_omit_c(curqys.col(i - 1)); arma::vec cz = (curqzs.col(i - 1)); @@ -104,22 +91,21 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, arma::mat k2(nq, nq); // Rcpp::Rcout << i << " "; // Rcpp::Rcout << cq; - for (arma::uword j = 0; j < nq; j++) - { + for (arma::uword j = 0; j < nq; j++) { double temppj = ondx(cq(j) - 1); - k1(j) = full_const(temppi - 1, temppj - 1) * - matern1(dist1(olocs.row(i - 1), olocs.row(cq(j) - 1)), - vijs(temppi - 1, temppj - 1), aijs(temppi - 1, temppj - 1)); - for (arma::uword k = j; k < nq; k++) - { + k1(j) = + full_const(temppi - 1, temppj - 1) * + matern1(dist1(olocs.row(i - 1), olocs.row(cq(j) - 1)), + vijs(temppi - 1, temppj - 1), aijs(temppi - 1, temppj - 1)); + for (arma::uword k = j; k < nq; k++) { double temppk = ondx(cq(k) - 1); - k2(j, k) = full_const(temppj - 1, temppk - 1) * - matern1(dist1(olocs.row(cq(j) - 1), olocs.row(cq(k) - 1)), - vijs(temppj - 1, temppk - 1), aijs(temppj - 1, temppk - 1)); + k2(j, k) = + full_const(temppj - 1, temppk - 1) * + matern1(dist1(olocs.row(cq(j) - 1), olocs.row(cq(k) - 1)), + vijs(temppj - 1, temppk - 1), aijs(temppj - 1, temppk - 1)); } // // add nugget if neighbor is not latent - if (j >= cy.n_elem) - { + if (j >= cy.n_elem) { k2(j, j) += nugget(temppj - 1); } } @@ -131,16 +117,14 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, arma::vec bi(nq); bool status = arma::solve(bi, k2, k1); double epsilon = 1e-9; - while (!status) - { + while (!status) { epsilon *= 2; k2.diag() += epsilon; status = arma::solve(bi, k2, k1); } double ritemp = arma::as_scalar(sig2(temppi - 1) - bi.t() * k1); epsilon = 1e-9; - while (ritemp <= 0) - { + while (ritemp <= 0) { epsilon *= 2; k2.diag() += epsilon; arma::solve(bi, k2, k1); @@ -154,11 +138,9 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, // uhat(2*i - 2, 2*i - 1) = -1.0 * uhat(2*i - 1, 2*i - 1); // uhat(2*i - 2, 2*i - 2) = 1.0 / ri; - for (arma::uword j = 0; j < cq.n_elem; j++) - { + for (arma::uword j = 0; j < cq.n_elem; j++) { arma::uword xind = 2 * cq(j) - 2; - if (j >= cy.n_elem) - { + if (j >= cy.n_elem) { xind += 1; } // feeder.col(ind) = {(double)xind,(double) 2*i - 2, bi(j)}; @@ -169,7 +151,9 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, // uhat(xind, 2.0*i - 2.0) = bi(j); } // vals[ind] = 1.0 / ri; - double feed_temp = pow(nugget(temppi - 1), -0.5); // used twice so temp storage to avoid recomputation + double feed_temp = + pow(nugget(temppi - 1), + -0.5); // used twice so temp storage to avoid recomputation // feeder.col(ind) = {(double)2*i - 2,(double)2*i - 2, 1.0 / ri}; // feeder.col(ind + 2) = {(double)2*i - 1,(double) 2*i - 1, feed_temp}; ndx.col(ind) = {2 * i - 2, 2 * i - 2}; @@ -178,7 +162,8 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, vals(ind + 2) = feed_temp; // vals[ind + 2] = pow(nugget(temppi - 1), -0.5); // vals[ind + 1] = -1.0 * vals[ind + 2]; - // feeder.col(ind + 1) = {(double)2*i - 2,(double) 2*i - 1, -1.0 * feed_temp}; + // feeder.col(ind + 1) = {(double)2*i - 2,(double) 2*i - 1, -1.0 * + // feed_temp}; ndx.col(ind + 1) = {2 * i - 2, 2 * i - 1}; vals(ind + 1) = -1.0 * feed_temp; ind += 3; diff --git a/tests/testthat/sim_multivariate.R b/tests/testthat/sim_multivariate.R index f33725b..6763784 100644 --- a/tests/testthat/sim_multivariate.R +++ b/tests/testthat/sim_multivariate.R @@ -13,35 +13,35 @@ set.seed(1234) # To keep the simulation small and tractable to quickly run, we # simulate 25 spatial locations and 5 temporal locations -n.spatial=25 -n.spatial.xy = sqrt(n.spatial) -n.temporal = 5 +n.spatial <- 25 +n.spatial.xy <- sqrt(n.spatial) +n.temporal <- 5 # Mean trend, X(s): spatiotemporal data, but the trend is spatial only. # Each outcome has a different combination of true betas # p = number of potential covariates -p=10 +p <- 10 -p.nz=4 -beta1= c(rep(0,p-p.nz),rep(1,p.nz)) -p.nz=2 -beta2= c(rep(1,p.nz),rep(0,p-p.nz)) -p.nz=3 -beta3= c(rep(1,p.nz),rep(0,p-p.nz)) +p.nz <- 4 +beta1 <- c(rep(0, p - p.nz), rep(1, p.nz)) +p.nz <- 2 +beta2 <- c(rep(1, p.nz), rep(0, p - p.nz)) +p.nz <- 3 +beta3 <- c(rep(1, p.nz), rep(0, p - p.nz)) # Combine all of the betas into 1 vector -beta.all <- c(beta1,beta2,beta3) +beta.all <- c(beta1, beta2, beta3) # correlated predictors -Sigma.X=exp(-rdist(sample(1:p))/3) -X=mvrnorm(n.spatial,rep(0,p),Sigma.X) +Sigma.X <- exp(-rdist(sample(1:p)) / 3) +X <- mvrnorm(n.spatial, rep(0, p), Sigma.X) # replicate for space-time for a single-outcomes -X.st <- rbind(X,X,X,X,X) +X.st <- rbind(X, X, X, X, X) # Create the multivariate X-matrix, for multiple outcomes # The rows are by outcome (i), space (j), and then time (k) -X.all <- psych::superMatrix(list(X.st,X.st,X.st)) +X.all <- psych::superMatrix(list(X.st, X.st, X.st)) # S-T design matrix with X's repeated is needed @@ -62,12 +62,12 @@ mean.trend.st <- X.all %*% beta.all # rho is the correlation between outcome S-T errors, 1 for with itself rho.vec <- c(0.8, 0.3, 0.5) -rho <- matrix(0,nrow=3,ncol=3) +rho <- matrix(0, nrow = 3, ncol = 3) rho[upper.tri(rho)] <- rho.vec -rho <- rho+t(rho)+diag(1,3) +rho <- rho + t(rho) + diag(1, 3) # nu, marginal smoothness: -marg.smoothness <- c(0.5,1,1.5) +marg.smoothness <- c(0.5, 1, 1.5) # Non-spatial error/noise/nugget # epsilon(s,t) will be simulated based on the nugget parameters, @@ -75,7 +75,7 @@ nuggets <- c(0.5, 0.7, 2) # marginal variances of the Matern # The marginal variances are scaled linearly with the nuggets -x.variance <- runif(3,1.5,4) +x.variance <- runif(3, 1.5, 4) marg.var <- nuggets + x.variance # ranges, the true ranges, we will scale coordinates by this, but then use a @@ -84,16 +84,16 @@ marg.var <- nuggets + x.variance # outcome1 space and time, outcome 2 space and time, outcome 3 space and time # Spatial domain is on the unit [0,1] # Temporal domain is a year [0,365] -ranges <- c(0.8,0.1,0.5) +ranges <- c(0.8, 0.1, 0.5) # set up the true spatial and temporal dimensions -space.dim1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -space.dim2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -space.dim3 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -time.dim1 <- seq(0,3,length.out = n.temporal)+rnorm(n.temporal,0,0.0001) -time.dim2 <- seq(0,3,length.out = n.temporal)+rnorm(n.temporal,0,0.0001) -time.dim3 <- seq(0,3,length.out = n.temporal)+rnorm(n.temporal,0,0.0001) +space.dim1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +space.dim2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +space.dim3 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +time.dim1 <- seq(0, 3, length.out = n.temporal) + rnorm(n.temporal, 0, 0.0001) +time.dim2 <- seq(0, 3, length.out = n.temporal) + rnorm(n.temporal, 0, 0.0001) +time.dim3 <- seq(0, 3, length.out = n.temporal) + rnorm(n.temporal, 0, 0.0001) # scale the space and time dimensions according to the true covariance ranges @@ -101,62 +101,69 @@ d11 <- fields::rdist(expand.grid(space.dim1, space.dim1, time.dim1)) d22 <- fields::rdist(expand.grid(space.dim2, space.dim2, time.dim2)) d33 <- fields::rdist(expand.grid(space.dim3, space.dim3, time.dim3)) -d12 <- fields::rdist(expand.grid(space.dim1, space.dim1, time.dim1), - expand.grid(space.dim2, space.dim2, time.dim2)) +d12 <- fields::rdist( + expand.grid(space.dim1, space.dim1, time.dim1), + expand.grid(space.dim2, space.dim2, time.dim2) +) d21 <- t(d12) -d13 <- fields::rdist(expand.grid(space.dim1, space.dim1, time.dim1), - expand.grid(space.dim3, space.dim3, time.dim3)) +d13 <- fields::rdist( + expand.grid(space.dim1, space.dim1, time.dim1), + expand.grid(space.dim3, space.dim3, time.dim3) +) d31 <- t(d13) -d23 <- fields::rdist(expand.grid(space.dim2, space.dim2, time.dim2), - expand.grid(space.dim3, space.dim3, time.dim3)) +d23 <- fields::rdist( + expand.grid(space.dim2, space.dim2, time.dim2), + expand.grid(space.dim3, space.dim3, time.dim3) +) d32 <- t(d23) ## Create the correlation matrices -Sigma11 <- marg.var[1] * fields::Matern(d11,range = ranges[1], smoothness = marg.smoothness[1]) -Sigma22 <- marg.var[2] * fields::Matern(d22,range = ranges[2], smoothness = marg.smoothness[2]) -Sigma33 <- marg.var[3] * fields::Matern(d33,range = ranges[3], smoothness = marg.smoothness[3]) +Sigma11 <- marg.var[1] * fields::Matern(d11, range = ranges[1], smoothness = marg.smoothness[1]) +Sigma22 <- marg.var[2] * fields::Matern(d22, range = ranges[2], smoothness = marg.smoothness[2]) +Sigma33 <- marg.var[3] * fields::Matern(d33, range = ranges[3], smoothness = marg.smoothness[3]) vii <- marg.smoothness[1] vjj <- marg.smoothness[2] -vij <- (vii+vjj)/2 -aii <- 1/ranges[1] -ajj <- 1/ranges[2] -aij <- sqrt((aii^2+ajj^2)/2) -Sigma12 <- rho[1,2] * sqrt(marg.var[1]) * sqrt(marg.var[2]) * aii^vii * - ajj^vjj * gamma(vij) / (aij^(2*vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(d12, smoothness=vij, alpha=aij) +vij <- (vii + vjj) / 2 +aii <- 1 / ranges[1] +ajj <- 1 / ranges[2] +aij <- sqrt((aii^2 + ajj^2) / 2) +Sigma12 <- rho[1, 2] * sqrt(marg.var[1]) * sqrt(marg.var[2]) * aii^vii * + ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * + fields::Matern(d12, smoothness = vij, alpha = aij) Sigma21 <- t(Sigma12) vii <- marg.smoothness[1] vjj <- marg.smoothness[3] -vij <- (vii+vjj)/2 -aii <- 1/ranges[1] -ajj <- 1/ranges[3] -aij <- sqrt((aii^2+ajj^2)/2) -Sigma13 <- rho[1,3] * sqrt(marg.var[1]) * sqrt(marg.var[3]) * aii^vii * - ajj^vjj * gamma(vij) / (aij^(2*vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(d13, smoothness=vij, alpha=aij) +vij <- (vii + vjj) / 2 +aii <- 1 / ranges[1] +ajj <- 1 / ranges[3] +aij <- sqrt((aii^2 + ajj^2) / 2) +Sigma13 <- rho[1, 3] * sqrt(marg.var[1]) * sqrt(marg.var[3]) * aii^vii * + ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * + fields::Matern(d13, smoothness = vij, alpha = aij) Sigma31 <- t(Sigma13) vii <- marg.smoothness[2] vjj <- marg.smoothness[3] -vij <- (vii+vjj)/2 -aii <- 1/ranges[2] -ajj <- 1/ranges[3] -aij <- sqrt((aii^2+ajj^2)/2) -Sigma23 <- rho[2,3] * sqrt(marg.var[2]) * sqrt(marg.var[3]) * aii^vii * - ajj^vjj * gamma(vij) / (aij^(2*vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(d23, smoothness=vij, alpha=aij) +vij <- (vii + vjj) / 2 +aii <- 1 / ranges[2] +ajj <- 1 / ranges[3] +aij <- sqrt((aii^2 + ajj^2) / 2) +Sigma23 <- rho[2, 3] * sqrt(marg.var[2]) * sqrt(marg.var[3]) * aii^vii * + ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * + fields::Matern(d23, smoothness = vij, alpha = aij) Sigma32 <- t(Sigma23) -#Combine into the super Multivariate covariance matrix -Sigma.All <- rbind( cbind(Sigma11,Sigma12,Sigma13), - cbind(Sigma21,Sigma22,Sigma23), - cbind(Sigma31,Sigma32,Sigma33) - ) +# Combine into the super Multivariate covariance matrix +Sigma.All <- rbind( + cbind(Sigma11, Sigma12, Sigma13), + cbind(Sigma21, Sigma22, Sigma23), + cbind(Sigma31, Sigma32, Sigma33) +) # Cholesky decomposition L.C <- chol(Sigma.All) diff --git a/tests/testthat/sim_multivariate_big.R b/tests/testthat/sim_multivariate_big.R index ce9b9da..aed9017 100755 --- a/tests/testthat/sim_multivariate_big.R +++ b/tests/testthat/sim_multivariate_big.R @@ -12,9 +12,9 @@ beta.all <- rep(beta1, ny) X.st <- list() for (i in 1:ny) { - set.seed(1212 + i) - Sigma.X <- exp(-rdist(sample(1:p)) / 3) - X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) + set.seed(1212 + i) + Sigma.X <- exp(-rdist(sample(1:p)) / 3) + X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) } X.all <- superMatrix(X.st) mean.trend.st <- X.all %*% beta.all @@ -35,39 +35,33 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() for (i in 1:ny) { - loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) } Sigma.All <- matrix(nrow = ny * n.spatial.xy^2, ncol = ny * n.spatial.xy^2) for (i in 1:ny) { - for (j in i:ny) { - if (i == j) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - dij <- fields::rdist(locs.list[[i]]) - Sigma.All[ndx1, ndx1] <- marg.var[i] * - fields::Matern(dij, - range = ranges[i], - smoothness = marg.smoothness[i] - ) - } else { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) - dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) - vii <- marg.smoothness[i] - vjj <- marg.smoothness[j] - vij <- (vii + vjj) / 2 - aii <- 1 / ranges[i] - ajj <- 1 / ranges[j] - aij <- sqrt((aii^2 + ajj^2) / 2) - Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * - sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness = vij, alpha = aij) - Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) - } + for (j in i:ny) { + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + dij <- fields::rdist(locs.list[[i]]) + Sigma.All[ndx1, ndx1] <- marg.var[i] * fields::Matern(dij, range = ranges[i], smoothness = marg.smoothness[i]) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) + dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) + vii <- marg.smoothness[i] + vjj <- marg.smoothness[j] + vij <- (vii + vjj) / 2 + aii <- 1 / ranges[i] + ajj <- 1 / ranges[j] + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } + } } L.C <- chol(Sigma.All) @@ -76,18 +70,18 @@ st.error <- rnorm(n.spatial.xy^2 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) } y.list <- list() for (i in 1:ny) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] } rm( - ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var + ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var ) diff --git a/tests/testthat/sim_multivariate_big_pred.R b/tests/testthat/sim_multivariate_big_pred.R index 6b6a972..2df4cb6 100755 --- a/tests/testthat/sim_multivariate_big_pred.R +++ b/tests/testthat/sim_multivariate_big_pred.R @@ -14,8 +14,8 @@ beta.all <- rep(beta1, ny) X.st <- list() for (i in 1:ny) { - Sigma.X <- exp(-rdist(sample(1:p)) / 3) - X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) + Sigma.X <- exp(-rdist(sample(1:p)) / 3) + X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) } X.all <- superMatrix(X.st) mean.trend.st <- X.all %*% beta.all @@ -36,39 +36,37 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() for (i in 1:ny) { - loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) } Sigma.All <- matrix(nrow = ny * n.spatial.xy^2, ncol = ny * n.spatial.xy^2) for (i in 1:ny) { - for (j in i:ny) { - if (i == j) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - dij <- fields::rdist(locs.list[[i]]) - Sigma.All[ndx1, ndx1] <- marg.var[i] * - fields::Matern(dij, - range = ranges[i], - smoothness = marg.smoothness[i] - ) - } else { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) - dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) - vii <- marg.smoothness[i] - vjj <- marg.smoothness[j] - vij <- (vii + vjj) / 2 - aii <- 1 / ranges[i] - ajj <- 1 / ranges[j] - aij <- sqrt((aii^2 + ajj^2) / 2) - Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * - sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness = vij, alpha = aij) - Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) - } + for (j in i:ny) { + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + dij <- fields::rdist(locs.list[[i]]) + Sigma.All[ndx1, ndx1] <- marg.var[i] * fields::Matern(dij, + range = ranges[i], + smoothness = marg.smoothness[i] + ) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) + dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) + vii <- marg.smoothness[i] + vjj <- marg.smoothness[j] + vij <- (vii + vjj) / 2 + aii <- 1 / ranges[i] + ajj <- 1 / ranges[j] + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * + fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } + } } L.C <- chol(Sigma.All) @@ -77,7 +75,7 @@ st.error <- rnorm(n.spatial.xy^2 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) } y.list <- list() @@ -88,22 +86,22 @@ y.list.otst <- list() X.st.otr <- list() X.st.otst <- list() for (i in 1:ny) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] - nn <- n.spatial.xy^2 - otr <- rep(FALSE, nn) - otr[sample(1:nn, size = floor(nn / 2))] <- TRUE - locs.list.otr[[i]] <- locs.list[[i]][otr, ] - locs.list.otst[[i]] <- locs.list[[i]][!otr, ] - y.list.otr[[i]] <- y.list[[i]][otr] - y.list.otst[[i]] <- y.list[[i]][!otr] - X.st.otr[[i]] <- X.st[[i]][otr, ] - X.st.otst[[i]] <- X.st[[i]][!otr, ] + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] + nn <- n.spatial.xy^2 + otr <- rep(FALSE, nn) + otr[sample(1:nn, size = floor(nn / 2))] <- TRUE + locs.list.otr[[i]] <- locs.list[[i]][otr, ] + locs.list.otst[[i]] <- locs.list[[i]][!otr, ] + y.list.otr[[i]] <- y.list[[i]][otr] + y.list.otst[[i]] <- y.list[[i]][!otr] + X.st.otr[[i]] <- X.st[[i]][otr, ] + X.st.otst[[i]] <- X.st[[i]][!otr, ] } rm( - ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var, x.variance, nn, otr + ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var, x.variance, nn, otr ) diff --git a/tests/testthat/sim_multivariate_big_st.R b/tests/testthat/sim_multivariate_big_st.R index 2408cb8..b0adc6c 100755 --- a/tests/testthat/sim_multivariate_big_st.R +++ b/tests/testthat/sim_multivariate_big_st.R @@ -14,8 +14,8 @@ beta.all <- rep(beta1, ny) X.st <- list() for (i in 1:ny) { - Sigma.X <- exp(-rdist(sample(1:p)) / 3) - X.st[[i]] <- mvrnorm(n.spatial.xy^3, rep(0, p), Sigma.X) + Sigma.X <- exp(-rdist(sample(1:p)) / 3) + X.st[[i]] <- mvrnorm(n.spatial.xy^3, rep(0, p), Sigma.X) } X.all <- superMatrix(X.st) mean.trend.st <- X.all %*% beta.all @@ -39,44 +39,38 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() locs.listb <- list() for (i in 1:ny) { - loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc3 <- seq(0, 365, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc1b <- loc1 / ranges[2 * i - 1] - loc2b <- loc2 / ranges[2 * i - 1] - loc3b <- loc3 / ranges[2 * i] - locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2, loc3)) - locs.listb[[i]] <- as.matrix(expand.grid(loc1b, loc2b, loc3b)) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc3 <- seq(0, 365, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc1b <- loc1 / ranges[2 * i - 1] + loc2b <- loc2 / ranges[2 * i - 1] + loc3b <- loc3 / ranges[2 * i] + locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2, loc3)) + locs.listb[[i]] <- as.matrix(expand.grid(loc1b, loc2b, loc3b)) } Sigma.All <- matrix(nrow = ny * n.spatial.xy^3, ncol = ny * n.spatial.xy^3) for (i in 1:ny) { - for (j in i:ny) { - if (i == j) { - ndx1 <- ((i - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * i) - dij <- fields::rdist(locs.listb[[i]]) - Sigma.All[ndx1, ndx1] <- marg.var[i] * - fields::Matern(dij, - range = 1, - smoothness = marg.smoothness[i] - ) - } else { - ndx1 <- ((i - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * i) - ndx2 <- ((j - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * j) - dij <- fields::rdist(locs.listb[[i]], locs.listb[[j]]) - vii <- marg.smoothness[i] - vjj <- marg.smoothness[j] - vij <- (vii + vjj) / 2 - aii <- 1 - ajj <- 1 - aij <- sqrt((aii^2 + ajj^2) / 2) - Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * - sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness = vij, alpha = aij) - Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) - } + for (j in i:ny) { + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * i) + dij <- fields::rdist(locs.listb[[i]]) + Sigma.All[ndx1, ndx1] <- marg.var[i] * fields::Matern(dij, range = 1, smoothness = marg.smoothness[i]) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * i) + ndx2 <- ((j - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * j) + dij <- fields::rdist(locs.listb[[i]], locs.listb[[j]]) + vii <- marg.smoothness[i] + vjj <- marg.smoothness[j] + vij <- (vii + vjj) / 2 + aii <- 1 + ajj <- 1 + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } + } } L.C <- chol(Sigma.All) @@ -85,18 +79,18 @@ st.error <- rnorm(n.spatial.xy^3 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^3)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^3)) } y.list <- list() for (i in 1:ny) { - ndx1 <- ((i - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * i) - y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] + ndx1 <- ((i - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * i) + y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] } rm( - ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var, locs.listb, loc3, loc1b, loc2b, loc3b + ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var, locs.listb, loc3, loc1b, loc2b, loc3b ) diff --git a/tests/testthat/sim_multivariate_lik.R b/tests/testthat/sim_multivariate_lik.R index 7bac2f5..d14dd05 100755 --- a/tests/testthat/sim_multivariate_lik.R +++ b/tests/testthat/sim_multivariate_lik.R @@ -23,39 +23,33 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() for (i in 1:ny) { - loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) } Sigma.All <- matrix(nrow = ny * n.spatial.xy^2, ncol = ny * n.spatial.xy^2) for (i in 1:ny) { - for (j in i:ny) { - if (i == j) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - dij <- fields::rdist(locs.list[[i]]) - Sigma.All[ndx1, ndx1] <- marg.var[i] * - fields::Matern(dij, - range = ranges[i], - smoothness = marg.smoothness[i] - ) - } else { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) - dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) - vii <- marg.smoothness[i] - vjj <- marg.smoothness[j] - vij <- (vii + vjj) / 2 - aii <- 1 / ranges[i] - ajj <- 1 / ranges[j] - aij <- sqrt((aii^2 + ajj^2) / 2) - Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * - sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness = vij, alpha = aij) - Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) - } + for (j in i:ny) { + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + dij <- fields::rdist(locs.list[[i]]) + Sigma.All[ndx1, ndx1] <- marg.var[i] * fields::Matern(dij, range = ranges[i], smoothness = marg.smoothness[i]) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) + dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) + vii <- marg.smoothness[i] + vjj <- marg.smoothness[j] + vij <- (vii + vjj) / 2 + aii <- 1 / ranges[i] + ajj <- 1 / ranges[j] + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } + } } L.C <- chol(Sigma.All) @@ -64,18 +58,18 @@ st.error <- rnorm(n.spatial.xy^2 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) } y.list <- list() for (i in 1:ny) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - y.list[[i]] <- st.error[ndx1] + nug.error[ndx1] + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + y.list[[i]] <- st.error[ndx1] + nug.error[ndx1] } rm( - ny, n.spatial.xy, i, j, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var + ny, n.spatial.xy, i, j, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var ) diff --git a/tests/testthat/sim_multivariate_small.R b/tests/testthat/sim_multivariate_small.R index 0fbcd93..aefab32 100755 --- a/tests/testthat/sim_multivariate_small.R +++ b/tests/testthat/sim_multivariate_small.R @@ -12,9 +12,9 @@ beta.all <- rep(beta1, ny) X.st <- list() for (i in 1:ny) { - set.seed(1212 + i) - Sigma.X <- exp(-rdist(sample(1:p)) / 3) - X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) + set.seed(1212 + i) + Sigma.X <- exp(-rdist(sample(1:p)) / 3) + X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) } X.all <- superMatrix(X.st) mean.trend.st <- X.all %*% beta.all @@ -35,39 +35,33 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() for (i in 1:ny) { - loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) } Sigma.All <- matrix(nrow = ny * n.spatial.xy^2, ncol = ny * n.spatial.xy^2) for (i in 1:ny) { - for (j in i:ny) { - if (i == j) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - dij <- fields::rdist(locs.list[[i]]) - Sigma.All[ndx1, ndx1] <- marg.var[i] * - fields::Matern(dij, - range = ranges[i], - smoothness = marg.smoothness[i] - ) - } else { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) - dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) - vii <- marg.smoothness[i] - vjj <- marg.smoothness[j] - vij <- (vii + vjj) / 2 - aii <- 1 / ranges[i] - ajj <- 1 / ranges[j] - aij <- sqrt((aii^2 + ajj^2) / 2) - Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * - sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness = vij, alpha = aij) - Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) - } + for (j in i:ny) { + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + dij <- fields::rdist(locs.list[[i]]) + Sigma.All[ndx1, ndx1] <- marg.var[i] * fields::Matern(dij, range = ranges[i], smoothness = marg.smoothness[i]) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) + dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) + vii <- marg.smoothness[i] + vjj <- marg.smoothness[j] + vij <- (vii + vjj) / 2 + aii <- 1 / ranges[i] + ajj <- 1 / ranges[j] + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } + } } L.C <- chol(Sigma.All) @@ -76,18 +70,18 @@ st.error <- rnorm(n.spatial.xy^2 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) } y.list <- list() for (i in 1:ny) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] } rm( - ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var + ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var ) diff --git a/tests/testthat/sim_multivariate_small_pred.R b/tests/testthat/sim_multivariate_small_pred.R index 0b5101e..961b63b 100755 --- a/tests/testthat/sim_multivariate_small_pred.R +++ b/tests/testthat/sim_multivariate_small_pred.R @@ -9,20 +9,20 @@ library(MASS) library(fields) library(psych) -beta1 <- c(rep(1,p.nz), rep(0,p-p.nz)) +beta1 <- c(rep(1, p.nz), rep(0, p - p.nz)) beta.all <- rep(beta1, ny) X.st <- list() for (i in 1:ny) { - Sigma.X <- exp(-rdist(sample(1:p))/3) - X.st[[i]] <- mvrnorm(n.spatial.xy^2,rep(0,p),Sigma.X) + Sigma.X <- exp(-rdist(sample(1:p)) / 3) + X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) } X.all <- superMatrix(X.st) mean.trend.st <- X.all %*% beta.all n.rho <- choose(ny, 2) rho.vec <- runif(n.rho, 0.2, 0.8) -rho <- matrix(0, nrow=ny, ncol=ny) +rho <- matrix(0, nrow = ny, ncol = ny) rho[upper.tri(rho)] <- rho.vec rho <- rho + t(rho) + diag(1, ny) @@ -36,38 +36,33 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() for (i in 1:ny) { - loc1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) - loc2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) - locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) } -Sigma.All <- matrix(nrow=ny*n.spatial.xy^2, ncol=ny*n.spatial.xy^2) +Sigma.All <- matrix(nrow = ny * n.spatial.xy^2, ncol = ny * n.spatial.xy^2) for (i in 1:ny) { - for (j in i:ny) { - if (i==j) { - ndx1 <- ((i-1)*n.spatial.xy^2+1):(n.spatial.xy^2*i) - dij <- fields::rdist(locs.list[[i]]) - Sigma.All[ndx1,ndx1] <- marg.var[i] * - fields::Matern(dij, range = ranges[i], - smoothness = marg.smoothness[i]) - } - else { - ndx1 <- ((i-1)*n.spatial.xy^2+1):(n.spatial.xy^2*i) - ndx2 <- ((j-1)*n.spatial.xy^2+1):(n.spatial.xy^2*j) - dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) - vii <- marg.smoothness[i] - vjj <- marg.smoothness[j] - vij <- (vii+vjj)/2 - aii <- 1/ranges[i] - ajj <- 1/ranges[j] - aij <- sqrt((aii^2+ajj^2)/2) - Sigma.All[ndx1,ndx2] <- rho[i,j] * sqrt(marg.var[i]) * - sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2*vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness=vij, alpha=aij) - Sigma.All[ndx2,ndx1] <- t(Sigma.All[ndx1,ndx2]) - } + for (j in i:ny) { + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + dij <- fields::rdist(locs.list[[i]]) + Sigma.All[ndx1, ndx1] <- marg.var[i] * fields::Matern(dij, range = ranges[i], smoothness = marg.smoothness[i]) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) + dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) + vii <- marg.smoothness[i] + vjj <- marg.smoothness[j] + vij <- (vii + vjj) / 2 + aii <- 1 / ranges[i] + ajj <- 1 / ranges[j] + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } + } } L.C <- chol(Sigma.All) @@ -76,7 +71,7 @@ st.error <- rnorm(n.spatial.xy^2 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i]*rnorm(n.spatial.xy^2)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) } y.list <- list() @@ -87,20 +82,22 @@ y.list.otst <- list() X.st.otr <- list() X.st.otst <- list() for (i in 1:ny) { - ndx1 <- ((i-1)*n.spatial.xy^2+1):(n.spatial.xy^2*i) - y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] - nn <- n.spatial.xy^2 - otr <- rep(FALSE, nn) - otr[sample(1:nn, size=floor(nn/2))] <- TRUE - locs.list.otr[[i]] <- locs.list[[i]][otr,] - locs.list.otst[[i]] <- locs.list[[i]][!otr,] - y.list.otr[[i]] <- y.list[[i]][otr] - y.list.otst[[i]] <- y.list[[i]][!otr] - X.st.otr[[i]] <- X.st[[i]][otr,] - X.st.otst[[i]] <- X.st[[i]][!otr,] + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] + nn <- n.spatial.xy^2 + otr <- rep(FALSE, nn) + otr[sample(1:nn, size = floor(nn / 2))] <- TRUE + locs.list.otr[[i]] <- locs.list[[i]][otr, ] + locs.list.otst[[i]] <- locs.list[[i]][!otr, ] + y.list.otr[[i]] <- y.list[[i]][otr] + y.list.otst[[i]] <- y.list[[i]][!otr] + X.st.otr[[i]] <- X.st[[i]][otr, ] + X.st.otst[[i]] <- X.st[[i]][!otr, ] } -rm(ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var, x.variance, nn, otr) +rm( + ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var, x.variance, nn, otr +) diff --git a/tests/testthat/sim_vecchia.R b/tests/testthat/sim_vecchia.R index 992e3f0..ecfb02c 100755 --- a/tests/testthat/sim_vecchia.R +++ b/tests/testthat/sim_vecchia.R @@ -7,10 +7,10 @@ n.spatial.xy <- 50 # number of spatial coordinates per dimension library(MASS) library(fields) -beta1 <- c(rep(1,p.nz), rep(0,p-p.nz)) +beta1 <- c(rep(1, p.nz), rep(0, p - p.nz)) -Sigma.X <- exp(-rdist(sample(1:p))/3) -X <- mvrnorm(n.spatial.xy^2,rep(0,p),Sigma.X) +Sigma.X <- exp(-rdist(sample(1:p)) / 3) +X <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) mean.trend.st <- as.vector(X %*% beta1) marg.smoothness <- 0.5 + rnorm(1, 0.1, 0.05) @@ -22,18 +22,22 @@ ranges <- runif(1, 0.5, 1.2) params.all <- c(x.variance, ranges, marg.smoothness, nuggets) -loc1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -loc2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) +loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2)) -Sigma.All <- marg.var*Matern(rdist(locs), range=ranges, - smoothness=marg.smoothness) +Sigma.All <- marg.var * Matern(rdist(locs), + range = ranges, + smoothness = marg.smoothness +) L.C <- chol(Sigma.All) st.error <- as.vector(rnorm(n.spatial.xy^2) %*% L.C) -nug.error <- nuggets*rnorm(n.spatial.xy^2) +nug.error <- nuggets * rnorm(n.spatial.xy^2) y <- mean.trend.st + st.error + nug.error -rm(p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, - st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance) +rm( + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, + st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, + x.variance +) diff --git a/tests/testthat/sim_vecchia_pred.R b/tests/testthat/sim_vecchia_pred.R index ffcb97e..566ff72 100755 --- a/tests/testthat/sim_vecchia_pred.R +++ b/tests/testthat/sim_vecchia_pred.R @@ -7,10 +7,10 @@ n.spatial.xy <- 50 # number of spatial coordinates per dimension library(MASS) library(fields) -beta1 <- c(rep(1,p.nz), rep(0,p-p.nz)) +beta1 <- c(rep(1, p.nz), rep(0, p - p.nz)) -Sigma.X <- exp(-rdist(sample(1:p))/3) -X <- mvrnorm(n.spatial.xy^2,rep(0,p),Sigma.X) +Sigma.X <- exp(-rdist(sample(1:p)) / 3) +X <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) mean.trend.st <- as.vector(X %*% beta1) marg.smoothness <- 0.5 + rnorm(1, 0.1, 0.05) @@ -22,29 +22,33 @@ ranges <- runif(1, 0.5, 1.2) params.all <- c(x.variance, ranges, marg.smoothness, nuggets) -loc1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -loc2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) +loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2)) -Sigma.All <- marg.var*Matern(rdist(locs), range=ranges, - smoothness=marg.smoothness) +Sigma.All <- marg.var * Matern(rdist(locs), + range = ranges, + smoothness = marg.smoothness +) L.C <- chol(Sigma.All) st.error <- as.vector(rnorm(n.spatial.xy^2) %*% L.C) -nug.error <- nuggets*rnorm(n.spatial.xy^2) +nug.error <- nuggets * rnorm(n.spatial.xy^2) y <- mean.trend.st + st.error + nug.error nn <- n.spatial.xy^2 otr <- rep(FALSE, nn) -otr[sample(1:nn, size=floor(nn/2))] <- TRUE +otr[sample(1:nn, size = floor(nn / 2))] <- TRUE -locs.otr <- locs[otr,] -locs.otst <- locs[!otr,] +locs.otr <- locs[otr, ] +locs.otst <- locs[!otr, ] y.otr <- y[otr] y.otst <- y[!otr] -X.otr <- X[otr,] -X.otst <- X[!otr,] - -rm(p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, - st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance, nn, otr) +X.otr <- X[otr, ] +X.otst <- X[!otr, ] + +rm( + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, + st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, + x.variance, nn, otr +) diff --git a/tests/testthat/sim_vecchia_small.R b/tests/testthat/sim_vecchia_small.R index 9c884eb..1bb9325 100755 --- a/tests/testthat/sim_vecchia_small.R +++ b/tests/testthat/sim_vecchia_small.R @@ -7,10 +7,10 @@ n.spatial.xy <- 10 # number of spatial coordinates per dimension library(MASS) library(fields) -beta1 <- c(rep(1,p.nz), rep(0,p-p.nz)) +beta1 <- c(rep(1, p.nz), rep(0, p - p.nz)) -Sigma.X <- exp(-rdist(sample(1:p))/3) -X <- mvrnorm(n.spatial.xy^2,rep(0,p),Sigma.X) +Sigma.X <- exp(-rdist(sample(1:p)) / 3) +X <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) mean.trend.st <- as.vector(X %*% beta1) marg.smoothness <- 0.5 + rnorm(1, 0.1, 0.05) @@ -22,18 +22,22 @@ ranges <- runif(1, 0.5, 1.2) params.all <- c(x.variance, ranges, marg.smoothness, nuggets) -loc1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -loc2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) +loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2)) -Sigma.All <- marg.var*Matern(rdist(locs), range=ranges, - smoothness=marg.smoothness) +Sigma.All <- marg.var * Matern(rdist(locs), + range = ranges, + smoothness = marg.smoothness +) L.C <- chol(Sigma.All) st.error <- as.vector(rnorm(n.spatial.xy^2) %*% L.C) -nug.error <- nuggets*rnorm(n.spatial.xy^2) +nug.error <- nuggets * rnorm(n.spatial.xy^2) y <- mean.trend.st + st.error + nug.error -rm(p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, - st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance) +rm( + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, + st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, + x.variance +) diff --git a/tests/testthat/sim_vecchia_small_pred.R b/tests/testthat/sim_vecchia_small_pred.R index 7f25507..c082605 100755 --- a/tests/testthat/sim_vecchia_small_pred.R +++ b/tests/testthat/sim_vecchia_small_pred.R @@ -7,10 +7,10 @@ n.spatial.xy <- 20 # number of spatial coordinates per dimension library(MASS) library(fields) -beta1 <- c(rep(1,p.nz), rep(0,p-p.nz)) +beta1 <- c(rep(1, p.nz), rep(0, p - p.nz)) -Sigma.X <- exp(-rdist(sample(1:p))/3) -X <- mvrnorm(n.spatial.xy^2,rep(0,p),Sigma.X) +Sigma.X <- exp(-rdist(sample(1:p)) / 3) +X <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) mean.trend.st <- as.vector(X %*% beta1) marg.smoothness <- 0.5 + rnorm(1, 0.1, 0.05) @@ -22,29 +22,33 @@ ranges <- runif(1, 0.5, 1.2) params.all <- c(x.variance, ranges, marg.smoothness, nuggets) -loc1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -loc2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) +loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2)) -Sigma.All <- marg.var*Matern(rdist(locs), range=ranges, - smoothness=marg.smoothness) +Sigma.All <- marg.var * Matern(rdist(locs), + range = ranges, + smoothness = marg.smoothness +) L.C <- chol(Sigma.All) st.error <- as.vector(rnorm(n.spatial.xy^2) %*% L.C) -nug.error <- nuggets*rnorm(n.spatial.xy^2) +nug.error <- nuggets * rnorm(n.spatial.xy^2) y <- mean.trend.st + st.error + nug.error nn <- n.spatial.xy^2 otr <- rep(FALSE, nn) -otr[sample(1:nn, size=floor(nn/2))] <- TRUE +otr[sample(1:nn, size = floor(nn / 2))] <- TRUE -locs.otr <- locs[otr,] -locs.otst <- locs[!otr,] +locs.otr <- locs[otr, ] +locs.otst <- locs[!otr, ] y.otr <- y[otr] y.otst <- y[!otr] -X.otr <- X[otr,] -X.otst <- X[!otr,] - -rm(p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, - st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance, nn, otr) +X.otr <- X[otr, ] +X.otst <- X[!otr, ] + +rm( + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, + st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, + x.variance, nn, otr +) diff --git a/tests/testthat/sim_vecchia_st.R b/tests/testthat/sim_vecchia_st.R index 9438ee3..401bbee 100755 --- a/tests/testthat/sim_vecchia_st.R +++ b/tests/testthat/sim_vecchia_st.R @@ -7,10 +7,10 @@ n.spatial.xy <- 15 # number of spatial coordinates per dimension library(MASS) library(fields) -beta1 <- c(rep(1,p.nz), rep(0,p-p.nz)) +beta1 <- c(rep(1, p.nz), rep(0, p - p.nz)) -Sigma.X <- exp(-rdist(sample(1:p))/3) -X <- mvrnorm(n.spatial.xy^3,rep(0,p),Sigma.X) +Sigma.X <- exp(-rdist(sample(1:p)) / 3) +X <- mvrnorm(n.spatial.xy^3, rep(0, p), Sigma.X) mean.trend.st <- as.vector(X %*% beta1) marg.smoothness <- 0.5 + rnorm(1, 0.1, 0.05) @@ -22,23 +22,27 @@ ranges <- runif(2, 0.5, 1.2) params.all <- c(x.variance, ranges, marg.smoothness, nuggets) -loc1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -loc2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -loc3 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) +loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +loc3 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2, loc3)) locss <- locs -locss[,1:2] <- locss[,1:2] / ranges[1] -locss[,3] <- locs[,3] / ranges[2] +locss[, 1:2] <- locss[, 1:2] / ranges[1] +locss[, 3] <- locs[, 3] / ranges[2] -Sigma.All <- marg.var*Matern(rdist(locss), range=1, - smoothness=marg.smoothness) +Sigma.All <- marg.var * Matern(rdist(locss), + range = 1, + smoothness = marg.smoothness +) L.C <- chol(Sigma.All) st.error <- as.vector(rnorm(n.spatial.xy^3) %*% L.C) -nug.error <- nuggets*rnorm(n.spatial.xy^3) +nug.error <- nuggets * rnorm(n.spatial.xy^3) y <- mean.trend.st + st.error + nug.error -rm(p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, loc3, L.C, - st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance, locss) +rm( + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, loc3, L.C, + st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, + x.variance, locss +) diff --git a/tests/testthat/test-Log_Likelihood.R b/tests/testthat/test-Log_Likelihood.R index d1ddb03..9326acb 100755 --- a/tests/testthat/test-Log_Likelihood.R +++ b/tests/testthat/test-Log_Likelihood.R @@ -1,359 +1,366 @@ test_that("negloglik_vecchia", { - set.seed(1234) - - y.orig <- rnorm(100) - - locs.dim1 <- seq(1, 10, length = 10) + rnorm(10, 0, 0.1) - locs.dim2 <- seq(1, 10, length = 10) + rnorm(10, 0, 0.1) - locs <- as.matrix(expand.grid(locs.dim1, locs.dim2)) - - covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(100) - - y <- y.orig %*% chol(covmat.true) - y <- as.vector(y) - - logparams <- create.initial.values.flex(0.9, 0.5, 0.5, 0.1, 1, 1) - pseq <- create.param.sequence(1) - vec.approx <- vecchia_specify(locs, 5) - LL.vecchia <- negloglik_vecchia(logparams, y, vec.approx, pseq) - expect_equal(194.75, LL.vecchia, tolerance=1e-2) - - vec.mapprox <- vecchia_Mspecify(list(locs), 5) - LL.mvecchia <- mvnegloglik(logparams, vec.mapprox, y, pseq, 1) - # Univariate likelihood should equal multivariate likelihood - expect_equal(LL.vecchia, LL.mvecchia, tolerance=1e-1) - - logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) - vec.approx2 <- vec.approx - vec.approx2$locsord <- vec.approx$locsord / 0.5 - LL.vecchia2 <- negloglik_vecchia(logparams2, y, vec.approx2, pseq) - LL.vecchia.st <- negloglik_vecchia_ST(logparams, y, vec.approx, pseq, 1, 1) - vec.mapprox2 <- vec.mapprox - vec.mapprox2$locsord <- vec.mapprox$locsord / 0.5 - LL.mvecchia2 <- mvnegloglik_ST(logparams2, vec.mapprox2, y, pseq, 1, 1, 1) - - # Likelihood should equal spatiotemporal likelihood after scaling locs - expect_equal(LL.vecchia, LL.vecchia2, tolerance=1e-3) - expect_equal(LL.vecchia, LL.vecchia.st, tolerance=1e-3) - expect_equal(LL.vecchia, LL.mvecchia2, tolerance=1e-1) + set.seed(1234) + + y.orig <- rnorm(100) + locs.dim1 <- seq(1, 10, length = 10) + rnorm(10, 0, 0.1) + locs.dim2 <- seq(1, 10, length = 10) + rnorm(10, 0, 0.1) + locs <- as.matrix(expand.grid(locs.dim1, locs.dim2)) + covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(100) + + y <- y.orig %*% chol(covmat.true) + y <- as.vector(y) + + logparams <- create.initial.values.flex(0.9, 0.5, 0.5, 0.1, 1, 1) + pseq <- create.param.sequence(1) + vec.approx <- vecchia_specify(locs, 5) + LL.vecchia <- negloglik_vecchia(logparams, y, vec.approx, pseq) + expect_equal(194.75, LL.vecchia, tolerance = 1e-2) + + vec.mapprox <- vecchia_Mspecify(list(locs), 5) + LL.mvecchia <- mvnegloglik(logparams, vec.mapprox, y, pseq, 1) + # Univariate likelihood should equal multivariate likelihood + expect_equal(LL.vecchia, LL.mvecchia, tolerance = 1e-1) + + logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) + vec.approx2 <- vec.approx + vec.approx2$locsord <- vec.approx$locsord / 0.5 + LL.vecchia2 <- negloglik_vecchia(logparams2, y, vec.approx2, pseq) + LL.vecchia.st <- negloglik_vecchia_ST(logparams, y, vec.approx, pseq, 1, 1) + vec.mapprox2 <- vec.mapprox + vec.mapprox2$locsord <- vec.mapprox$locsord / 0.5 + LL.mvecchia2 <- mvnegloglik_ST(logparams2, vec.mapprox2, y, pseq, 1, 1, 1) + + # Likelihood should equal spatiotemporal likelihood after scaling locs + expect_equal(LL.vecchia, LL.vecchia2, tolerance = 1e-3) + expect_equal(LL.vecchia, LL.vecchia.st, tolerance = 1e-3) + expect_equal(LL.vecchia, LL.mvecchia2, tolerance = 1e-1) }) test_that("negloglik_vecchia_ST", { - set.seed(1234) - - y.orig <- rnorm(125) - - locs.dim1 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) - locs.dim2 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) - locs.dim3 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) - locs <- as.matrix(expand.grid(locs.dim1, locs.dim2, locs.dim3)) - - covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(125) - - y <- y.orig %*% chol(covmat.true) - y <- as.vector(y) - - logparams <- create.initial.values.flex(0.9, c(2, 3), 0.5, 0.1, 1, 1) - pseq <- create.param.sequence(1, 2) - vec.approx <- vecchia_specify(locs, 5) - LL.vecchia <- negloglik_vecchia_ST(logparams, y, vec.approx, pseq, - c(1, 1, 2), 2) - expect_equal(284.73, LL.vecchia, tolerance=1e-2) - - logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) - pseq2 <- create.param.sequence(1) - vec.approx2 <- vec.approx - vec.approx2$locsord[,1:2] <- vec.approx$locsord[,1:2] / 2 - vec.approx2$locsord[,3] <- vec.approx$locsord[,3] / 3 - LL.vecchia2 <- negloglik_vecchia(logparams2, y, vec.approx2, pseq2) - expect_equal(LL.vecchia, LL.vecchia2, tolerance=1e-3) + set.seed(1234) + + y.orig <- rnorm(125) + + locs.dim1 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs.dim2 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs.dim3 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs <- as.matrix(expand.grid(locs.dim1, locs.dim2, locs.dim3)) + + covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(125) + + y <- y.orig %*% chol(covmat.true) + y <- as.vector(y) + + logparams <- create.initial.values.flex(0.9, c(2, 3), 0.5, 0.1, 1, 1) + pseq <- create.param.sequence(1, 2) + vec.approx <- vecchia_specify(locs, 5) + LL.vecchia <- negloglik_vecchia_ST( + logparams, y, vec.approx, pseq, + c(1, 1, 2), 2 + ) + expect_equal(284.73, LL.vecchia, tolerance = 1e-2) + + logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) + pseq2 <- create.param.sequence(1) + vec.approx2 <- vec.approx + vec.approx2$locsord[, 1:2] <- vec.approx$locsord[, 1:2] / 2 + vec.approx2$locsord[, 3] <- vec.approx$locsord[, 3] / 3 + LL.vecchia2 <- negloglik_vecchia(logparams2, y, vec.approx2, pseq2) + expect_equal(LL.vecchia, LL.vecchia2, tolerance = 1e-3) }) test_that("negloglik.full", { - set.seed(1234) + set.seed(1234) - y.orig <- rnorm(100) + y.orig <- rnorm(100) - locs.dim <- seq(1, 10, length = sqrt(length(y.orig))) - locs <- as.matrix(expand.grid(locs.dim, locs.dim)) + locs.dim <- seq(1, 10, length = sqrt(length(y.orig))) + locs <- as.matrix(expand.grid(locs.dim, locs.dim)) - covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(100) + covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(100) - y <- y.orig %*% chol(covmat.true) - y <- as.vector(y) + y <- y.orig %*% chol(covmat.true) + y <- as.vector(y) - params.init <- rep(NA, 4) - params.init[1] <- 0.9 * var(y) - params.init[2] <- 1 - params.init[3] <- 0.5 - params.init[4] <- 0.1 * var(y) + params.init <- rep(NA, 4) + params.init[1] <- 0.9 * var(y) + params.init[2] <- 1 + params.init[3] <- 0.5 + params.init[4] <- 0.1 * var(y) - params.init <- create.initial.values.flex( - params.init[1], - params.init[2], - params.init[3], - params.init[4], - 1, 1 - ) + params.init <- create.initial.values.flex( + params.init[1], + params.init[2], + params.init[3], + params.init[4], + 1, 1 + ) - d <- rdist(locs) - pseq <- create.param.sequence(1) + d <- rdist(locs) + pseq <- create.param.sequence(1) - res.optim.NM <- optim(par=params.init, fn=negloglik.full, d=d, y=y, - param.seq=pseq, control=list(maxit=5000)) + res.optim.NM <- optim( + par = params.init, fn = negloglik.full, d = d, y = y, + param.seq = pseq, control = list(maxit = 5000) + ) - LL.full <- negloglik.full(res.optim.NM$par, d, y, pseq) + LL.full <- negloglik.full(res.optim.NM$par, d, y, pseq) - params.final <- c(exp(res.optim.NM$par[1:2]), - gtools::inv.logit(res.optim.NM$par[3], 0, 2.5), - exp(res.optim.NM$par[4])) + params.final <- c( + exp(res.optim.NM$par[1:2]), + gtools::inv.logit(res.optim.NM$par[3], 0, 2.5), + exp(res.optim.NM$par[4]) + ) - pgp.params <- create.initial.values.flex(params.final[1], - params.final[2], - params.final[3], - params.final[4], - 1, 1) + pgp.params <- create.initial.values.flex( + params.final[1], + params.final[2], + params.final[3], + params.final[4], + 1, 1 + ) - LL.full.pgp <- mvnegloglik.full(pgp.params, list(locs), y, pseq) + LL.full.pgp <- mvnegloglik.full(pgp.params, list(locs), y, pseq) - vec.approx <- vecchia_specify(locs, 99) + vec.approx <- vecchia_specify(locs, 99) - LL.vecchia <- -1 * vecchia_likelihood( - y, vec.approx, params.final[1:3], - params.final[4] - ) + LL.vecchia <- -1 * vecchia_likelihood( + y, vec.approx, params.final[1:3], + params.final[4] + ) - vec.approx.pgp <- vecchia_Mspecify(list(locs), 99) - vec.U.pgp <- createUMultivariate(vec.approx.pgp, c(params.final, 1)) + vec.approx.pgp <- vecchia_Mspecify(list(locs), 99) + vec.U.pgp <- createUMultivariate(vec.approx.pgp, c(params.final, 1)) - LL.pgp <- -1 * GPvecchia:::vecchia_likelihood_U(y, vec.U.pgp) + LL.pgp <- -1 * GPvecchia:::vecchia_likelihood_U(y, vec.U.pgp) - expect_equal(173.315, LL.full, tolerance = 1e-3) - # Univariate likelihood should equal the multivariate likelihood - expect_equal(LL.full, LL.full.pgp, tolerance = 1e-3) - # Full likelihood should equal both the univariate and multivariate - # Vecchia approximations - expect_equal(LL.full, LL.vecchia, tolerance = 1e-3) - expect_equal(LL.full, LL.pgp, tolerance = 1e-3) + expect_equal(173.315, LL.full, tolerance = 1e-3) + # Univariate likelihood should equal the multivariate likelihood + expect_equal(LL.full, LL.full.pgp, tolerance = 1e-3) + # Full likelihood should equal both the univariate and multivariate + # Vecchia approximations + expect_equal(LL.full, LL.vecchia, tolerance = 1e-3) + expect_equal(LL.full, LL.pgp, tolerance = 1e-3) }) test_that("negloglik_full_ST", { - set.seed(1234) - - y.orig <- rnorm(125) - - locs.dim1 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) - locs.dim2 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) - locs.dim3 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) - locs <- as.matrix(expand.grid(locs.dim1, locs.dim2, locs.dim3)) - - covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(125) - - y <- y.orig %*% chol(covmat.true) - y <- as.vector(y) - - logparams <- create.initial.values.flex(0.9, c(2, 3), 0.5, 0.1, 1, 1) - pseq <- create.param.sequence(1, 2) - LL.full <- negloglik_full_ST(logparams, locs, y, pseq, - c(1, 1, 2), 2) - expect_equal(286.84, LL.full, tolerance=1e-2) - - logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) - pseq2 <- create.param.sequence(1) - locs2 <- locs - locs2[,1:2] <- locs[,1:2] / 2 - locs2[,3] <- locs[,3] / 3 - LL.full2 <- negloglik.full(logparams2, rdist(locs2), y, pseq2) - expect_equal(LL.full, LL.full2, tolerance=1e-3) + set.seed(1234) + + y.orig <- rnorm(125) + + locs.dim1 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs.dim2 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs.dim3 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs <- as.matrix(expand.grid(locs.dim1, locs.dim2, locs.dim3)) + + covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(125) + + y <- y.orig %*% chol(covmat.true) + y <- as.vector(y) + + logparams <- create.initial.values.flex(0.9, c(2, 3), 0.5, 0.1, 1, 1) + pseq <- create.param.sequence(1, 2) + LL.full <- negloglik_full_ST( + logparams, locs, y, pseq, + c(1, 1, 2), 2 + ) + expect_equal(286.84, LL.full, tolerance = 1e-2) + + logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) + pseq2 <- create.param.sequence(1) + locs2 <- locs + locs2[, 1:2] <- locs[, 1:2] / 2 + locs2[, 3] <- locs[, 3] / 3 + LL.full2 <- negloglik.full(logparams2, rdist(locs2), y, pseq2) + expect_equal(LL.full, LL.full2, tolerance = 1e-3) }) test_that("mvnegloglik", { - load("sim_multivariate_big.RData") - set.seed(1234) - P <- 3 - Y <- cbind(runif(10), runif(10), runif(10)) - cor.matrix <- cor(Y) - cov_mat <- c(cor.matrix[upper.tri(cor.matrix)]) - logparams <- create.initial.values.flex( - rep(0.9, P), # marginal variance - rep(0.5, P), # range - rep(0.5, P), # smoothness - rep(0.1, P), # nuggets - cov_mat, - P - ) - pseq <- create.param.sequence(P) - vec.approx <- vecchia_Mspecify(locs.list, 25) - neg_likelihood <- mvnegloglik(logparams, vec.approx, - unlist(y.list), pseq, P) - expect_equal(34384.23, neg_likelihood, tolerance=1e-2) + load("sim_multivariate_big.RData") + set.seed(1234) + P <- 3 + Y <- cbind(runif(10), runif(10), runif(10)) + cor.matrix <- cor(Y) + cov_mat <- c(cor.matrix[upper.tri(cor.matrix)]) + logparams <- create.initial.values.flex( + rep(0.9, P), # marginal variance + rep(0.5, P), # range + rep(0.5, P), # smoothness + rep(0.1, P), # nuggets + cov_mat, + P + ) + pseq <- create.param.sequence(P) + vec.approx <- vecchia_Mspecify(locs.list, 25) + neg_likelihood <- mvnegloglik( + logparams, vec.approx, + unlist(y.list), pseq, P + ) + expect_equal(34384.23, neg_likelihood, tolerance = 1e-2) }) test_that("mvnegloglik_ST", { - source("sim_multivariate_big_st.R") - P <- 3 - Y <- cbind(runif(10), runif(10), runif(10)) - cor.matrix <- cor(Y) - cov_mat <- c(cor.matrix[upper.tri(cor.matrix)]) - logparams <- create.initial.values.flex( - rep(0.9, P), # marginal variance - rep(c(2, 3), P), # range - rep(0.5, P), # smoothness - rep(0.1, P), # nuggets - cov_mat, - P - ) - pseq <- create.param.sequence(P, 2) - vec.approx <- vecchia_Mspecify(locs.list, 25) - neg_likelihood <- mvnegloglik_ST( - logparams, vec.approx, - unlist(y.list), pseq, P, c(1, 1, 2), 2 - ) - expect_equal(35106.73, neg_likelihood, tolerance = 1e-2) - - vec.approx2 <- vec.approx - for (i in 1:P) { - vec.approx2$locsord[, 1:2] <- vec.approx$locsord[, 1:2] / 2 - vec.approx2$locsord[, 3] <- vec.approx$locsord[, 3] / 3 - } - logparams2 <- logparams - logparams2[pseq[2, 1]:pseq[2, 2]] <- 0 - neg_likelihood2 <- mvnegloglik_ST( - logparams2, vec.approx2, - unlist(y.list), pseq, P, c(1, 1, 2), 2 - ) - expect_equal(neg_likelihood, neg_likelihood2, tolerance = 1e-3) - - logparams <- create.initial.values.flex( - rep(0.9, P), # marginal variance - 2:7, # range - rep(0.5, P), # smoothness - rep(0.1, P), # nuggets - cov_mat, - P - ) - neg_likelihood <- mvnegloglik_ST( - logparams, vec.approx, - unlist(y.list), pseq, P, c(1, 1, 2), 2 - ) - expect_equal(36354.9, neg_likelihood, tolerance = 1e-2) - - vec.approx2 <- vec.approx - vec.approx2$locsord[vec.approx$ondx == 1, 1:2] <- - vec.approx$locsord[vec.approx$ondx == 1, 1:2] / 2 - vec.approx2$locsord[vec.approx$ondx == 1, 3] <- - vec.approx$locsord[vec.approx$ondx == 1, 3] / 3 - vec.approx2$locsord[vec.approx$ondx == 2, 1:2] <- - vec.approx$locsord[vec.approx$ondx == 2, 1:2] / 4 - vec.approx2$locsord[vec.approx$ondx == 2, 3] <- - vec.approx$locsord[vec.approx$ondx == 2, 3] / 5 - vec.approx2$locsord[vec.approx$ondx == 3, 1:2] <- - vec.approx$locsord[vec.approx$ondx == 3, 1:2] / 6 - vec.approx2$locsord[vec.approx$ondx == 3, 3] <- - vec.approx$locsord[vec.approx$ondx == 3, 3] / 7 - - neg_likelihood2 <- mvnegloglik_ST( - logparams2, vec.approx2, - unlist(y.list), pseq, P, c(1, 1, 2), 2 - ) - expect_equal(neg_likelihood, neg_likelihood2, tolerance = 1e-3) + source("sim_multivariate_big_st.R") + P <- 3 + Y <- cbind(runif(10), runif(10), runif(10)) + cor.matrix <- cor(Y) + cov_mat <- c(cor.matrix[upper.tri(cor.matrix)]) + logparams <- create.initial.values.flex( + rep(0.9, P), # marginal variance + rep(c(2, 3), P), # range + rep(0.5, P), # smoothness + rep(0.1, P), # nuggets + cov_mat, + P + ) + pseq <- create.param.sequence(P, 2) + vec.approx <- vecchia_Mspecify(locs.list, 25) + neg_likelihood <- mvnegloglik_ST( + logparams, vec.approx, + unlist(y.list), pseq, P, c(1, 1, 2), 2 + ) + expect_equal(35106.73, neg_likelihood, tolerance = 1e-2) + + vec.approx2 <- vec.approx + for (i in 1:P) { + vec.approx2$locsord[, 1:2] <- vec.approx$locsord[, 1:2] / 2 + vec.approx2$locsord[, 3] <- vec.approx$locsord[, 3] / 3 + } + logparams2 <- logparams + logparams2[pseq[2, 1]:pseq[2, 2]] <- 0 + neg_likelihood2 <- mvnegloglik_ST( + logparams2, vec.approx2, + unlist(y.list), pseq, P, c(1, 1, 2), 2 + ) + expect_equal(neg_likelihood, neg_likelihood2, tolerance = 1e-3) + + logparams <- create.initial.values.flex( + rep(0.9, P), # marginal variance + 2:7, # range + rep(0.5, P), # smoothness + rep(0.1, P), # nuggets + cov_mat, + P + ) + neg_likelihood <- mvnegloglik_ST( + logparams, vec.approx, + unlist(y.list), pseq, P, c(1, 1, 2), 2 + ) + expect_equal(36354.9, neg_likelihood, tolerance = 1e-2) + + vec.approx2 <- vec.approx + vec.approx2$locsord[vec.approx$ondx == 1, 1:2] <- + vec.approx$locsord[vec.approx$ondx == 1, 1:2] / 2 + vec.approx2$locsord[vec.approx$ondx == 1, 3] <- + vec.approx$locsord[vec.approx$ondx == 1, 3] / 3 + vec.approx2$locsord[vec.approx$ondx == 2, 1:2] <- + vec.approx$locsord[vec.approx$ondx == 2, 1:2] / 4 + vec.approx2$locsord[vec.approx$ondx == 2, 3] <- + vec.approx$locsord[vec.approx$ondx == 2, 3] / 5 + vec.approx2$locsord[vec.approx$ondx == 3, 1:2] <- + vec.approx$locsord[vec.approx$ondx == 3, 1:2] / 6 + vec.approx2$locsord[vec.approx$ondx == 3, 3] <- + vec.approx$locsord[vec.approx$ondx == 3, 3] / 7 + + neg_likelihood2 <- mvnegloglik_ST( + logparams2, vec.approx2, + unlist(y.list), pseq, P, c(1, 1, 2), 2 + ) + expect_equal(neg_likelihood, neg_likelihood2, tolerance = 1e-3) }) test_that("mvnegloglik.full", { - source("sim_multivariate_lik.R") - - pseq <- create.param.sequence(3) - - param.marg.var <- 0.9 * unlist(lapply(y.list, var)) - param.marg.scale <- rep(1, 3) - param.marg.smooth <- rep(0.5, 3) - param.marg.nugget <- 0.1 * unlist(lapply(y.list, var)) - param.rho <- rep(0.5, 3) - params.init <- create.initial.values.flex( - param.marg.var, - param.marg.scale, - param.marg.smooth, - param.marg.nugget, - param.rho, 3 - ) - - res.optim.NM <- optim( - par = params.init, fn = mvnegloglik.full, - locs = locs.list, y = unlist(y.list), - param.seq = pseq, - method = "Nelder-Mead", - control = list(trace = 0, maxit = 10000, reltol = 1e-4) - ) - - LL.full.mv <- mvnegloglik.full( - res.optim.NM$par, locs.list, - unlist(y.list), pseq - ) - - param.seq.begin <- pseq[, 1] - param.seq.end <- pseq[, 2] - params.init.final <- res.optim.NM$par - params.init.final.t <- c( - exp(params.init.final[param.seq.begin[1]: - param.seq.end[2]]), - gtools::inv.logit(params.init.final[ - param.seq.begin[3]: - param.seq.end[3] - ], 0, 2.5), - exp(params.init.final[param.seq.begin[4]: - param.seq.end[4]]), - tanh(params.init.final[ - param.seq.begin[5]:param.seq.end[5] - ]) - ) - - params.final.flex.test <- create.initial.values.flex( - params.init.final.t[param.seq.begin[1]:param.seq.end[1]], - params.init.final.t[param.seq.begin[2]:param.seq.end[2]], - params.init.final.t[param.seq.begin[3]:param.seq.end[3]], - params.init.final.t[param.seq.begin[4]:param.seq.end[4]], - params.init.final.t[param.seq.begin[5]:param.seq.end[5]], 3 - ) - - cov.list <- create.cov.upper.flex( - 3, - params.init.final.t[param.seq.begin[1]:param.seq.end[1]], - params.init.final.t[param.seq.begin[2]:param.seq.end[2]], - params.init.final.t[param.seq.begin[3]:param.seq.end[3]], - params.init.final.t[param.seq.begin[4]:param.seq.end[4]], - params.init.final.t[param.seq.begin[5]:param.seq.end[5]] - ) - - cov.mat <- cat.covariances( - locs.list, cov.list$variance, - cov.list$range, - cov.list$smoothness, cov.list$nugget - ) - - LL.full.calc <- -1 * mvtnorm::dmvnorm(unlist(y.list), - rep(0, length(unlist(y.list))), - cov.mat, - log = TRUE - ) - - vec.mapprox <- vecchia_Mspecify(locs.list, length(unlist(y.list)) - 1) - U.mobj <- createUMultivariate(vec.mapprox, params.init.final.t) - - LL.vecchia.mv <- -1 * GPvecchia:::vecchia_likelihood_U(unlist(y.list), U.mobj) - - expect_equal(541.31, LL.full.mv, tolerance = 1e-3) - expect_equal(LL.full.calc, LL.full.mv, tolerance = 1e-3) - # Full likelihood should equal the Vecchia likelihood - expect_equal(LL.full.mv, LL.vecchia.mv, tolerance = 1e-3) - # Check create.initial.values.flex: - expect_equal(params.init.final, params.final.flex.test, tolerance = 1e-3) - # Check create.cov.upper.flex: - cov.list.true <- readRDS("covlist.rds") - expect_equal(cov.list$variance, cov.list.true$variance, tolerance = 1e-3) - expect_equal(cov.list$range, cov.list.true$range, tolerance = 1e-3) - expect_equal(cov.list$smoothness, cov.list.true$smoothness, tolerance = 1e-3) - expect_equal(cov.list$nugget, cov.list.true$nugget, tolerance = 1e-3) - # Check cat.covariances: - cov.mat.true <- readRDS("covmat.rds") - expect_equal(cov.mat, cov.mat.true, tolerance = 1e-3) + source("sim_multivariate_lik.R") + + pseq <- create.param.sequence(3) + + param.marg.var <- 0.9 * unlist(lapply(y.list, var)) + param.marg.scale <- rep(1, 3) + param.marg.smooth <- rep(0.5, 3) + param.marg.nugget <- 0.1 * unlist(lapply(y.list, var)) + param.rho <- rep(0.5, 3) + params.init <- create.initial.values.flex( + param.marg.var, + param.marg.scale, + param.marg.smooth, + param.marg.nugget, + param.rho, 3 + ) + + res.optim.NM <- optim( + par = params.init, fn = mvnegloglik.full, + locs = locs.list, y = unlist(y.list), + param.seq = pseq, + method = "Nelder-Mead", + control = list(trace = 0, maxit = 10000, reltol = 1e-4) + ) + + LL.full.mv <- mvnegloglik.full( + res.optim.NM$par, locs.list, + unlist(y.list), pseq + ) + + param.seq.begin <- pseq[, 1] + param.seq.end <- pseq[, 2] + params.init.final <- res.optim.NM$par + params.init.final.t <- c( + exp(params.init.final[param.seq.begin[1]:param.seq.end[2]]), + gtools::inv.logit(params.init.final[ + param.seq.begin[3]:param.seq.end[3] + ], 0, 2.5), + exp(params.init.final[param.seq.begin[4]:param.seq.end[4]]), + tanh(params.init.final[ + param.seq.begin[5]:param.seq.end[5] + ]) + ) + + params.final.flex.test <- create.initial.values.flex( + params.init.final.t[param.seq.begin[1]:param.seq.end[1]], + params.init.final.t[param.seq.begin[2]:param.seq.end[2]], + params.init.final.t[param.seq.begin[3]:param.seq.end[3]], + params.init.final.t[param.seq.begin[4]:param.seq.end[4]], + params.init.final.t[param.seq.begin[5]:param.seq.end[5]], 3 + ) + + cov.list <- create.cov.upper.flex( + 3, + params.init.final.t[param.seq.begin[1]:param.seq.end[1]], + params.init.final.t[param.seq.begin[2]:param.seq.end[2]], + params.init.final.t[param.seq.begin[3]:param.seq.end[3]], + params.init.final.t[param.seq.begin[4]:param.seq.end[4]], + params.init.final.t[param.seq.begin[5]:param.seq.end[5]] + ) + + cov.mat <- cat.covariances( + locs.list, cov.list$variance, + cov.list$range, + cov.list$smoothness, cov.list$nugget + ) + + LL.full.calc <- -1 * mvtnorm::dmvnorm(unlist(y.list), + rep(0, length(unlist(y.list))), + cov.mat, + log = TRUE + ) + + vec.mapprox <- vecchia_Mspecify(locs.list, length(unlist(y.list)) - 1) + U.mobj <- createUMultivariate(vec.mapprox, params.init.final.t) + + LL.vecchia.mv <- -1 * GPvecchia:::vecchia_likelihood_U(unlist(y.list), U.mobj) + + expect_equal(541.31, LL.full.mv, tolerance = 1e-3) + expect_equal(LL.full.calc, LL.full.mv, tolerance = 1e-3) + # Full likelihood should equal the Vecchia likelihood + expect_equal(LL.full.mv, LL.vecchia.mv, tolerance = 1e-3) + # Check create.initial.values.flex: + expect_equal(params.init.final, params.final.flex.test, tolerance = 1e-3) + # Check create.cov.upper.flex: + cov.list.true <- readRDS("covlist.rds") + expect_equal(cov.list$variance, cov.list.true$variance, tolerance = 1e-3) + expect_equal(cov.list$range, cov.list.true$range, tolerance = 1e-3) + expect_equal(cov.list$smoothness, cov.list.true$smoothness, tolerance = 1e-3) + expect_equal(cov.list$nugget, cov.list.true$nugget, tolerance = 1e-3) + # Check cat.covariances: + cov.mat.true <- readRDS("covmat.rds") + expect_equal(cov.mat, cov.mat.true, tolerance = 1e-3) }) diff --git a/tests/testthat/test-PrestoGP_CreateU_Multivariate.R b/tests/testthat/test-PrestoGP_CreateU_Multivariate.R index cd0fb65..9c98a09 100755 --- a/tests/testthat/test-PrestoGP_CreateU_Multivariate.R +++ b/tests/testthat/test-PrestoGP_CreateU_Multivariate.R @@ -1,144 +1,148 @@ test_that("create.param.sequence", { - seq <- create.param.sequence(1) - colnames(seq) <- NULL - expect_equal(2, ncol(seq)) - expect_equal(5, nrow(seq)) - expect_equal(c(1, 1), seq[1, ]) - expect_equal(c(2, 2), seq[2, ]) - expect_equal(c(3, 3), seq[3, ]) - expect_equal(c(4, 4), seq[4, ]) - expect_equal(c(5, 5), seq[5, ]) - - seq <- create.param.sequence(3) - colnames(seq) <- NULL - expect_equal(2, ncol(seq)) - expect_equal(5, nrow(seq)) - expect_equal(c(1, 3), seq[1, ]) - expect_equal(c(4, 6), seq[2, ]) - expect_equal(c(7, 9), seq[3, ]) - expect_equal(c(10, 12), seq[4, ]) - expect_equal(c(13, 15), seq[5, ]) - - seq <- create.param.sequence(3, 2) - colnames(seq) <- NULL - expect_equal(2, ncol(seq)) - expect_equal(5, nrow(seq)) - expect_equal(c(1, 3), seq[1, ]) - expect_equal(c(4, 9), seq[2, ]) - expect_equal(c(10, 12), seq[3, ]) - expect_equal(c(13, 15), seq[4, ]) - expect_equal(c(16, 18), seq[5, ]) + seq <- create.param.sequence(1) + colnames(seq) <- NULL + expect_equal(2, ncol(seq)) + expect_equal(5, nrow(seq)) + expect_equal(c(1, 1), seq[1, ]) + expect_equal(c(2, 2), seq[2, ]) + expect_equal(c(3, 3), seq[3, ]) + expect_equal(c(4, 4), seq[4, ]) + expect_equal(c(5, 5), seq[5, ]) + + seq <- create.param.sequence(3) + colnames(seq) <- NULL + expect_equal(2, ncol(seq)) + expect_equal(5, nrow(seq)) + expect_equal(c(1, 3), seq[1, ]) + expect_equal(c(4, 6), seq[2, ]) + expect_equal(c(7, 9), seq[3, ]) + expect_equal(c(10, 12), seq[4, ]) + expect_equal(c(13, 15), seq[5, ]) + + seq <- create.param.sequence(3, 2) + colnames(seq) <- NULL + expect_equal(2, ncol(seq)) + expect_equal(5, nrow(seq)) + expect_equal(c(1, 3), seq[1, ]) + expect_equal(c(4, 9), seq[2, ]) + expect_equal(c(10, 12), seq[3, ]) + expect_equal(c(13, 15), seq[4, ]) + expect_equal(c(16, 18), seq[5, ]) }) test_that("max_min_ordering", { - set.seed(7919) - load("multivariate_sim_spatial3.Rdata") - order <- max_min_ordering(locs_train, fields::rdist) - exp_order <- c( - 63, 21, 50, 30, 76, 78, 40, 36, 23, 69, 9, 67, 32, 2, 20, 62, - 22, 31, 74, 39, 35, 58, 68, 54, 41, 3, 80, 46, 6, 88, 12, 47, - 72, 42, 13, 83, 25, 52, 11, 60, 24, 1, 28, 84, 29, 64, 66, 81, - 82, 55, 61, 87, 17, 33, 43, 45, 10, 79, 53, 75, 89, 51, 73, 27, - 26, 77, 44, 38, 65, 16, 19, 37, 57, 70, 15, 4, 5, 86, 14, 49, - 85, 34, 48, 59, 18, 8, 71, 7, 56, 90 - ) - order.gpv <- order_maxmin_exact(locs_train) - expect_equal(exp_order, order) - expect_equal(order, order.gpv) + set.seed(7919) + load("multivariate_sim_spatial3.Rdata") + order <- max_min_ordering(locs_train, fields::rdist) + exp_order <- c( + 63, 21, 50, 30, 76, 78, 40, 36, 23, 69, 9, 67, 32, 2, 20, 62, + 22, 31, 74, 39, 35, 58, 68, 54, 41, 3, 80, 46, 6, 88, 12, 47, + 72, 42, 13, 83, 25, 52, 11, 60, 24, 1, 28, 84, 29, 64, 66, 81, + 82, 55, 61, 87, 17, 33, 43, 45, 10, 79, 53, 75, 89, 51, 73, 27, + 26, 77, 44, 38, 65, 16, 19, 37, 57, 70, 15, 4, 5, 86, 14, 49, + 85, 34, 48, 59, 18, 8, 71, 7, 56, 90 + ) + order.gpv <- order_maxmin_exact(locs_train) + expect_equal(exp_order, order) + expect_equal(order, order.gpv) }) test_that("sparseNN", { - set.seed(1212) - - locs <- matrix(nrow = 100, ncol = 2) - locs[1, ] <- rep(0, 2) - for (i in 2:nrow(locs)) { - cur.r <- rnorm(1, 5) - cur.t <- runif(1, 0, 2 * pi) - locs[i, ] <- locs[i - 1, ] + c(cur.r * cos(cur.t), cur.r * sin(cur.t)) + set.seed(1212) + + locs <- matrix(nrow = 100, ncol = 2) + locs[1, ] <- rep(0, 2) + for (i in 2:nrow(locs)) { + cur.r <- rnorm(1, 5) + cur.t <- runif(1, 0, 2 * pi) + locs[i, ] <- locs[i - 1, ] + c(cur.r * cos(cur.t), cur.r * sin(cur.t)) + } + + mm.order <- order_maxmin_exact(locs) + olocs <- locs[mm.order, ] + pgp.nn <- sparseNN(olocs, 5, fields::rdist, "rdist") + gpv.nn <- GpGp:::find_ordered_nn(olocs, 5) + + indices <- matrix(nrow = nrow(olocs), ncol = 5) + distances <- indices + for (i in seq_len(nrow(olocs))) { + if (i <= 5) { + cur.dist <- fields::rdist( + olocs[(1:(5 + 1)), ][-i, ], + olocs[i, , drop = FALSE] + ) + indices[i, ] <- order(cur.dist) + } else { + cur.dist <- fields::rdist(olocs[(1:(i - 1)), ], olocs[i, , drop = FALSE]) + indices[i, ] <- order(cur.dist)[1:5] } - - mm.order <- order_maxmin_exact(locs) - olocs <- locs[mm.order, ] - pgp.nn <- sparseNN(olocs, 5, fields::rdist, "rdist") - gpv.nn <- GpGp:::find_ordered_nn(olocs, 5) - - indices <- matrix(nrow = nrow(olocs), ncol = 5) - distances <- indices - for (i in 1:nrow(olocs)) { - if (i <= 5) { - cur.dist <- fields::rdist( - olocs[(1:(5 + 1)), ][-i, ], - olocs[i, , drop = FALSE] - ) - indices[i, ] <- order(cur.dist) - } else { - cur.dist <- fields::rdist(olocs[(1:(i - 1)), ], olocs[i, , drop = FALSE]) - indices[i, ] <- order(cur.dist)[1:5] - } - distances[i, ] <- cur.dist[indices[i, ]] - } - - # Should produce the same nearest neighbors as GPvecchia - expect_equal(pgp.nn$indices[-(1:5), ], gpv.nn[-(1:5), -1]) - # Should obtain the same nearest neighbors and distances when we calculate - # the neighbors by brute force. - expect_equal(pgp.nn$indices[-(1:5), ], indices[-(1:5), ]) - expect_equal(pgp.nn$distances[-(1:5), ], distances[-(1:5), ], tolerance = 1e-2) + distances[i, ] <- cur.dist[indices[i, ]] + } + + # Should produce the same nearest neighbors as GPvecchia + expect_equal(pgp.nn$indices[-(1:5), ], gpv.nn[-(1:5), -1]) + # Should obtain the same nearest neighbors and distances when we calculate + # the neighbors by brute force. + expect_equal(pgp.nn$indices[-(1:5), ], indices[-(1:5), ]) + expect_equal(pgp.nn$distances[-(1:5), ], distances[-(1:5), ], tolerance = 1e-2) }) test_that("createUMultivariate", { - set.seed(1212) - - dim1 <- (0:9)^2 + rnorm(10, 0, 1e-2) - dim2 <- (1:10)^2 + rnorm(10, 0, 1e-2) - - locs <- as.matrix(expand.grid(dim1, dim2)) - - params <- c(3,1.5,0.6,2) - - vec.approx <- vecchia_specify(locs, m=5) - U.obj <- createU(vec.approx, covparms=params[1:3], nuggets=params[4]) - - vec.mapprox <- vecchia_Mspecify(list(locs), 5) - U.mobj <- createUMultivariate(vec.mapprox, c(params,1)) -# Should produce the same ordered locs as GPvecchia in the univariate case - expect_equal(vec.approx$locsord, vec.mapprox$locsord) - expect_equal(vec.approx$ord, vec.mapprox$ord) -# Should produce the same U matrix as GPvecchia in the univariate case - expect_equal(sum(abs(U.obj$U-U.mobj$U)), 0, tolerance=1e-3) - -# Should identify the same nearest neighbors and latents as GPvecchia in the -# univariate case - NNarray <- vec.approx$U.prep$revNNarray[,-6] - cond <- vec.approx$U.prep$revCond[,-6] - - for (i in 6:100) { - cur.y <- NNarray[i,cond[i,]] - expect_equal(sort(cur.y), sort(vec.mapprox$q.list$q.y[[i]])) - if (sum(!cond[i,])>0) { - cur.z <- NNarray[i,!cond[i,]] - expect_equal(sort(cur.z), sort(vec.mapprox$q.list$q.z[[i]])) - } - } + set.seed(1212) + + dim1 <- (0:9)^2 + rnorm(10, 0, 1e-2) + dim2 <- (1:10)^2 + rnorm(10, 0, 1e-2) + + locs <- as.matrix(expand.grid(dim1, dim2)) - source("sim_multivariate.R") - - vec.mapprox <- vecchia_Mspecify(locs.list, 25) - U.mobj <- createUMultivariate(vec.mapprox, c(marg.var, ranges, - marg.smoothness, - nuggets, rho.vec)) - Umat.tst <- readRDS("Umat.rds") - expect_equal(sum(abs(U.mobj$U-Umat.tst)), 0, tolerance=1e-4) - - vec.mapprox.full <- vecchia_Mspecify(locs.list, 374) - U.mobj.full <- createUMultivariate(vec.mapprox.full, c(marg.var, ranges, - marg.smoothness, - nuggets, rho.vec)) - - Sigma.hat <- solve(U.mobj.full$U %*% t(U.mobj.full$U)) - Sigma.hat <- Sigma.hat[U.mobj.full$latent,U.mobj.full$latent] - Sigma.hat <- Sigma.hat[order(U.mobj.full$ord),order(U.mobj.full$ord)] -# Vecchia approximation of covariance should equal true covariance when m=n-1 - expect_equal(sum(abs(Sigma.All-Sigma.hat)), 0, tolerance=1e-4) + params <- c(3, 1.5, 0.6, 2) + + vec.approx <- vecchia_specify(locs, m = 5) + U.obj <- createU(vec.approx, covparms = params[1:3], nuggets = params[4]) + + vec.mapprox <- vecchia_Mspecify(list(locs), 5) + U.mobj <- createUMultivariate(vec.mapprox, c(params, 1)) + # Should produce the same ordered locs as GPvecchia in the univariate case + expect_equal(vec.approx$locsord, vec.mapprox$locsord) + expect_equal(vec.approx$ord, vec.mapprox$ord) + # Should produce the same U matrix as GPvecchia in the univariate case + expect_equal(sum(abs(U.obj$U - U.mobj$U)), 0, tolerance = 1e-3) + + # Should identify the same nearest neighbors and latents as GPvecchia in the + # univariate case + NNarray <- vec.approx$U.prep$revNNarray[, -6] + cond <- vec.approx$U.prep$revCond[, -6] + + for (i in 6:100) { + cur.y <- NNarray[i, cond[i, ]] + expect_equal(sort(cur.y), sort(vec.mapprox$q.list$q.y[[i]])) + if (sum(!cond[i, ]) > 0) { + cur.z <- NNarray[i, !cond[i, ]] + expect_equal(sort(cur.z), sort(vec.mapprox$q.list$q.z[[i]])) + } + } + + source("sim_multivariate.R") + + vec.mapprox <- vecchia_Mspecify(locs.list, 25) + U.mobj <- createUMultivariate(vec.mapprox, c( + marg.var, ranges, + marg.smoothness, + nuggets, rho.vec + )) + Umat.tst <- readRDS("Umat.rds") + expect_equal(sum(abs(U.mobj$U - Umat.tst)), 0, tolerance = 1e-4) + + vec.mapprox.full <- vecchia_Mspecify(locs.list, 374) + U.mobj.full <- createUMultivariate(vec.mapprox.full, c( + marg.var, ranges, + marg.smoothness, + nuggets, rho.vec + )) + + Sigma.hat <- solve(U.mobj.full$U %*% t(U.mobj.full$U)) + Sigma.hat <- Sigma.hat[U.mobj.full$latent, U.mobj.full$latent] + Sigma.hat <- Sigma.hat[order(U.mobj.full$ord), order(U.mobj.full$ord)] + # Vecchia approximation of covariance should equal true covariance when m=n-1 + expect_equal(sum(abs(Sigma.All - Sigma.hat)), 0, tolerance = 1e-4) }) diff --git a/tests/testthat/test-PrestoGP_Model.R b/tests/testthat/test-PrestoGP_Model.R index a13cc6c..a9f7a16 100755 --- a/tests/testthat/test-PrestoGP_Model.R +++ b/tests/testthat/test-PrestoGP_Model.R @@ -1,155 +1,209 @@ test_that("beta.hat not numeric", { source("sim_vecchia_small.R") pgp.model1 <- new("VecchiaModel") - expect_error(prestogp_fit(pgp.model1, y, X, locs, beta.hat="foo"), - "beta.hat parameter must be a numeric vector") + expect_error( + prestogp_fit(pgp.model1, y, X, locs, beta.hat = "foo"), + "beta.hat parameter must be a numeric vector" + ) }) test_that("beta.hat not a vector", { source("sim_vecchia_small.R") pgp.model1 <- new("VecchiaModel") - expect_error(prestogp_fit(pgp.model1, y, X, locs, beta.hat=matrix(1:3)), - "beta.hat parameter must be a numeric vector") + expect_error( + prestogp_fit(pgp.model1, y, X, locs, beta.hat = matrix(1:3)), + "beta.hat parameter must be a numeric vector" + ) }) test_that("beta.hat incorrect length", { source("sim_vecchia_small.R") pgp.model1 <- new("VecchiaModel") - expect_error(prestogp_fit(pgp.model1, y, X, locs, beta.hat=1:3), - "Length of beta.hat must match the number of predictors") + expect_error( + prestogp_fit(pgp.model1, y, X, locs, beta.hat = 1:3), + "Length of beta.hat must match the number of predictors" + ) }) test_that("tol not numeric", { source("sim_vecchia_small.R") pgp.model1 <- new("VecchiaModel") - expect_error(prestogp_fit(pgp.model1, y, X, locs, tol="foo"), - "tol must be numeric") + expect_error( + prestogp_fit(pgp.model1, y, X, locs, tol = "foo"), + "tol must be numeric" + ) }) test_that("tol not a scalar", { source("sim_vecchia_small.R") pgp.model1 <- new("VecchiaModel") - expect_error(prestogp_fit(pgp.model1, y, X, locs, tol=1:2), - "tol must be a scalar") + expect_error( + prestogp_fit(pgp.model1, y, X, locs, tol = 1:2), + "tol must be a scalar" + ) }) test_that("tol out of range", { source("sim_vecchia_small.R") pgp.model1 <- new("VecchiaModel") - expect_error(prestogp_fit(pgp.model1, y, X, locs, tol=1.1), - "tol must satisfy 0=n. Changing to m=n-1") + pgp.model1 <- new("VecchiaModel", n_neighbors = 101) + expect_warning( + prestogp_fit(pgp.model1, y, X, locs), + "Conditioning set size m chosen to be >=n. Changing to m=n-1" + ) }) diff --git a/tests/testthat/test-PrestoGP_Multivariate_Vecchia_Model.R b/tests/testthat/test-PrestoGP_Multivariate_Vecchia_Model.R index 7a2a58a..52f8e38 100755 --- a/tests/testthat/test-PrestoGP_Multivariate_Vecchia_Model.R +++ b/tests/testthat/test-PrestoGP_Multivariate_Vecchia_Model.R @@ -1,306 +1,435 @@ test_that("Invalid locs input", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(as.matrix(1:3)), list(as.matrix(1:3)), - as.matrix(1:3)), - "locs must be a list for multivariate models") + expect_error( + prestogp_fit( + model, list(as.matrix(1:3)), list(as.matrix(1:3)), + as.matrix(1:3) + ), + "locs must be a list for multivariate models" + ) }) test_that("Invalid Y input", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, as.matrix(1:3), list(as.matrix(1:3)), - list(as.matrix(1:3))), - "Y must be a list for multivariate models") + expect_error( + prestogp_fit( + model, as.matrix(1:3), list(as.matrix(1:3)), + list(as.matrix(1:3)) + ), + "Y must be a list for multivariate models" + ) }) test_that("Invalid X input", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(as.matrix(1:3)), as.matrix(1:3), - list(as.matrix(1:3))), - "X must be a list for multivariate models") + expect_error( + prestogp_fit( + model, list(as.matrix(1:3)), as.matrix(1:3), + list(as.matrix(1:3)) + ), + "X must be a list for multivariate models" + ) }) test_that("locs/Y length mismatch", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(1:3, 2:4), list(as.matrix(1:3)), - list(as.matrix(1:3))), - "locs and Y must have the same length") + expect_error( + prestogp_fit( + model, list(1:3, 2:4), list(as.matrix(1:3)), + list(as.matrix(1:3)) + ), + "locs and Y must have the same length" + ) }) test_that("locs/X length mismatch", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(1:3), list(as.matrix(1:3), - as.matrix(2:4)), - list(as.matrix(1:3))), - "locs and X must have the same length") + expect_error( + prestogp_fit( + model, list(1:3), list( + as.matrix(1:3), + as.matrix(2:4) + ), + list(as.matrix(1:3)) + ), + "locs and X must have the same length" + ) }) test_that("locs not a matrix", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(1:3), list(as.matrix(1:3)), - list(1:3)), - "Each locs must be a matrix") + expect_error( + prestogp_fit( + model, list(1:3), list(as.matrix(1:3)), + list(1:3) + ), + "Each locs must be a matrix" + ) }) test_that("X not a matrix", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(1:3), list(1:3), - list(as.matrix(1:3))), - "Each X must be a matrix") + expect_error( + prestogp_fit( + model, list(1:3), list(1:3), + list(as.matrix(1:3)) + ), + "Each X must be a matrix" + ) }) test_that("Y not a numeric matrix/vector", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list("foo"), list(as.matrix(1:3)), - list(as.matrix(1:3))), - "Each Y must be a numeric vector or matrix") + expect_error( + prestogp_fit( + model, list("foo"), list(as.matrix(1:3)), + list(as.matrix(1:3)) + ), + "Each Y must be a numeric vector or matrix" + ) }) test_that("locs with differing numbers of columns", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(1:4, 1:4), - list(as.matrix(1:4), as.matrix(1:4)), - list(as.matrix(1:4), matrix(1:4, nrow=2))), - "All locs must have the same number of columns") + expect_error( + prestogp_fit( + model, list(1:4, 1:4), + list(as.matrix(1:4), as.matrix(1:4)), + list(as.matrix(1:4), matrix(1:4, nrow = 2)) + ), + "All locs must have the same number of columns" + ) }) test_that("Y has multiple columns", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(matrix(1:4, nrow=2)), - list(as.matrix(1:4)), list(as.matrix(1:4))), - "Each Y must have only 1 column") + expect_error( + prestogp_fit( + model, list(matrix(1:4, nrow = 2)), + list(as.matrix(1:4)), list(as.matrix(1:4)) + ), + "Each Y must have only 1 column" + ) }) test_that("nrow(Y) != nrow(locs)", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(as.matrix(1:4)), - list(as.matrix(1:4)), list(as.matrix(1:3))), - "Each Y must have the same number of rows as locs") + expect_error( + prestogp_fit( + model, list(as.matrix(1:4)), + list(as.matrix(1:4)), list(as.matrix(1:3)) + ), + "Each Y must have the same number of rows as locs" + ) }) test_that("nrow(Y) != nrow(X)", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(as.matrix(1:4)), - list(as.matrix(1:3)), list(as.matrix(1:4))), - "Each Y must have the same number of rows as X") + expect_error( + prestogp_fit( + model, list(as.matrix(1:4)), + list(as.matrix(1:3)), list(as.matrix(1:4)) + ), + "Each Y must have the same number of rows as X" + ) }) test_that("Simulated dataset multivariate spatial", { - load("sim_multivariate_big.RData") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=25) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list, X.st, locs.list, - scaling = c(1, 1), apanasovich = TRUE, verbose = FALSE, - optim.control = list( - trace = 0, maxit = 5000, - reltol = 1e-3 - ) + load("sim_multivariate_big.RData") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 25) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list, X.st, locs.list, + scaling = c(1, 1), apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 ) - beta.out <- as.vector(pgp.mmodel1@beta) - params.out <- pgp.mmodel1@covparams + ) + beta.out <- as.vector(pgp.mmodel1@beta) + params.out <- pgp.mmodel1@covparams - expect_length(beta.out, 31) - expect_length(params.out, 15) - expect_equal(beta.out, c(0, 0.96, 0.93, 0.92, 0.89, rep(0, 2), 0.03, - rep(0,3), 0.57, 0.72, 1.11, 1, 0, 0.06, rep(0,4), - 0.93, 0.87, 1.03, 0.92, 0.05, rep(0, 2), 0.01, - 0.05, 0), tolerance=0.05) - expect_equal(params.out[1], 1.7, tolerance=2.5) - expect_equal(params.out[2], 2.9, tolerance=3.3) - expect_equal(params.out[3], 3.2, tolerance=5) - expect_equal(params.out[4], 0.71, tolerance=10) - expect_equal(params.out[5], 0.56, tolerance=2.8) - expect_equal(params.out[6], 0.4, tolerance=1.7) - expect_equal(params.out[7], 0.63, tolerance=1.3) - expect_equal(params.out[8], 0.47, tolerance=1.2) - expect_equal(params.out[9], 0.74, tolerance=1.1) - expect_equal(params.out[10], 1.8, tolerance=1.4) - expect_equal(params.out[11], 2.6, tolerance=2.4) - expect_equal(params.out[12], 1.4, tolerance=0.8) - expect_equal(params.out[13], 0.15, tolerance=0.6) - expect_equal(params.out[14], 0.32, tolerance=0.4) - expect_equal(params.out[15], 0.17, tolerance=0.3) + expect_length(beta.out, 31) + expect_length(params.out, 15) + expect_equal(beta.out, c( + 0, 0.96, 0.93, 0.92, 0.89, rep(0, 2), 0.03, + rep(0, 3), 0.57, 0.72, 1.11, 1, 0, 0.06, rep(0, 4), + 0.93, 0.87, 1.03, 0.92, 0.05, rep(0, 2), 0.01, + 0.05, 0 + ), tolerance = 0.05) + expect_equal(params.out[1], 1.7, tolerance = 2.5) + expect_equal(params.out[2], 2.9, tolerance = 3.3) + expect_equal(params.out[3], 3.2, tolerance = 5) + expect_equal(params.out[4], 0.71, tolerance = 10) + expect_equal(params.out[5], 0.56, tolerance = 2.8) + expect_equal(params.out[6], 0.4, tolerance = 1.7) + expect_equal(params.out[7], 0.63, tolerance = 1.3) + expect_equal(params.out[8], 0.47, tolerance = 1.2) + expect_equal(params.out[9], 0.74, tolerance = 1.1) + expect_equal(params.out[10], 1.8, tolerance = 1.4) + expect_equal(params.out[11], 2.6, tolerance = 2.4) + expect_equal(params.out[12], 1.4, tolerance = 0.8) + expect_equal(params.out[13], 0.15, tolerance = 0.6) + expect_equal(params.out[14], 0.32, tolerance = 0.4) + expect_equal(params.out[15], 0.17, tolerance = 0.3) }) test_that("Simulated dataset multivariate spatiotemporal", { - source("sim_multivariate_big_st.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 25) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list, X.st, locs.list, - scaling = c(1, 1, 2), verbose = FALSE, - optim.control = list( - trace = 0, maxit = 5000, - reltol = 1e-3 - ) + source("sim_multivariate_big_st.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 25) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list, X.st, locs.list, + scaling = c(1, 1, 2), verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 ) - beta.out <- as.vector(pgp.mmodel1@beta) - params.out <- pgp.mmodel1@covparams + ) + beta.out <- as.vector(pgp.mmodel1@beta) + params.out <- pgp.mmodel1@covparams - expect_length(beta.out, 31) - expect_length(params.out, 18) - expect_equal(beta.out, c( - 0, 0.91, 0.86, 0.82, 0.97, rep(0, 6), 0.95, - 0.97, 0.92, 0.78, rep(0, 6), 0.8, 0.97, - 1.04, 0.81, rep(0, 6) - ), tolerance = 1.1) + expect_length(beta.out, 31) + expect_length(params.out, 18) + expect_equal(beta.out, c( + 0, 0.91, 0.86, 0.82, 0.97, rep(0, 6), 0.95, + 0.97, 0.92, 0.78, rep(0, 6), 0.8, 0.97, + 1.04, 0.81, rep(0, 6) + ), tolerance = 1.1) }) test_that("Invalid locs input for prediction", { - source("sim_multivariate_small_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=5) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - expect_error(prestogp_predict(pgp.mmodel1, X.st.otst, as.matrix(1:3)), - "locs must be a list for multivariate models") + source("sim_multivariate_small_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 5) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + expect_error( + prestogp_predict(pgp.mmodel1, X.st.otst, as.matrix(1:3)), + "locs must be a list for multivariate models" + ) }) test_that("Invalid X input for prediction", { - source("sim_multivariate_small_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=5) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - expect_error(prestogp_predict(pgp.mmodel1, as.matrix(1:3), locs.list.otst), - "X must be a list for multivariate models") + source("sim_multivariate_small_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 5) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + expect_error( + prestogp_predict(pgp.mmodel1, as.matrix(1:3), locs.list.otst), + "X must be a list for multivariate models" + ) }) test_that("locs/X length mismatch for prediction", { - source("sim_multivariate_small_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=5) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - expect_error(prestogp_predict(pgp.mmodel1, X.st.otst, list(as.matrix(1:3))), - "locs and X must have the same length") + source("sim_multivariate_small_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 5) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + expect_error( + prestogp_predict(pgp.mmodel1, X.st.otst, list(as.matrix(1:3))), + "locs and X must have the same length" + ) }) test_that("locs/locs_train length mismatch", { - source("sim_multivariate_small_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=5) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - expect_error(prestogp_predict(pgp.mmodel1, list(1:3, 1:3, 1:3), - list(1:3, 1:3, 1:3)), - "Training and test set locs must have the same length") + source("sim_multivariate_small_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 5) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + expect_error( + prestogp_predict( + pgp.mmodel1, list(1:3, 1:3, 1:3), + list(1:3, 1:3, 1:3) + ), + "Training and test set locs must have the same length" + ) }) test_that("locs not a matrix for prediction", { - source("sim_multivariate_small_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=5) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - expect_error(prestogp_predict(pgp.mmodel1, list(as.matrix(1:3), - as.matrix(1:3)), - list(1:3, as.matrix(1:3))), - "Each locs must be a matrix") + source("sim_multivariate_small_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 5) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + expect_error( + prestogp_predict( + pgp.mmodel1, list( + as.matrix(1:3), + as.matrix(1:3) + ), + list(1:3, as.matrix(1:3)) + ), + "Each locs must be a matrix" + ) }) test_that("ncol(locs) != ncol(locs_train)", { - source("sim_multivariate_small_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=5) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - expect_error(prestogp_predict(pgp.mmodel1, list(as.matrix(1:3), - as.matrix(1:3)), - list(as.matrix(1:3), as.matrix(1:3))), - "All locs must have the same number of columns as locs_train") + source("sim_multivariate_small_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 5) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + expect_error( + prestogp_predict( + pgp.mmodel1, list( + as.matrix(1:3), + as.matrix(1:3) + ), + list(as.matrix(1:3), as.matrix(1:3)) + ), + "All locs must have the same number of columns as locs_train" + ) }) test_that("X not a matrix for prediction", { - source("sim_multivariate_small_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=5) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - X.st.otst[[1]] <- as.vector(X.st.otst[[1]]) - expect_error(prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst), - "Each X must be a matrix") + source("sim_multivariate_small_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 5) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + X.st.otst[[1]] <- as.vector(X.st.otst[[1]]) + expect_error( + prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst), + "Each X must be a matrix" + ) }) test_that("nrow(X) != nrow(locs) for prediction", { - source("sim_multivariate_small_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=5) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - locs.list.otst[[1]] <- rbind(locs.list.otst[[1]], rep(0, 2)) - expect_error(prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst), - "Each X must have the same number of rows as locs") + source("sim_multivariate_small_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 5) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + locs.list.otst[[1]] <- rbind(locs.list.otst[[1]], rep(0, 2)) + expect_error( + prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst), + "Each X must have the same number of rows as locs" + ) }) test_that("ncol(X) != ncol(X_train) for prediction", { - source("sim_multivariate_small_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=5) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - X.st.otst[[1]] <- cbind(X.st.otst[[1]], rep(0, nrow(X.st.otst[[1]]))) - expect_error(prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst), - "X and X_train must have the same number of predictors") + source("sim_multivariate_small_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 5) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + X.st.otst[[1]] <- cbind(X.st.otst[[1]], rep(0, nrow(X.st.otst[[1]]))) + expect_error( + prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst), + "X and X_train must have the same number of predictors" + ) }) test_that("m too small for prediction", { - source("sim_multivariate_small_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=5) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - expect_error(prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst, m=2), - "m must be at least 3") + source("sim_multivariate_small_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 5) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + expect_error( + prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst, m = 2), + "m must be at least 3" + ) }) test_that("m too large for prediction", { - source("sim_multivariate_small_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=5) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - expect_warning(prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst, - m=51), - "Conditioning set size m chosen to be >=n. Changing to m=n-1") + source("sim_multivariate_small_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 5) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + expect_warning( + prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst, + m = 51 + ), + "Conditioning set size m chosen to be >=n. Changing to m=n-1" + ) }) test_that("Simulated dataset multivariate spatial prediction", { - source("sim_multivariate_big_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 25) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, - scaling = c(1, 1), - apanasovich = TRUE, verbose = FALSE, - optim.control = list( - trace = 0, maxit = 5000, - reltol = 1e-3 - ) + source("sim_multivariate_big_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 25) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 ) + ) - pgp.mmodel1.pred <- prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst) + pgp.mmodel1.pred <- prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst) - mse <- mean((pgp.mmodel1.pred$means - unlist(y.list.otst))^2) - me <- mean(pgp.mmodel1.pred$means - unlist(y.list.otst)) + mse <- mean((pgp.mmodel1.pred$means - unlist(y.list.otst))^2) + me <- mean(pgp.mmodel1.pred$means - unlist(y.list.otst)) - expect_equal(mse, 1.99, tolerance = 0.1) - expect_equal(me, -0.04, tolerance = 0.03) + expect_equal(mse, 1.99, tolerance = 0.1) + expect_equal(me + 0.04, 0, tolerance = 0.03) }) diff --git a/tests/testthat/test-PrestoGP_Univariate.R b/tests/testthat/test-PrestoGP_Univariate.R index 3ee7435..89747c6 100755 --- a/tests/testthat/test-PrestoGP_Univariate.R +++ b/tests/testthat/test-PrestoGP_Univariate.R @@ -1,242 +1,326 @@ test_that("Invalid locs input", { model <- new("VecchiaModel") - expect_error(prestogp_fit(model, as.matrix(1:3), as.matrix(1:3), - 1:3), - "locs must be a matrix") + expect_error( + prestogp_fit( + model, as.matrix(1:3), as.matrix(1:3), + 1:3 + ), + "locs must be a matrix" + ) }) test_that("Invalid Y input", { model <- new("VecchiaModel") - expect_error(prestogp_fit(model, "y", as.matrix(1:3), - as.matrix(1:3)), - "Y must be a numeric vector or matrix") + expect_error( + prestogp_fit( + model, "y", as.matrix(1:3), + as.matrix(1:3) + ), + "Y must be a numeric vector or matrix" + ) }) test_that("Invalid X input", { model <- new("VecchiaModel") - expect_error(prestogp_fit(model, as.matrix(1:3), 1:3, - as.matrix(1:3)), - "X must be a matrix") + expect_error( + prestogp_fit( + model, as.matrix(1:3), 1:3, + as.matrix(1:3) + ), + "X must be a matrix" + ) }) test_that("Y has multiple columns", { model <- new("VecchiaModel") - expect_error(prestogp_fit(model, matrix(1:4, nrow=2), as.matrix(1:4), - as.matrix(1:4)), - "Y must have only 1 column") + expect_error( + prestogp_fit( + model, matrix(1:4, nrow = 2), as.matrix(1:4), + as.matrix(1:4) + ), + "Y must have only 1 column" + ) }) test_that("nrow(Y) != nrow(locs)", { model <- new("VecchiaModel") - expect_error(prestogp_fit(model, as.matrix(1:4), as.matrix(1:4), - as.matrix(1:3)), - "Y must have the same number of rows as locs") + expect_error( + prestogp_fit( + model, as.matrix(1:4), as.matrix(1:4), + as.matrix(1:3) + ), + "Y must have the same number of rows as locs" + ) }) test_that("nrow(Y) != nrow(X)", { model <- new("VecchiaModel") - expect_error(prestogp_fit(model, as.matrix(1:4), as.matrix(1:3), - as.matrix(1:4)), - "Y must have the same number of rows as X") + expect_error( + prestogp_fit( + model, as.matrix(1:4), as.matrix(1:3), + as.matrix(1:4) + ), + "Y must have the same number of rows as X" + ) }) test_that("Simulated dataset spatial", { - load("sim_vecchia.RData") - pgp.model1 <- new("VecchiaModel", n_neighbors=25) - pgp.model1 <- prestogp_fit(pgp.model1, y, X, locs, - scaling=c(1,1), apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - beta.out <- as.vector(pgp.model1@beta) - params.out <- pgp.model1@covparams - - expect_length(beta.out, 11) - expect_length(params.out, 5) - expect_equal(beta.out, c(0.01, 0.86, 0.98, 0.94, 0.9, rep(0, 6)), - tolerance=0.03) - expect_equal(params.out[1], 1.6, tolerance=0.5) - expect_equal(params.out[2], 0.4, tolerance=0.2) - expect_equal(params.out[3], 0.59, tolerance=0.2) - expect_equal(params.out[4], 2.0, tolerance=0.15) - - pgp.model2 <- new("FullModel") - pgp.model2 <- prestogp_fit(pgp.model2, y, X, locs, - scaling=c(1,1), apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - beta.out2 <- as.vector(pgp.model2@beta) - params.out2 <- pgp.model2@covparams - - expect_length(beta.out2, 11) - expect_length(params.out2, 5) - expect_equal(beta.out2, c(0.01, 0.86, 0.98, 0.95, 0.9, rep(0, 6)), - tolerance=0.03) - expect_equal(params.out2[1], 1.5, tolerance=0.6) - expect_equal(params.out2[2], 0.4, tolerance=0.15) - expect_equal(params.out2[3], 0.62, tolerance=0.2) - expect_equal(params.out2[4], 2.0, tolerance=0.15) - - # Vecchia and full models should be approximately equal - expect_equal(beta.out[1], beta.out2[1], tolerance=0.07) - expect_equal(beta.out[-1], beta.out2[-1], tolerance=0.04) - expect_equal(params.out[1], params.out2[1], tolerance=1) - expect_equal(params.out[2]-params.out2[2], 0, tolerance=0.2) - expect_equal(params.out[3], params.out2[3], tolerance=0.3) - expect_equal(params.out[4], params.out2[4], tolerance=0.2) + load("sim_vecchia.RData") + pgp.model1 <- new("VecchiaModel", n_neighbors = 25) + pgp.model1 <- prestogp_fit(pgp.model1, y, X, locs, + scaling = c(1, 1), apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + beta.out <- as.vector(pgp.model1@beta) + params.out <- pgp.model1@covparams + + expect_length(beta.out, 11) + expect_length(params.out, 5) + expect_equal(beta.out, c(0.01, 0.86, 0.98, 0.94, 0.9, rep(0, 6)), + tolerance = 0.03 + ) + expect_equal(params.out[1], 1.6, tolerance = 0.5) + expect_equal(params.out[2], 0.4, tolerance = 0.2) + expect_equal(params.out[3], 0.59, tolerance = 0.2) + expect_equal(params.out[4], 2.0, tolerance = 0.15) + + pgp.model2 <- new("FullModel") + pgp.model2 <- prestogp_fit(pgp.model2, y, X, locs, + scaling = c(1, 1), apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + beta.out2 <- as.vector(pgp.model2@beta) + params.out2 <- pgp.model2@covparams + + expect_length(beta.out2, 11) + expect_length(params.out2, 5) + expect_equal(beta.out2, c(0.01, 0.86, 0.98, 0.95, 0.9, rep(0, 6)), + tolerance = 0.03 + ) + expect_equal(params.out2[1], 1.5, tolerance = 0.6) + expect_equal(params.out2[2], 0.4, tolerance = 0.15) + expect_equal(params.out2[3], 0.62, tolerance = 0.2) + expect_equal(params.out2[4], 2.0, tolerance = 0.15) + + # Vecchia and full models should be approximately equal + expect_equal(beta.out[1], beta.out2[1], tolerance = 0.07) + expect_equal(beta.out[-1], beta.out2[-1], tolerance = 0.04) + expect_equal(params.out[1], params.out2[1], tolerance = 1) + expect_equal(params.out[2] - params.out2[2], 0, tolerance = 0.2) + expect_equal(params.out[3], params.out2[3], tolerance = 0.3) + expect_equal(params.out[4], params.out2[4], tolerance = 0.2) }) test_that("Simulated dataset spatiotemporal", { - source("sim_vecchia_st.R") - pgp.model1 <- new("VecchiaModel", n_neighbors=25) - pgp.model1 <- prestogp_fit(pgp.model1, y, X, locs, scaling=c(1,1,2), - apanasovich=FALSE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - beta.out <- as.vector(pgp.model1@beta) - params.out <- pgp.model1@covparams - - expect_length(beta.out, 11) - expect_length(params.out, 6) - expect_equal(beta.out, c(0.01, 0.93, 1.01, 0.91, 0.99, rep(0, 6)), - tolerance=0.015) - expect_equal(params.out[1], 1.7, tolerance=0.55) - expect_equal(params.out[2]-0.19, 0, tolerance=0.05) - expect_equal(params.out[3], 0.22, tolerance=0.05) - expect_equal(params.out[4], 1.12, tolerance=0.3) - expect_equal(params.out[5], 0.62, tolerance=0.05) - - pgp.model2 <- new("FullModel", n_neighbors=25) - pgp.model2 <- prestogp_fit(pgp.model2, y, X, locs, scaling=c(1,1,2), - apanasovich=FALSE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - beta.out2 <- as.vector(pgp.model2@beta) - params.out2 <- pgp.model2@covparams - - expect_length(beta.out2, 11) - expect_length(params.out2, 6) - expect_equal(beta.out2, c(-0.03, 0.93, 1, 0.91, 0.98, rep(0, 6)), - tolerance=0.02) - expect_equal(params.out2[1], 1.6, tolerance=0.5) - expect_equal(params.out2[2]-0.19, 0, tolerance=0.05) - expect_equal(params.out2[3]-0.22, 0, tolerance=0.05) - expect_equal(params.out2[4], 1.18, tolerance=0.15) - expect_equal(params.out2[5], 0.64, tolerance=0.05) - - # Vecchia and full models should be approximately equal - expect_equal(beta.out[1], beta.out2[1], tolerance=0.08) - expect_equal(beta.out[-1], beta.out2[-1], tolerance=0.03) - expect_equal(params.out[1], params.out2[1], tolerance=0.6) - expect_equal(params.out[2], params.out2[2], tolerance=0.06) - expect_equal(params.out[3]-params.out2[3], 0, tolerance=0.06) - expect_equal(params.out[4], params.out2[4], tolerance=0.3) - expect_equal(params.out[5], params.out2[5], tolerance=0.1) + source("sim_vecchia_st.R") + pgp.model1 <- new("VecchiaModel", n_neighbors = 25) + pgp.model1 <- prestogp_fit(pgp.model1, y, X, locs, + scaling = c(1, 1, 2), + apanasovich = FALSE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + beta.out <- as.vector(pgp.model1@beta) + params.out <- pgp.model1@covparams + + expect_length(beta.out, 11) + expect_length(params.out, 6) + expect_equal(beta.out, c(0.01, 0.93, 1.01, 0.91, 0.99, rep(0, 6)), + tolerance = 0.015 + ) + expect_equal(params.out[1], 1.7, tolerance = 0.55) + expect_equal(params.out[2] - 0.19, 0, tolerance = 0.05) + expect_equal(params.out[3], 0.22, tolerance = 0.05) + expect_equal(params.out[4], 1.12, tolerance = 0.3) + expect_equal(params.out[5], 0.62, tolerance = 0.05) + + pgp.model2 <- new("FullModel", n_neighbors = 25) + pgp.model2 <- prestogp_fit(pgp.model2, y, X, locs, + scaling = c(1, 1, 2), + apanasovich = FALSE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + beta.out2 <- as.vector(pgp.model2@beta) + params.out2 <- pgp.model2@covparams + + expect_length(beta.out2, 11) + expect_length(params.out2, 6) + expect_equal(beta.out2, c(-0.03, 0.93, 1, 0.91, 0.98, rep(0, 6)), + tolerance = 0.02 + ) + expect_equal(params.out2[1], 1.6, tolerance = 0.5) + expect_equal(params.out2[2] - 0.19, 0, tolerance = 0.05) + expect_equal(params.out2[3] - 0.22, 0, tolerance = 0.05) + expect_equal(params.out2[4], 1.18, tolerance = 0.15) + expect_equal(params.out2[5], 0.64, tolerance = 0.05) + + # Vecchia and full models should be approximately equal + expect_equal(beta.out[1], beta.out2[1], tolerance = 0.08) + expect_equal(beta.out[-1], beta.out2[-1], tolerance = 0.03) + expect_equal(params.out[1], params.out2[1], tolerance = 0.6) + expect_equal(params.out[2], params.out2[2], tolerance = 0.06) + expect_equal(params.out[3] - params.out2[3], 0, tolerance = 0.06) + expect_equal(params.out[4], params.out2[4], tolerance = 0.3) + expect_equal(params.out[5], params.out2[5], tolerance = 0.1) }) test_that("Invalid locs input for prediction", { - source("sim_vecchia_small_pred.R") - pgp.model1 <- new("VecchiaModel", n_neighbors=5) - pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, - locs.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - expect_error(prestogp_predict(pgp.model1, X.otst, 1:3), - "locs must be a matrix") + source("sim_vecchia_small_pred.R") + pgp.model1 <- new("VecchiaModel", n_neighbors = 5) + pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, + locs.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + expect_error( + prestogp_predict(pgp.model1, X.otst, 1:3), + "locs must be a matrix" + ) }) test_that("Invalid X input for prediction", { - source("sim_vecchia_small_pred.R") - pgp.model1 <- new("VecchiaModel", n_neighbors=5) - pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, - locs.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - expect_error(prestogp_predict(pgp.model1, 1:3, locs.otst), - "X must be a matrix") + source("sim_vecchia_small_pred.R") + pgp.model1 <- new("VecchiaModel", n_neighbors = 5) + pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, + locs.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + expect_error( + prestogp_predict(pgp.model1, 1:3, locs.otst), + "X must be a matrix" + ) }) test_that("ncol(locs) != ncol(locs_train)", { - source("sim_vecchia_small_pred.R") - pgp.model1 <- new("VecchiaModel", n_neighbors=5) - pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, - locs.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - locs.otst <- cbind(locs.otst, rep(1, nrow(locs.otst))) - expect_error(prestogp_predict(pgp.model1, X.otst, locs.otst), - "locs must have the same number of columns as locs_train") + source("sim_vecchia_small_pred.R") + pgp.model1 <- new("VecchiaModel", n_neighbors = 5) + pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, + locs.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + locs.otst <- cbind(locs.otst, rep(1, nrow(locs.otst))) + expect_error( + prestogp_predict(pgp.model1, X.otst, locs.otst), + "locs must have the same number of columns as locs_train" + ) }) test_that("nrow(X) != nrow(locs) for prediction", { - source("sim_vecchia_small_pred.R") - pgp.model1 <- new("VecchiaModel", n_neighbors=5) - pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, - locs.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - X.otst <- rbind(X.otst, rep(1, ncol(X.otst))) - expect_error(prestogp_predict(pgp.model1, X.otst, locs.otst), - "X must have the same number of rows as locs") + source("sim_vecchia_small_pred.R") + pgp.model1 <- new("VecchiaModel", n_neighbors = 5) + pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, + locs.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + X.otst <- rbind(X.otst, rep(1, ncol(X.otst))) + expect_error( + prestogp_predict(pgp.model1, X.otst, locs.otst), + "X must have the same number of rows as locs" + ) }) test_that("ncol(X) != ncol(X_train) for prediction", { - source("sim_vecchia_small_pred.R") - pgp.model1 <- new("VecchiaModel", n_neighbors=5) - pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, - locs.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - X.otst <- cbind(X.otst, rep(1, nrow(X.otst))) - expect_error(prestogp_predict(pgp.model1, X.otst, locs.otst), - "X and X_train must have the same number of predictors") + source("sim_vecchia_small_pred.R") + pgp.model1 <- new("VecchiaModel", n_neighbors = 5) + pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, + locs.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + X.otst <- cbind(X.otst, rep(1, nrow(X.otst))) + expect_error( + prestogp_predict(pgp.model1, X.otst, locs.otst), + "X and X_train must have the same number of predictors" + ) }) test_that("m too small for prediction", { - source("sim_vecchia_small_pred.R") - pgp.model1 <- new("VecchiaModel", n_neighbors=5) - pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, - locs.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - expect_error(prestogp_predict(pgp.model1, X.otst, locs.otst, m=2), - "m must be at least 3") + source("sim_vecchia_small_pred.R") + pgp.model1 <- new("VecchiaModel", n_neighbors = 5) + pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, + locs.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + expect_error( + prestogp_predict(pgp.model1, X.otst, locs.otst, m = 2), + "m must be at least 3" + ) }) test_that("m too large for prediction", { - source("sim_vecchia_small_pred.R") - pgp.model1 <- new("VecchiaModel", n_neighbors=5) - pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, - locs.otr, scaling=c(1,1), - apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - expect_warning(prestogp_predict(pgp.model1, X.otst, locs.otst, m=201), - "Conditioning set size m chosen to be >=n. Changing to m=n-1") + source("sim_vecchia_small_pred.R") + pgp.model1 <- new("VecchiaModel", n_neighbors = 5) + pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, + locs.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + expect_warning( + prestogp_predict(pgp.model1, X.otst, locs.otst, m = 201), + "Conditioning set size m chosen to be >=n. Changing to m=n-1" + ) }) test_that("Simulated spatial prediction", { - source("sim_vecchia_pred.R") - pgp.model1 <- new("VecchiaModel", n_neighbors = 25) - pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, locs.otr, - scaling = c(1, 1), - apanasovich = TRUE, verbose = FALSE, - optim.control = list( - trace = 0, maxit = 5000, - reltol = 1e-3 - ) + source("sim_vecchia_pred.R") + pgp.model1 <- new("VecchiaModel", n_neighbors = 25) + pgp.model1 <- prestogp_fit(pgp.model1, y.otr, X.otr, locs.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 ) + ) - pgp.model1.pred <- prestogp_predict(pgp.model1, X.otst, locs.otst) + pgp.model1.pred <- prestogp_predict(pgp.model1, X.otst, locs.otst) - mse <- mean((pgp.model1.pred$means - y.otst)^2) - me <- mean(pgp.model1.pred$means - y.otst) + mse <- mean((pgp.model1.pred$means - y.otst)^2) + me <- mean(pgp.model1.pred$means - y.otst) - expect_equal(mse, 2.24, tolerance = 0.07) - expect_equal(me, 0.03, tolerance = 0.05) + expect_equal(mse, 2.24, tolerance = 0.07) + expect_equal(me, 0.03, tolerance = 0.05) })