diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION
index cddd719..795e00a 100644
--- a/CRAN-SUBMISSION
+++ b/CRAN-SUBMISSION
@@ -1,3 +1,3 @@
-Version: 2.1.0
-Date: 2024-05-25 23:39:06 UTC
-SHA: 2d3961f6300b8e585c8d859e3cfba49dcbe12d96
+Version: 2.1.1
+Date: 2024-05-25 23:54:54 UTC
+SHA: 4119733f949e30b97c6a83fe737cdcc45235e9b9
diff --git a/DESCRIPTION b/DESCRIPTION
index eed1cfb..297a63e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: rxode2random
 Title: Random Number Generation Functions for 'rxode2'
-Version: 2.1.1
+Version: 2.1.1.9000
 Authors@R: c(
     person("Matthew L.", "Fidler", , "matthew.fidler@gmail.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0001-8538-6691")),
diff --git a/NEWS.md b/NEWS.md
index a429f11..05b5d9d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+# rxode2random (development version)
+
+- Fix a bug when simulating nested variables (#25)
+
 # rxode2random 2.1.0
 
 - **Breaking Change** changed distributions from the standard C++
diff --git a/src/cvPost.cpp b/src/cvPost.cpp
index b169378..e0532af 100644
--- a/src/cvPost.cpp
+++ b/src/cvPost.cpp
@@ -873,7 +873,7 @@ SEXP expandPars_(SEXP objectS, SEXP paramsS, SEXP eventsS, SEXP controlS) {
       ni = nestingInfo_(omegaLotri, eventsL);
       IntegerVector idIV =  as<IntegerVector>(ni["id"]);//length(levels(.ni$id))
       nid = Rf_length(Rf_getAttrib(idIV, R_LevelsSymbol));
-      if (nid <= 1){
+      if (nid <= 1) {
       } else if (nSub <= 1) {
         // control[Rxc_nSub] = nid;
         nSub = nid;
@@ -934,11 +934,20 @@ SEXP expandPars_(SEXP objectS, SEXP paramsS, SEXP eventsS, SEXP controlS) {
     methodStr = as<std::string>(control[Rxc_omegaSeparation]);
     methodInt = getMethodInt(methodStr, allNames, et);
 
+    if (Rf_isNull(lotriAbove)) {
+      aboveSEXP = R_NilValue;
+      rxModelsAssign(".nestEvents", events);
+      rxModelsAssign(".nestEta",    R_NilValue);
+      rxModelsAssign(".nestTheta",  R_NilValue);
+    }
+
     if (!Rf_isNull(aboveSEXP)) {
       // Create an extra theta matrix list
       //
       // Note this is for between study variability and the
-      // method/specification comes from omega, so methodInt comes from omegaSeparation
+      // method/specification comes from omega, so methodInt comes from
+      // omegaSeparation
+
       SEXP thetaList = PROTECT(cvPost_(et, lotriAbove,
                                        nStudS,
                                        LogicalVector::create(false),
@@ -955,16 +964,16 @@ SEXP expandPars_(SEXP objectS, SEXP paramsS, SEXP eventsS, SEXP controlS) {
       NumericVector upper = bounds[0];
       NumericVector lower = bounds[1];
       // With
-           NumericMatrix aboveMat = rxRmvnSEXP(IntegerVector::create(1),
-                                               R_NilValue, thetaList,
-                                               upper, lower, // lower upper
-                                               control[Rxc_nCoresRV],
-                                               LogicalVector::create(false), // isChol
-                                               LogicalVector::create(true), // keepNames
-                                               NumericVector::create(0.4), // a
-                                               NumericVector::create(2.05), // tol
-                                               NumericVector::create(1e-10), // nlTol
-                                               IntegerVector::create(100)); // nlMaxiter
+      NumericMatrix aboveMat = rxRmvnSEXP(IntegerVector::create(1),
+                                          R_NilValue, thetaList,
+                                          upper, lower, // lower upper
+                                          control[Rxc_nCoresRV],
+                                          LogicalVector::create(false), // isChol
+                                          LogicalVector::create(true), // keepNames
+                                          NumericVector::create(0.4), // a
+                                          NumericVector::create(2.05), // tol
+                                          NumericVector::create(1e-10), // nlTol
+                                          IntegerVector::create(100)); // nlMaxiter
       DataFrame newLst = as<DataFrame>(aboveMat);
       if (!Rf_isNull(et) && Rf_length(et) != 0){
         CharacterVector etListNames = asCv(Rf_getAttrib(et, R_NamesSymbol), "names(et)");
@@ -990,6 +999,9 @@ SEXP expandPars_(SEXP objectS, SEXP paramsS, SEXP eventsS, SEXP controlS) {
     } else {
       rxModelsAssign(".thetaL", R_NilValue);
     }
+    if (Rf_isNull(lotriBelow)) {
+      belowSEXP = R_NilValue;
+    }
     if (!Rf_isNull(belowSEXP)) {
       // below to sample matrix
       SEXP omegaList = PROTECT(cvPost_(et, // In case needed