diff --git a/DESCRIPTION b/DESCRIPTION
index d3072f1..edeff54 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: PrestoGP
 Type: Package
 Title: Penalized Regression for Spatio-Temporal Outcomes via Gaussian Processes 
-Version: 0.2.0.9037
+Version: 0.2.0.9039
 Authors@R: c(
     person(given = "Eric",
            family = "Bair",
@@ -64,4 +64,4 @@ Suggests:
     geoR,
     testthat (>= 3.0.0)
 Config/testthat/edition: 3
-URL: https://niehs.github.io/PrestoGP/
+URL: https://niehs.github.io/PrestoGP/, https://github.com/NIEHS/PrestoGP
diff --git a/NAMESPACE b/NAMESPACE
index 3a7f54f..d581ac2 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,6 +12,7 @@ export(create.param.sequence)
 export(createUMultivariate)
 export(enquo)
 export(enquos)
+export(get_Y)
 export(max_min_ordering)
 export(prestogp_predict)
 export(vecchia_Mlikelihood)
@@ -19,10 +20,12 @@ export(vecchia_Mprediction)
 export(vecchia_Mspecify)
 exportMethods(get_beta)
 exportMethods(get_converged)
+exportMethods(get_linear_model)
 exportMethods(get_neighbors)
 exportMethods(get_pen_loglik)
 exportMethods(get_scaling)
 exportMethods(get_theta)
+exportMethods(plot_beta)
 exportMethods(prestogp_fit)
 exportMethods(show)
 import(GPvecchia)
diff --git a/NEWS.md b/NEWS.md
new file mode 100644
index 0000000..f6e86bd
--- /dev/null
+++ b/NEWS.md
@@ -0,0 +1,42 @@
+# PrestoGP 0.2.0.9039 (2024-12-26)
+
+## MINOR IMPROVEMENTS
+
+  * Rewrote some tests that were failing on some verions of R
+
+# PrestoGP 0.2.0.9038 (2024-12-23)
+
+## NEW FEATURES
+
+  * Added the method `get_Y` to extract the Y values (including the imputed
+    Y values) from a PrestoGP model
+
+  * Added the method `get_linear_model` to extract the fitted cv.glmnet object
+    from a PrestoGP model
+	
+  * Added the method `plot_beta` to plot the glide path of the
+    regression coefficients
+
+## MINOR IMPROVEMENTS
+
+  * Added a NEWS.md file to track changes in each update
+
+  * Modified the `show` method to not print zero betas
+  
+  * Modified `prestogp_predict` to retun predictions in the form of a
+    list for multivariate models
+
+## BUG FIXES
+
+  * Fixed a bug in the intercept calculation in `get_beta`
+
+  * Fixed a bug where multivariate models with no variable names for X were
+    being assigned arbitrary variable names
+
+## DOCUMENTATION FIXES
+
+  * Updated the documention in `PrestoGPModel-class` to reflect that the
+    Y_train slot is updated to contain imputed values during model
+    fitting
+	
+  * Added documentation for all new methods
diff --git a/R/PrestoGP_Model.R b/R/PrestoGP_Model.R
index 6ea423b..9ded66e 100644
--- a/R/PrestoGP_Model.R
+++ b/R/PrestoGP_Model.R
@@ -21,6 +21,8 @@ setOldClass("cv.glmnet")
 #' "super matrix" for multivariate models. See
 #' \code{\link[psych]{superMatrix}}.
 #' @slot Y_train A column matrix containing the original response values.
+#' After fitting the model, missing values will be replaced with the imputed
+#' values if impute.y is TRUE. See \code{link{prestogp_fit}}.
 #' @slot X_ndx A vector used to find the elements of beta corresponding to
 #' specific outcomes. The \emph{i}th element of X_ndx is the index of the
 #' last element of beta corresponding to predictors for outcome \emph{i}.
@@ -117,12 +119,49 @@ setMethod("initialize", "PrestoGPModel", function(.Object, ...) {
   .Object
 })
 
+#' Extract the Y values for a PrestoGP model.
+#'
+#' This method extracts the values of Y for a PrestoGP model. If imputation
+#' was used when fitting the model, the missing Y's will be replaced with
+#' their imputed values.
+#'
+#' @param model The PrestoGP model object
+#'
+#' @return A vector or list containing the values of Y. Any missing Y's will
+#' be replaced with their imputed values.
+#'
+#' @seealso \code{\link{PrestoGPModel-class}}, \code{\link{prestogp_fit}}
+#'
+#' @references
+#' \itemize{
+#' \item Messier, K.P. and Katzfuss, M. "Scalable penalized spatiotemporal
+#' land-use regression for ground-level nitrogen dioxide", The Annals of
+#' Applied Statistics (2021) 15(2):688-710.
+#' }
+#'
+#' @rdname get_Y
+#' @export
+#'
+#' @examples
+#' data(soil)
+#' soil <- soil[!is.na(soil[,5]),] # remove rows with NA's
+#' y <- soil[,4]                   # predict moisture content
+#' X <- as.matrix(soil[,5:9])
+#' locs <- as.matrix(soil[,1:2])
+#'
+#' soil.vm <- new("VecchiaModel", n_neighbors = 10)
+#' soil.vm <- prestogp_fit(soil.vm, y, X, locs)
+#' get_Y(soil.vm)
+setGeneric("get_Y", function(model) standardGeneric("get_Y"))
 setGeneric("get_theta", function(model) standardGeneric("get_theta"))
 setGeneric("get_beta", function(model) standardGeneric("get_beta"))
+setGeneric(
+  "get_linear_model", function(model) standardGeneric("get_linear_model"))
 setGeneric("get_neighbors", function(model) standardGeneric("get_neighbors"))
 setGeneric("get_scaling", function(model) standardGeneric("get_scaling"))
 setGeneric("get_converged", function(model) standardGeneric("get_converged"))
 setGeneric("get_pen_loglik", function(model) standardGeneric("get_pen_loglik"))
+setGeneric("plot_beta", function(model, ...) standardGeneric("plot_beta"))
 setGeneric(
   "prestogp_fit",
   function(model, Y, X, locs, Y.names = NULL, X.names = NULL, scaling = NULL,
@@ -262,7 +301,11 @@ setMethod(
     cat("Matern covariance parameters (theta):", "\n")
     print(get_theta(object))
     cat("Regression coefficients (beta):", "\n")
-    print(get_beta(object))
+    cur.coef <- get_beta(object)
+    for (i in seq_along(cur.coef)) {
+      cur.coef[[i]] <- cur.coef[[i]][cur.coef[[i]] != 0]
+    }
+    print(cur.coef)
     cat("Model type:", is(object)[1], "\n")
     if (object@n_neighbors > 0) {
       cat("Nearest neighbors:", object@n_neighbors, "\n")
@@ -411,13 +454,14 @@ setMethod(
     for (i in seq_len(P)) {
       if (i == 1) {
         indx <- 1
-        junk[[P + 1]] <- model@beta[1]
+        junk[[P + 1]] <- model@beta[1] + model@Y_bar[1]
         names(junk[[P + 1]]) <- "(Intercept)"
       } else {
         indx <- model@X_ndx[i - 1] + 1
         new.names <- c(names(junk[[P + 1]]),
           names(model@locs_train)[i])
-        junk[[P + 1]] <- c(junk[[P + 1]], model@beta[indx])
+        junk[[P + 1]] <- c(junk[[P + 1]], model@beta[indx] + model@Y_bar[i] -
+            model@Y_bar[1])
         names(junk[[P + 1]]) <- new.names
       }
       pred.ndx <- (indx + 1):model@X_ndx[i]
@@ -433,6 +477,51 @@ setMethod(
   }
 )
 
+#' Extract the fitted linear model for a PrestoGP model
+#'
+#' This method return the fitted linear model (of class cv.glmnet) for a
+#' PrestoGP model.
+#'
+#' @param model The PrestoGP model object
+#'
+#' @details It is important to note that the model is fit to the
+#' transformed data. The CV error rate and predicted values of Y will not be
+#' correct for the original (untransformed) data. This method should be
+#' used primarily for examining the coefficient path and generating plots.
+#'
+#' @return The fitted linear model (of class cv.glmnet).
+#'
+#' @seealso \code{\link{PrestoGPModel-class}}, \code{\link{prestogp_fit}},
+#' \code{\link[glmnet]{cv.glmnet}}
+#'
+#' @references
+#' \itemize{
+#' \item Messier, K.P. and Katzfuss, M. "Scalable penalized spatiotemporal
+#' land-use regression for ground-level nitrogen dioxide", The Annals of
+#' Applied Statistics (2021) 15(2):688-710.
+#' }
+#'
+#' @aliases get_linear_model
+#' @export
+#'
+#' @examples
+#' data(soil)
+#' soil <- soil[!is.na(soil[,5]),] # remove rows with NA's
+#' y <- soil[,4]                   # predict moisture content
+#' X <- as.matrix(soil[,5:9])
+#' locs <- as.matrix(soil[,1:2])
+#'
+#' soil.vm <- new("VecchiaModel", n_neighbors = 10)
+#' soil.vm <- prestogp_fit(soil.vm, y, X, locs)
+#' get_linear_model(soil.vm)
+
+setMethod(
+  "get_linear_model", "PrestoGPModel",
+  function(model) {
+    return(model@linear_model)
+  }
+)
+
 #' Extract the number of neighbors for a PrestoGP model.
 #'
 #' This method extracts the number of nearest neighbors that were used to
@@ -590,6 +679,47 @@ setMethod(
   }
 )
 
+#' Plots the glide path for the coefficients for a PrestoGP model.
+#'
+#' This method generates a plot showing the coefficients of the model for
+#' different values of the tuning parameter. It is a wrapper for
+#' \code{\link[glmnet]{plot.glmnet}}.
+#'
+#' @param model The PrestoGP model object
+#'
+#' @param ... Additional parameters to \code{\link[glmnet]{plot.glmnet}}
+#'
+#' @seealso \code{\link{PrestoGPModel-class}}, \code{\link{prestogp_fit}},
+#' \code{\link[glmnet]{plot.glmnet}}
+#'
+#' @references
+#' \itemize{
+#' \item Messier, K.P. and Katzfuss, M. "Scalable penalized spatiotemporal
+#' land-use regression for ground-level nitrogen dioxide", The Annals of
+#' Applied Statistics (2021) 15(2):688-710.
+#' }
+#'
+#' @aliases plot_beta
+#' @export
+#'
+#' @examples
+#' data(soil)
+#' soil <- soil[!is.na(soil[,5]),] # remove rows with NA's
+#' y <- soil[,4]                   # predict moisture content
+#' X <- as.matrix(soil[,5:9])
+#' locs <- as.matrix(soil[,1:2])
+#'
+#' soil.vm <- new("VecchiaModel", n_neighbors = 10)
+#' soil.vm <- prestogp_fit(soil.vm, y, X, locs)
+#' plot_beta(soil.vm)
+
+setMethod(
+  "plot_beta", "PrestoGPModel",
+  function(model, ...) {
+    plot(model@linear_model$glmnet.fit, ...)
+  }
+)
+
 #' Train a PrestoGP model.
 #'
 #' This method fits a PrestoGP model given a set of locations and predictor
diff --git a/R/PrestoGP_Multivariate_Vecchia.R b/R/PrestoGP_Multivariate_Vecchia.R
index df13c49..f8106b3 100644
--- a/R/PrestoGP_Multivariate_Vecchia.R
+++ b/R/PrestoGP_Multivariate_Vecchia.R
@@ -31,6 +31,20 @@ setMethod("initialize", "MultivariateVecchiaModel", function(.Object, n_neighbor
   .Object
 })
 
+#' @rdname get_Y
+setMethod("get_Y", "MultivariateVecchiaModel",
+  function(model) {
+    Y.all <- model@Y_train
+    Y.out <- list()
+    for (i in seq_along(model@locs_train)) {
+      cur.y <- seq_len(nrow(model@locs_train[[i]]))
+      Y.out[[i]] <- Y.all[cur.y] + model@Y_bar[i]
+      Y.all <- Y.all[-cur.y]
+    }
+    return(Y.out)
+  }
+)
+
 #' @rdname prestogp_predict
 setMethod("prestogp_predict", "MultivariateVecchiaModel",
   function(model, X, locs, m = NULL, ordering.pred = c("obspred", "general"), pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) {
@@ -93,18 +107,23 @@ setMethod("prestogp_predict", "MultivariateVecchiaModel",
 
     # prediction function can return both mean and sds
     # returns a list with elements mu.pred,mu.obs,var.pred,var.obs,V.ord
-    Vec.mean <- Y_bar + pred$mu.pred + Vecchia.Pred # residual + mean trend
+    Vec.mean.all <- Y_bar + pred$mu.pred + Vecchia.Pred # residual + mean trend
+    ndx.out <- NULL
+    for (i in seq_along(locs)) {
+      ndx.out <- c(ndx.out, rep(i, nrow(locs[[i]])))
+    }
     if (return.values == "mean") {
+      Vec.mean <- list()
+      for (i in seq_along(locs)) {
+        Vec.mean[[i]] <- Vec.mean.all[ndx.out == i]
+      }
+      names(Vec.mean) <- names(model@locs_train)
       return.list <- list(means = Vec.mean)
     } else {
       warning("Variance estimates do not include model fitting variance and are anticonservative. Use with caution.")
-      Vec.sds <- pred$var.pred
-      ndx.out <- NULL
+      Vec.sds <- list()
       for (i in seq_along(locs)) {
-        ndx.out <- c(ndx.out, rep(i, nrow(locs[[i]])))
-      }
-      for (i in seq_along(locs)) {
-        Vec.sds[ndx.out == i] <- sqrt(Vec.sds[ndx.out == i] +
+        Vec.sds[[i]] <- sqrt(pred$var.pred[ndx.out == i] +
             model@covparams[model@param_sequence[4, i]])
       }
       return.list <- list(means = Vec.mean, sds = Vec.sds)
@@ -141,9 +160,7 @@ setMethod("check_input", "MultivariateVecchiaModel", function(model, Y, X, locs,
     }
   }
   if (is.null(names(locs))) {
-    for (i in seq_along(locs)) {
-      names(locs)[i] <- paste0("Y", i)
-    }
+    names(locs) <- paste0("Y", seq_along(locs))
   }
   if (!is.null(lod)) {
     if (!is.list(lod)) {
@@ -241,8 +258,15 @@ setMethod("check_input", "MultivariateVecchiaModel", function(model, Y, X, locs,
   if (length(X) == 1) {
     model@X_train <- X[[1]]
   } else {
+    if (is.null(colnames(X[[1]]))) {
+      colnames(X[[1]]) <- paste0(names(locs)[1], "_", seq_len(ncol(X[[1]])))
+    }
     for (i in 2:length(X)) {
       X[[i]] <- cbind(rep(1, nrow(X[[i]])), X[[i]])
+      if (is.null(colnames(X[[i]]))) {
+        colnames(X[[i]]) <- paste0(names(locs)[i], "_",
+          (seq_len(ncol(X[[i]])) - 1))
+      }
       model@X_ndx <- c(model@X_ndx, model@X_ndx[i - 1] + ncol(X[[i]]))
     }
     model@X_train <- psych::superMatrix(X)
diff --git a/R/PrestoGP_Vecchia.R b/R/PrestoGP_Vecchia.R
index 3b04d1b..2de181e 100644
--- a/R/PrestoGP_Vecchia.R
+++ b/R/PrestoGP_Vecchia.R
@@ -28,6 +28,13 @@ setMethod("initialize", "VecchiaModel", function(.Object, n_neighbors = 25, ...)
   .Object
 })
 
+#' @rdname get_Y
+setMethod("get_Y", "VecchiaModel",
+  function(model) {
+    return(as.vector(model@Y_train + model@Y_bar))
+  }
+)
+
 #' @rdname prestogp_predict
 setMethod("prestogp_predict", "VecchiaModel",
   function(model, X, locs, m = NULL, ordering.pred = c("obspred", "general"), pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) {
diff --git a/man/PrestoGP-package.Rd b/man/PrestoGP-package.Rd
index 674f053..84e2edb 100644
--- a/man/PrestoGP-package.Rd
+++ b/man/PrestoGP-package.Rd
@@ -14,6 +14,7 @@ Simultaneous variable seletion and estimation of LUR models with spatiotemporall
 Useful links:
 \itemize{
   \item \url{https://niehs.github.io/PrestoGP/}
+  \item \url{https://github.com/NIEHS/PrestoGP}
 }
 
 }
diff --git a/man/PrestoGPModel-class.Rd b/man/PrestoGPModel-class.Rd
index e1a9abc..e1ba7e4 100644
--- a/man/PrestoGPModel-class.Rd
+++ b/man/PrestoGPModel-class.Rd
@@ -33,7 +33,9 @@ the glmnet model. See \code{\link[glmnet]{glmnet}}.}
 "super matrix" for multivariate models. See
 \code{\link[psych]{superMatrix}}.}
 
-\item{\code{Y_train}}{A column matrix containing the original response values.}
+\item{\code{Y_train}}{A column matrix containing the original response values.
+After fitting the model, missing values will be replaced with the imputed
+values if impute.y is TRUE. See \code{link{prestogp_fit}}.}
 
 \item{\code{X_ndx}}{A vector used to find the elements of beta corresponding to
 specific outcomes. The \emph{i}th element of X_ndx is the index of the
diff --git a/man/get_Y.Rd b/man/get_Y.Rd
new file mode 100644
index 0000000..309dcd9
--- /dev/null
+++ b/man/get_Y.Rd
@@ -0,0 +1,48 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PrestoGP_Model.R, R/PrestoGP_Vecchia.R,
+%   R/PrestoGP_Multivariate_Vecchia.R
+\name{get_Y}
+\alias{get_Y}
+\alias{get_Y,VecchiaModel-method}
+\alias{get_Y,MultivariateVecchiaModel-method}
+\title{Extract the Y values for a PrestoGP model.}
+\usage{
+get_Y(model)
+
+\S4method{get_Y}{VecchiaModel}(model)
+
+\S4method{get_Y}{MultivariateVecchiaModel}(model)
+}
+\arguments{
+\item{model}{The PrestoGP model object}
+}
+\value{
+A vector or list containing the values of Y. Any missing Y's will
+be replaced with their imputed values.
+}
+\description{
+This method extracts the values of Y for a PrestoGP model. If imputation
+was used when fitting the model, the missing Y's will be replaced with
+their imputed values.
+}
+\examples{
+data(soil)
+soil <- soil[!is.na(soil[,5]),] # remove rows with NA's
+y <- soil[,4]                   # predict moisture content
+X <- as.matrix(soil[,5:9])
+locs <- as.matrix(soil[,1:2])
+
+soil.vm <- new("VecchiaModel", n_neighbors = 10)
+soil.vm <- prestogp_fit(soil.vm, y, X, locs)
+get_Y(soil.vm)
+}
+\references{
+\itemize{
+\item Messier, K.P. and Katzfuss, M. "Scalable penalized spatiotemporal
+land-use regression for ground-level nitrogen dioxide", The Annals of
+Applied Statistics (2021) 15(2):688-710.
+}
+}
+\seealso{
+\code{\link{PrestoGPModel-class}}, \code{\link{prestogp_fit}}
+}
diff --git a/man/get_linear_model-PrestoGPModel-method.Rd b/man/get_linear_model-PrestoGPModel-method.Rd
new file mode 100644
index 0000000..658c7db
--- /dev/null
+++ b/man/get_linear_model-PrestoGPModel-method.Rd
@@ -0,0 +1,47 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PrestoGP_Model.R
+\name{get_linear_model,PrestoGPModel-method}
+\alias{get_linear_model,PrestoGPModel-method}
+\alias{get_linear_model}
+\title{Extract the fitted linear model for a PrestoGP model}
+\usage{
+\S4method{get_linear_model}{PrestoGPModel}(model)
+}
+\arguments{
+\item{model}{The PrestoGP model object}
+}
+\value{
+The fitted linear model (of class cv.glmnet).
+}
+\description{
+This method return the fitted linear model (of class cv.glmnet) for a
+PrestoGP model.
+}
+\details{
+It is important to note that the model is fit to the
+transformed data. The CV error rate and predicted values of Y will not be
+correct for the original (untransformed) data. This method should be
+used primarily for examining the coefficient path and generating plots.
+}
+\examples{
+data(soil)
+soil <- soil[!is.na(soil[,5]),] # remove rows with NA's
+y <- soil[,4]                   # predict moisture content
+X <- as.matrix(soil[,5:9])
+locs <- as.matrix(soil[,1:2])
+
+soil.vm <- new("VecchiaModel", n_neighbors = 10)
+soil.vm <- prestogp_fit(soil.vm, y, X, locs)
+get_linear_model(soil.vm)
+}
+\references{
+\itemize{
+\item Messier, K.P. and Katzfuss, M. "Scalable penalized spatiotemporal
+land-use regression for ground-level nitrogen dioxide", The Annals of
+Applied Statistics (2021) 15(2):688-710.
+}
+}
+\seealso{
+\code{\link{PrestoGPModel-class}}, \code{\link{prestogp_fit}},
+\code{\link[glmnet]{cv.glmnet}}
+}
diff --git a/man/plot_beta-PrestoGPModel-method.Rd b/man/plot_beta-PrestoGPModel-method.Rd
new file mode 100644
index 0000000..2c5bafc
--- /dev/null
+++ b/man/plot_beta-PrestoGPModel-method.Rd
@@ -0,0 +1,41 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PrestoGP_Model.R
+\name{plot_beta,PrestoGPModel-method}
+\alias{plot_beta,PrestoGPModel-method}
+\alias{plot_beta}
+\title{Plots the glide path for the coefficients for a PrestoGP model.}
+\usage{
+\S4method{plot_beta}{PrestoGPModel}(model, ...)
+}
+\arguments{
+\item{model}{The PrestoGP model object}
+
+\item{...}{Additional parameters to \code{\link[glmnet]{plot.glmnet}}}
+}
+\description{
+This method generates a plot showing the coefficients of the model for
+different values of the tuning parameter. It is a wrapper for
+\code{\link[glmnet]{plot.glmnet}}.
+}
+\examples{
+data(soil)
+soil <- soil[!is.na(soil[,5]),] # remove rows with NA's
+y <- soil[,4]                   # predict moisture content
+X <- as.matrix(soil[,5:9])
+locs <- as.matrix(soil[,1:2])
+
+soil.vm <- new("VecchiaModel", n_neighbors = 10)
+soil.vm <- prestogp_fit(soil.vm, y, X, locs)
+plot_beta(soil.vm)
+}
+\references{
+\itemize{
+\item Messier, K.P. and Katzfuss, M. "Scalable penalized spatiotemporal
+land-use regression for ground-level nitrogen dioxide", The Annals of
+Applied Statistics (2021) 15(2):688-710.
+}
+}
+\seealso{
+\code{\link{PrestoGPModel-class}}, \code{\link{prestogp_fit}},
+\code{\link[glmnet]{plot.glmnet}}
+}
diff --git a/tests/testthat/sim_multivariate_big_st.RData b/tests/testthat/sim_multivariate_big_st.RData
new file mode 100644
index 0000000..88e27d7
Binary files /dev/null and b/tests/testthat/sim_multivariate_big_st.RData differ
diff --git a/tests/testthat/test-Log_Likelihood.R b/tests/testthat/test-Log_Likelihood.R
index ff6bb26..061485d 100644
--- a/tests/testthat/test-Log_Likelihood.R
+++ b/tests/testthat/test-Log_Likelihood.R
@@ -200,9 +200,8 @@ test_that("mvnegloglik", {
 })
 
 test_that("mvnegloglik_ST", {
-  source("sim_multivariate_big_st.R")
+  load("sim_multivariate_big_st.RData")
   P <- 3
-  Y <- cbind(runif(10), runif(10), runif(10))
   cor.matrix <- cor(Y)
   cov_mat <- c(cor.matrix[upper.tri(cor.matrix)])
   logparams <- create.initial.values.flex(
diff --git a/tests/testthat/test-PrestoGP_Multivariate_Vecchia_Model.R b/tests/testthat/test-PrestoGP_Multivariate_Vecchia_Model.R
index ffa25cb..54a41b3 100644
--- a/tests/testthat/test-PrestoGP_Multivariate_Vecchia_Model.R
+++ b/tests/testthat/test-PrestoGP_Multivariate_Vecchia_Model.R
@@ -271,8 +271,13 @@ test_that("Simulated dataset multivariate spatial", {
     )
   )
 
+  expect_silent(plot_beta(pgp.mmodel1))
+  dev.off()
+
   expect_true(validObject(pgp.mmodel1))
 
+  expect_equal(get_Y(pgp.mmodel1), y.list)
+  expect_equal(get_linear_model(pgp.mmodel1), pgp.mmodel1@linear_model)
   expect_equal(get_neighbors(pgp.mmodel1), pgp.mmodel1@n_neighbors)
   expect_equal(get_scaling(pgp.mmodel1), pgp.mmodel1@scaling)
   expect_equal(get_converged(pgp.mmodel1), pgp.mmodel1@converged)
@@ -326,12 +331,12 @@ test_that("Simulated dataset multivariate spatial", {
   expect_identical(as.vector(beta.out[[1]]), as.vector(pgp.mmodel1@beta[2:11]))
   expect_identical(as.vector(beta.out[[2]]), as.vector(pgp.mmodel1@beta[13:22]))
   expect_identical(as.vector(beta.out[[3]]), as.vector(pgp.mmodel1@beta[24:33]))
-  expect_identical(as.numeric(beta.out[[4]][1]),
+  expect_equal(as.numeric(beta.out[[4]][1] - mean(y.list[[1]])),
     as.numeric(pgp.mmodel1@beta[1]))
-  expect_identical(as.numeric(beta.out[[4]][2]),
-    as.numeric(pgp.mmodel1@beta[12]))
-  expect_identical(as.numeric(beta.out[[4]][3]),
-    as.numeric(pgp.mmodel1@beta[23]))
+  expect_identical(as.numeric(beta.out[[4]][2] + mean(y.list[[1]]) -
+        mean(y.list[[2]])), as.numeric(pgp.mmodel1@beta[12]))
+  expect_identical(as.numeric(beta.out[[4]][3] + mean(y.list[[1]]) -
+        mean(y.list[[3]])), as.numeric(pgp.mmodel1@beta[23]))
 
   expect_equal(as.numeric(beta.out[[1]]), c(0.98, 0.97, 0.92, 0.96, rep(0, 3),
       0.02, rep(0, 2)), tolerance = 0.04)
@@ -339,7 +344,9 @@ test_that("Simulated dataset multivariate spatial", {
       0.04, rep(0, 2)), tolerance = 0.04)
   expect_equal(as.numeric(beta.out[[3]]), c(1.0, 0.99, 0.96, 0.94, rep(0, 6)),
     tolerance = 0.04)
-  expect_equal(as.numeric(beta.out[[4]]), c(-0.01, 0, 0), tolerance = 0.04)
+  expect_equal(as.numeric(beta.out[[4]]), c(-0.01 + mean(y.list[[1]]),
+      mean(y.list[[2]]) - mean(y.list[[1]]),
+      mean(y.list[[3]]) - mean(y.list[[1]])), tolerance = 0.04)
 
   expect_equal(params.out[1], 2.3, tolerance = 5.1)
   expect_equal(params.out[2], 4.4, tolerance = 4.6)
@@ -399,6 +406,12 @@ test_that("Simulated dataset multivariate spatial", {
   expect_equal(rownames(theta.out2[[5]]), c("Y1", "Y2", "Y3"))
 
   expect_named(beta.out2, c("Y1", "Y2", "Y3", "(Intercept)"))
+  expect_named(beta.out2[[1]], c("Y1_1", "Y1_2", "Y1_3", "Y1_4", "Y1_5",
+      "Y1_6", "Y1_7", "Y1_8", "Y1_9", "Y1_10"))
+  expect_named(beta.out2[[2]], c("Y2_1", "Y2_2", "Y2_3", "Y2_4", "Y2_5",
+      "Y2_6", "Y2_7", "Y2_8", "Y2_9", "Y2_10"))
+  expect_named(beta.out2[[3]], c("Y3_1", "Y3_2", "Y3_3", "Y3_4", "Y3_5",
+      "Y3_6", "Y3_7", "Y3_8", "Y3_9", "Y3_10"))
   expect_named(beta.out2[[4]], c("(Intercept)", "Y2", "Y3"))
 
   # Results should be the same after imputation
@@ -442,8 +455,14 @@ test_that("Simulated dataset multivariate spatial", {
     tolerance = 0.05)
   expect_equal(as.numeric(beta.out[[3]]), as.numeric(beta.out3[[3]]),
     tolerance = 0.05)
-  expect_equal(as.numeric(beta.out[[4]]), as.numeric(beta.out3[[4]]),
-    tolerance = 0.05)
+  expect_equal(as.numeric(beta.out[[4]][1] - pgp.mmodel1@Y_bar[1]),
+    as.numeric(beta.out3[[4]][1] - pgp.mmodel3@Y_bar[1]), tolerance = 0.05)
+  expect_equal(as.numeric(beta.out[[4]][2] - pgp.mmodel1@Y_bar[2] +
+        pgp.mmodel1@Y_bar[1]), as.numeric(beta.out3[[4]][2] -
+        pgp.mmodel3@Y_bar[2] + pgp.mmodel3@Y_bar[1]), tolerance = 0.05)
+  expect_equal(as.numeric(beta.out[[4]][3] - pgp.mmodel1@Y_bar[3] +
+        pgp.mmodel1@Y_bar[1]), as.numeric(beta.out3[[4]][3] -
+        pgp.mmodel3@Y_bar[3] + pgp.mmodel3@Y_bar[1]), tolerance = 0.05)
   expect_equal(params.out[1], params.out3[1], tolerance = 7.5)
   expect_equal(params.out[2], params.out3[2], tolerance = 7.2)
   expect_equal(params.out[3], params.out3[3], tolerance = 7.1)
@@ -783,11 +802,31 @@ test_that("Simulated dataset multivariate spatial prediction", {
 
   pgp.mmodel1.pred <- prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst)
 
-  mse <- mean((pgp.mmodel1.pred$means - unlist(y.list.otst))^2)
-  me <- mean(pgp.mmodel1.pred$means - unlist(y.list.otst))
-
+  expect_named(pgp.mmodel1.pred$means, c("Y1", "Y2", "Y3"))
+
+  mse1 <- mean((pgp.mmodel1.pred$means[[1]] - y.list.otst[[1]])^2)
+  me1 <- mean(pgp.mmodel1.pred$means[[1]] - y.list.otst[[1]])
+  mse2 <- mean((pgp.mmodel1.pred$means[[2]] - y.list.otst[[2]])^2)
+  me2 <- mean(pgp.mmodel1.pred$means[[2]] - y.list.otst[[2]])
+  mse3 <- mean((pgp.mmodel1.pred$means[[3]] - y.list.otst[[3]])^2)
+  me3 <- mean(pgp.mmodel1.pred$means[[3]] - y.list.otst[[3]])
+  mse <- mean((unlist(pgp.mmodel1.pred$means) - unlist(y.list.otst))^2)
+  me <- mean(unlist(pgp.mmodel1.pred$means) - unlist(y.list.otst))
+
+  expect_gt(mse1, 1.7)
+  expect_lt(mse1, 2.1)
+  expect_gt(me1, -0.12)
+  expect_lt(me1, -0.01)
+  expect_gt(mse2, 1.1)
+  expect_lt(mse2, 1.4)
+  expect_gt(me2, -0.29)
+  expect_lt(me2, -0.15)
+  expect_gt(mse3, 2.7)
+  expect_lt(mse3, 3.0)
+  expect_gt(me3, 0.26)
+  expect_lt(me3, 0.4)
   expect_gt(mse, 1.9)
-  expect_lt(mse, 2.3)
-  expect_gt(me, -0.02)
-  expect_lt(me, 0.035)
+  expect_lt(mse, 3.0)
+  expect_gt(me, -0.012)
+  expect_lt(me, 0.042)
 })
diff --git a/tests/testthat/test-PrestoGP_Univariate.R b/tests/testthat/test-PrestoGP_Univariate.R
index 3b39251..ef5ccf3 100644
--- a/tests/testthat/test-PrestoGP_Univariate.R
+++ b/tests/testthat/test-PrestoGP_Univariate.R
@@ -153,9 +153,14 @@ test_that("Simulated dataset spatial", {
     )
   )
 
+  expect_silent(plot_beta(pgp.model1))
+  dev.off()
+
   expect_true(validObject(pgp.model1))
   show(pgp.model1)
 
+  expect_equal(get_Y(pgp.model1), y)
+  expect_equal(get_linear_model(pgp.model1), pgp.model1@linear_model)
   expect_equal(get_neighbors(pgp.model1), pgp.model1@n_neighbors)
   expect_equal(get_scaling(pgp.model1), pgp.model1@scaling)
   expect_equal(get_converged(pgp.model1), pgp.model1@converged)
@@ -188,12 +193,13 @@ test_that("Simulated dataset spatial", {
   expect_identical(as.numeric(theta.out[[2]]), params.out[2])
   expect_identical(as.numeric(theta.out[[3]]), params.out[3])
   expect_identical(as.numeric(theta.out[[4]]), params.out[4])
-  expect_identical(as.numeric(beta.out[[2]]), as.numeric(pgp.model1@beta[1]))
+  expect_equal(as.numeric(beta.out[[2]] - mean(y)),
+    as.numeric(pgp.model1@beta[1]))
   expect_identical(as.vector(beta.out[[1]]), as.vector(pgp.model1@beta[2:11]))
 
   expect_equal(as.numeric(beta.out[[1]]), c(0.86, 0.98, 0.94, 0.9, rep(0, 6)),
     tolerance = 0.03)
-  expect_equal(as.numeric(beta.out[[2]]), 0.01, tolerance = 0.03)
+  expect_equal(as.numeric(beta.out[[2]] - mean(y)), 0.01, tolerance = 0.03)
   expect_equal(params.out[1], 1.6, tolerance = 0.5)
   expect_equal(params.out[2] - 0.4, 0, tolerance = 0.2)
   expect_equal(params.out[3], 0.59, tolerance = 0.2)
@@ -222,7 +228,7 @@ test_that("Simulated dataset spatial", {
 
   expect_equal(beta.out2[[1]], c(0.86, 0.98, 0.95, 0.9, rep(0, 6)),
     tolerance = 0.03)
-  expect_equal(as.numeric(beta.out2[[2]]), 0.01, tolerance = 0.03)
+  expect_equal(as.numeric(beta.out2[[2]] - mean(y)), 0.01, tolerance = 0.03)
   expect_equal(params.out2[1], 1.5, tolerance = 0.6)
   expect_equal(params.out2[2] - 0.4, 0, tolerance = 0.15)
   expect_equal(params.out2[3], 0.62, tolerance = 0.2)
@@ -256,11 +262,12 @@ test_that("Simulated dataset spatial", {
   params.out3 <- pgp.model3@covparams
 
   # Results should be the same after imputation
-  expect_equal(as.numeric(beta.out[[2]]), beta.out3[1], tolerance = 0.07)
+  expect_equal(as.numeric(beta.out[[2]] - mean(y)), beta.out3[1],
+    tolerance = 0.07)
   expect_equal(as.numeric(beta.out[[1]]), beta.out3[-1], tolerance = 0.07)
   expect_equal(params.out[1], params.out3[1], tolerance = 1.1)
   expect_equal(params.out[2] - params.out3[2], 0, tolerance = 0.3)
-  expect_equal(params.out[3], params.out3[3], tolerance = 0.3)
+  expect_equal(params.out[3] - params.out3[3], 0, tolerance = 0.3)
   expect_equal(params.out[4], params.out3[4], tolerance = 0.4)
 
   # Missing data with lod
@@ -282,7 +289,8 @@ test_that("Simulated dataset spatial", {
   params.out4 <- pgp.model4@covparams
 
   # Results should be the same after imputation
-  expect_equal(as.numeric(beta.out[[2]]), beta.out4[1], tolerance = 0.09)
+  expect_equal(as.numeric(beta.out[[2]] - mean(y)), beta.out4[1],
+    tolerance = 0.09)
   expect_equal(as.numeric(beta.out[[1]]), beta.out4[-1], tolerance = 0.09)
   expect_equal(params.out[1], params.out3[1], tolerance = 1.3)
   expect_equal(params.out[2], params.out3[2], tolerance = 0.8)