diff --git a/DESCRIPTION b/DESCRIPTION
index dbc7b8b..e13aeb3 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -23,7 +23,7 @@ Description: 'Monolix' is a tool for running mixed effects model using
 License: MIT + file LICENSE
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 LinkingTo:
     dparser,
     rxode2parse,
diff --git a/NAMESPACE b/NAMESPACE
index 72b0a87..99c7cff 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -69,6 +69,7 @@ export(logit)
 export(lotri)
 export(mlxTxt)
 export(mlxtran)
+export(mlxtranGetMutate)
 export(model)
 export(monolix2rx)
 export(monolixDataImport)
diff --git a/R/rxSolve.R b/R/rxSolve.R
index 708d610..5d6eb09 100644
--- a/R/rxSolve.R
+++ b/R/rxSolve.R
@@ -5,33 +5,36 @@ rxSolve.monolix2rx <- function(object, params = NULL, events = NULL,
         "dop853", "indLin"), sigdig = NULL, atol = 1e-08, rtol = 1e-06, 
     maxsteps = 70000L, hmin = 0, hmax = NA_real_, hmaxSd = 0, 
     hini = 0, maxordn = 12L, maxords = 5L, ..., cores, covsInterpolation = c("locf", 
-        "linear", "nocb", "midpoint"), addCov = TRUE, sigma = NULL, 
-    sigmaDf = NULL, sigmaLower = -Inf, sigmaUpper = Inf, nCoresRV = 1L, 
-    sigmaIsChol = FALSE, sigmaSeparation = c("auto", "lkj", "separation"), 
-    sigmaXform = c("identity", "variance", "log", "nlmixrSqrt", 
-        "nlmixrLog", "nlmixrIdentity"), nDisplayProgress = 10000L, 
-    amountUnits = NA_character_, timeUnits = "hours", theta = NULL, 
-    thetaLower = -Inf, thetaUpper = Inf, eta = NULL, addDosing = FALSE, 
-    stateTrim = Inf, updateObject = FALSE, omega = NULL, omegaDf = NULL, 
-    omegaIsChol = FALSE, omegaSeparation = c("auto", "lkj", "separation"), 
-    omegaXform = c("variance", "identity", "log", "nlmixrSqrt", 
-        "nlmixrLog", "nlmixrIdentity"), omegaLower = -Inf, omegaUpper = Inf, 
-    nSub = 1L, thetaMat = NULL, thetaDf = NULL, thetaIsChol = FALSE, 
-    nStud = 1L, dfSub = 0, dfObs = 0, returnType = c("rxSolve", 
-        "matrix", "data.frame", "data.frame.TBS", "data.table", 
-        "tbl", "tibble"), seed = NULL, nsim = NULL, minSS = 10L, 
-    maxSS = 1000L, infSSstep = 12, strictSS = TRUE, istateReset = TRUE, 
-    subsetNonmem = TRUE, maxAtolRtolFactor = 0.1, from = NULL, 
-    to = NULL, by = NULL, length.out = NULL, iCov = NULL, keep = NULL, 
-    indLinPhiTol = 1e-07, indLinPhiM = 0L, indLinMatExpType = c("expokit", 
-        "Al-Mohy", "arma"), indLinMatExpOrder = 6L, drop = NULL, 
-    idFactor = TRUE, mxhnil = 0, hmxi = 0, warnIdSort = TRUE, 
-    warnDrop = TRUE, ssAtol = 1e-08, ssRtol = 1e-06, safeZero = TRUE, 
-    sumType = c("pairwise", "fsum", "kahan", "neumaier", "c"), 
-    prodType = c("long double", "double", "logify"), sensType = c("advan", 
-        "autodiff", "forward", "central"), linDiff = c(tlag = 1.5e-05, 
-        f = 1.5e-05, rate = 1.5e-05, dur = 1.5e-05, tlag2 = 1.5e-05, 
-        f2 = 1.5e-05, rate2 = 1.5e-05, dur2 = 1.5e-05), linDiffCentral = c(tlag = TRUE, 
+        "linear", "nocb", "midpoint"), naInterpolation = c("locf", 
+        "nocb"), keepInterpolation = c("locf", "nocb", "na"), 
+    addCov = TRUE, sigma = NULL, sigmaDf = NULL, sigmaLower = -Inf, 
+    sigmaUpper = Inf, nCoresRV = 1L, sigmaIsChol = FALSE, sigmaSeparation = c("auto", 
+        "lkj", "separation"), sigmaXform = c("identity", "variance", 
+        "log", "nlmixrSqrt", "nlmixrLog", "nlmixrIdentity"), 
+    nDisplayProgress = 10000L, amountUnits = NA_character_, timeUnits = "hours", 
+    theta = NULL, thetaLower = -Inf, thetaUpper = Inf, eta = NULL, 
+    addDosing = FALSE, stateTrim = Inf, updateObject = FALSE, 
+    omega = NULL, omegaDf = NULL, omegaIsChol = FALSE, omegaSeparation = c("auto", 
+        "lkj", "separation"), omegaXform = c("variance", "identity", 
+        "log", "nlmixrSqrt", "nlmixrLog", "nlmixrIdentity"), 
+    omegaLower = -Inf, omegaUpper = Inf, nSub = 1L, thetaMat = NULL, 
+    thetaDf = NULL, thetaIsChol = FALSE, nStud = 1L, dfSub = 0, 
+    dfObs = 0, returnType = c("rxSolve", "matrix", "data.frame", 
+        "data.frame.TBS", "data.table", "tbl", "tibble"), seed = NULL, 
+    nsim = NULL, minSS = 10L, maxSS = 1000L, infSSstep = 12, 
+    strictSS = TRUE, istateReset = TRUE, subsetNonmem = TRUE, 
+    maxAtolRtolFactor = 0.1, from = NULL, to = NULL, by = NULL, 
+    length.out = NULL, iCov = NULL, keep = NULL, indLinPhiTol = 1e-07, 
+    indLinPhiM = 0L, indLinMatExpType = c("expokit", "Al-Mohy", 
+        "arma"), indLinMatExpOrder = 6L, drop = NULL, idFactor = TRUE, 
+    mxhnil = 0, hmxi = 0, warnIdSort = TRUE, warnDrop = TRUE, 
+    ssAtol = 1e-08, ssRtol = 1e-06, safeZero = TRUE, safeLog = TRUE, 
+    safePow = TRUE, sumType = c("pairwise", "fsum", "kahan", 
+        "neumaier", "c"), prodType = c("long double", "double", 
+        "logify"), sensType = c("advan", "autodiff", "forward", 
+        "central"), linDiff = c(tlag = 1.5e-05, f = 1.5e-05, 
+        rate = 1.5e-05, dur = 1.5e-05, tlag2 = 1.5e-05, f2 = 1.5e-05, 
+        rate2 = 1.5e-05, dur2 = 1.5e-05), linDiffCentral = c(tlag = TRUE, 
         f = TRUE, rate = TRUE, dur = TRUE, tlag2 = TRUE, f2 = TRUE, 
         rate2 = TRUE, dur2 = TRUE), resample = NULL, resampleID = TRUE, 
     maxwhile = 1e+05, atolSens = 1e-08, rtolSens = 1e-06, ssAtolSens = 1e-08, 
@@ -118,6 +121,7 @@ rxSolve.monolix2rx <- function(object, params = NULL, events = NULL,
         atol = atol, rtol = rtol, maxsteps = maxsteps, hmin = hmin, 
         hmax = hmax, hmaxSd = hmaxSd, hini = hini, maxordn = maxordn, 
         maxords = maxords, ..., cores = cores, covsInterpolation = covsInterpolation, 
+        naInterpolation = naInterpolation, keepInterpolation = keepInterpolation, 
         addCov = addCov, sigma = sigma, sigmaDf = sigmaDf, sigmaLower = sigmaLower, 
         sigmaUpper = sigmaUpper, nCoresRV = nCoresRV, sigmaIsChol = sigmaIsChol, 
         sigmaSeparation = sigmaSeparation, sigmaXform = sigmaXform, 
@@ -138,13 +142,14 @@ rxSolve.monolix2rx <- function(object, params = NULL, events = NULL,
         indLinMatExpType = indLinMatExpType, indLinMatExpOrder = indLinMatExpOrder, 
         drop = drop, idFactor = idFactor, mxhnil = mxhnil, hmxi = hmxi, 
         warnIdSort = warnIdSort, warnDrop = warnDrop, ssAtol = ssAtol, 
-        ssRtol = ssRtol, safeZero = safeZero, sumType = sumType, 
-        prodType = prodType, sensType = sensType, linDiff = linDiff, 
-        linDiffCentral = linDiffCentral, resample = resample, 
-        resampleID = resampleID, maxwhile = maxwhile, atolSens = atolSens, 
-        rtolSens = rtolSens, ssAtolSens = ssAtolSens, ssRtolSens = ssRtolSens, 
-        simVariability = simVariability, nLlikAlloc = nLlikAlloc, 
-        useStdPow = useStdPow, naTimeHandle = naTimeHandle, addlKeepsCov = addlKeepsCov, 
-        addlDropSs = addlDropSs, ssAtDoseTime = ssAtDoseTime, 
-        ss2cancelAllPending = ss2cancelAllPending, envir = envir)
+        ssRtol = ssRtol, safeZero = safeZero, safeLog = safeLog, 
+        safePow = safePow, sumType = sumType, prodType = prodType, 
+        sensType = sensType, linDiff = linDiff, linDiffCentral = linDiffCentral, 
+        resample = resample, resampleID = resampleID, maxwhile = maxwhile, 
+        atolSens = atolSens, rtolSens = rtolSens, ssAtolSens = ssAtolSens, 
+        ssRtolSens = ssRtolSens, simVariability = simVariability, 
+        nLlikAlloc = nLlikAlloc, useStdPow = useStdPow, naTimeHandle = naTimeHandle, 
+        addlKeepsCov = addlKeepsCov, addlDropSs = addlDropSs, 
+        ssAtDoseTime = ssAtDoseTime, ss2cancelAllPending = ss2cancelAllPending, 
+        envir = envir)
 }
diff --git a/man/mlxtranGetMutate.Rd b/man/mlxtranGetMutate.Rd
new file mode 100644
index 0000000..feaf1c3
--- /dev/null
+++ b/man/mlxtranGetMutate.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/dataCov.R
+\name{mlxtranGetMutate}
+\alias{mlxtranGetMutate}
+\title{Get the string of the mutate statement for input dataset based on
+mlxtran}
+\usage{
+mlxtranGetMutate(mlxtran)
+}
+\arguments{
+\item{mlxtran}{input mlxtran file}
+}
+\value{
+mlxtran string that can be applied to a model (by evaluating it)
+}
+\description{
+Get the string of the mutate statement for input dataset based on
+mlxtran
+}
+\examples{
+
+covD <- system.file("cov", package="monolix2rx")
+m <- mlxtran(file.path(covD, "phenobarbital_project.mlxtran"))
+message(mlxtranGetMutate(m))
+}
+\author{
+Matthew L. Fidler
+}
diff --git a/man/monolix2rx.Rd b/man/monolix2rx.Rd
index 43502c2..5117585 100644
--- a/man/monolix2rx.Rd
+++ b/man/monolix2rx.Rd
@@ -60,6 +60,8 @@ pkgTheo <- system.file("theo", package="monolix2rx")
 
 rx <- monolix2rx(file.path(pkgTheo, "theophylline_project.mlxtran"))
 
+rx <- monolix2rx(file.path(system.file("cov", package="monolix2rx"), "warfarin_covariate3_project.mlxtran"))
+
 rx
 }
 \author{