diff --git a/NAMESPACE b/NAMESPACE
index 2bd09cb8e..ddf06cac1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -215,5 +215,4 @@ importFrom(utils,sessionInfo)
 importFrom(utils,stack)
 importFrom(utils,str)
 importFrom(utils,tail)
-importFrom(vpc,vpc)
 useDynLib(nlmixr, .registration=TRUE)
diff --git a/R/vpc.R b/R/vpc.R
index 4b5369307..5d36e8d22 100644
--- a/R/vpc.R
+++ b/R/vpc.R
@@ -56,10 +56,13 @@ vpc <- function (sim, ...)
 vpc.default <- function(sim, ...){
     ns <- loadNamespace("vpc");
     if (exists("vpc_vpc",ns)){
-        vpc::vpc_vpc(sim, ...)
+        vpcn <- "vpc_vpc"
     } else {
-        vpc::vpc(sim, ...)
+        vpcn <- "vpc"
     }
+    call <- as.list(match.call(expand.dots=TRUE))[-1];
+    call <- call[names(call) %in% methods::formalArgs(getFromNamespace(vpcn,"vpc"))]
+    p = do.call(getFromNamespace(vpcn,"vpc"), call, envir = parent.frame(1))
 }
 
 #' Visual predictive check (VPC) for nlmixr nlme objects
@@ -69,6 +72,7 @@ vpc.default <- function(sim, ...){
 #' @param fit nlme fit object
 #' @param nsim number of simulations
 #' @param condition conditional variable
+#' @param ... Additional arguments
 #' @inheritParams vpc::vpc
 #' @return NULL
 #' @examples
@@ -76,7 +80,7 @@ vpc.default <- function(sim, ...){
 #' fit <- nlme_lin_cmpt(theo_md, par_model=specs, ncmt=1, verbose=TRUE)
 #' vpc(fit, nsim = 100, condition = NULL)
 #' @export
-vpc_nlmixr_nlme = function(fit, nsim=100, condition=NULL)
+vpc_nlmixr_nlme = function(fit, nsim=100, condition=NULL, ...)
 {
     nlmeModList(fit$env);
     on.exit({nlmeModList(new.env(parent=emptyenv()))})
diff --git a/R/vpc.nlme.R b/R/vpc.nlme.R
index 64a054ae1..b8af61a79 100644
--- a/R/vpc.nlme.R
+++ b/R/vpc.nlme.R
@@ -1,4 +1,4 @@
-vpc.nlme = function (fit, nsim = 100, by = NULL, ...) 
+vpc.nlme = function (fit, nsim = 100, by = NULL, ...)
 {
     if (class(fit) == "nlmixr.ui.nlme") fit = as.nlme(fit)
     dat = getData(fit)
@@ -30,11 +30,19 @@ vpc.nlme = function (fit, nsim = 100, by = NULL, ...)
     xs = do.call("cbind", s)
 
     df = cbind(xd[ord, c("ID", "TIME", "grp")], DV = as.vector(xs), SIM = sim)
-    if (!is.null(by)) {
-        p = vpc::vpc(sim = df, obs = dat, strat = c("grp"), facet = "wrap", ...)
+    ns <- loadNamespace("vpc");
+    if (exists("vpc_vpc",ns)){
+        vpcn <- "vpc_vpc"
     } else {
-        p = vpc::vpc(sim = df, obs = dat, ...)
+        vpcn <- "vpc"
+    }
+    call <- as.list(match.call(expand.dots=TRUE))[-1];
+    if (!is.null(by)) {
+        call$strat <- c("grp")
+        call$facet <- "wrap";
     }
+    call <- call[names(call) %in% methods::formalArgs(getFromNamespace(vpcn,"vpc"))]
+    p = do.call(getFromNamespace(vpcn,"vpc"), c(list(sim=df, obs=dat), call), envir = parent.frame(1))
     print(p)
     invisible(df)
 }
diff --git a/R/vpc.saemFit.R b/R/vpc.saemFit.R
index 117f2c069..265e838a4 100644
--- a/R/vpc.saemFit.R
+++ b/R/vpc.saemFit.R
@@ -57,11 +57,19 @@ vpc_saemFit = function(fit, dat, nsim = 100, by=NULL, ...) {
   xs = do.call("cbind",s)
   df = cbind(xd[ord, c("ID", "TIME", "grp")], DV=as.vector(xs), SIM=sim)
 
-  if (!is.null(by)) {
-      p = vpc::vpc(sim = df, obs = dat, strat=c("grp"), facet="wrap", ...)
+  ns <- loadNamespace("vpc");
+  if (exists("vpc_vpc",ns)){
+      vpcn <- "vpc_vpc"
   } else {
-      p = vpc::vpc(sim = df, obs = dat, ...)
+      vpcn <- "vpc"
+  }
+  call <- as.list(match.call(expand.dots=TRUE))[-1];
+  if (!is.null(by)) {
+      call$strat <- c("grp")
+      call$facet <- "wrap";
   }
+  call <- call[names(call) %in% methods::formalArgs(getFromNamespace(vpcn,"vpc"))]
+  p = do.call(getFromNamespace(vpcn,"vpc"), c(list(sim=df, obs=dat), call), envir = parent.frame(1))
   print(p)
 
   invisible(df)
diff --git a/R/vpc.ui.R b/R/vpc.ui.R
index 92e6d19db..5282b11a9 100644
--- a/R/vpc.ui.R
+++ b/R/vpc.ui.R
@@ -9,7 +9,7 @@
 ##' @inheritParams vpc::vpc
 ##' @inheritParams RxODE::rxSolve
 ##' @param ... Args sent to \code{\link[RxODE]{rxSolve}}
-##' @inheritParams vpc::vpc
+##' @return Simulated dataset (invisibly)
 ##' @author Matthew L. Fidler
 ##' @export
 vpc_ui <- function(fit, data=NULL, n=100, bins = "jenks",
@@ -100,7 +100,9 @@ vpc_ui <- function(fit, data=NULL, n=100, bins = "jenks",
     call$obs_cols = list(id="id", dv="dv", idv="time")
     call$sim_cols = list(id="id", dv="dv", idv="time")
     call$stratify = stratify
-    do.call(getFromNamespace(vpcn,"vpc"), c(list(sim=sim, obs=dat), call), envir = parent.frame(1))
+    p = do.call(getFromNamespace(vpcn,"vpc"), c(list(sim=sim, obs=dat), call), envir = parent.frame(1))
+    print(p);
+    return(invisible(sim));
 }
 
 
@@ -128,7 +130,3 @@ vpc.nlmixr.ui.nlme <- function(sim, ...){
 vpc.ui <- function(sim, ...){
     vpc_ui(fit=sim, ...);
 }
-
-##' @importFrom vpc vpc
-##' @export
-vpc::vpc
diff --git a/man/reexports.Rd b/man/reexports.Rd
index ca37ee98d..07c61b740 100644
--- a/man/reexports.Rd
+++ b/man/reexports.Rd
@@ -1,11 +1,9 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/simulate.R, R/vpc.ui.R
+% Please edit documentation in R/simulate.R
 \docType{import}
 \name{reexports}
 \alias{reexports}
 \alias{rxSolve}
-\alias{reexports}
-\alias{vpc}
 \title{Objects exported from other packages}
 \keyword{internal}
 \description{
@@ -14,7 +12,5 @@ below to see their documentation.
 
 \describe{
   \item{RxODE}{\code{\link[RxODE]{rxSolve}}}
-
-  \item{vpc}{\code{\link[vpc]{vpc}}}
 }}
 
diff --git a/man/vpc.Rd b/man/vpc.Rd
new file mode 100644
index 000000000..ef0be59c4
--- /dev/null
+++ b/man/vpc.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/vpc.R
+\name{vpc}
+\alias{vpc}
+\alias{vpc.default}
+\title{Vpc function for nlmixr}
+\usage{
+vpc(sim, ...)
+
+\method{vpc}{default}(sim, ...)
+}
+\arguments{
+\item{sim}{Observed data frame or fit object}
+
+\item{...}{Other parameters}
+}
+\description{
+Vpc function for nlmixr
+}
diff --git a/man/vpc_nlmixr_nlme.Rd b/man/vpc_nlmixr_nlme.Rd
index 087b100ef..681255deb 100644
--- a/man/vpc_nlmixr_nlme.Rd
+++ b/man/vpc_nlmixr_nlme.Rd
@@ -5,7 +5,7 @@
 \alias{vpc.nlmixr_nlme}
 \title{Visual predictive check (VPC) for nlmixr nlme objects}
 \usage{
-vpc_nlmixr_nlme(fit, nsim = 100, condition = NULL)
+vpc_nlmixr_nlme(fit, nsim = 100, condition = NULL, ...)
 
 \method{vpc}{nlmixr_nlme}(sim, ...)
 }
@@ -16,6 +16,8 @@ vpc_nlmixr_nlme(fit, nsim = 100, condition = NULL)
 
 \item{condition}{conditional variable}
 
+\item{...}{Additional arguments}
+
 \item{sim}{a data.frame with observed data, containing the independent and dependent variable, a column indicating the individual, and possibly covariates. E.g. load in from NONMEM using \link{read_table_nm}}
 }
 \description{
diff --git a/man/vpc_ui.Rd b/man/vpc_ui.Rd
index e06a49406..895ca00cb 100644
--- a/man/vpc_ui.Rd
+++ b/man/vpc_ui.Rd
@@ -96,6 +96,9 @@ parameters.}
 
 \item{sim}{a data.frame with observed data, containing the independent and dependent variable, a column indicating the individual, and possibly covariates. E.g. load in from NONMEM using \link{read_table_nm}}
 }
+\value{
+Simulated dataset (invisibly)
+}
 \description{
 VPC based on ui model
 }