From 9b1020106f819a0ba290964a3f037cdfa4b19447 Mon Sep 17 00:00:00 2001 From: Shail Choksi Date: Sat, 13 Jan 2024 12:14:22 -0500 Subject: [PATCH 1/3] Pull request #31: Fixed merge issues caused by linting. Fixed tests Merge in STAT/prestogp from master to to-git Squashed commit of the following: commit db6cd4f6f86e760ee8393adf224cd0231ccde8cf Merge: 55607e2 7da2588 Author: Eric Bair Date: Sat Jan 13 12:11:56 2024 -0500 Pull request #29: Fixed some bugs caused by the linter Merge in STAT/prestogp from eb-dev to master * commit '7da2588abba4491c1ff490de95e834d00f5df88e': Fixed some bugs caused by the linter commit 7da2588abba4491c1ff490de95e834d00f5df88e Author: Eric Bair Date: Sat Jan 13 01:37:58 2024 -0500 Fixed some bugs caused by the linter commit 55607e291477d174b8323bf258a8e9aa5c4db84c Merge: c7ed481 81b675d Author: Shail Choksi Date: Fri Jan 12 19:37:55 2024 -0500 Pull request #28: Final linting fixes Merge in STAT/prestogp from build-workflow to master * commit '81b675dda9b57f2a362c34876364adf5bcf45560': Fix remaining linting issues Fix all indentation errors Add indentation_linter configuration to lintr config file commit 81b675dda9b57f2a362c34876364adf5bcf45560 Author: sciome-bot Date: Fri Jan 12 19:36:05 2024 -0500 Fix remaining linting issues commit 06673a48f68b5be9e6ba480e8482db0dbd9d0c3a Author: sciome-bot Date: Fri Jan 12 19:23:11 2024 -0500 Fix all indentation errors commit c7ed481479c2ca5a4fc111f429affd01b6c3ce53 Merge: 723c3eb 7bcfa1a Author: Shail Choksi Date: Fri Jan 12 19:02:59 2024 -0500 Pull request #27: Lintr fixes. Added release and sanitizer actions Merge in STAT/prestogp from build-workflow to master * commit '7bcfa1aca78cc6b994cc0227c083357c58d50ca0': (21 commits) Add R_LINTR_LINTER_FILE to lint action to point to global .lintr file More lintr fixes for 1:length, 1:nrow and line length. Disabled object_length_linter WIP: Fix linter warnings for 1:nrow, 1:ncol and 1:length. Increase line length to 160 chars Don't build vignettes during the release action Move release action file to correct directory Add release action Reorder imports in RcppExports file Add missing dependency in Namespace/Description Remove unneeded exports from NAMESPACE Add missing comma in Imports section of DESCRIPTION Add missing comma in imports section of DESCRIPTION Rerun auto-formatter WIP - linting Add ignore rules for .lintr and .github for R build remove linter options from lint action as we have added .lintr project file. Fix all vector_logic_linter warnings Remove line length linter to see the remaining errors/warnings workaround for lintr bug: https://github.com/REditorSupport/languageserver/issues/89 Run auto-lint on vscode enable linting on build-workflow Add missing imports in DESCRIPTION ... commit 29d82cd628a7662ab2fe3fb8f8b03088e3d08bd0 Author: sciome-bot Date: Fri Jan 12 18:40:17 2024 -0500 Add indentation_linter configuration to lintr config file commit 7bcfa1aca78cc6b994cc0227c083357c58d50ca0 Author: sciome-bot Date: Fri Jan 12 18:16:01 2024 -0500 Add R_LINTR_LINTER_FILE to lint action to point to global .lintr file commit b6df473e4109722884c70262d504a5d9ed3ca6da Author: sciome-bot Date: Fri Jan 12 17:54:50 2024 -0500 More lintr fixes for 1:length, 1:nrow and line length. Disabled object_length_linter commit 6f6f91dc6a81372e1ff2bfb106d5528760485f01 Author: sciome-bot Date: Fri Jan 12 17:32:12 2024 -0500 WIP: Fix linter warnings for 1:nrow, 1:ncol and 1:length. Increase line length to 160 chars commit 2db2e1946bc0804dbb616fe54001e89e0b70b233 Author: sciome-bot Date: Fri Jan 12 16:48:59 2024 -0500 Don't build vignettes during the release action commit c53922341968e64d1b25920cff2175f91d871126 Author: sciome-bot Date: Fri Jan 12 16:39:52 2024 -0500 Move release action file to correct directory commit f83df2bdc9dc12aca807252596a67717e4a85253 Author: sciome-bot Date: Fri Jan 12 16:36:03 2024 -0500 Add release action commit c84e15959f0c5073ee6ff05ec2e4ed9fa3564fa3 Author: sciome-bot Date: Fri Jan 12 16:04:36 2024 -0500 Reorder imports in RcppExports file commit 0460202bf717b996cec224ec12c54f5bd2c66c32 Merge: 245eaa3 43ae272 Author: sciome-bot Date: Fri Jan 12 15:48:30 2024 -0500 Merge branch 'build-workflow' of sciome-bot-git:Spatiotemporal-Exposures-and-Toxicology/PrestoGP into build-workflow commit 245eaa33e9891125cdd6b3d2b8368868d4e73ad6 Merge: 5355e24 c060815 Author: sciome-bot Date: Fri Jan 12 15:47:42 2024 -0500 Merge branch 'main' of sciome-bot-git:Spatiotemporal-Exposures-and-Toxicology/PrestoGP into build-workflow commit c060815eefd29ecf38e2ad1ba25ae25dbfe46b21 Merge: 731da76 6d69b69 Author: {SET}group <127860447+Spatiotemporal-Exposures-and-Toxicology@users.noreply.github.com> Date: Thu Jan 11 16:02:52 2024 -0500 Merge pull request #44 from Spatiotemporal-Exposures-and-Toxicology/main-sciome Sciome Update 1/10/2024 commit 6d69b692cc54c815cfd76f8c844ea92b8f2811a3 Author: sciome-bot Date: Thu Jan 11 14:39:39 2024 -0500 Add missing dependency in Namespace/Description commit bdf52c697b24a5b7916a9fa8ec6263c874c2b24a Author: sciome-bot Date: Thu Jan 11 14:30:38 2024 -0500 Remove unneeded exports from NAMESPACE commit 06918b038775fb8374409f92750915b3c2483c93 Author: sciome-bot Date: Thu Jan 11 14:15:45 2024 -0500 Add missing comma in Imports section of DESCRIPTION ... and 13 more commits --- .Rbuildignore | 3 +- .github/workflows/lint.yaml | 6 +- .github/workflows/release.yaml | 109 +++ .github/workflows/sanitizers.yaml | 24 +- .lintr | 13 + DESCRIPTION | 8 +- NAMESPACE | 4 - R/Log_Likelihood.R | 59 +- R/PrestoGP_CreateU_Multivariate.R | 71 +- R/PrestoGP_Model.R | 225 +++--- R/PrestoGP_Multivariate_Vecchia.R | 162 +++-- R/PrestoGP_Util_Functions.R | 34 +- R/PrestoGP_Vecchia.R | 97 +-- R/RcppExports.R | 0 inst/doc/test.R | 1 - man/compute_error-PrestoGPModel-method.Rd | 3 +- man/estimate_betas-PrestoGPModel-method.Rd | 2 +- man/prestogp_fit-PrestoGPModel-method.Rd | 17 +- man/prestogp_predict-VecchiaModel-method.Rd | 10 +- man/show_theta-PrestoGPModel-method.Rd | 3 +- man/sparseNN.Rd | 3 +- src/RcppExports.cpp | 0 src/createUMC.cpp | 105 ++- tests/testthat/sim_multivariate.R | 129 ++-- tests/testthat/sim_multivariate_big.R | 70 +- tests/testthat/sim_multivariate_big_pred.R | 90 ++- tests/testthat/sim_multivariate_big_st.R | 78 +-- tests/testthat/sim_multivariate_lik.R | 64 +- tests/testthat/sim_multivariate_small.R | 70 +- tests/testthat/sim_multivariate_small_pred.R | 93 ++- tests/testthat/sim_vecchia.R | 26 +- tests/testthat/sim_vecchia_pred.R | 38 +- tests/testthat/sim_vecchia_small.R | 26 +- tests/testthat/sim_vecchia_small_pred.R | 38 +- tests/testthat/sim_vecchia_st.R | 32 +- tests/testthat/test-Log_Likelihood.R | 651 +++++++++--------- .../test-PrestoGP_CreateU_Multivariate.R | 266 +++---- tests/testthat/test-PrestoGP_Model.R | 160 +++-- ...test-PrestoGP_Multivariate_Vecchia_Model.R | 547 +++++++++------ tests/testthat/test-PrestoGP_Univariate.R | 454 +++++++----- 40 files changed, 2108 insertions(+), 1683 deletions(-) create mode 100644 .github/workflows/release.yaml create mode 100644 .lintr mode change 100644 => 100755 R/PrestoGP_Util_Functions.R mode change 100644 => 100755 R/RcppExports.R mode change 100644 => 100755 man/compute_error-PrestoGPModel-method.Rd mode change 100644 => 100755 man/estimate_betas-PrestoGPModel-method.Rd mode change 100644 => 100755 man/prestogp_predict-VecchiaModel-method.Rd mode change 100644 => 100755 man/show_theta-PrestoGPModel-method.Rd mode change 100644 => 100755 man/sparseNN.Rd mode change 100644 => 100755 src/RcppExports.cpp 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) }) From ee193de709a301ce570c70945402639e30540fac Mon Sep 17 00:00:00 2001 From: sciome-bot Date: Sat, 13 Jan 2024 12:29:53 -0500 Subject: [PATCH 2/3] Fix lingering linting issues --- R/PrestoGP_Multivariate_Vecchia.R | 4 +--- R/RcppExports.R | 7 +++---- 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/R/PrestoGP_Multivariate_Vecchia.R b/R/PrestoGP_Multivariate_Vecchia.R index 642c125..65810e2 100755 --- a/R/PrestoGP_Multivariate_Vecchia.R +++ b/R/PrestoGP_Multivariate_Vecchia.R @@ -110,9 +110,7 @@ setMethod("prestogp_predict", "MultivariateVecchiaModel", 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]]) + 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) } diff --git a/R/RcppExports.R b/R/RcppExports.R index dd16aaf..c50230c 100755 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -2,10 +2,9 @@ # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 na_omit_c <- function(x) { - .Call('_PrestoGP_na_omit_c', PACKAGE = 'PrestoGP', x) + .Call('_PrestoGP_na_omit_c', PACKAGE = 'PrestoGP', x) } createU_helper_mat <- function(olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning) { - .Call('_PrestoGP_createU_helper_mat', PACKAGE = 'PrestoGP', olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning) -} - + .Call('_PrestoGP_createU_helper_mat', PACKAGE = 'PrestoGP', olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning) +} \ No newline at end of file From 8ffb20ba76ae39eeeb80de5330f87c266bddb5e9 Mon Sep 17 00:00:00 2001 From: Shail Choksi Date: Sun, 14 Jan 2024 11:00:33 -0500 Subject: [PATCH 3/3] Pull request #34: Fix tolerance and expected values for tests Merge in STAT/prestogp from master to to-git Squashed commit of the following: commit 439790cb0d96d85e80a3782ad533eee78afe41dc Merge: f4af0d1 56f881f Author: Eric Bair Date: Sun Jan 14 10:58:49 2024 -0500 Pull request #32: Fixed some testthat issues Merge in STAT/prestogp from eb-dev to master * commit '56f881f96a0e55e60bca3c8a840ec31291665594': Fixed some testthat issues commit 56f881f96a0e55e60bca3c8a840ec31291665594 Author: Eric Bair Date: Sat Jan 13 14:05:12 2024 -0500 Fixed some testthat issues commit f4af0d1be3490dd2edd4f4fd0ae97bf729464748 Merge: db6cd4f ee193de Author: sciome-bot Date: Sat Jan 13 12:37:28 2024 -0500 Merge branch 'main-sciome' of sciome-bot-git:Spatiotemporal-Exposures-and-Toxicology/PrestoGP commit db6cd4f6f86e760ee8393adf224cd0231ccde8cf Merge: 55607e2 7da2588 Author: Eric Bair Date: Sat Jan 13 12:11:56 2024 -0500 Pull request #29: Fixed some bugs caused by the linter Merge in STAT/prestogp from eb-dev to master * commit '7da2588abba4491c1ff490de95e834d00f5df88e': Fixed some bugs caused by the linter commit 7da2588abba4491c1ff490de95e834d00f5df88e Author: Eric Bair Date: Sat Jan 13 01:37:58 2024 -0500 Fixed some bugs caused by the linter commit 55607e291477d174b8323bf258a8e9aa5c4db84c Merge: c7ed481 81b675d Author: Shail Choksi Date: Fri Jan 12 19:37:55 2024 -0500 Pull request #28: Final linting fixes Merge in STAT/prestogp from build-workflow to master * commit '81b675dda9b57f2a362c34876364adf5bcf45560': Fix remaining linting issues Fix all indentation errors Add indentation_linter configuration to lintr config file commit 81b675dda9b57f2a362c34876364adf5bcf45560 Author: sciome-bot Date: Fri Jan 12 19:36:05 2024 -0500 Fix remaining linting issues commit 06673a48f68b5be9e6ba480e8482db0dbd9d0c3a Author: sciome-bot Date: Fri Jan 12 19:23:11 2024 -0500 Fix all indentation errors commit c7ed481479c2ca5a4fc111f429affd01b6c3ce53 Merge: 723c3eb 7bcfa1a Author: Shail Choksi Date: Fri Jan 12 19:02:59 2024 -0500 Pull request #27: Lintr fixes. Added release and sanitizer actions Merge in STAT/prestogp from build-workflow to master * commit '7bcfa1aca78cc6b994cc0227c083357c58d50ca0': (21 commits) Add R_LINTR_LINTER_FILE to lint action to point to global .lintr file More lintr fixes for 1:length, 1:nrow and line length. Disabled object_length_linter WIP: Fix linter warnings for 1:nrow, 1:ncol and 1:length. Increase line length to 160 chars Don't build vignettes during the release action Move release action file to correct directory Add release action Reorder imports in RcppExports file Add missing dependency in Namespace/Description Remove unneeded exports from NAMESPACE Add missing comma in Imports section of DESCRIPTION Add missing comma in imports section of DESCRIPTION Rerun auto-formatter WIP - linting Add ignore rules for .lintr and .github for R build remove linter options from lint action as we have added .lintr project file. Fix all vector_logic_linter warnings Remove line length linter to see the remaining errors/warnings workaround for lintr bug: https://github.com/REditorSupport/languageserver/issues/89 Run auto-lint on vscode enable linting on build-workflow Add missing imports in DESCRIPTION ... commit 29d82cd628a7662ab2fe3fb8f8b03088e3d08bd0 Author: sciome-bot Date: Fri Jan 12 18:40:17 2024 -0500 Add indentation_linter configuration to lintr config file commit 7bcfa1aca78cc6b994cc0227c083357c58d50ca0 Author: sciome-bot Date: Fri Jan 12 18:16:01 2024 -0500 Add R_LINTR_LINTER_FILE to lint action to point to global .lintr file commit b6df473e4109722884c70262d504a5d9ed3ca6da Author: sciome-bot Date: Fri Jan 12 17:54:50 2024 -0500 More lintr fixes for 1:length, 1:nrow and line length. Disabled object_length_linter commit 6f6f91dc6a81372e1ff2bfb106d5528760485f01 Author: sciome-bot Date: Fri Jan 12 17:32:12 2024 -0500 WIP: Fix linter warnings for 1:nrow, 1:ncol and 1:length. Increase line length to 160 chars commit 2db2e1946bc0804dbb616fe54001e89e0b70b233 Author: sciome-bot Date: Fri Jan 12 16:48:59 2024 -0500 Don't build vignettes during the release action commit c53922341968e64d1b25920cff2175f91d871126 Author: sciome-bot Date: Fri Jan 12 16:39:52 2024 -0500 Move release action file to correct directory commit f83df2bdc9dc12aca807252596a67717e4a85253 Author: sciome-bot Date: Fri Jan 12 16:36:03 2024 -0500 Add release action commit c84e15959f0c5073ee6ff05ec2e4ed9fa3564fa3 Author: sciome-bot Date: Fri Jan 12 16:04:36 2024 -0500 Reorder imports in RcppExports file commit 0460202bf717b996cec224ec12c54f5bd2c66c32 Merge: 245eaa3 43ae272 Author: sciome-bot Date: Fri Jan 12 15:48:30 2024 -0500 Merge branch 'build-workflow' of sciome-bot-git:Spatiotemporal-Exposures-and-Toxicology/PrestoGP into build-workflow commit 245eaa33e9891125cdd6b3d2b8368868d4e73ad6 Merge: 5355e24 c060815 Author: sciome-bot Date: Fri Jan 12 15:47:42 2024 -0500 Merge branch 'main' of sciome-bot-git:Spatiotemporal-Exposures-and-Toxicology/PrestoGP into build-workflow commit c060815eefd29ecf38e2ad1ba25ae25dbfe46b21 Merge: 731da76 6d69b69 Author: {SET}group <127860447+Spatiotemporal-Exposures-and-Toxicology@users.noreply.github.com> Date: Thu Jan 11 16:02:52 2024 -0500 Merge pull request #44 from Spatiotemporal-Exposures-and-Toxicology/main-sciome Sciome Update 1/10/2024 ... and 13 more commits --- DESCRIPTION | 2 +- NAMESPACE | 1 - tests/testthat/test-PrestoGP_Univariate.R | 4 ++-- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6b3e42e..db498eb 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.9022 +Version: 0.2.0.9023 Authors@R: c( person(given = "Eric", family = "Bair", diff --git a/NAMESPACE b/NAMESPACE index c303063..961e489 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -29,7 +29,6 @@ import(glmnet) import(ncvreg) import(readxl) import(scoringRules) -import(covr) importFrom(aod,wald.test) importFrom(dplyr,"%>%") importFrom(foreach,"%dopar%") diff --git a/tests/testthat/test-PrestoGP_Univariate.R b/tests/testthat/test-PrestoGP_Univariate.R index 89747c6..179128a 100755 --- a/tests/testthat/test-PrestoGP_Univariate.R +++ b/tests/testthat/test-PrestoGP_Univariate.R @@ -138,7 +138,7 @@ test_that("Simulated dataset spatiotemporal", { ) 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[3] - 0.22, 0, tolerance = 0.05) expect_equal(params.out[4], 1.12, tolerance = 0.3) expect_equal(params.out[5], 0.62, tolerance = 0.05) @@ -169,7 +169,7 @@ test_that("Simulated dataset spatiotemporal", { 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[2] - params.out2[2], 0, 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)