context("GPModel_non_Gaussian_data") # Avoid being tested on CRAN if(Sys.getenv("GPBOOST_ALL_TESTS") == "GPBOOST_ALL_TESTS"){ TOLERANCE_LOOSE <- 1e-3 TOLERANCE_STRICT <- 1E-6 TOLERANCE_ITERATIVE <- 1e-1 DEFAULT_OPTIM_PARAMS <- list(optimizer_cov = "gradient_descent", optimizer_coef = "gradient_descent", use_nesterov_acc = TRUE, lr_cov=0.1, lr_coef = 0.1, maxit = 1000) DEFAULT_OPTIM_PARAMS_STD <- c(DEFAULT_OPTIM_PARAMS, list(std_dev = TRUE)) # Function that simulates uniform random variables sim_rand_unif <- function(n, init_c=0.1){ mod_lcg <- 2^32 # modulus for linear congruential generator (random0 used) sim <- rep(NA, n) sim[1] <- floor(init_c * mod_lcg) for(i in 2:n) sim[i] <- (22695477 * sim[i-1] + 1) %% mod_lcg return(sim / mod_lcg) } # Simulate data n <- 100 # number of samples # Simulate locations / features of GP d <- 2 # dimension of GP locations coords <- matrix(sim_rand_unif(n=n*d, init_c=0.1), ncol=d) D <- as.matrix(dist(coords)) # Simulate GP sigma2_1 <- 1^2 # marginal variance of GP rho <- 0.1 # range parameter Sigma <- sigma2_1 * exp(-D/rho) + diag(1E-20,n) L <- t(chol(Sigma)) b_1 <- qnorm(sim_rand_unif(n=n, init_c=0.8)) # GP random coefficients Z_SVC <- matrix(sim_rand_unif(n=n*2, init_c=0.6), ncol=2) # covariate data for random coefficients colnames(Z_SVC) <- c("var1","var2") b_2 <- qnorm(sim_rand_unif(n=n, init_c=0.17)) b_3 <- qnorm(sim_rand_unif(n=n, init_c=0.42)) # First grouped random effects model m <- 10 # number of categories / levels for grouping variable group <- rep(1,n) # grouping variable for(i in 1:m) group[((i-1)*n/m+1):(i*n/m)] <- i Z1 <- model.matrix(rep(1,n) ~ factor(group) - 1) b_gr_1 <- qnorm(sim_rand_unif(n=m, init_c=0.565)) # Second grouped random effect n_obs_gr <- n/m # number of samples per group group2 <- rep(1,n) # grouping variable for(i in 1:m) group2[(1:n_obs_gr)+n_obs_gr*(i-1)] <- 1:n_obs_gr Z2 <- model.matrix(rep(1,n)~factor(group2)-1) b_gr_2 <- qnorm(sim_rand_unif(n=n_obs_gr, init_c=0.36)) # Grouped random slope / coefficient x <- cos((1:n-n/2)^2*5.5*pi/n) # covariate data for random slope Z3 <- diag(x) %*% Z1 b_gr_3 <- qnorm(sim_rand_unif(n=m, init_c=0.5678)) # Data for linear mixed effects model X <- cbind(rep(1,n),sin((1:n-n/2)^2*2*pi/n)) # design matrix / covariate data for fixed effect beta <- c(0.1,2) # regression coefficients # cluster_ids cluster_ids <- c(rep(1,0.4*n),rep(2,0.6*n)) # GP with multiple observations at the same locations coords_multiple <- matrix(sim_rand_unif(n=n*d/4, init_c=0.1), ncol=d) coords_multiple <- rbind(coords_multiple,coords_multiple,coords_multiple,coords_multiple) D_multiple <- as.matrix(dist(coords_multiple)) Sigma_multiple <- sigma2_1*exp(-D_multiple/rho)+diag(1E-10,n) L_multiple <- t(chol(Sigma_multiple)) b_multiple <- qnorm(sim_rand_unif(n=n, init_c=0.8)) test_that("Binary classification with Gaussian process model ", { probs <- pnorm(L %*% b_1) y <- as.numeric(sim_rand_unif(n=n, init_c=0.2341) < probs) init_cov_pars <- c(1,mean(dist(coords))/3) # Label needs to have correct format expect_error(fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", y = b_1, params = list(optimizer_cov = "gradient_descent"))) yw <- y yw[3] <- yw[3] + 1E-6 expect_error(fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", y = yw, params = list(optimizer_cov = "gradient_descent"))) # Only gradient descent can be used expect_error(fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", y = y, params = list(optimizer_cov = "fisher_scoring"))) # Estimation using gradient descent gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit") capture.output( fit(gp_model, y = y, params = list(optimizer_cov = "gradient_descent", lr_cov = 0.1, use_nesterov_acc = FALSE, convergence_criterion = "relative_change_in_parameters", init_cov_pars = init_cov_pars)), file='NUL') cov_pars <- c(0.9419234, 0.1866877) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 40) # Can switch between likelihoods gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit") gp_model$set_likelihood("gaussian") gp_model$set_likelihood("bernoulli_probit") capture.output( fit(gp_model, y = y, params = list(optimizer_cov = "gradient_descent", lr_cov = 0.1, use_nesterov_acc = FALSE, convergence_criterion = "relative_change_in_parameters", init_cov_pars = init_cov_pars)), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) # Estimation using gradient descent and Nesterov acceleration gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit") capture.output( fit(gp_model, y = y, params = list(optimizer_cov = "gradient_descent", lr_cov = 0.01, use_nesterov_acc = TRUE, acc_rate_cov = 0.5, init_cov_pars = init_cov_pars)), file='NUL') cov_pars2 <- c(0.9646422, 0.1844797) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars2)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 26) # Estimation using Nelder-Mead gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit") capture.output( fit(gp_model, y = y, params = list(optimizer_cov = "nelder_mead", delta_rel_conv=1e-6, init_cov_pars = init_cov_pars)) , file='NUL') cov_pars3 <- c(0.9998047, 0.1855072) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars3)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 6) # Estimation using BFGS gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit") capture.output( fit(gp_model, y = y, params = list(optimizer_cov = "bfgs", init_cov_pars = init_cov_pars)), file='NUL') cov_pars3 <- c(0.9419084, 0.1866882) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars3)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 4) # Estimation using Adam gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit") capture.output( fit(gp_model, y = y, params = list(optimizer_cov = "adam", init_cov_pars = init_cov_pars)), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars3)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 200) # Prediction capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", y = y, params = list(optimizer_cov = "gradient_descent", lr_cov=0.01, use_nesterov_acc=FALSE, init_cov_pars = init_cov_pars)) , file='NUL') coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(-0.6595663, -0.6638940, 0.4997690) expected_cov <- c(0.6482224576, 0.5765285950, -0.0001030520, 0.5765285950, 0.6478191338, -0.0001163496, -0.0001030520, -0.0001163496, 0.4435551436) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Predict variances pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$var)-expected_cov[c(1,5,9)])),TOLERANCE_STRICT) # Predict response pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_var=TRUE, predict_response = TRUE) expected_mu <- c(0.3037139, 0.3025143, 0.6612807) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(pred$var-expected_mu*(1-expected_mu))),TOLERANCE_STRICT) # Predict training data random effects training_data_random_effects <- predict_training_data_random_effects(gp_model, predict_var = TRUE) preds <- predict(gp_model, gp_coords_pred = coords, predict_response = FALSE, predict_var = TRUE) expect_lt(sum(abs(training_data_random_effects[,1] - preds$mu)),TOLERANCE_STRICT) expect_lt(sum(abs(training_data_random_effects[,2] - preds$var)),TOLERANCE_STRICT) # Evaluate approximate negative marginal log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.9,0.2),y=y) expect_lt(abs(nll-63.6205917),TOLERANCE_STRICT) # Do optimization using optim and e.g. Nelder-Mead gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit") opt <- optim(par=c(1,0.1), fn=gp_model$neg_log_likelihood, y=y, method="Nelder-Mead") cov_pars <- c(0.9419234, 0.1866877) expect_lt(sum(abs(opt$par-cov_pars)),TOLERANCE_LOOSE) expect_lt(abs(opt$value-(63.6126363)),TOLERANCE_LOOSE) expect_equal(as.integer(opt$counts[1]), 47) ################### ## Random coefficient GPs ################### probs <- pnorm(as.vector(L %*% b_1 + Z_SVC[,1] * L %*% b_2 + Z_SVC[,2] * L %*% b_3)) y <- as.numeric(sim_rand_unif(n=n, init_c=0.543) < probs) init_cov_pars_RC <- rep(init_cov_pars, 3) # Fit model capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", gp_rand_coef_data = Z_SVC, y = y, likelihood = "bernoulli_probit", params = list(optimizer_cov = "gradient_descent", lr_cov = 1, use_nesterov_acc = TRUE, acc_rate_cov=0.5, maxit=1000, init_cov_pars=init_cov_pars_RC)) , file='NUL') expected_values <- c(0.3701097, 0.2846740, 2.1160325, 0.3305266, 0.1241462, 0.1846456) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-expected_values)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 39) # Prediction gp_model <- GPModel(gp_coords = coords, gp_rand_coef_data = Z_SVC, cov_function = "exponential", likelihood = "bernoulli_probit") coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) Z_SVC_test <- cbind(c(0.1,0.3,0.7),c(0.5,0.2,0.4)) pred <- gp_model$predict(y = y, gp_coords_pred = coord_test, gp_rand_coef_data_pred=Z_SVC_test, cov_pars = c(1,0.1,0.8,0.15,1.1,0.08), predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.18346008, 0.03479258, -0.17247579) expected_cov <- c(1.039879e+00, 7.521981e-01, -3.256500e-04, 7.521981e-01, 8.907289e-01, -6.719282e-05, -3.256500e-04, -6.719282e-05, 9.147899e-01) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Evaluate negative log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(1,0.1,0.8,0.15,1.1,0.08),y=y) expect_lt(abs(nll-65.1768199),TOLERANCE_LOOSE) ################### ## Multiple cluster IDs ################### probs <- pnorm(L %*% b_1) y <- as.numeric(sim_rand_unif(n=n, init_c=0.2978341) < probs) capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", y = y, cluster_ids = cluster_ids,likelihood = "bernoulli_probit", params = list(optimizer_cov = "gradient_descent", lr_cov=0.2, use_nesterov_acc=FALSE, init_cov_pars=init_cov_pars)) , file='NUL') cov_pars <- c(0.5085134, 0.2011667) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 20) # Prediction coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) cluster_ids_pred = c(1,3,1) gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", cluster_ids = cluster_ids,likelihood = "bernoulli_probit") pred <- gp_model$predict(y = y, gp_coords_pred = coord_test, cluster_ids_pred = cluster_ids_pred, cov_pars = c(1.5,0.15), predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.1509569, 0.0000000, 0.9574946) expected_cov <- c(1.2225959453, 0.0000000000, 0.0003074858, 0.0000000000, 1.5000000000, 0.0000000000, 0.0003074858, 0.0000000000, 1.0761874845) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) }) test_that("Binary classification with Gaussian process model with multiple observations at the same location", { eps_multiple <- as.vector(L_multiple %*% b_multiple) probs <- pnorm(eps_multiple) y <- as.numeric(sim_rand_unif(n=n, init_c=0.9341) < probs) capture.output( gp_model <- fitGPModel(gp_coords = coords_multiple, cov_function = "exponential", y = y,likelihood = "bernoulli_probit", params = list(optimizer_cov = "gradient_descent", lr_cov=0.1, use_nesterov_acc = TRUE)), file='NUL') cov_pars <- c(0.6857065, 0.2363754) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 8) # Prediction coord_test <- cbind(c(0.1,0.11,0.11),c(0.9,0.91,0.91)) pred <- gp_model$predict(y = y, gp_coords_pred = coord_test, cov_pars = c(1.5,0.15), predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(-0.2633282, -0.2637633, -0.2637633) expected_cov <- c(0.9561355, 0.8535206, 0.8535206, 0.8535206, 1.0180227, 1.0180227, 0.8535206, 1.0180227, 1.0180227) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) pred_resp <- gp_model$predict(y = y, gp_coords_pred = coord_test, cov_pars = c(1.5,0.15), predict_var = TRUE, predict_response = TRUE) expect_lt(sum(abs(pred_resp$mu-c(0.4253296, 0.4263502, 0.4263502))),TOLERANCE_STRICT) expect_lt(sum(abs(pred_resp$var-c(0.2444243, 0.2445757, 0.2445757))),TOLERANCE_STRICT) # Predict training data random effects training_data_random_effects <- predict_training_data_random_effects(gp_model, predict_var = TRUE) preds <- predict(gp_model, gp_coords_pred = coords_multiple, predict_response = FALSE, predict_var = TRUE) expect_lt(sum(abs(training_data_random_effects[,1] - preds$mu)),TOLERANCE_STRICT) expect_lt(sum(abs(training_data_random_effects[,2] - preds$var)),TOLERANCE_STRICT) # Multiple cluster IDs and multiple observations coord_test <- cbind(c(0.1,0.11,0.11),c(0.9,0.91,0.91)) cluster_ids_pred = c(0L,3L,3L) pred <- gp_model$predict(y = y, gp_coords_pred = coord_test, cluster_ids_pred = cluster_ids_pred, cov_pars = c(1.5,0.15), predict_cov_mat = TRUE, predict_response = FALSE) expected_cov <- c(0.9561355, 0.0000000, 0.0000000, 0.0000000, 1.5000000, 1.5000000, 0.0000000, 1.5000000, 1.5000000) expect_lt(sum(abs(pred$mu-c(-0.2633282, rep(0,2)))),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) pred_resp <- gp_model$predict(y = y, gp_coords_pred = coord_test, cluster_ids_pred = cluster_ids_pred, cov_pars = c(1.5,0.15), predict_var = TRUE, predict_response = TRUE) expect_lt(sum(abs(pred_resp$mu-c(0.4253296, 0.5000000, 0.5000000))),TOLERANCE_STRICT) expect_lt(sum(abs(pred_resp$var-c(0.2444243, 0.2500000, 0.2500000))),TOLERANCE_STRICT) }) test_that("Binary classification with one grouped random effects ", { probs <- pnorm(Z1 %*% b_gr_1) y <- as.numeric(sim_rand_unif(n=n, init_c=0.823431) < probs) init_cov_pars <- c(1) # Estimation using gradient descent gp_model <- GPModel(group_data = group, likelihood = "bernoulli_probit") fit(gp_model, y = y, params = list(optimizer_cov = "gradient_descent", lr_cov = 0.1, use_nesterov_acc = FALSE, convergence_criterion = "relative_change_in_parameters", init_cov_pars=init_cov_pars)) cov_pars <- c(0.40255) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 62) # Can switch between likelihoods gp_model <- GPModel(group_data = group, likelihood = "bernoulli_probit") gp_model$set_likelihood("gaussian") gp_model$set_likelihood("bernoulli_probit") fit(gp_model, y = y, params = list(optimizer_cov = "gradient_descent", lr_cov = 0.1, use_nesterov_acc = FALSE, convergence_criterion = "relative_change_in_parameters", init_cov_pars=init_cov_pars)) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) # Estimation using gradient descent and Nesterov acceleration gp_model <- GPModel(group_data = group, likelihood = "bernoulli_probit") fit(gp_model, y = y, params = list(optimizer_cov = "gradient_descent", lr_cov = 0.1, use_nesterov_acc = TRUE, acc_rate_cov = 0.5, init_cov_pars=init_cov_pars)) cov_pars2 <- c(0.4012595) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars2)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 10) # Estimation using gradient descent and too large learning rate gp_model <- GPModel(group_data = group, likelihood = "bernoulli_probit") fit(gp_model, y = y, params = list(optimizer_cov = "gradient_descent", lr_cov = 10, use_nesterov_acc = FALSE, init_cov_pars=init_cov_pars)) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 5) # Prediction gp_model <- fitGPModel(group_data = group, likelihood = "bernoulli_probit", y = y, params = list(optimizer_cov = "gradient_descent", use_nesterov_acc = FALSE, lr_cov = 0.1, init_cov_pars=init_cov_pars)) group_test <- c(1,3,3,9999) pred <- predict(gp_model, y=y, group_data_pred = group_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.000000, -0.796538, -0.796538, 0.000000) expected_cov <- c(0.1133436, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1407783, 0.1407783, 0.0000000, 0.0000000, 0.1407783, 0.1407783, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.4070775) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Predict variances pred <- predict(gp_model, y=y, group_data_pred = group_test, predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$var)-expected_cov[c(1,6,11,16)])),TOLERANCE_STRICT) # Predict response pred <- predict(gp_model, y=y, group_data_pred = group_test, predict_response = TRUE) expected_mu <- c(0.5000000, 0.2279027, 0.2279027, 0.5000000) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) # Prediction for only new groups group_test <- c(-1,-1,-2,-2) pred <- predict(gp_model, y=y, group_data_pred = group_test, predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-rep(0,4))),TOLERANCE_STRICT) expect_lt(sum(abs(pred$var-rep(0,0.4070775))),TOLERANCE_STRICT) pred <- predict(gp_model, y=y, group_data_pred = group_test, predict_response = TRUE) expect_lt(sum(abs(pred$mu-rep(0.5,4))),TOLERANCE_STRICT) # Prediction for only new cluster_ids cluster_ids_pred <- c(-1L,-1L,-2L,-2L) group_test <- c(1,99999,3,3) pred <- predict(gp_model, y=y, group_data_pred = group_test, cluster_ids_pred = cluster_ids_pred, predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-rep(0,4))),TOLERANCE_STRICT) expect_lt(sum(abs(pred$var-rep(0.4070771,4))),TOLERANCE_STRICT) pred <- predict(gp_model, y=y, group_data_pred = group_test, cluster_ids_pred = cluster_ids_pred, predict_response = TRUE) expect_lt(sum(abs(pred$mu-rep(0.5,4))),TOLERANCE_STRICT) # Predict training data random effects all_training_data_random_effects <- predict_training_data_random_effects(gp_model, predict_var = TRUE) first_occurences <- match(unique(group), group) training_data_random_effects <- all_training_data_random_effects[first_occurences,] group_unique <- unique(group) preds <- predict(gp_model, group_data_pred = group_unique, predict_response = FALSE, predict_var = TRUE) expect_lt(sum(abs(training_data_random_effects[,1] - preds$mu)),1E-6) expect_lt(sum(abs(training_data_random_effects[,2] - preds$var)),1E-6) # Estimation using Nelder-Mead gp_model <- GPModel(group_data = group, likelihood = "bernoulli_probit") fit(gp_model, y = y, params = list(optimizer_cov = "nelder_mead", delta_rel_conv=1e-6, init_cov_pars=init_cov_pars)) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-0.4027452)),TOLERANCE_STRICT) # Prediction group_test <- c(1,3,3,9999) pred <- predict(gp_model, y=y, group_data_pred = group_test, predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-c(0.0000000, -0.7935873, -0.7935873, 0.0000000))),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$var)-c(0.1130051, 0.1401125, 0.1401125, 0.4027452))),TOLERANCE_STRICT) # Estimation using BFGS gp_model <- GPModel(group_data = group, likelihood = "bernoulli_probit") fit(gp_model, y = y, params = list(optimizer_cov = "bfgs", init_cov_pars=init_cov_pars)) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-0.4025483)),TOLERANCE_STRICT) # Prediction group_test <- c(1,3,3,9999) pred <- predict(gp_model, y=y, group_data_pred = group_test, predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-c(0.0000000, -0.7934523, -0.7934523, 0.0000000))),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$var)-c(0.1129896, 0.1400821, 0.1400821, 0.4025483))),TOLERANCE_STRICT) # Evaluate approximate negative marginal log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.9),y=y) expect_lt(abs(nll-65.8590638),TOLERANCE_STRICT) # Do optimization using optim gp_model <- GPModel(group_data = group, likelihood = "bernoulli_probit") opt <- optim(par=c(2), fn=gp_model$neg_log_likelihood, y=y, method="Brent", lower=0, upper=1E9) cov_pars <- c(0.40255) expect_lt(sum(abs(opt$par-cov_pars)),TOLERANCE_LOOSE) expect_lt(abs(opt$value-(65.2599674)),TOLERANCE_LOOSE) }) test_that("Binary classification with multiple grouped random effects ", { probs <- pnorm(Z1 %*% b_gr_1 + Z2 %*% b_gr_2 + Z3 %*% b_gr_3) y <- as.numeric(sim_rand_unif(n=n, init_c=0.57341) < probs) init_cov_pars <- rep(1,3) capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, y = y, likelihood = "bernoulli_probit", params = list(optimizer_cov = "gradient_descent", init_cov_pars=init_cov_pars, lr_cov = 0.2, use_nesterov_acc = FALSE, maxit=100)) , file='NUL') expected_values <- c(0.3060671, 0.9328884, 0.3146682) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-expected_values)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 37) # Predict training data random effects cov_pars <- gp_model$get_cov_pars() all_training_data_random_effects <- predict_training_data_random_effects(gp_model, predict_var = TRUE) first_occurences_1 <- match(unique(group), group) first_occurences_2 <- match(unique(group2), group2) pred_random_effects <- all_training_data_random_effects[first_occurences_1,c(1,4)] pred_random_slopes <- all_training_data_random_effects[first_occurences_1,c(3,6)] head(pred_random_slopes) pred_random_effects_crossed <- all_training_data_random_effects[first_occurences_2,c(2,5)] group_unique <- unique(group) group_data_pred = cbind(group_unique,rep(-1,length(group_unique))) x_pr = rep(0,length(group_unique)) preds <- predict(gp_model, group_data_pred=group_data_pred, group_rand_coef_data_pred=x_pr, predict_response = FALSE, predict_var = TRUE) expect_lt(sum(abs(pred_random_effects[,1] - preds$mu)),TOLERANCE_STRICT) expect_lt(sum(abs(pred_random_effects[,2] - (preds$var-cov_pars[2]))),TOLERANCE_STRICT) # Check whether random slopes are correct x_pr = rep(1,length(group_unique)) preds2 <- predict(gp_model, group_data_pred=group_data_pred, group_rand_coef_data_pred=x_pr, predict_response = FALSE) expect_lt(sum(abs(pred_random_slopes[,1] - (preds2$mu-preds$mu))),TOLERANCE_STRICT) # Check whether crossed random effects are correct group_unique <- unique(group2) group_data_pred = cbind(rep(-1,length(group_unique)),group_unique) x_pr = rep(0,length(group_unique)) preds <- predict(gp_model, group_data_pred=group_data_pred, group_rand_coef_data_pred=x_pr, predict_response = FALSE, predict_var = TRUE) expect_lt(sum(abs(pred_random_effects_crossed[,1] - preds$mu)),TOLERANCE_LOOSE) expect_lt(sum(abs(pred_random_effects_crossed[,2] - (preds$var-cov_pars[1]))),TOLERANCE_STRICT) # Prediction gp_model <- GPModel(likelihood = "bernoulli_probit", group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1) group_data_pred = cbind(c(1,1,77),c(2,1,98)) group_rand_coef_data_pred = c(0,0.1,0.3) pred <- gp_model$predict(y = y, group_data_pred=group_data_pred, group_rand_coef_data_pred=group_rand_coef_data_pred, cov_pars = c(0.9,0.8,1.2), predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.5195889, -0.6411954, 0.0000000) expected_cov <- c(0.3422367, 0.1554011, 0.0000000, 0.1554011, 0.3457334, 0.0000000, 0.0000000, 0.0000000, 1.8080000) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Predict variances pred <- gp_model$predict(y = y, group_data_pred=group_data_pred, group_rand_coef_data_pred=group_rand_coef_data_pred, cov_pars = c(0.9,0.8,1.2), predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$var)-expected_cov[c(1,5,9)])),TOLERANCE_STRICT) # Multiple random effects: training with Nelder-Mead capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, y = y, likelihood = "bernoulli_probit", params = list(optimizer_cov = "nelder_mead", delta_rel_conv=1e-6, init_cov_pars=init_cov_pars)) , file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-c(0.3055487, 0.9300562, 0.3048811))),TOLERANCE_STRICT) # Multiple random effects: training with BFGS capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, y = y, likelihood = "bernoulli_probit", params = list(optimizer_cov = "bfgs", init_cov_pars=init_cov_pars)), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-c(0.3030693, 0.9293106, 0.3037503))),TOLERANCE_STRICT) # Evaluate negative log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.9,0.8,1.2),y=y) expect_lt(abs(nll-60.6422359),TOLERANCE_LOOSE) # Multiple cluster_ids capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, y = y, cluster_ids = cluster_ids, likelihood = "bernoulli_probit", params = list(optimizer_cov = "gradient_descent", init_cov_pars=init_cov_pars, lr_cov = 0.2, use_nesterov_acc = FALSE, maxit=100)), file='NUL') expected_values <- c(0.1634433, 0.8952201, 0.3219087) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-expected_values)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 42) # Prediction cluster_ids_pred = c(1,3,1) gp_model <- GPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, cluster_ids = cluster_ids, likelihood = "bernoulli_probit") pred <- gp_model$predict(y = y, group_data_pred=group_data_pred, group_rand_coef_data_pred=group_rand_coef_data_pred, cov_pars = c(0.9,0.8,1.2), cluster_ids_pred = cluster_ids_pred, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(-0.2159939, 0.0000000, 0.0000000) expected_cov <- c(0.4547941, 0.0000000, 0.0000000, 0.0000000, 1.7120000, 0.0000000, 0.0000000, 0.0000000, 1.8080000) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Only one RE and random coefficient probs <- pnorm(Z1 %*% b_gr_1 + Z3 %*% b_gr_3) y <- as.numeric(sim_rand_unif(n=n, init_c=0.957341) < probs) init_cov_pars <- c(1,1) capture.output( gp_model <- fitGPModel(group_data = group, group_rand_coef_data = x, ind_effect_group_rand_coef = 1, y = y, likelihood = "bernoulli_probit", params = list(optimizer_cov = "gradient_descent", init_cov_pars=init_cov_pars, lr_cov = 0.1, use_nesterov_acc = TRUE, maxit=100)) , file='NUL') expected_values <- c(1.00742383, 0.02612587) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-expected_values)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 100) # Random coefficients with intercept random effect dropped probs <- pnorm(Z2 %*% b_gr_2 + Z3 %*% b_gr_3) y <- as.numeric(sim_rand_unif(n=n, init_c=0.8341) < probs) capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, drop_intercept_group_rand_effect = c(TRUE,FALSE), y = y, likelihood = "bernoulli_probit", params = list(optimizer_cov = "gradient_descent", init_cov_pars=init_cov_pars)), file='NUL') expected_values <- c(1.0044712, 0.6549656) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-expected_values)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 18) # Predict training data random effects all_training_data_random_effects <- predict_training_data_random_effects(gp_model) first_occurences_1 <- match(unique(group), group) first_occurences_2 <- match(unique(group2), group2) pred_random_slopes <- all_training_data_random_effects[first_occurences_1,2] pred_random_effects_crossed <- all_training_data_random_effects[first_occurences_2,1] group_unique <- unique(group) group_data_pred = cbind(group_unique,rep(-1,length(group_unique))) # Check whether random slopes are correct x_pr = rep(1,length(group_unique)) preds <- predict(gp_model, group_data_pred=group_data_pred, group_rand_coef_data_pred=x_pr, predict_response = FALSE) expect_lt(sum(abs(pred_random_slopes - preds$mu)),TOLERANCE_LOOSE) # Check whether crossed random effects are correct group_unique <- unique(group2) group_data_pred = cbind(rep(-1,length(group_unique)),group_unique) x_pr = rep(0,length(group_unique)) preds <- predict(gp_model, group_data_pred=group_data_pred, group_rand_coef_data_pred=x_pr, predict_response = FALSE) expect_lt(sum(abs(pred_random_effects_crossed - preds$mu)),TOLERANCE_LOOSE) # Prediction gp_model <- GPModel(likelihood = "bernoulli_probit", group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, drop_intercept_group_rand_effect = c(TRUE,FALSE)) group_data_pred = cbind(c(1,1,77),c(2,1,98)) group_rand_coef_data_pred = c(0,0.1,0.3) pred <- gp_model$predict(y = y, group_data_pred=group_data_pred, group_rand_coef_data_pred=group_rand_coef_data_pred, cov_pars = c(0.8,1.2), predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.8493404, -0.2338359, 0.0000000) expected_cov <- c(0.206019606, -0.001276366, 0.0000000, -0.001276366, 0.155209578, 0.0000000, 0.0000000, 0.0000000, 0.908000000) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Predict variances pred <- gp_model$predict(y = y, group_data_pred=group_data_pred, group_rand_coef_data_pred=group_rand_coef_data_pred, cov_pars = c(0.8,1.2), predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$var)-expected_cov[c(1,5,9)])),TOLERANCE_STRICT) # Including linear fixed effects probs <- pnorm(Z1 %*% b_gr_1 + Z2 %*% b_gr_2 + Z3 %*% b_gr_3 + X%*%beta) y_lin <- as.numeric(sim_rand_unif(n=n, init_c=0.41) < probs) capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, y = y_lin, X=X, likelihood = "bernoulli_probit", params = DEFAULT_OPTIM_PARAMS) , file='NUL') cov_pars <- c(0.7135209, 1.4289386, 1.6037208) coef <- c(-0.382657, 2.413484) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coef)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 19) # Prediction group_data_pred = cbind(c(1,1,77),c(2,1,98)) group_rand_coef_data_pred = c(0,0.1,0.3) X_test <- cbind(rep(1,3),c(-0.5,0.4,1)) pred <- gp_model$predict(group_data_pred=group_data_pred, group_rand_coef_data_pred=group_rand_coef_data_pred, X_pred = X_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(-0.5136463, 0.8644346, 2.0308268) expected_cov <- c(0.5546899, 0.1847860, 0.0000000, 0.1847860, 0.5615866, 0.0000000, 0.0000000, 0.0000000, 2.2867943) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) }) test_that("Binary classification for combined Gaussian process and grouped random effects ", { probs <- pnorm(L %*% b_1 + Z1 %*% b_gr_1) y <- as.numeric(sim_rand_unif(n=n, init_c=0.67341) < probs) init_cov_pars <- c(1,1,mean(dist(coords))/3) # Estimation using gradient descent gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", group_data = group, likelihood = "bernoulli_probit") capture.output( fit(gp_model, y = y, params = list(optimizer_cov = "gradient_descent", init_cov_pars=init_cov_pars, lr_cov = 0.2, use_nesterov_acc = FALSE, convergence_criterion = "relative_change_in_parameters")) , file='NUL') cov_pars <- c(0.3181509, 1.2788456, 0.1218680) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 55) # Prediction capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", group_data = group, y = y, params = list(optimizer_cov = "gradient_descent", init_cov_pars=init_cov_pars, use_nesterov_acc = FALSE, lr_cov = 0.2)) , file='NUL') coord_test <- cbind(c(0.1,0.21,0.7),c(0.9,0.91,0.55)) group_test <- c(1,3,9999) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, group_data_pred = group_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.1217634, -0.9592585, -0.2694489) expected_cov <- c(1.0745455607, 0.2190063794, 0.0040797451, 0.2190063794, 1.0089298170, 0.0000629706, 0.0040797451, 0.0000629706, 1.0449941968) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Predict variances pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, group_data_pred = group_test, predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$var)-expected_cov[c(1,5,9)])),TOLERANCE_STRICT) # Predict response pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, group_data_pred = group_test, predict_response = TRUE) expected_mu <- c(0.5336859, 0.2492699, 0.4252731) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) # Predict training data random effects training_data_random_effects <- predict_training_data_random_effects(gp_model) pred_GP <- predict(gp_model, gp_coords_pred = coords, group_data_pred=rep(-1,dim(coords)[1]), predict_response = FALSE) expect_lt(sum(abs(training_data_random_effects[,2] - pred_GP$mu)),1E-6) # Grouped REs preds <- predict(gp_model, group_data_pred = group, gp_coords_pred = coords, predict_response = FALSE) pred_RE <- preds$mu - pred_GP$mu expect_lt(sum(abs(training_data_random_effects[,1] - pred_RE)),1E-6) # Estimation using Nelder-Mead gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", group_data = group, likelihood = "bernoulli_probit") capture.output( fit(gp_model, y = y, params = list(optimizer_cov = "nelder_mead", delta_rel_conv=1E-8, init_cov_pars=init_cov_pars)), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-c(0.3181320, 1.2795124, 0.1218866))),TOLERANCE_STRICT) # Evaluate approximate negative marginal log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(1.1,0.9,0.2),y=y) expect_lt(abs(nll-65.7219266),TOLERANCE_STRICT) # Do optimization using optim and e.g. Nelder-Mead gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", group_data = group, likelihood = "bernoulli_probit") capture.output( opt <- optim(par=c(1.5,1,0.1), fn=gp_model$neg_log_likelihood, y=y, method="Nelder-Mead"), file='NUL') cov_pars <- c(0.3181509, 1.2788456, 0.1218680) expect_lt(sum(abs(opt$par-cov_pars)),TOLERANCE_LOOSE) expect_lt(abs(opt$value-(63.7432077)),TOLERANCE_LOOSE) expect_equal(as.integer(opt$counts[1]), 164) }) test_that("Combined GP and grouped random effects model with random coefficients ", { probs <- pnorm(as.vector(L %*% b_1 + Z_SVC[,1] * L %*% b_2 + Z_SVC[,2] * L %*% b_3) + Z1 %*% b_gr_1 + Z2 %*% b_gr_2 + Z3 %*% b_gr_3) y <- as.numeric(sim_rand_unif(n=n, init_c=0.9867234) < probs) init_cov_pars <- c(rep(1,3),rep(c(1,mean(dist(coords))/3),3)) # Fit model capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", gp_rand_coef_data = Z_SVC, group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, y = y, likelihood = "bernoulli_probit", params = list(optimizer_cov = "gradient_descent", init_cov_pars=init_cov_pars, lr_cov = 0.2, use_nesterov_acc = FALSE, maxit=10)) , file='NUL') expected_values <- c(0.09859312, 0.35813763, 0.50164573, 0.67372019, 0.08825524, 0.77807532, 0.10896128, 1.03921290, 0.09538707) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-expected_values)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 10) # Prediction gp_model <- GPModel(gp_coords = coords, gp_rand_coef_data = Z_SVC, cov_function = "exponential", likelihood = "bernoulli_probit", group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1) coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) Z_SVC_test <- cbind(c(0.1,0.3,0.7),c(0.5,0.2,0.4)) group_data_pred = cbind(c(1,1,7),c(2,1,3)) group_rand_coef_data_pred = c(0,0.1,0.3) pred <- gp_model$predict(y = y, gp_coords_pred = coord_test, gp_rand_coef_data_pred=Z_SVC_test, group_data_pred=group_data_pred, group_rand_coef_data_pred=group_rand_coef_data_pred, cov_pars = c(0.9,0.8,1.2,1,0.1,0.8,0.15,1.1,0.08), predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(1.612451, 1.147407, -1.227187) expected_cov <- c(1.63468526, 1.02982815, -0.01916993, 1.02982815, 1.43601348, -0.03404720, -0.01916993, -0.03404720, 1.55017397) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Evaluate negative log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.9,0.8,1.2,1,0.1,0.8,0.15,1.1,0.08),y=y) expect_lt(abs(nll-71.4286594),TOLERANCE_LOOSE) }) test_that("Combined GP and grouped random effects model with cluster_id's not constant ", { probs <- pnorm(L %*% b_1 + Z1 %*% b_gr_1) y <- as.numeric(sim_rand_unif(n=n, init_c=0.2341) < probs) init_cov_pars <- c(1,1,mean(dist(coords[cluster_ids==1,]))/3) capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", group_data = group, y = y, cluster_ids = cluster_ids,likelihood = "bernoulli_probit", params = list(optimizer_cov = "gradient_descent", lr_cov=0.2, use_nesterov_acc = FALSE, init_cov_pars=init_cov_pars)) , file='NUL') cov_pars <- c(0.276476226, 0.007278016, 0.132195703) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 261) # Prediction coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) group_data_pred = c(1,1,9999) cluster_ids_pred = c(1,3,1) gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", group_data = group, cluster_ids = cluster_ids,likelihood = "bernoulli_probit") pred <- gp_model$predict(y = y, gp_coords_pred = coord_test, group_data_pred = group_data_pred, cluster_ids_pred = cluster_ids_pred, cov_pars = c(1.5,1,0.15), predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.1074035, 0.0000000, 0.2945508) expected_cov <- c(0.98609786, 0.00000000, -0.02013244, 0.00000000, 2.50000000, 0.00000000, -0.02013244, 0.00000000, 2.28927616) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) }) test_that("Binary classification Gaussian process model with Vecchia approximation", { probs <- pnorm(L %*% b_1) y <- as.numeric(sim_rand_unif(n=n, init_c=0.2341) < probs) init_cov_pars <- c(1,mean(dist(coords))/3) params_vecchia <- list(optimizer_cov = "gradient_descent", init_cov_pars=init_cov_pars, lr_cov = 0.1, use_nesterov_acc = FALSE, convergence_criterion = "relative_change_in_parameters", cg_delta_conv = sqrt(1e-6), num_rand_vec_trace = 500, cg_preconditioner_type = "piv_chol_on_Sigma") for(inv_method in c("cholesky", "iterative")){ if(inv_method == "iterative"){ tolerance_loc_1 <- TOLERANCE_ITERATIVE tolerance_loc_2 <- TOLERANCE_ITERATIVE } else{ tolerance_loc_1 <- TOLERANCE_STRICT tolerance_loc_2 <- TOLERANCE_LOOSE } # Estimation using gradient descent capture.output( gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "vecchia", num_neighbors = n-1, vecchia_ordering = "none", matrix_inversion_method = inv_method), file='NUL') capture.output( fit(gp_model, y = y, params = params_vecchia) , file='NUL') cov_pars <- c(0.9419234, 0.1866877) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),tolerance_loc_1) if(inv_method != "iterative") { expect_equal(gp_model$get_num_optim_iter(), 40) } # Vecchia approximation with random ordering capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", vecchia_ordering="random", likelihood = "bernoulli_probit", gp_approx = "vecchia", num_neighbors = n-1, y = y, params = params_vecchia, matrix_inversion_method = inv_method), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),tolerance_loc_1) if(inv_method != "iterative") { expect_equal(gp_model$get_num_optim_iter(), 40) } # init_cov_pars_ not given params_vecchia_no_init <- params_vecchia params_vecchia_no_init$init_cov_pars <- NULL capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", vecchia_ordering="random", likelihood = "bernoulli_probit", gp_approx = "vecchia", num_neighbors = n-1, y = y, params = params_vecchia_no_init, matrix_inversion_method = inv_method), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),tolerance_loc_1) if(inv_method != "iterative") { expect_equal(gp_model$get_num_optim_iter(), 40) } if(inv_method != "iterative"){ # Same estimation without Vecchia approximation capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "none", y = y, params = params_vecchia), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),tolerance_loc_1) expect_equal(gp_model$get_num_optim_iter(), 40) } # Prediction capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "vecchia", num_neighbors = n-1, vecchia_ordering = "none", y = y, matrix_inversion_method = inv_method, params = list(optimizer_cov = "gradient_descent", use_nesterov_acc = FALSE, lr_cov = 0.01, init_cov_pars = init_cov_pars, cg_preconditioner_type = "piv_chol_on_Sigma")), file='NUL') coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) gp_model$set_prediction_data(vecchia_pred_type = "latent_order_obs_first_cond_all", num_neighbors_pred = n+2, nsim_var_pred = 100000) capture.output( pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_cov_mat = TRUE, predict_response = FALSE), file='NUL') expected_mu <- c(-0.6595662, -0.6638940, 0.4997690) expected_cov <- c(0.6482223800, 0.5765285294, -0.0001030520, 0.5765285294, 0.6478190560, -0.0001163496, -0.0001030520, -0.0001163496, 0.4435550921) expect_lt(sum(abs(pred$mu-expected_mu)),tolerance_loc_1) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),tolerance_loc_1) # Predict variances capture.output( pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_var = TRUE, predict_response = FALSE), file='NUL') expect_lt(sum(abs(pred$mu-expected_mu)),tolerance_loc_1) expect_lt(sum(abs(as.vector(pred$var)-expected_cov[c(1,5,9)])),tolerance_loc_1) # Predict response capture.output( pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_response = TRUE, predict_var = TRUE), file='NUL') expected_mu_resp <- c(0.3037139, 0.3025143, 0.6612807) expected_var_resp <- c(0.2114718, 0.2109994, 0.2239885) expect_lt(sum(abs(pred$mu-expected_mu_resp)),tolerance_loc_1) expect_lt(sum(abs(pred$var-expected_var_resp)),tolerance_loc_1) if(inv_method != "iterative"){ # Same prediction without Vecchia approximation capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "none", y = y, params = list(optimizer_cov = "gradient_descent", use_nesterov_acc = FALSE, lr_cov=0.01, init_cov_pars=init_cov_pars)), file='NUL') pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_cov_mat = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),tolerance_loc_1) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),tolerance_loc_1) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),tolerance_loc_1) expect_lt(sum(abs(as.vector(pred$var)-expected_cov[c(1,5,9)])),tolerance_loc_1) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_response = TRUE, predict_var = TRUE) expect_lt(sum(abs(pred$mu-expected_mu_resp)),tolerance_loc_1) expect_lt(sum(abs(pred$var-expected_var_resp)),tolerance_loc_1) # Evaluate approximate negative marginal log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.9,0.2),y=y) expect_lt(abs(nll-63.6205917),TOLERANCE_STRICT) } ####################### ## Less neighbors than observations ####################### capture.output( gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "vecchia", num_neighbors=30, vecchia_ordering = "none", matrix_inversion_method = inv_method), file='NUL') capture.output( fit(gp_model, y = y, params = list(optimizer_cov = "gradient_descent", init_cov_pars=init_cov_pars, lr_cov = 0.1, use_nesterov_acc = FALSE, convergence_criterion = "relative_change_in_parameters", cg_delta_conv = sqrt(1e-6), num_rand_vec_trace = 500, cg_preconditioner_type = "piv_chol_on_Sigma")), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),tolerance_loc_2) if(inv_method == "cholesky"){ expect_equal(gp_model$get_num_optim_iter(), 40) } else { expect_equal(gp_model$get_num_optim_iter(), 5) } if(inv_method == "iterative"){ ## Cannot change cg_preconditioner_type after a model has been fitted expect_error( capture.output( fit(gp_model, y = y, params = list(optimizer_cov = "gradient_descent", init_cov_pars=init_cov_pars, lr_cov = 0.1, use_nesterov_acc = FALSE, convergence_criterion = "relative_change_in_parameters", cg_delta_conv = 1e-6, num_rand_vec_trace = 500, cg_preconditioner_type = "Sigma_inv_plus_BtWB")), file='NUL')) capture.output( gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "vecchia", num_neighbors=30, vecchia_ordering = "none", matrix_inversion_method = inv_method), file='NUL') capture.output( fit(gp_model, y = y, params = list(optimizer_cov = "gradient_descent", init_cov_pars=init_cov_pars, lr_cov = 0.1, use_nesterov_acc = FALSE, convergence_criterion = "relative_change_in_parameters", cg_delta_conv = 1e-6, num_rand_vec_trace = 500, cg_preconditioner_type = "Sigma_inv_plus_BtWB")), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),tolerance_loc_2) } # Random ordering capture.output( gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", vecchia_ordering="random", likelihood = "bernoulli_probit", gp_approx = "vecchia", num_neighbors=30, matrix_inversion_method = inv_method), file='NUL') capture.output( fit(gp_model, y = y, params = list(optimizer_cov = "gradient_descent", init_cov_pars=init_cov_pars, lr_cov = 0.1, use_nesterov_acc = FALSE, convergence_criterion = "relative_change_in_parameters", cg_delta_conv = 1e-6, num_rand_vec_trace = 500)) , file='NUL') if(inv_method != "iterative"){ expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),0.01) # Note (02.03.2023): the results differ among compilers (e.g., gcc) # as pseudo random numbers from the C++ RNG 'mt19937' can differ across compilers expect_equal(gp_model$get_num_optim_iter(), 40) } else{ expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),tolerance_loc_1) } # Prediction capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "vecchia", num_neighbors=30, vecchia_ordering = "none", matrix_inversion_method = inv_method, y = y, params = list(optimizer_cov = "gradient_descent", use_nesterov_acc = FALSE, lr_cov=0.01, init_cov_pars=init_cov_pars, cg_delta_conv = 1e-6, num_rand_vec_trace = 500)), file='NUL') gp_model$set_prediction_data(vecchia_pred_type = "latent_order_obs_first_cond_all", num_neighbors_pred = 30, nsim_var_pred=5000) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(-0.6602899, -0.6646044, 0.5004378) expected_cov <- c(0.6476328648, 0.5760722706, -0.0001037986, 0.5760722706, 0.6472246191, -0.0001150096, -0.0001037986, -0.0001150096, 0.4430490207) expect_lt(sum(abs(pred$mu-expected_mu)),tolerance_loc_1) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),tolerance_loc_1) # Predict variances pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),tolerance_loc_1) expect_lt(sum(abs(as.vector(pred$var)-expected_cov[c(1,5,9)])),tolerance_loc_1) # Use vecchia_pred_type = "order_obs_first_cond_all" gp_model$set_prediction_data(vecchia_pred_type = "order_obs_first_cond_all") pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_cov_mat = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),tolerance_loc_1) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),tolerance_loc_1) # Predict response pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_response = TRUE, predict_var = TRUE) expected_mu_lat <- c(0.3034847, 0.3022886, 0.6615111) expected_var_lat <- c(0.2113818, 0.2109102, 0.2239142) expect_lt(sum(abs(pred$mu-expected_mu_lat)),tolerance_loc_1) expect_lt(sum(abs(pred$var-expected_var_lat)), tolerance_loc_1) # Use vecchia_pred_type = "latent_order_obs_first_cond_obs_only" gp_model$set_prediction_data(vecchia_pred_type = "latent_order_obs_first_cond_obs_only", nsim_var_pred=2000) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_cov <- c(0.6476328648, 0.2954937370, -0.0001037986, 0.2954937370, 0.6472245580, -0.0001133886, -0.0001037986, -0.0001133886, 0.4430490207) expect_lt(sum(abs(pred$mu-expected_mu)),tolerance_loc_2) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),tolerance_loc_1) # Use vecchia_pred_type = "order_obs_first_cond_obs_only" gp_model$set_prediction_data(vecchia_pred_type = "order_obs_first_cond_obs_only") pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_cov_mat = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),tolerance_loc_2) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),tolerance_loc_1) # Predict training data random effects training_data_random_effects <- predict_training_data_random_effects(gp_model, predict_var = TRUE) gp_model$set_prediction_data(vecchia_pred_type = "order_obs_first_cond_obs_only") preds <- predict(gp_model, gp_coords_pred = coords, predict_response = FALSE, predict_var = TRUE) expect_lt(sum(abs(training_data_random_effects[,1] - preds$mu)),tolerance_loc_1) if(inv_method == "iterative"){ expect_lt(mean(abs(training_data_random_effects[,2] - preds$var)), TOLERANCE_ITERATIVE) #Different RNG-Status }else{ expect_lt(sum(abs(training_data_random_effects[,2] - preds$var)),TOLERANCE_STRICT) } # Evaluate approximate negative marginal log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.9,0.2),y=y) expect_lt(abs(nll-63.6211118),tolerance_loc_1) } ################### ## Random coefficient GPs ################### probs <- pnorm(as.vector(L %*% b_1 + Z_SVC[,1] * L %*% b_2 + Z_SVC[,2] * L %*% b_3)) y <- as.numeric(sim_rand_unif(n=n, init_c=0.543) < probs) init_cov_pars_RC <- rep(init_cov_pars, 3) # Estimation capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", gp_rand_coef_data = Z_SVC, y = y, likelihood = "bernoulli_probit", gp_approx = "vecchia", num_neighbors = n-1, vecchia_ordering = "none", params = list(optimizer_cov = "gradient_descent", lr_cov = 1, use_nesterov_acc = TRUE, acc_rate_cov=0.5, maxit=1000, init_cov_pars=init_cov_pars_RC)), file='NUL') expected_values <- c(0.3701097, 0.2846740, 2.1160323, 0.3305266, 0.1241462, 0.1846456) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-expected_values)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 39) # Same estimation without Vecchia approximation capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", gp_rand_coef_data = Z_SVC, y = y, likelihood = "bernoulli_probit", gp_approx = "none", params = list(optimizer_cov = "gradient_descent", lr_cov = 1, use_nesterov_acc = TRUE, acc_rate_cov=0.5, maxit=1000, init_cov_pars=init_cov_pars_RC)), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-expected_values)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 39) # Prediction capture.output( gp_model <- GPModel(gp_coords = coords, gp_rand_coef_data = Z_SVC, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "vecchia", num_neighbors = n-1, vecchia_ordering = "none"), file='NUL') coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) Z_SVC_test <- cbind(c(0.1,0.3,0.7),c(0.5,0.2,0.4)) gp_model$set_prediction_data(vecchia_pred_type = "latent_order_obs_first_cond_all", num_neighbors_pred=n+2) pred <- gp_model$predict(y = y, gp_coords_pred = coord_test, gp_rand_coef_data_pred=Z_SVC_test, cov_pars = c(1,0.1,0.8,0.15,1.1,0.08), predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.18346009, 0.03479259, -0.17247579) expected_cov <- c(1.039879e+00, 7.521981e-01, -3.256500e-04, 7.521981e-01, 8.907289e-01, -6.719282e-05, -3.256500e-04, -6.719282e-05, 9.147899e-01) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Same prediction without Veccchia approximation capture.output( gp_model <- GPModel(gp_coords = coords, gp_rand_coef_data = Z_SVC, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "none"), file='NUL') pred <- gp_model$predict(y = y, gp_coords_pred = coord_test, gp_rand_coef_data_pred=Z_SVC_test, cov_pars = c(1,0.1,0.8,0.15,1.1,0.08), predict_cov_mat = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Evaluate negative log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(1,0.1,0.8,0.15,1.1,0.08),y=y) expect_lt(abs(nll-65.1768199),TOLERANCE_LOOSE) ################### ## Multiple cluster IDs ################### probs <- pnorm(L %*% b_1) y <- as.numeric(sim_rand_unif(n=n, init_c=0.2978341) < probs) init_cov_pars <- c(1,mean(dist(coords[cluster_ids==1,]))/3) capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", y = y, cluster_ids = cluster_ids, likelihood = "bernoulli_probit", gp_approx = "vecchia", num_neighbors = n-1, vecchia_ordering = "none", params = list(optimizer_cov = "gradient_descent", lr_cov=0.2, use_nesterov_acc = FALSE, init_cov_pars=init_cov_pars)), file='NUL') cov_pars <- c(0.5085134, 0.2011667) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 20) # Prediction coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) cluster_ids_pred = c(1,3,1) capture.output( gp_model <- GPModel(gp_coords = coords, cov_function = "exponential", cluster_ids = cluster_ids,likelihood = "bernoulli_probit"), file='NUL') pred <- gp_model$predict(y = y, gp_coords_pred = coord_test, cluster_ids_pred = cluster_ids_pred, cov_pars = c(1.5,0.15), predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.1509569, 0.0000000, 0.9574946) expected_cov <- c(1.2225959453, 0.0000000000, 0.0003074858, 0.0000000000, 1.5000000000, 0.0000000000, 0.0003074858, 0.0000000000, 1.0761874845) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) ## Cannot have duplicate coordinates expect_error( capture.output( gp_model <- GPModel(gp_coords = coords_multiple, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "vecchia", vecchia_ordering = "none"), file='NUL') ) }) test_that("Binary classification Gaussian process model with Wendland covariance function", { probs <- pnorm(L %*% b_1) y <- as.numeric(sim_rand_unif(n=n, init_c=0.2341) < probs) init_cov_pars <- c(mean(dist(coords))/3) # Estimation using gradient descent capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "wendland", cov_fct_taper_shape = 0, cov_fct_taper_range = 0.1, y = y, likelihood = "bernoulli_probit", params = list(optimizer_cov = "gradient_descent", std_dev = FALSE, lr_cov = 0.1, use_nesterov_acc = TRUE, acc_rate_cov = 0.5, init_cov_pars=init_cov_pars)), file='NUL') cov_pars <- c(0.5553221) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 33) # Prediction coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(-0.05440076, -0.05767809, 0.05060592) expected_cov <- c(0.5539199, 0.4080647, 0.0000000, 0.4080647, 0.5533222, 0.0000000, 0.0000000, 0.0000000, 0.5146614) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Predict variances pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$var)-expected_cov[c(1,5,9)])),TOLERANCE_STRICT) # Predict response pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_response = TRUE) expected_mu <- c(0.4825954, 0.4815441, 0.5163995) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) }) test_that("Binary classification with linear predictor and grouped random effects model ", { probs <- pnorm(Z1 %*% b_gr_1 + X%*%beta) y <- as.numeric(sim_rand_unif(n=n, init_c=0.542) < probs) init_cov_pars = c(1) # Estimation using gradient descent and Nesterov acceleration gp_model <- fitGPModel(group_data = group, likelihood = "bernoulli_probit", y = y, X=X, params = list(optimizer_cov = "gradient_descent", optimizer_coef = "gradient_descent", lr_cov = 0.05, lr_coef = 1, use_nesterov_acc = TRUE, acc_rate_cov = 0.2, acc_rate_coef = 0.1, init_cov_pars=init_cov_pars)) cov_pars <- c(0.4072025) coef <- c(-0.1113238, 1.5178339) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coef)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 43) # Estimation using Nelder-Mead gp_model <- fitGPModel(group_data = group, likelihood = "bernoulli_probit", y = y, X=X, params = list(optimizer_cov = "nelder_mead", optimizer_coef = "nelder_mead", delta_rel_conv=1e-12, init_cov_pars=init_cov_pars)) cov_pars <- c(0.399973) coef <- c(-0.1109516, 1.5149596) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coef)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 188) # init_cov_pars not given gp_model <- fitGPModel(group_data = group, likelihood = "bernoulli_probit", y = y, X=X, params = list(optimizer_cov = "nelder_mead", optimizer_coef = "nelder_mead", delta_rel_conv=1e-12)) cov_pars <- c(0.399973) coef <- c(-0.1109516, 1.5149596) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coef)),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 194) # # Estimation using BFGS # Does not converge to correct solution (version 0.7.2, 21.02.2022) # gp_model <- fitGPModel(group_data = group, likelihood = "bernoulli_probit", # y = y, X=X, params = list(optimizer_cov = "bfgs", # optimizer_coef = "bfgs")) # cov_pars <- c(0.3999729) # coef <- c(-0.1109532, 1.5149592) # expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) # expect_lt(sum(abs(as.vector(gp_model$get_coef())-coef)),TOLERANCE_LOOSE) # expect_equal(gp_model$get_num_optim_iter(), 17) # Prediction gp_model <- fitGPModel(group_data = group, likelihood = "bernoulli_probit", y = y, X=X, params = list(optimizer_cov = "gradient_descent", optimizer_coef = "gradient_descent", use_nesterov_acc=FALSE, lr_coef=1, init_cov_pars=init_cov_pars)) X_test <- cbind(rep(1,4),c(-0.5,0.2,0.4,1)) group_test <- c(1,3,3,9999) pred <- predict(gp_model, y=y, group_data_pred = group_test, X_pred = X_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(-0.81132150, -0.08574588, 0.21768684, 1.40591430) expected_cov <- c(0.1380238, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1688248, 0.1688248, 0.0000000, 0.0000000, 0.1688248, 0.1688248, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.4051185) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Predict response pred <- predict(gp_model, y=y, group_data_pred = group_test, X_pred = X_test, predict_response = TRUE) expected_mu <- c(0.2234684, 0.4683923, 0.5797886, 0.8821984) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) # Predict training data random effects all_training_data_random_effects <- predict_training_data_random_effects(gp_model, predict_var = TRUE) first_occurences <- match(unique(group), group) training_data_random_effects <- all_training_data_random_effects[first_occurences,] group_unique <- unique(group) X_zero <- cbind(rep(0,length(group_unique)),rep(0,length(group_unique))) preds <- predict(gp_model, group_data_pred = group_unique, X_pred = X_zero, predict_response = FALSE, predict_var = TRUE) expect_lt(sum(abs(training_data_random_effects[,1] - preds$mu)),TOLERANCE_STRICT) expect_lt(sum(abs(training_data_random_effects[,2] - preds$var)),TOLERANCE_STRICT) # Standard deviations capture.output( gp_model <- fitGPModel(group_data = group, likelihood = "bernoulli_probit", y = y, X=X, params = list(std_dev = TRUE, optimizer_cov = "gradient_descent", optimizer_coef = "gradient_descent", init_cov_pars=init_cov_pars, use_nesterov_acc = TRUE, lr_cov = 0.1, lr_coef = 1)), file='NUL') cov_pars <- c(0.4006247) coef <- c(-0.1112284, 0.2566224, 1.5160319, 0.2636918) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coef)),TOLERANCE_STRICT) # Providing initial covariance parameters and coefficients cov_pars <- c(1) coef <- c(2,5) gp_model <- fitGPModel(group_data = group, likelihood = "bernoulli_probit", y = y, X=X, params = list(maxit=0, init_cov_pars=cov_pars, init_coef=coef)) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coef)),TOLERANCE_LOOSE) # Large data n_L <- 1e6 # number of samples m_L <- n_L/10 # number of categories / levels for grouping variable group_L <- rep(1,n_L) # grouping variable for(i in 1:m_L) group_L[((i-1)*n_L/m_L+1):(i*n_L/m_L)] <- i keps <- 1E-10 b1_L <- qnorm(sim_rand_unif(n=m_L, init_c=0.671)*(1-keps) + keps/2) X_L <- cbind(rep(1,n_L),sim_rand_unif(n=n_L, init_c=0.8671)-0.5) # design matrix / covariate data for fixed effect probs_L <- pnorm(b1_L[group_L] + X_L%*%beta) y_L <- as.numeric(sim_rand_unif(n=n_L, init_c=0.12378)*(1-keps) + keps/2 < probs_L) # Estimation using gradient descent and Nesterov acceleration gp_model <- fitGPModel(group_data = group_L, likelihood = "bernoulli_probit", y = y_L, X=X_L, params = list(optimizer_cov = "gradient_descent", optimizer_coef = "gradient_descent", lr_cov = 0.05, lr_coef = 0.1, use_nesterov_acc = TRUE, init_cov_pars=init_cov_pars)) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-0.9759329)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(gp_model$get_coef())-c(0.0971936, 1.9950664))),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 6) # Estimation using Nelder-Mead gp_model <- fitGPModel(group_data = group_L, likelihood = "bernoulli_probit", y = y_L, X=X_L, params = list(optimizer_cov = "nelder_mead", optimizer_coef = "nelder_mead", delta_rel_conv=1e-6, init_cov_pars=init_cov_pars)) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-0.9712287)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(gp_model$get_coef())-c(0.09758098, 1.99473498))),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 109) }) test_that("Binary classification with linear predictor and Gaussian process model ", { probs <- pnorm(L %*% b_1 + X%*%beta) y <- as.numeric(sim_rand_unif(n=n, init_c=0.199) < probs) # Estimation gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", y = y, X=X, params = DEFAULT_OPTIM_PARAMS) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-c(1.3992407, 0.261976))),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(gp_model$get_coef())-c(0.2764603, 1.5556477))),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 11) # Prediction coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) X_test <- cbind(rep(1,3),c(-0.5,0.2,1)) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, X_pred = X_test, predict_var = TRUE, predict_response = FALSE) expected_mu <- c(-0.7480014, 0.3297389, 2.7039005) expected_var <- c(0.8596074, 0.8574038, 0.5016189) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(pred$var)-expected_var)),TOLERANCE_LOOSE) # Estimation using Nelder-Mead gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", y = y, X=X, params = list(optimizer_cov = "nelder_mead", optimizer_coef = "nelder_mead", maxit=1000, delta_rel_conv=1e-12)) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-c(1.2717516, 0.2875537))),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(gp_model$get_coef())-c(0.1999365, 1.4666199))),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 231) # Standard deviations capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", y = y, X=X, params = DEFAULT_OPTIM_PARAMS_STD), file='NUL') coef <- c(0.2764603, 0.5420554, 1.5556477, 0.3146670) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coef)),TOLERANCE_STRICT) }) test_that("Tapering for binary classification", { probs <- pnorm(L %*% b_1 + X%*%beta) y <- as.numeric(sim_rand_unif(n=n, init_c=0.199) < probs) # No tapering gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", y = y, X=X, params = DEFAULT_OPTIM_PARAMS) cov_pars <- c(1.3992407, 0.261976) coefs <- c(0.2764603, 1.5556477) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coefs)),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 11) # Prediction coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) X_test <- cbind(rep(1,3),c(-0.5,0.2,1)) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, X_pred = X_test, predict_var = TRUE, predict_response = FALSE) expected_mu <- c(-0.7480014, 0.3297389, 2.7039005) expected_var <- c(0.8596074, 0.8574038, 0.5016189) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(pred$var)-expected_var)),TOLERANCE_LOOSE) # With tapering and very large tapering range capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "tapering", cov_fct_taper_shape = 0, cov_fct_taper_range = 1e6, y = y, X=X, params = DEFAULT_OPTIM_PARAMS), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coefs)),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 11) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, X_pred = X_test, predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(pred$var)-expected_var)),TOLERANCE_LOOSE) # With tapering and small tapering range capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "matern", likelihood = "bernoulli_probit", cov_fct_shape = 2.5, gp_approx = "tapering", cov_fct_taper_shape = 1, cov_fct_taper_range = 0.5, y = y, X=X, params = DEFAULT_OPTIM_PARAMS), file='NUL') cov_pars <- c(0.9397216, 0.1193397) coefs <- c(0.4465991, 1.4398973) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coefs)),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 16) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, X_pred = X_test, predict_var = TRUE, predict_response = FALSE) expected_mu <- c(-0.3887768, 0.6451925, 2.5398402) expected_var <- c(0.7738142, 0.7868789, 0.4606071) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(pred$var)-expected_var)),TOLERANCE_LOOSE) # Multiple observations at the same location eps_multiple <- as.vector(L_multiple %*% b_multiple) probs <- pnorm(eps_multiple + X%*%beta) y <- as.numeric(sim_rand_unif(n=n, init_c=0.41) < probs) #No tapering capture.output( gp_model <- fitGPModel(gp_coords = coords_multiple, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "none", cov_fct_taper_shape = 0, cov_fct_taper_range = 1e6, y = y, X=X, params = DEFAULT_OPTIM_PARAMS), file='NUL') cov_pars <- c(1.09933629, 0.08239163) coefs <- c(0.5652697, 1.7696475) num_it <- 15 expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coefs)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), num_it) # Prediction coord_test <- cbind(c(0.1,0.11,0.11),c(0.9,0.91,0.91)) X_test <- cbind(rep(1,3),c(-0.5,0.2,1)) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, X_pred = X_test, predict_var = TRUE, predict_response = FALSE) expected_mu <- c(-0.4170239, 0.8119773, 2.2276953) expected_var <- c(0.9473460, 0.9779394, 0.9779394) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(pred$var)-expected_var)),TOLERANCE_LOOSE) # With tapering and very large tapering range capture.output( gp_model <- fitGPModel(gp_coords = coords_multiple, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "tapering", cov_fct_taper_shape = 0, cov_fct_taper_range = 1e6, y = y, X=X, params = DEFAULT_OPTIM_PARAMS), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coefs)),1e-2) expect_equal(gp_model$get_num_optim_iter(), num_it) # Prediction pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, X_pred = X_test, predict_var = TRUE, predict_response = FALSE) expect_lt(sum(abs(pred$mu-expected_mu)),1e-2) expect_lt(sum(abs(as.vector(pred$var)-expected_var)),1e-2) # With tapering and smalltapering range capture.output( gp_model <- fitGPModel(gp_coords = coords_multiple, cov_function = "exponential", likelihood = "bernoulli_probit", gp_approx = "tapering", cov_fct_taper_shape = 0, cov_fct_taper_range = 0.5, y = y, X=X, params = DEFAULT_OPTIM_PARAMS), file='NUL') cov_pars <- c(1.1062329, 0.1479239) coefs <- c(0.5562694, 1.7715742) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coefs)),1e-2) expect_equal(gp_model$get_num_optim_iter(), 48) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, X_pred = X_test, predict_var = TRUE, predict_response = FALSE) expected_mu <- c(-0.4272906, 0.8014820, 2.2187414) expected_var <- c(0.9295132, 0.9637862, 0.9637862) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(pred$var)-expected_var)),TOLERANCE_LOOSE) }) test_that("Binary classification with Gaussian process model and logit link function", { probs <- 1/(1+exp(- L %*% b_1)) y <- as.numeric(sim_rand_unif(n=n, init_c=0.2341) < probs) # Estimation capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "bernoulli_logit", y = y, params = list(optimizer_cov = "gradient_descent", use_nesterov_acc = TRUE, lr_cov=0.01)) , file='NUL') cov_pars <- c(1.4300136, 0.1891952) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 85) # Prediction coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(-0.7792960, -0.7876208, 0.5476390) expected_cov <- c(1.024267e+00, 9.215206e-01, 5.561435e-05, 9.215206e-01, 1.022897e+00, 2.028618e-05, 5.561435e-05, 2.028618e-05, 7.395747e-01) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Predict response pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_var=TRUE, predict_response = TRUE) expected_mu <- c(0.3442815, 0.3426873, 0.6159933) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(pred$var-expected_mu*(1-expected_mu))),TOLERANCE_STRICT) # Evaluate approximate negative marginal log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.9,0.2),y=y) expect_lt(abs(nll-66.299571),TOLERANCE_STRICT) }) test_that("Poisson regression ", { # Single level grouped random effects mu <- exp(Z1 %*% b_gr_1) y <- qpois(sim_rand_unif(n=n, init_c=0.04532), lambda = mu) # Estimation capture.output( gp_model <- fitGPModel(group_data = group, likelihood = "poisson", y = y, params = list(optimizer_cov = "gradient_descent", use_nesterov_acc = TRUE, lr_cov=0.1)) , file='NUL') cov_pars <- c(0.4033406) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 8) # Prediction group_test <- c(1,3,3,9999) pred <- predict(gp_model, y=y, group_data_pred = group_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.07765297, -0.87488533, -0.87488533, 0.00000000) expected_cov <- c(0.07526284, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.15041230, 0.15041230, 0.00000000, 0.00000000, 0.15041230, 0.15041230, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.40334058) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Predict response pred <- predict(gp_model, y=y, group_data_pred = group_test, predict_var=TRUE, predict_response = TRUE) expected_mu <- c(1.1221925, 0.4494731, 0.4494731, 1.2234446) expected_var <- c(1.2206301, 0.4822647, 0.4822647, 1.9670879) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(pred$var-expected_var)),TOLERANCE_STRICT) # Evaluate negative log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.9),y=y) expect_lt(abs(nll-140.4554806),TOLERANCE_LOOSE) # Multiple random effects mu <- exp(Z1 %*% b_gr_1 + Z2 %*% b_gr_2 + Z3 %*% b_gr_3) y <- qpois(sim_rand_unif(n=n, init_c=0.74532), lambda = mu) # Estimation capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, likelihood = "poisson", y = y, params = list(optimizer_cov = "gradient_descent", use_nesterov_acc = TRUE, lr_cov=0.1)) , file='NUL') cov_pars <- c(0.4069344, 1.6988978, 1.3415016) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 7) # Prediction group_data_pred = cbind(c(1,1,77),c(2,1,98)) group_rand_coef_data_pred = c(0,0.1,0.3) pred <- gp_model$predict(y = y, group_data_pred=group_data_pred, group_rand_coef_data_pred=group_rand_coef_data_pred, cov_pars = c(0.9,0.8,1.2), predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.92620057, -0.08200469, 0.00000000) expected_cov <- c(0.07730896, 0.04403442, 0.00000000, 0.04403442, 0.11600469, 0.00000000, 0.00000000, 0.00000000, 1.80800000) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Gaussian process model mu <- exp(L %*% b_1) y <- qpois(sim_rand_unif(n=n, init_c=0.435), lambda = mu) # Estimation capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "poisson", y = y, params = list(optimizer_cov = "gradient_descent", use_nesterov_acc = TRUE)) , file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-c(1.1853922, 0.1500197))),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 6) # Prediction coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.4329068, 0.4042531, 0.6833738) expected_cov <- c(6.550626e-01, 5.553938e-01, -8.406290e-06, 5.553938e-01, 6.631295e-01, -7.658261e-06, -8.406290e-06, -7.658261e-06, 4.170417e-01) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Predict response pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_var=TRUE, predict_response = TRUE) expected_mu <- c(2.139213, 2.087188, 2.439748) expected_var <- c(6.373433, 6.185895, 5.519896) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(pred$var-expected_var)),TOLERANCE_LOOSE) # Evaluate approximate negative marginal log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.9,0.2),y=y) expect_lt(abs(nll-195.03708036),TOLERANCE_STRICT) ## Grouped random effects model with a linear predictor mu_lin <- exp(Z1 %*% b_gr_1 + X%*%beta) y_lin <- qpois(sim_rand_unif(n=n, init_c=0.84532), lambda = mu_lin) gp_model <- fitGPModel(group_data = group, likelihood = "poisson", y = y_lin, X=X, params = list(optimizer_cov = "gradient_descent", optimizer_coef = "gradient_descent", lr_cov = 0.1, lr_coef = 0.1, use_nesterov_acc = TRUE, acc_rate_cov = 0.5)) cov_pars <- c(0.2993371) coef <- c(-0.1526089, 2.1267601) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_LOOSE) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coef)),TOLERANCE_LOOSE) expect_equal(gp_model$get_num_optim_iter(), 29) }) test_that("Gamma regression ", { params <- list(optimizer_cov = "gradient_descent", optimizer_coef = "gradient_descent", estimate_aux_pars = FALSE, init_aux_pars = 1., lr_cov = 0.1, lr_coef = 0.1, use_nesterov_acc = TRUE, acc_rate_cov = 0.5) params_shape <- params params_shape$estimate_aux_pars <- TRUE shape <- 1 # Single level grouped random effects mu <- exp(Z1 %*% b_gr_1) y <- qgamma(sim_rand_unif(n=n, init_c=0.04532), scale = mu/shape, shape = shape) # Estimation capture.output( gp_model <- fitGPModel(group_data = group, likelihood = "gamma", y = y, params = params) , file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-0.5174554)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 6) # Prediction group_test <- c(1,3,3,9999) pred <- predict(gp_model, y=y, group_data_pred = group_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.2095341, -0.9170767, -0.9170767, 0.0000000) expected_cov <- c(0.08105393, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.09842279, 0.09842279, 0.00000000, 0.00000000, 0.09842279, 0.09842279, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.51745540) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Predict response pred <- predict(gp_model, y=y, group_data_pred = group_test, predict_var=TRUE, predict_response = TRUE) expected_mu <- c(1.2841038, 0.4198468, 0.4198468, 1.2952811) expected_var <- c(1.9273575, 0.2127346, 0.2127346, 3.9519573) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(pred$var-expected_var)),TOLERANCE_STRICT) # Evaluate negative log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.9),y=y) expect_lt(abs(nll-105.676137),TOLERANCE_LOOSE) # Also estimate shape parameter params_shape$optimizer_cov <- "nelder_mead" capture.output( gp_model <- fitGPModel(group_data = group, likelihood = "gamma", y = y, params = params_shape), file='NUL') cov_pars <- c(0.5141632) aux_pars <- c(0.9719373) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 45) # Also estimate shape parameter with gradient descent params_shape$optimizer_cov <- "gradient_descent" capture.output( gp_model <- fitGPModel(group_data = group, likelihood = "gamma", y = y, params = params_shape), file='NUL') cov_pars <- c(0.5198431) aux_pars <- c(0.970999) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 23) # Can set learning rate for auxiliary parameters via lr_cov params_temp <- params_shape params_temp$maxit = 1 capture.output( gp_model <- fitGPModel(group_data = group, likelihood = "gamma", y = y, params = params_temp), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-0.9058829)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-0.9297985)),TOLERANCE_STRICT) params_temp$lr_cov = 0.001 capture.output( gp_model <- fitGPModel(group_data = group, likelihood = "gamma", y = y, params = params_temp), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-0.998025)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-0.9985453)),TOLERANCE_STRICT) # Multiple random effects mu <- exp(Z1 %*% b_gr_1 + Z2 %*% b_gr_2 + Z3 %*% b_gr_3) y <- qgamma(sim_rand_unif(n=n, init_c=0.04532), scale = mu/shape, shape = shape) # Estimation capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, likelihood = "gamma", y = y, params = params) , file='NUL') cov_pars <- c(0.5050690, 1.2043329, 0.5280103) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 10) # Prediction group_data_pred = cbind(c(1,1,77),c(2,1,98)) group_rand_coef_data_pred = c(0,0.1,0.3) pred <- gp_model$predict(y = y, group_data_pred=group_data_pred, group_rand_coef_data_pred=group_rand_coef_data_pred, cov_pars = c(0.9,0.8,1.2), predict_var = TRUE, predict_response = FALSE) expected_mu <- c(0.1121777, 0.1972216, 0.0000000) expected_var <- c(0.2405621, 0.2259258, 1.8080000) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$var)-expected_var)),TOLERANCE_STRICT) # Also estimate shape parameter params_shape$optimizer_cov <- "nelder_mead" capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, likelihood = "gamma", y = y, params = params_shape), file='NUL') cov_pars <- c(0.5050897, 1.2026241, 0.5232070) aux_pars <- c(0.9819755) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 215) # Also estimate shape parameter with gradient descent params_shape$optimizer_cov <- "gradient_descent" capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, likelihood = "gamma", y = y, params = params_shape), file='NUL') cov_pars <- c(0.5065183, 1.2028488, 0.5360939) aux_pars <- c(0.9827199) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 31) # Also estimate shape parameter with adam params_shape$optimizer_cov <- "adam" capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, likelihood = "gamma", y = y, params = params_shape), file='NUL') cov_pars <- c(0.5052794, 1.2018843, 0.5230190) aux_pars <- c(0.9820493) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 279) # Also estimate shape parameter with bfgs params_shape$optimizer_cov <- "bfgs" capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, likelihood = "gamma", y = y, params = params_shape), file='NUL') cov_pars <- c(0.5052794, 1.2018842, 0.5230190) aux_pars <- c(0.9820493) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 8) # Also estimate shape parameter with gradient descent using internal initialization params_shape_no_init <- params_shape params_shape_no_init$init_aux_pars <- NULL params_shape_no_init$optimizer_cov <- "gradient_descent" capture.output( gp_model <- fitGPModel(group_data = cbind(group,group2), group_rand_coef_data = x, ind_effect_group_rand_coef = 1, likelihood = "gamma", y = y, params = params_shape_no_init), file='NUL') cov_pars <- c(0.5064068, 1.2028118, 0.5355322) aux_pars <- c(0.9826897) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 34) # Gaussian process model mu <- exp(L %*% b_1) y <- qgamma(sim_rand_unif(n=n, init_c=0.435), scale = mu/shape, shape = shape) # Estimation capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "gamma", y = y, params = params) , file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-c(1.0649094, 0.2738999))),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 8) # Prediction coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.3376250, 0.3023855, 0.7810425) expected_cov <- c(0.4567916157, 0.4033257822, -0.0002256179, 0.4033257822, 0.4540419202, -0.0002258048, -0.0002256179, -0.0002258048, 0.3368598330) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Evaluate approximate negative marginal log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.9,0.2),y=y) expect_lt(abs(nll-154.4561783),TOLERANCE_STRICT) # Also estimate shape parameter params_shape$optimizer_cov <- "nelder_mead" capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "gamma", y = y, params = params_shape) , file='NUL') cov_pars <- c(1.0447764, 0.2973086) aux_pars <- c(0.9402551) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 104) # Also estimate shape parameter with gradient descent params_shape$optimizer_cov <- "gradient_descent" capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "gamma", y = y, params = params_shape) , file='NUL') cov_pars <- c(1.0323441, 0.2898717) aux_pars <- c(0.9413081) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 26) ## Grouped random effects model with a linear predictor mu_lin <- exp(Z1 %*% b_gr_1 + X%*%beta) y_lin <- qgamma(sim_rand_unif(n=n, init_c=0.532), scale = mu_lin/shape, shape = shape) gp_model <- fitGPModel(group_data = group, likelihood = "gamma", y = y_lin, X=X, params = params) cov_pars <- c(0.474273) coef <- c(-0.07802971, 1.89766436) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coef)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 11) # Also estimate shape parameter params_shape$optimizer_cov <- "nelder_mead" gp_model <- fitGPModel(group_data = group, likelihood = "gamma", y = y_lin, X=X, params = params_shape) cov_pars <- c(0.5097316) coef <- c(-0.08623548, 1.90033132) aux_pars <- c(1.350364) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coef)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 258) # Also estimate shape parameter with gradient descent params_shape$optimizer_cov <- "gradient_descent" gp_model <- fitGPModel(group_data = group, likelihood = "gamma", y = y_lin, X=X, params = params_shape) cov_pars <- c(0.5179147) coef <- c(-0.08646342, 1.90053164) aux_pars <- c(1.351008) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_coef())-coef)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 42) ## Combined grouped random effects and Gaussian process model mu <- exp(L %*% b_1 + Z1 %*% b_gr_1) y <- qgamma(sim_rand_unif(n=n, init_c=0.987), scale = mu/shape, shape = shape) # Estimation capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", group_data = group, likelihood = "gamma", y = y, params = params) , file='NUL') cov_pars <- c(0.56585185, 0.62507125, 0.08278787) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 9) # Prediction coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) group_test <- c(1,3,3) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, group_data_pred=group_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(0.28574903, -0.67562130, 0.08821624) expected_cov <- c(0.649420831, 0.448952853, 0.007415143, 0.448952853, 0.683363103, 0.126645556, 0.007415143, 0.126645556, 0.531015480) expect_lt(sum(abs(pred$mu-expected_mu)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),TOLERANCE_STRICT) # Evaluate approximate negative marginal log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.6,0.9,0.2),y=y) expect_lt(abs(nll-123.3965571),TOLERANCE_STRICT) # Also estimate shape parameter params_shape$optimizer_cov <- "nelder_mead" capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", group_data = group, likelihood = "gamma", y = y, params = params_shape) , file='NUL') cov_pars <- c(0.62184856, 0.98925230, 0.07441182) aux_pars <- c(1.702741) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 222) # Also estimate shape parameter with gradient descent params_shape$optimizer_cov <- "gradient_descent" capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", group_data = group, likelihood = "gamma", y = y, params = params_shape) , file='NUL') cov_pars <- c(0.62143448, 0.98703748, 0.07443428) aux_pars <- c(1.707991) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),TOLERANCE_STRICT) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),TOLERANCE_STRICT) expect_equal(gp_model$get_num_optim_iter(), 27) # Gaussian process model with Vecchia approximation for(inv_method in c("cholesky", "iterative")){ if(inv_method == "iterative"){ tolerance_loc_1 <- TOLERANCE_ITERATIVE tolerance_loc_2 <- TOLERANCE_ITERATIVE } else{ tolerance_loc_1 <- 0.01 tolerance_loc_2 <- TOLERANCE_LOOSE } mu <- exp(0.75 * L %*% b_1) y <- qgamma(sim_rand_unif(n=n, init_c=0.7654), scale = mu/shape, shape = shape) # Estimation if(inv_method=="iterative"){ params$cg_delta_conv = 1e-6 params$num_rand_vec_trace=500 } capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "gamma", y = y, params = params, gp_approx = "vecchia", num_neighbors = 30, vecchia_ordering = "random", matrix_inversion_method = inv_method), file='NUL') expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-c(0.9484890, 0.0731435))),tolerance_loc_2) if(inv_method != "iterative"){ expect_lt(gp_model$get_num_optim_iter(), 11) expect_gt(gp_model$get_num_optim_iter(), 8) } # Prediction coord_test <- cbind(c(0.1,0.11,0.7),c(0.9,0.91,0.55)) gp_model$set_prediction_data(nsim_var_pred = 10000) pred <- predict(gp_model, y=y, gp_coords_pred = coord_test, predict_cov_mat = TRUE, predict_response = FALSE) expected_mu <- c(-0.1159426, -0.1028064, -0.3223582) expected_cov <- c(8.091398e-01, 1.079958e-01, -4.403387e-07, 1.079958e-01, 8.055727e-01, -4.442709e-07, -4.403387e-07, -4.442709e-07, 6.957873e-01) expect_lt(sum(abs(pred$mu-expected_mu)),tolerance_loc_1) adjust_tol <- 1 if (inv_method == "iterative") adjust_tol <- 1.5 expect_lt(sum(abs(as.vector(pred$cov)-expected_cov)),adjust_tol*tolerance_loc_1) # Evaluate approximate negative marginal log-likelihood nll <- gp_model$neg_log_likelihood(cov_pars=c(0.9,0.2),y=y) if(inv_method=="iterative"){ expect_lt(abs(nll-159.9221359),TOLERANCE_ITERATIVE) } else{ expect_lt(abs(nll-159.9221359),0.05) } # Also estimate shape parameter params_shape$optimizer_cov <- "nelder_mead" if(inv_method=="iterative"){ params_shape$cg_delta_conv = 1e-6 params_shape$num_rand_vec_trace=500 } capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "gamma", y = y, params = params_shape, gp_approx = "vecchia", num_neighbors = 30, vecchia_ordering = "random") , file='NUL') cov_pars <- c(1.14184253, 0.03605877) aux_pars <- c(1.328749) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),tolerance_loc_2) if(inv_method!="iterative"){ expect_gt(gp_model$get_num_optim_iter(), 117) expect_lt(gp_model$get_num_optim_iter(), 131) } # Also estimate shape parameter with gradient descent params_shape$optimizer_cov <- "gradient_descent" capture.output( gp_model <- fitGPModel(gp_coords = coords, cov_function = "exponential", likelihood = "gamma", y = y, params = params_shape, gp_approx = "vecchia", num_neighbors = 30, vecchia_ordering = "random") , file='NUL') cov_pars <- c(1.13722505, 0.03706853) aux_pars <- c(1.321834) expect_lt(sum(abs(as.vector(gp_model$get_cov_pars())-cov_pars)),tolerance_loc_2) expect_lt(sum(abs(as.vector(gp_model$get_aux_pars())-aux_pars)),tolerance_loc_2) if(inv_method!="iterative"){ expect_equal(gp_model$get_num_optim_iter(), 58) } } }) test_that("Saving a GPModel and loading from file works for non-Gaussian data", { probs <- pnorm(Z1 %*% b_gr_1 + X%*%beta) y <- as.numeric(sim_rand_unif(n=n, init_c=0.542) < probs) # Train model gp_model <- fitGPModel(group_data = group, likelihood = "bernoulli_probit", y = y, X=X, params = list(optimizer_cov = "gradient_descent", optimizer_coef = "gradient_descent", use_nesterov_acc = TRUE)) # Make predictions X_test <- cbind(rep(1,4),c(-0.5,0.2,0.4,1)) group_test <- c(1,3,3,9999) pred <- predict(gp_model, y=y, group_data_pred = group_test, X_pred = X_test, predict_cov_mat = TRUE, predict_response = FALSE) # Predict response pred_resp <- predict(gp_model, y=y, group_data_pred = group_test, X_pred = X_test, predict_var = TRUE, predict_response = TRUE) # Save model to file filename <- tempfile(fileext = ".json") saveGPModel(gp_model,filename = filename) # Delete model rm(gp_model) # Load from file and make predictions again gp_model_loaded <- loadGPModel(filename = filename) pred_loaded <- predict(gp_model_loaded, group_data_pred = group_test, X_pred = X_test, predict_cov_mat = TRUE, predict_response = FALSE) pred_resp_loaded <- predict(gp_model_loaded, y=y, group_data_pred = group_test, X_pred = X_test, predict_var = TRUE, predict_response = TRUE) expect_equal(pred$mu, pred_loaded$mu) expect_equal(pred$cov, pred_loaded$cov) expect_equal(pred_resp$mu, pred_resp_loaded$mu) expect_equal(pred_resp$var, pred_resp_loaded$var) }) }