diff --git a/.Rbuildignore b/.Rbuildignore
index b9a994a..8c79af2 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -4,3 +4,4 @@ MASTER.R
 TODO.*
 cran-comments.md
 ^CRAN-SUBMISSION$
+revdep
diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION
index b8194e5..47fc187 100644
--- a/CRAN-SUBMISSION
+++ b/CRAN-SUBMISSION
@@ -1,3 +1,3 @@
-Version: 1.0.0
-Date: 2022-07-07 17:31:28 UTC
-SHA: f25c091cb14623d70c2ce89f349d8e338ae4e71d
+Version: 1.0.1
+Date: 2024-03-29 15:33:13 UTC
+SHA: a6166b74bbd812f1e28652612228db1a837c6e8a
diff --git a/DESCRIPTION b/DESCRIPTION
index 4913d0c..7967b43 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: cgaim
 Title: Constrained Groupwise Additive Index Models
-Version: 1.0.0
+Version: 1.0.1
 Authors@R: 
     person(given = "Pierre",
            family = "Masselot",
@@ -13,7 +13,7 @@ Encoding: UTF-8
 LazyData: true
 Suggests: 
     testthat (>= 2.1.0)
-RoxygenNote: 7.2.0
+RoxygenNote: 7.3.1
 Imports: 
     graphics,
     stats,
diff --git a/NAMESPACE b/NAMESPACE
index 7c63ecc..21a9eff 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,6 +5,8 @@ S3method(confint,cgaim)
 S3method(plot,cgaim)
 S3method(predict,cgaim)
 S3method(print,cgaim)
+S3method(vcov,boot.cgaim)
+S3method(vcov,cgaim)
 export(boot.cgaim)
 export(build_constraints)
 export(cgaim)
diff --git a/NEWS.md b/NEWS.md
index 25a4491..e19d8b4 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,12 @@
+## 1.0.1
+
+### New features
+- New `vcov` method
+
+### Bug fixes
+- Adapts `smooth_scam` function to new version of package `gratia`
+- Fixes names in confidence intervals from `confint.cgaim`
+
 ## 1.0.0
 
 ### Changes
diff --git a/R/confint.boot.cgaim.R b/R/confint.boot.cgaim.R
index 34a4d66..fe6c1a0 100644
--- a/R/confint.boot.cgaim.R
+++ b/R/confint.boot.cgaim.R
@@ -36,7 +36,7 @@ confint.boot.cgaim <- function(object, parm, level = 0.95, ...)
       alims, na.rm = TRUE))
     
     # Names
-    rownames(res$alpha) <- names(unlist(object$alpha))
+    rownames(res$alpha) <- names(unlist(object$obs$alpha))
     colnames(res$alpha) <- level.labels
   }
   if (2 %in% parm){
diff --git a/R/smooth_scam.R b/R/smooth_scam.R
index 52cf43b..65ed61f 100644
--- a/R/smooth_scam.R
+++ b/R/smooth_scam.R
@@ -17,7 +17,7 @@ smooth_scam <- function(x, y, formula, Xcov, ...)
     # print("scam"); flush.console()
   } else {
     gfit <- mgcv::gam(formula, data = scam_data, ...)
-    derivs <- as.data.frame(gratia::derivatives(gfit, newdata = scam_data))
+    derivs <- as.data.frame(gratia::derivatives(gfit, data = scam_data))
     # print("gam"); flush.console()
   }
   # Extract estimated terms
@@ -43,8 +43,8 @@ smooth_scam <- function(x, y, formula, Xcov, ...)
       if (iscons){
         dgx[,j] <- scam::derivative.scam(gfit, jind)$d
       } else {
-        dgx[,j] <- derivs[derivs$smooth == gfit$smooth[[jind]]$label, 
-          "derivative"]
+        dgx[,j] <- derivs[derivs$.smooth == gfit$smooth[[jind]]$label, 
+          ".derivative"]
       }
     } else {
       if (is.numeric(scam_data[j])){
diff --git a/R/vcov.boot.cgaim.R b/R/vcov.boot.cgaim.R
new file mode 100644
index 0000000..11fa738
--- /dev/null
+++ b/R/vcov.boot.cgaim.R
@@ -0,0 +1,39 @@
+#' @rdname vcov.cgaim
+#' @order 2
+#' 
+#' @export
+vcov.boot.cgaim <- function(object, parm = c("alpha", "beta"), complete = TRUE, 
+  ...)
+{
+  #----- Header
+  
+  # Check parm
+  parm <- match.arg(parm)
+  
+  #----- Create variance covariane matrix
+  
+  # Compute vcov
+  vres <- stats::var(t(object$boot[[parm]]))
+    
+  # Set names
+  rownames(vres) <- colnames(vres) <- names(unlist(object$obs[[parm]]))
+ 
+  #----- Complete for potentially aliased coefficients
+  if (isTRUE(complete)){
+    
+    # Check if any aliased coefficients
+    aliased <- is.na(unlist(object$obs[[parm]]))
+    
+    # Complete
+    vresa <- matrix(NA, length(aliased), length(aliased),
+      dimnames = list(names(unlist(object$obs[[parm]])), 
+        names(unlist(object$obs[[parm]]))))
+    vresa[!aliased, !aliased] <- vres
+  
+    #----- Return 
+    return(vresa)
+  } else {
+    return(vres)  
+  }
+  
+}
\ No newline at end of file
diff --git a/R/vcov.cgaim.R b/R/vcov.cgaim.R
new file mode 100644
index 0000000..0da823c
--- /dev/null
+++ b/R/vcov.cgaim.R
@@ -0,0 +1,124 @@
+#' Calculate Variance-Covariance Matrix for a Fitted CGAIM Object
+#'
+#' Returns the variance covariance matrix of the main parameters of a fitted \code{cgaim} object. These parameters correspond to the index weights \code{alpha} and the scaling coefficients \code{beta}.
+#'    
+#' @param object A \code{cgaim} or \code{boot.cgaim} object.
+#' @param parm The model components for which to get confidence intervals. 
+#'    Either \code{"alpha"} (the default) for index weights or \code{"beta"} for scaling coefficients.
+#' @param type The type of confidence intervals. Either \code{"normal"} (the default)
+#'    or \code{"bootstrap"}. See details.
+#' @param B The number of samples to be simulated.
+#' @param complete Indicates whether the full variance-covariance matrix should be returned when some of the parameters could not be estimated. If so, the matrix is padded with \code{NA}s.
+#' @param ... Additional parameters to be passed to \code{\link{boot.cgaim}} for bootstrap replications.
+#'
+#' @details 
+#' Two types of computation are currently implemented in the function.
+#' When \code{type = "normal"}, variance-covariance matrices are computed assuming 
+#' components are normally distributed. Beta coefficients are treated as
+#' regular linear regression coefficients and alpha
+#' coefficients are assumed to follow a Truncated Multivariate Normal distribution. 
+#' The latter is obtained by simulating from TMVN (see \code{\link[TruncatedNormal]{tmvnorm}}) 
+#' and computing the empirical variance covariance matrix from these simulations. The parameter \code{B} controls the number of simulations from the TMVN (and is not used when \code{parm = "beta"}).
+#' 
+#' When \code{type = "bootstrap"}, the variance-covariance matrix is computed on Bootstrap replications. In this case \code{\link{boot.cgaim}} is called internally and \code{B} corresponds to the number of replications. Alternatively, the user can directly call \code{\link{boot.cgaim}} and feed the result into \code{vcov.boot.cgaim} (see examples).
+#' 
+#' @returns A variance-covariance matrix object. 
+#' 
+#' @seealso \code{\link{boot.cgaim}} for bootstrapping and \code{\link{confint.cgaim}} for confidence intervals. 
+#'
+#' @references
+#'   Masselot, P. and others, 2022. Constrained groupwise additive index models.
+#'     Biostatistics.
+#'     
+#'   Pya, N., Wood, S.N., 2015. Shape constrained additive models. 
+#'    Stat. Comput. 25, 543–559. 
+#'    
+#'   Wood, S.N., 2017. Generalized Additive Models: An Introduction with R, 
+#'     2nd ed, Texts in Statistical Science. Chapman and Hall/CRC.
+#'
+#' @examples 
+#' # A simple CGAIM
+#' n <- 200
+#' x1 <- rnorm(n)
+#' x2 <- x1 + rnorm(n)
+#' z <- x1 + x2
+#' y <- z + rnorm(n)
+#' df1 <- data.frame(y, x1, x2) 
+#' ans <- cgaim(y ~ g(x1, x2, acons = list(monotone = 1)), data = df1)
+#' 
+#' # (Truncated) Normal variance-covariance matrix
+#' set.seed(1)
+#' vcov(ans, B = 1000)
+#' set.seed(1)
+#' vcov(ans, parm = "alpha", B = 1000) # Same result
+#' vcov(ans, parm = "beta", B = 1000)
+#' 
+#' # Confidence intervals by bootstrap (more computationally intensive, B should be increased)
+#' set.seed(2)
+#' vcov(ans, type = "boot", B = 10)
+#' 
+#' # Alternatively, bootstrap samples can be performed beforehand
+#' set.seed(2) 
+#' boot1 <- boot.cgaim(ans, B = 10)
+#' vcov(boot1)
+#'
+#' @order 1
+#' @export
+vcov.cgaim <- function(object, parm = c("alpha", "beta"), 
+  type = c("normal", "bootstrap"), B = 100, complete = TRUE, ...)
+{
+  #----- Header
+  
+  # Check parm
+  parm <- match.arg(parm)
+  
+  # Check type
+  type <- match.arg(type)
+  
+  #----- Compute vcov if "normal"
+  if (type == "normal"){
+    # For alpha: simulate from truncated normal
+    if(parm == "alpha"){
+      
+      # Simulate alpha from truncated multivariate normal
+      simures <- simul_tmvnorm(object, B = B)
+      
+      # Compute vcov
+      vres <- stats::var(simures)
+    }
+    
+    # For beta
+    if (parm == "beta"){
+      
+      # use the usual vcov matrix
+      vres <- vcov_beta(object)
+    }
+    
+    # Set names
+    rownames(vres) <- colnames(vres) <- names(unlist(object[[parm]]))
+    
+    # Fill potentially for aliased coefficients
+    if (isTRUE(complete)){
+      
+      # Check if any aliased coefficients
+      aliased <- is.na(unlist(object[[parm]]))
+      
+      # Complete
+      vresa <- matrix(NA, length(aliased), length(aliased),
+        dimnames = list(names(unlist(object[[parm]])), 
+          names(unlist(object[[parm]]))))
+      vresa[!aliased, !aliased] <- vres
+      vres <- vresa
+    }
+  } else {
+    
+    # Simulate
+    simures <- boot.cgaim(object, B = B, ...)
+    
+    # Compute CI
+    vres <- vcov.boot.cgaim(simures, parm = parm, complete = complete)
+  }
+  
+  #----- Return
+  return(vres)
+}
\ No newline at end of file
diff --git a/README.md b/README.md
index 37c0e27..cd04777 100644
--- a/README.md
+++ b/README.md
@@ -28,5 +28,5 @@ The package is available for installation from the usual CRAN repository. Altern
 
 ## Functions
 
-The main function of the package is the eponymous `cgaim` that fits the model. Then, the `print.cgaim` method allows displaying the results. The `confint.cgaim` method allows computing confidence intervals, either using normal approximation or through bootstrap. Ridge functions can be displayed using the `plot.cgaim` method, with the possibility to add confidence intervals. Finally, `predict.cgaim` allows computing the indices and perform prediction on new observations.
+The main function of the package is the eponymous `cgaim` that fits the model. Then, the `print.cgaim` method allows displaying the results. The `confint.cgaim` and `vcov.cgaim` methods allow computing confidence intervals and variance-covariance matrix using various approaches. Ridge functions can be displayed using the `plot.cgaim` method, with the possibility to add confidence intervals. Finally, `predict.cgaim` allows computing the indices and perform prediction on new observations.
 See the help of each function.
diff --git a/cran-comments.md b/cran-comments.md
index a6a4ff6..75a6517 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -1,35 +1,9 @@
-# 1.0.0
+# 1.0.1
 
-## Resubmission
-This is a resubmission. In this version I have:
-
-* Added a \donttest statement for R CMD check to ignore the most computationally expensive example.
-* Added the reference to the methodological paper in the DESCRIPTION file.
-* Replaced the T and F that previously escaped me by TRUE and FALSE.
-* Added \value fields to the documentation of exported functions missing one.
+Package update after check errors from 27/04/2024. Fixes error induced by the update of the `gratia` dependency.
 
 ## R CMD check results
-There were no ERRORs or WARNINGs.
-
-There was one NOTE on winbuilder only:
-
-* checking CRAN incoming feasibility ... NOTE
-Maintainer: 'Pierre Masselot <pierre.masselot@lshtm.ac.uk>'
-
-New submission
-
-Possibly misspelled words in DESCRIPTION:
-  Campagna (10:182)
-  Chebana (10:173)
-  Gosselin (10:209)
-  Groupwise (2:20)
-  Lavigne (10:192)
-  Masselot (10:163)
-  Ouarda (10:201)
-  groupwise (10:31, 10:238)
-
-I believe these notes are just flagging the fact it is a new submission.
-Regarding misspelled words, this is how these words are spelled and used in associated publications and are therefore correct (groupwise) and the names of authors in the citation.
+0 errors | 0 warnings | 0 notes
 
 ## Downstream dependencies
-There are currently no downstream dependencies for this package.
+Only bug fix and addition of new convenience functions, no change in the specification nor output of the existing function.
diff --git a/inst/CITATION b/inst/CITATION
index 1ea7c95..62f6fde 100644
--- a/inst/CITATION
+++ b/inst/CITATION
@@ -1,8 +1,8 @@
 citHeader("To cite cgaim in publications use:")
 
-citEntry(entry = "Article",
+bibentry(bibtype = "Article",
   title        = "Constrained groupwise additive index models",
-  author       = personList(as.person("Pierre Masselot"),
+  author       = c(as.person("Pierre Masselot"),
                    as.person("Fateh Chebana"),
                    as.person("Céline Campagna"),
                    as.person("Éric Lavigne"),
@@ -10,14 +10,14 @@ citEntry(entry = "Article",
                    as.person("Pierre Gosselin")),
   journal      = "Biostatistics",
   year         = "2022",
-  volume       = "00",
-  number       = "00",
-  pages        = "1--19",
+  volume       = "24",
+  number       = "4",
+  pages        = "1066–1084",
   url          = "https://doi.org/10.1093/biostatistics/kxac023",
 
   textVersion  =
-  paste("Pierre Masselot, Fateh Chebana, Céline Campagna, Éric Lavigne, Taha B M J Ouarda, Pierre Gosselin (2022).",
+  paste("Pierre Masselot, Fateh Chebana, Céline Campagna, Éric Lavigne, Taha B.M.J. Ouarda, Pierre Gosselin (2022).",
         "Constrained groupwise additive index models.",
-        "Biostatistics, 00(00), 1-19.",
+        "Biostatistics, 24(4), 1066–1084.",
         "URL https://doi.org/10.1093/biostatistics/kxac023.")
 )
diff --git a/man/vcov.cgaim.Rd b/man/vcov.cgaim.Rd
new file mode 100644
index 0000000..cc16b94
--- /dev/null
+++ b/man/vcov.cgaim.Rd
@@ -0,0 +1,85 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/vcov.cgaim.R, R/vcov.boot.cgaim.R
+\name{vcov.cgaim}
+\alias{vcov.cgaim}
+\alias{vcov.boot.cgaim}
+\title{Calculate Variance-Covariance Matrix for a Fitted CGAIM Object}
+\usage{
+\method{vcov}{cgaim}(object, parm = c("alpha", "beta"), type = c("normal",
+  "bootstrap"), B = 100, complete = TRUE, ...)
+
+\method{vcov}{boot.cgaim}(object, parm = c("alpha", "beta"),
+  complete = TRUE, ...)
+}
+\arguments{
+\item{object}{A \code{cgaim} or \code{boot.cgaim} object.}
+
+\item{parm}{The model components for which to get confidence intervals. 
+Either \code{"alpha"} (the default) for index weights or \code{"beta"} for scaling coefficients.}
+
+\item{type}{The type of confidence intervals. Either \code{"normal"} (the default)
+or \code{"bootstrap"}. See details.}
+
+\item{B}{The number of samples to be simulated.}
+
+\item{complete}{Indicates whether the full variance-covariance matrix should be returned when some of the parameters could not be estimated. If so, the matrix is padded with \code{NA}s.}
+
+\item{...}{Additional parameters to be passed to \code{\link{boot.cgaim}} for bootstrap replications.}
+}
+\value{
+A variance-covariance matrix object.
+}
+\description{
+Returns the variance covariance matrix of the main parameters of a fitted \code{cgaim} object. These parameters correspond to the index weights \code{alpha} and the scaling coefficients \code{beta}.
+}
+\details{
+Two types of computation are currently implemented in the function.
+When \code{type = "normal"}, variance-covariance matrices are computed assuming 
+components are normally distributed. Beta coefficients are treated as
+regular linear regression coefficients and alpha
+coefficients are assumed to follow a Truncated Multivariate Normal distribution. 
+The latter is obtained by simulating from TMVN (see \code{\link[TruncatedNormal]{tmvnorm}}) 
+and computing the empirical variance covariance matrix from these simulations. The parameter \code{B} controls the number of simulations from the TMVN (and is not used when \code{parm = "beta"}).
+
+When \code{type = "bootstrap"}, the variance-covariance matrix is computed on Bootstrap replications. In this case \code{\link{boot.cgaim}} is called internally and \code{B} corresponds to the number of replications. Alternatively, the user can directly call \code{\link{boot.cgaim}} and feed the result into \code{vcov.boot.cgaim} (see examples).
+}
+\examples{
+# A simple CGAIM
+n <- 200
+x1 <- rnorm(n)
+x2 <- x1 + rnorm(n)
+z <- x1 + x2
+y <- z + rnorm(n)
+df1 <- data.frame(y, x1, x2) 
+ans <- cgaim(y ~ g(x1, x2, acons = list(monotone = 1)), data = df1)
+
+# (Truncated) Normal variance-covariance matrix
+set.seed(1)
+vcov(ans, B = 1000)
+set.seed(1)
+vcov(ans, parm = "alpha", B = 1000) # Same result
+vcov(ans, parm = "beta", B = 1000)
+
+# Confidence intervals by bootstrap (more computationally intensive, B should be increased)
+set.seed(2)
+vcov(ans, type = "boot", B = 10)
+
+# Alternatively, bootstrap samples can be performed beforehand
+set.seed(2) 
+boot1 <- boot.cgaim(ans, B = 10)
+vcov(boot1)
+
+}
+\references{
+Masselot, P. and others, 2022. Constrained groupwise additive index models.
+    Biostatistics.
+    
+  Pya, N., Wood, S.N., 2015. Shape constrained additive models. 
+   Stat. Comput. 25, 543–559. 
+   
+  Wood, S.N., 2017. Generalized Additive Models: An Introduction with R, 
+    2nd ed, Texts in Statistical Science. Chapman and Hall/CRC.
+}
+\seealso{
+\code{\link{boot.cgaim}} for bootstrapping and \code{\link{confint.cgaim}} for confidence intervals.
+}
diff --git a/tests/testthat/test-confint.R b/tests/testthat/test-confint.R
deleted file mode 100644
index 6a1c655..0000000
--- a/tests/testthat/test-confint.R
+++ /dev/null
@@ -1,63 +0,0 @@
-
-#----- A model
-set.seed(2020)
-n <- 200
-x1 <- rnorm(n)
-x2 <- rnorm(n)
-x3 <- rnorm(n)
-x4 <- rnorm(n)
-x5 <- rnorm(n)
-mu <- 4 * exp(8 * x1) / (1 + exp(8 * x1)) + exp(x3) + sin(x4 + x5)
-y <- mu + rnorm(n)
-df3 <- data.frame(y, x1, x2, x3, x4, x5)
-
-ans <- cgaim(y ~ g(x1, x2, label = "i1", s_opts = list(sp = 0)) + 
-    s(x3, s_opts = list(bs = "cr")) + g(x4, x5, label = "i2", fcons = "inc"),
-  data = df3)
-
-#----- Check confint
-
-test_that("both methods work", {
-  cin <- confint(ans, type = "norm", B = 10)
-  expect_equal(dim(cin$alpha), c(4, 2))
-  expect_equal(dim(cin$beta), c(4, 2))
-  
-  cib <- confint(ans, type = "boot", B = 10)
-  expect_equal(dim(cib$alpha), c(4, 2))
-  expect_equal(dim(cib$beta), c(4, 2))
-})
-
-test_that("parm argument works", {
-  cin <- confint(ans, parm = 1, type = "norm", B = 10)
-  expect_length(cin, 1)
-  expect_named(cin, "alpha")
-  
-  cin <- confint(ans, parm = 2, type = "norm", B = 10)
-  expect_length(cin, 1)
-  expect_named(cin, "beta")
-  
-  cin <- confint(ans, parm = c("beta", "alpha"), type = "norm", B = 10)
-  expect_length(cin, 2)
-  expect_named(cin, c("alpha", "beta"))
-  
-  cib <- confint(ans, parm = c("g", "alpha"), type = "boot", B = 10)
-  expect_length(cib, 2)
-  expect_named(cib, c("alpha", "g"))
-})
-
-test_that("both confint methods are equivalent", {
-  
-  # Calling confint.cgaim directly
-  set.seed(1)
-  cib1 <- confint(ans, type = "boot", B = 10)
-  
-  # Creating replications first
-  set.seed(1)
-  bans <- boot.cgaim(ans, B = 10)
-  cib2 <- confint(bans)
-  
-  # Check they are identical
-  expect_equal(cib1, cib2)
-})
-
-
diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R
new file mode 100644
index 0000000..a8da407
--- /dev/null
+++ b/tests/testthat/test-uncertainty.R
@@ -0,0 +1,115 @@
+
+#----- A model
+set.seed(2020)
+n <- 200
+x1 <- rnorm(n)
+x2 <- rnorm(n)
+x3 <- rnorm(n)
+x4 <- rnorm(n)
+x5 <- rnorm(n)
+mu <- 4 * exp(8 * x1) / (1 + exp(8 * x1)) + exp(x3) + sin(x4 + x5)
+y <- mu + rnorm(n)
+df3 <- data.frame(y, x1, x2, x3, x4, x5)
+
+ans <- cgaim(y ~ g(x1, x2, label = "i1", s_opts = list(sp = 0)) + 
+    s(x3, s_opts = list(bs = "cr")) + g(x4, x5, label = "i2", fcons = "inc"),
+  data = df3)
+
+#----- Check confint
+
+test_that("both methods work", {
+  cin <- confint(ans, type = "norm", B = 10)
+  expect_equal(dim(cin$alpha), c(4, 2))
+  expect_equal(dim(cin$beta), c(4, 2))
+  
+  cib <- confint(ans, type = "boot", B = 10)
+  expect_equal(dim(cib$alpha), c(4, 2))
+  expect_equal(dim(cib$beta), c(4, 2))
+})
+
+test_that("parm argument works", {
+  cin <- confint(ans, parm = 1, type = "norm", B = 10)
+  expect_length(cin, 1)
+  expect_named(cin, "alpha")
+  
+  cin <- confint(ans, parm = 2, type = "norm", B = 10)
+  expect_length(cin, 1)
+  expect_named(cin, "beta")
+  
+  cin <- confint(ans, parm = c("beta", "alpha"), type = "norm", B = 10)
+  expect_length(cin, 2)
+  expect_named(cin, c("alpha", "beta"))
+  
+  cib <- confint(ans, parm = c("g", "alpha"), type = "boot", B = 10)
+  expect_length(cib, 2)
+  expect_named(cib, c("alpha", "g"))
+})
+
+test_that("both confint methods are equivalent", {
+  
+  # Calling confint.cgaim directly
+  set.seed(1)
+  cib1 <- confint(ans, type = "boot", B = 10)
+  
+  # Creating replications first
+  set.seed(1)
+  bans <- boot.cgaim(ans, B = 10)
+  cib2 <- confint(bans)
+  
+  # Check they are identical
+  expect_equal(cib1, cib2)
+})
+
+#----- Check vcov
+
+test_that("vcov methods work", {
+  v <- vcov(ans, type = "norm", B = 10)
+  expect_equal(colnames(v), names(unlist(ans$alpha)))
+  
+  v <- vcov(ans, type = "boot", B = 10)
+  expect_equal(colnames(v), names(unlist(ans$alpha)))
+})
+
+test_that("vcov parameter specification works", {
+  v <- vcov(ans); expect_equal(colnames(v), names(unlist(ans$alpha)))
+  v <- vcov(ans, parm = "alpha"); expect_equal(colnames(v), 
+    names(unlist(ans$alpha)))
+  v <- vcov(ans, parm = "beta"); expect_equal(colnames(v), names(ans$beta))
+  expect_error(vcov(ans, parm = "kappa"))
+  v <- vcov(ans, parm = "a"); expect_equal(colnames(v), 
+    names(unlist(ans$alpha)))
+  v <- vcov(ans, parm = "b"); expect_equal(colnames(v), names(ans$beta))
+  v <- vcov(ans, parm = NULL); expect_equal(colnames(v), 
+    names(unlist(ans$alpha)))
+})
+
+test_that("both bootstrap vcov methods are equivalent", {
+  
+  # Calling vcov.cgaim directly
+  set.seed(1)
+  v1 <- vcov(ans, type = "boot", B = 10)
+  
+  # Creating replications first
+  set.seed(1)
+  bans <- boot.cgaim(ans, B = 10)
+  v2 <- vcov(bans)
+  
+  # Check they are identical
+  expect_equal(v1, v2)
+})
+
+test_that("vcov with constraints work", {
+  
+  # new model with constraints on alpha coefficients
+  ans2 <- cgaim(y ~ g(x1, x2, label = "i1", s_opts = list(sp = 0), 
+      acons = list(monotone = -1)) + s(x3, s_opts = list(bs = "cr")) + 
+      g(x4, x5, label = "i2", fcons = "inc", acons = list(sign = 1)),
+    data = df3)
+  
+  # Check we get results
+  v <- vcov(ans2, type = "norm", B = 10)
+  expect_equal(colnames(v), names(unlist(ans$alpha)))
+  
+  v <- vcov(ans2, type = "boot", B = 10)
+  expect_equal(colnames(v), names(unlist(ans$alpha)))
+})
\ No newline at end of file
diff --git a/tests/testthat/test_updateAlpha.R b/tests/testthat/test_updateAlpha.R
index 324c68a..dfaaf85 100644
--- a/tests/testthat/test_updateAlpha.R
+++ b/tests/testthat/test_updateAlpha.R
@@ -22,8 +22,8 @@ alpha[3:4] <- cgaim:::normalize(alpha[3:4], "1")
 dfgam <- data.frame(y, i1 = cbind(x1, x2) %*% alpha[1:2], 
   i2 = cbind(x3, x4) %*% alpha[3:4])
 gamres <- mgcv::gam(y ~ s(i1) + s(i2), data = dfgam)
-derivs <- gratia::derivatives(gamres, newdata = dfgam)
-dgz <- matrix(derivs[["derivative"]], ncol = 2)
+derivs <- gratia::derivatives(gamres, data = dfgam)
+dgz <- matrix(derivs[[".derivative"]], ncol = 2)
 gz <- stats::predict(gamres, type = "terms")
 
 # Test with constraints