diff --git a/DESCRIPTION b/DESCRIPTION
index 33a3c921..cdb1efc0 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: semlbci
 Title: Likelihood Based Confidence Interval for SEM
-Version: 0.0.0.6
+Version: 0.0.0.7
 Authors@R: 
     person(given = "Shu Fai",
            family = "Cheung",
diff --git a/NAMESPACE b/NAMESPACE
index d371e06b..b632e354 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,12 +1,15 @@
 # Generated by roxygen2: do not edit by hand
 
 export(ci_bound_i)
+export(ci_bound_nm_i)
 export(ci_i)
+export(find_dependent)
 export(find_free)
 export(find_variance_in_free)
 export(get_lhs_op_rhs)
 export(semlbci)
 export(set_constraint)
+export(set_constraint_nm)
 export(set_start)
 export(syntax_to_i)
 export(update_model)
diff --git a/NEWS.md b/NEWS.md
index ebbaf8b4..49228a7f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+# semlbci 0.0.0.7
+
+- Added diagnostic info to the output.
+
 # semlbci 0.0.0.6
 
 * Added functions to implement the Neale-Miller 1992 approach. Preliminary tests passed. Can be invoked by `method = "nm"`.
diff --git a/R/ci_bound_i.R b/R/ci_bound_i.R
index 25cd9dd2..27172c3f 100644
--- a/R/ci_bound_i.R
+++ b/R/ci_bound_i.R
@@ -25,6 +25,12 @@
 #'                       be used as the starting value for the target parameter if it is not a userd-defined paramter.
 #' @param standardized If TRUE, the LBCI is for the standardized estimate. 
 #' @param opts Options to be passed to \code{nloptr}
+#' @param ci_limit_ratio_tol The tolerance for the ratio of a to b, where a is the distance between an LBCI limit and the 
+#'               point estimate, and the b is the distance between the original [`lavaan`] CI limit (by default the Wald CI)
+#'               and the point estimate. If the ratio is larger than this value, a warning is set in the status code.
+#'               Default is 1.5.
+#' @param verbose If [`TRUE`], will store more diagnostic information in the attribute [`diag`].GlobalEnv
+#'                Default is [`FALSE`].
 #' @param ... Optional arguments. Not used.
 #' 
 #'@examples
@@ -51,6 +57,8 @@ ci_bound_i <- function(i = NULL,
                        wald_ci_start = TRUE,
                        standardized = FALSE,
                        opts = list(),
+                       ci_limit_ratio_tol = 1.5,
+                       verbose = FALSE,
                        ...) {
     k <- switch(which,
                 lbound = 1,
@@ -58,6 +66,34 @@ ci_bound_i <- function(i = NULL,
     # Check if the parameter is a user-defined parameter
     p_table <- lavaan::parameterTable(sem_out) 
     i_op <- p_table[i, "op"]
+
+    # Get original point estiamte and CI
+    if (standardized) {
+        p_est <- lavaan::standardizedSolution(sem_out,
+                    type = "std.all",
+                    se = TRUE,
+                    zstat = FALSE,
+                    pvalue = FALSE,
+                    ci = TRUE,
+                    remove.eq = FALSE,
+                    remove.ineq = FALSE,
+                    remove.def = FALSE,
+                    output = "data.frame")
+        i_est <- p_est[i, "est.std"]
+        i_org_ci_lower <- p_est[i, "ci.lower"]
+        i_org_ci_upper <- p_est[i, "ci.upper"]
+      } else {
+        p_est <- lavaan::parameterEstimates(sem_out,
+                                            se = TRUE,
+                                            zstat = FALSE,
+                                            fmi = FALSE,
+                                            rsquare = TRUE,
+                                            output = "data.frame")
+        i_est <- p_est[i, "est"]
+        i_org_ci_lower <- p_est[i, "ci.lower"]
+        i_org_ci_upper <- p_est[i, "ci.upper"]
+      }
+
     if (standardized && (i_op %in% c("=~", "~", "~~", ":="))) {
         p_std <- lavaan::standardizedSolution(sem_out,
                                               type = "std.all",
@@ -103,50 +139,6 @@ ci_bound_i <- function(i = NULL,
       }
     if (i_op == ":=") {
         if (standardized) {
-            # # To be deleted
-            # # Copied from the block for standardized free parameters.
-            # # If works, should merge the two blocks together
-            # p_std <- lavaan::standardizedSolution(sem_out,
-            #                                       type = "std.all",
-            #                                       se = FALSE,
-            #                                       zstat = FALSE,
-            #                                       pvalue = FALSE,
-            #                                       ci = FALSE,
-            #                                       remove.eq = FALSE,
-            #                                       remove.ineq = FALSE,
-            #                                       remove.def = FALSE,
-            #                                       output = "data.frame")
-            # p_std$id <- seq_len(nrow(p_std))
-            # i_lor <- get_lhs_op_rhs(i, sem_out)
-            # i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs"))$id
-            # start0 <- lavaan::parameterTable(sem_out)
-            # # The function to be minimized.
-            # lbci_b_f <- function(param, sem_out, debug, lav_warn) {
-            #     start1 <- start0
-            #     start1[start1$free > 0, "est"] <- param
-            #     sem_out2 <- sem_out
-            #     sem_out2@ParTable <- as.list(start1)
-            #     sem_model <- sem_out2@Model
-            #     sem_model <- update_model(sem_model, 
-            #                               start1[start1$free > 0, "est"] )
-            #     sem_out2@Model <- sem_model
-            #     std0 <- lavaan::standardizedSolution(sem_out2,
-            #                                     type = "std.all",
-            #                                     se = FALSE,
-            #                                     zstat = FALSE,
-            #                                     pvalue = FALSE,
-            #                                     ci = FALSE,
-            #                                     remove.eq = FALSE,
-            #                                     remove.ineq = FALSE,
-            #                                     remove.def = FALSE,
-            #                                     output = "data.frame")
-            #     k * std0[i_std, "est.std"]
-            #   }
-            # # The gradient of the function to be minimized
-            # lbci_b_grad <- function(param, sem_out, debug, lav_warn) {
-            #     numDeriv::grad(lbci_b_f, param, sem_out = sem_out, 
-            #                             debug = debug, lav_warn = lav_warn)
-            #   }
           } else {
             # Get the name of the defined parameter
             i_name <- p_table[i, "label"]
@@ -162,48 +154,6 @@ ci_bound_i <- function(i = NULL,
           }
       } else {
         if (standardized) {
-            # # To be deleted:
-            # p_std <- lavaan::standardizedSolution(sem_out,
-            #                                       type = "std.all",
-            #                                       se = FALSE,
-            #                                       zstat = FALSE,
-            #                                       pvalue = FALSE,
-            #                                       ci = FALSE,
-            #                                       remove.eq = FALSE,
-            #                                       remove.ineq = FALSE,
-            #                                       remove.def = FALSE,
-            #                                       output = "data.frame")
-            # p_std$id <- seq_len(nrow(p_std))
-            # i_lor <- get_lhs_op_rhs(i, sem_out)
-            # i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs"))$id
-            # start0 <- lavaan::parameterTable(sem_out)
-            # # The function to be minimized.
-            # lbci_b_f <- function(param, sem_out, debug, lav_warn) {
-            #     start1 <- start0
-            #     start1[start1$free > 0, "est"] <- param
-            #     sem_out2 <- sem_out
-            #     sem_out2@ParTable <- as.list(start1)
-            #     sem_model <- sem_out2@Model
-            #     sem_model <- update_model(sem_model, 
-            #                               start1[start1$free > 0, "est"] )
-            #     sem_out2@Model <- sem_model
-            #     std0 <- lavaan::standardizedSolution(sem_out2,
-            #                                     type = "std.all",
-            #                                     se = FALSE,
-            #                                     zstat = FALSE,
-            #                                     pvalue = FALSE,
-            #                                     ci = FALSE,
-            #                                     remove.eq = FALSE,
-            #                                     remove.ineq = FALSE,
-            #                                     remove.def = FALSE,
-            #                                     output = "data.frame")
-            #     k * std0[i_std, "est.std"]
-            #   }
-            # # The gradient of the function to be minimized
-            # lbci_b_grad <- function(param, sem_out, debug, lav_warn) {
-            #     numDeriv::grad(lbci_b_f, param, sem_out = sem_out, 
-            #                             debug = debug, lav_warn = lav_warn)
-            #   }
           } else {
             # The function to be minimized.
             lbci_b_f <- function(param, sem_out, debug, lav_warn) {
@@ -252,6 +202,31 @@ ci_bound_i <- function(i = NULL,
                         lav_warn = FALSE,
                         debug = FALSE)
     bound <- k * out$objective
+    bound_unchecked <- bound
+
+    # Initialize the status code
+    status <- 0
+
+    # Check convergence
+    if (out$status < 0) {
+        # Verify the range that denotes error
+        status <- 1
+        bound <- NA
+      }
+
+    # Check the limit
+    i_org_ci_limit <- switch(which,
+                        lbound = i_org_ci_lower,
+                        ubound = i_org_ci_upper)
+    if (!is.na(bound)) {
+        ci_limit_ratio <- abs((bound - i_est) / (i_org_ci_limit - i_est))
+        if (ci_limit_ratio > ci_limit_ratio_tol) {
+            status <- 1
+            # Do not set the bound to NA because the limit may still be valid.
+          }
+      } else {
+        ci_limit_ratio <- NA
+      }
 
     # Check whether admissible
     start0 <- lavaan::parameterTable(sem_out)
@@ -263,11 +238,30 @@ ci_bound_i <- function(i = NULL,
                                 check.vcov = TRUE)
     fit_post_check <- lavaan::lavInspect(fit_final, "post.check")
     if (!fit_post_check) {
+        status <- 0
         bound <- NA
-        warning("Optimization converged but the final solution is not admissible.")
+        # The warning should be raised by the calling function, not this one
+        # warning("Optimization converged but the final solution is not admissible.")
+      }
+    diag <- list(status = status,
+                 est_org = i_est,
+                 ci_org_limit = i_org_ci_limit,
+                 ci_limit_ratio = ci_limit_ratio,
+                 fit_post_check = fit_post_check,
+                 start_values = xstart,
+                 bound_unchecked = bound_unchecked
+                 )
+    if (verbose) {
+        diag$history <- out
+        diag$fit_final <- fit_final
+      } else {
+        diag$history <- NULL
+        diag$fit_final <- NULL
       }
-    if (history) {
-        attr(bound, "history") <- out
+    if (out$status < 0) {
+        # If convergence status < 0, override verbose
+        diag$history <- out        
       }
+    attr(bound, "diag") <- diag
     bound
   }
diff --git a/R/ci_bound_nm_i.R b/R/ci_bound_nm_i.R
index eedb88b6..5362292f 100644
--- a/R/ci_bound_nm_i.R
+++ b/R/ci_bound_nm_i.R
@@ -104,21 +104,6 @@ ci_bound_nm_i <- function(i = NULL,
         i_org_ci_lower <- p_est[i, "ci.lower"]
         i_org_ci_upper <- p_est[i, "ci.upper"]
       }
-    # TO DELETE
-    # if (standardized) {
-    #     i_est <- lavaan::standardizedSolution(sem_out,
-    #                 type = "std.all",
-    #                 se = TRUE,
-    #                 zstat = FALSE,
-    #                 pvalue = FALSE,
-    #                 ci = TRUE,
-    #                 remove.eq = FALSE,
-    #                 remove.ineq = FALSE,
-    #                 remove.def = FALSE,
-    #                 output = "data.frame")[i, "est.std"]
-    #   } else {
-    #     i_est <- p_table[i, "est"]
-    #   }
     if (wald_ci_start) {
         if (i_op == ":=") {
               xstart <- set_start(i, sem_out, which)
@@ -212,7 +197,6 @@ ci_bound_nm_i <- function(i = NULL,
 
             # Enforce linear equality constraints
             if (!is.null(eq_K)) {
-                # stop("Not yet work for models with equality constraints")
                 param_i_reduced <- rep(0, npar_reduced)
                 param_i_reduced[i_depend_reduced] <- param_depend
                 param_i <- as.numeric(eq_K %*% matrix(param_i_reduced, npar_reduced, 1) + eq_k0)
@@ -233,46 +217,13 @@ ci_bound_nm_i <- function(i = NULL,
                                        baseline = FALSE,
                                        h1 = FALSE)
                 assign("f_i_shared", fit2, envir0)
-                # p_table2 <- lavaan::parameterTable(fit2)
-                # p_table2$esty <- p_table2$est
               } else {
-                # p_table2 <- lavaan::parameterTable(envir0$f_i_shared)
-                # p_table2[i_depend, "est"] <- param_depend
-                # p_table2[i_exp_target, "est"] <- p_table2[i_exp_source, "est"] 
-                # p_table2$esty <- p_table2$est
               }
             if (standardized) {
-                # p_table2 <- lavaan::parameterTable(envir0$f_i_shared)
-                # p_table2 <- as.data.frame(envir0$f_i_shared@ParTable)
-                # p_table2[i_depend, "est"] <- param_depend
-                # p_table2[i_exp_target, "est"] <- p_table2[i_exp_source, "est"] 
-                # p_table2$esty <- p_table2$est
-
                 p_table2 <- envir0$f_i_shared@ParTable
                 p_table2$est[i_depend] <- param_depend
                 p_table2$est[i_exp_target] <- p_table2$est[i_exp_source] 
-
-                # sem_out2 <- sem_out
-                # p_table2 <- as.data.frame(envir0$f_i_shared@ParTable)
-                # k1 <- paste0(p_table$lhs, p_table$op, p_table$rhs)
-                # k2 <- paste0(p_table2$lhs, p_table2$op, p_table2$rhs)
-                # start3 <- p_table
-                # start3[match(k2, k1), "est"] <- p_table2[k2 %in% k1, "est"]
-                # sem_out2@ParTable <- as.list(start3)
-                # sem_model <- sem_out2@Model
-                # sem_model <- update_model(sem_model, 
-                #                           start3[start3$free > 0, "est"] )
-                # sem_out2@Model <- sem_model
-                # sem_out2@Model@GLIST <- envir0$f_i_shared@Model@GLIST
-
                 sem_out2 <- sem_out
-
-                # start3 <- merge(p_table, p_table2[, c("lhs", "op", "rhs", "esty")], 
-                #             by = c("lhs", "op", "rhs"), all.x = TRUE, all.y = FALSE,
-                #             sort = FALSE)
-                # start3$est <- start3$esty
-                # start3$esty <- NULL
-
                 start3 <- p_table
                 k1 <- paste0(p_table$lhs, p_table$op, p_table$rhs)
                 k2 <- paste0(p_table2$lhs, p_table2$op, p_table2$rhs)
@@ -293,24 +244,12 @@ ci_bound_nm_i <- function(i = NULL,
                                                 remove.ineq = FALSE,
                                                 remove.def = FALSE,
                                                 output = "data.frame")
-                # return(k * std0[i, "est.std"])
                 return(k * ( std0[i, "est.std"] - i_est))
               } else {
                 if (i_op == ":=") {
-                  # return(k * sem_out@Model@def.function(param_i)[i_name])
                   return(k * (sem_out@Model@def.function(param_i)[i_name] - i_est))
                 } else {
-                  # return(k * p_table2[i, "est"])
-                  # tmp <- mapply(function(x, y, z) {
-                  #     z[y[x == i]]
-                  #     },
-                  #       sem_out@Model@x.free.idx,
-                  #       sem_out@Model@m.free.idx, 
-                  #       envir0$f_i_shared@Model@GLIST)
-                  # tmp <- unlist(tmp)
-                  # return(k * (tmp - i_est))
                   return(k * (envir0$f_i_shared@ParTable$est[i] - i_est))
-                  # return(k * (p_table2[i, "est"] - i_est))
                 }
               }
           }
@@ -338,7 +277,6 @@ ci_bound_nm_i <- function(i = NULL,
 
             g3_tmp <- lavaan::lavTech(fit3, "gradient")
             if (!is.null(eq_K)) {
-                # stop("Not yet work for models with equality constraints")
                 g3_tmp <- (g3_tmp - eq_k0) %*% eq_K
               }
             list(constraints = rbind(f_obj),
@@ -360,149 +298,7 @@ ci_bound_nm_i <- function(i = NULL,
         fit_lb[find_variance_in_free(sem_out)[i_depend]] <- lb_var
         xstart <-  xstart[i_depend_reduced]
       }
-    if (standardized && (i_op %in% c("=~", "~", "~~", ":=")) && !test_generic) {
-        # TODO: Not ready
-        # Get the name of the defined parameter
-        i_name <- p_table[i, "label"]
-        # The function to be minimized.
-        i_depend <- find_dependent(i, sem_out = sem_out, standardized = TRUE)
-        # The function to be minimized.
-        param_i <- rep(NA, npar)
-        start0 <- lavaan::parameterTable(sem_out)
-        envir0 <- new.env()
-        assign("f_i_shared", sem_out, envir = envir0)
-        assign("f_i_free_shared", sem_out, envir = envir0)
-        # f_i_shared <- sem_out
-        # f_i_free_shared <- sem_out
-        lbci_b_f <- function(param_depend, sem_out, debug, lav_warn) {
-            force(i_depend)
-            force(param_i)
-            force(i_name)
-            force(start0)
-            force(envir0)
-            if (is.null(attr(param_depend, "refit"))) {
-                param_i[i_depend] <- param_depend
-                start1 <- start0[!(start0$op == ":="), ]
-                start1[i_depend, "free"] <- 0
-                start1[i_depend, "ustart"] <- 0
-                start1[i_depend, "est"] <- param_depend
-                fit2 <- lavaan::update(sem_out, start1)
-                # f_i_shared <<- fit2
-                assign("f_i_shared", fit2, envir0)
-                p_table2 <- lavaan::parameterTable(fit2)
-                p_table2$esty <- p_table2$est
-              } else {
-                p_table2 <- lavaan::parameterTable(envir0$f_i_shared)
-                p_table2[i_depend, "est"] <- param_depend
-                p_table2$esty <- p_table2$est
-              }
-            # start1 <- start0[start0$op != "==", ]
-            # start1[i_depend, "est"] <- param_i
-            sem_out2 <- sem_out
-            start3 <- merge(p_table, p_table2[, c("lhs", "op", "rhs", "esty")], 
-                        by = c("lhs", "op", "rhs"), all.x = TRUE, all.y = FALSE,
-                        sort = FALSE)
-            start3$est <- start3$esty
-            start3$esty <- NULL
-            sem_out2@ParTable <- as.list(start3)
-            sem_model <- sem_out2@Model
-            sem_model <- update_model(sem_model, 
-                                      start3[start3$free > 0, "est"] )
-            sem_out2@Model <- sem_model
-            std0 <- lavaan::standardizedSolution(sem_out2,
-                                            type = "std.all",
-                                            se = FALSE,
-                                            zstat = FALSE,
-                                            pvalue = FALSE,
-                                            ci = FALSE,
-                                            remove.eq = FALSE,
-                                            remove.ineq = FALSE,
-                                            remove.def = FALSE,
-                                            output = "data.frame")
-            k * std0[i, "est.std"]
-          }
-        lbci_b_grad <- function(param_depend, sem_out, debug, lav_warn) {
-            param_depend_g <- param_depend
-            attr(param_depend_g, "refit") <- FALSE 
-            lavaan::lav_func_gradient_simple(lbci_b_f, param_depend_g, sem_out = sem_out, 
-                                    debug = debug, lav_warn = lav_warn)
-          }
-        f_constr = set_constraint_nm(i_depend, sem_out, envir = envir0)
-        # Set shared variables
-        # lbci_b_grad <- NULL
-        fit_lb <- rep(-Inf, length(i_depend))
-        fit_ub <- rep( Inf, length(i_depend))
-        fit_lb[find_variance_in_free(sem_out)[i_depend]] <- lb_i_var
-        xstart <-  xstart[i_depend]
-      }
-    if ((i_op == ":=") && !test_generic) {
-        if (standardized) {
-          } else {
-            # Get the name of the defined parameter
-            i_name <- p_table[i, "label"]
-            # The function to be minimized.
-            i_depend <- find_dependent(i, sem_out = sem_out, standardized = FALSE)
-            # Set shared variables
-            envir0 <- new.env()
-            assign("f_i_shared", sem_out, envir = envir0)
-            assign("f_i_free_shared", sem_out, envir = envir0)
-            # f_i_shared <- sem_out
-            # f_i_free_shared <- sem_out
-            # The function to be minimized.
-            param_i <- rep(NA, npar)
-            lbci_b_f <- function(param_depend, sem_out, debug, lav_warn) {
-                force(i_depend)
-                force(param_i)
-                force(i_name)
-                param_i[i_depend] <- param_depend
-                k * sem_out@Model@def.function(param_i)[i_name]
-              }
-            # lbci_b_grad <- function(param_depend, sem_out, debug, lav_warn) {
-            #     numDeriv::grad(lbci_b_f, param_depend, sem_out = sem_out, 
-            #                             debug = debug, lav_warn = lav_warn)
-            #   }
-            lbci_b_grad <- function(param_depend, sem_out, debug, lav_warn) {
-                force(i)
-                force(envir0)
-                force(i_depend)
-                lavaan::lavTech(envir0$f_i_free_shared, "gradient")[i_depend]
-              }
-            # f_constr = set_constraint_nm(i_depend, sem_out)
-            f_constr = set_constraint_nm(i_depend, sem_out, get_fit_from_envir = FALSE,
-                                         envir = envir0)
-            # lbci_b_grad <- NULL
-            fit_lb <- rep(-Inf, length(i_depend))
-            fit_ub <- rep( Inf, length(i_depend))
-            fit_lb[find_variance_in_free(sem_out)[i_depend]] <- lb_var
-            xstart <-  xstart[i_depend]
-          }
-      } else {
-        if (!test_generic) {
-          if (standardized) {
-            } else {
-              # Set shared variables
-              envir0 <- new.env()
-              assign("f_i_shared", sem_out, envir = envir0)
-              assign("f_i_free_shared", sem_out, envir = envir0)
-              # The function to be minimized.
-              lbci_b_f <- function(p_f, sem_out, debug, lav_warn) {
-                  k * p_f
-                }
-              # The gradient of the function to be minimized
-              lbci_b_grad <- function(p_f, sem_out, debug, lav_warn) {
-                  force(i)
-                  force(envir0)
-                  lavaan::lavTech(envir0$f_i_free_shared, "gradient")[i]
-                }
-              f_constr = set_constraint_nm(i, sem_out, get_fit_from_envir = FALSE,
-                                          envir = envir0)
-              # lbci_b_grad <- NULL
-              fit_lb <- -Inf
-              fit_ub <-  Inf
-              xstart <-  xstart[i]
-            }
-        }
-      }
+
     # fit_lb <- rep(-Inf, npar)
     # fit_lb[find_variance_in_free(sem_out)] <- lb_var
     opts_final <- modifyList(list("algorithm" = "NLOPT_LD_SLSQP",
@@ -521,7 +317,6 @@ ci_bound_nm_i <- function(i = NULL,
                         sem_out = sem_out,
                         lav_warn = FALSE,
                         debug = FALSE)
-    # bound <- k * out$objective
     bound <- i_est + k * out$objective
     bound_unchecked <- bound
 
@@ -539,10 +334,14 @@ ci_bound_nm_i <- function(i = NULL,
     i_org_ci_limit <- switch(which,
                         lbound = i_org_ci_lower,
                         ubound = i_org_ci_upper)
-    ci_limit_ratio <- abs((bound - i_est) / (i_org_ci_limit - i_est))
-    if (ci_limit_ratio > ci_limit_ratio_tol) {
-        status <- 1
-        # Do not set the bound to NA because the limit may still be valid.
+    if (!is.na(bound)) {
+        ci_limit_ratio <- abs((bound - i_est) / (i_org_ci_limit - i_est))
+        if (ci_limit_ratio > ci_limit_ratio_tol) {
+            status <- 1
+            # Do not set the bound to NA because the limit may still be valid.
+          }
+      } else {
+        ci_limit_ratio <- NA
       }
 
     # Check whether the solution at the limit is admissible
@@ -569,13 +368,7 @@ ci_bound_nm_i <- function(i = NULL,
         # The warning should be raised by the calling function, not this one
         # warning("Optimization converged but the final solution is not admissible.")
       }
-    # TO DELETE
-    # if (history) {
-    #     attr(bound, "history") <- out
-    #     attr(bound, "fit") <- envir0$f_i_shared
-    #     attr(bound, "xstart") <- xstart
-    #     attr(bound, "i_est") <- i_est
-    #   }
+
     diag <- list(status = status,
                  est_org = i_est,
                  ci_org_limit = i_org_ci_limit,
diff --git a/R/ci_i.R b/R/ci_i.R
index b7dc0653..ad8b8915 100644
--- a/R/ci_i.R
+++ b/R/ci_i.R
@@ -30,19 +30,18 @@
 
 ci_i <- function(i, method = "wn", ...) {
     if (method == "wn") {
-        lb <- ci_bound_i(i, which = "lbound", ...)
-        ub <- ci_bound_i(i, which = "ubound", ...)
+        lb_time <- system.time(lb <- ci_bound_i(i, which = "lbound", ...))
+        ub_time <- system.time(ub <- ci_bound_i(i, which = "ubound", ...))
       }
     if (method == "nm") {
-        lb <- ci_bound_nm_i(i, which = "lbound", ...)
-        ub <- ci_bound_nm_i(i, which = "ubound", ...)
+        lb_time <- system.time(lb <- ci_bound_nm_i(i, which = "lbound", ...))
+        ub_time <- system.time(ub <- ci_bound_nm_i(i, which = "ubound", ...))
       }
     out <- c(lb, ub)
-    if (!is.null(attr(lb, "history"))) {
-        attr(out, "lb_history") <- attr(lb, "history")
-      }
-    if (!is.null(attr(ub, "history"))) {
-        attr(out, "ub_history") <- attr(ub, "history")
-      }
+    attr(out, "lb_diag") <- attr(lb, "diag")
+    attr(out, "ub_diag") <- attr(ub, "diag")
+    attr(out, "method") <- method
+    attr(out, "lb_time") <- lb_time[3]
+    attr(out, "ub_time") <- ub_time[3]
     out
   }
diff --git a/R/data.R b/R/data.R
index 6dc36ad6..82e0d6ed 100644
--- a/R/data.R
+++ b/R/data.R
@@ -1,8 +1,15 @@
-#'Sample dataset for testing
-#' 
 #' Generated from a two-factor model with six variables, n 500
 #' 
 "cfa_two_factors"
 
 #' Generated from a simple mediation model, n 500
-"simple_med"
\ No newline at end of file
+#' 
+"simple_med"
+
+#' Generated from a regression model, with one correlation close to one.
+#' 
+"reg_cor_near_one"
+
+#' Generated from a two-factor model, with one error variacne close to zero.
+#' 
+"cfa_evar_near_zero"
\ No newline at end of file
diff --git a/R/semlbci.R b/R/semlbci.R
index 78379a64..3649c7eb 100644
--- a/R/semlbci.R
+++ b/R/semlbci.R
@@ -19,7 +19,7 @@
 #'               is .95.
 #' @param standardized If TRUE, the LBCI is for the standardized estimate.
 #' @param method The approach to be used. Can be "wn" (Wu-Neale-2012) or "nm" 
-#'               (Neale-Miller-1997). Default is "wn".
+#'               (Neale-Miller-1997). Default is "nm".
 #' @param ... Arguments to be passed to \code{ci_bound_i}.
 #' @param parallel If \code{TRUE}, will use parallel. Currently disabled.
 #'                  Need to find out how to make \code{lavaan::udpate} works
@@ -43,14 +43,14 @@ semlbci <- function(sem_out,
                     pars = NULL,
                     ciperc = .95,
                     standardized = FALSE,
-                    method = "wn",
+                    method = "nm",
                     ...,
                     parallel = FALSE,
                     ncpu = 2) {
     if (!inherits(sem_out, "lavaan")) {
         stop("sem_out is not a supported object.")
       }
-    ptable <- lavaan::parameterTable(sem_out)
+    ptable <- as.data.frame(lavaan::parameterTable(sem_out))
     # Do not check for 
     i <- ptable$free > 0
     #i_id <- ptable$id[i]
@@ -104,7 +104,7 @@ semlbci <- function(sem_out,
         #     stop("Error occured in parallel mode. Try setting parallel = FALSE")
         #   }
       } else {
-        out <- lapply(pars, ci_i, 
+        out_raw <- lapply(pars, ci_i, 
                       npar = npar,
                       sem_out = sem_out,
                       standardized = standardized,
@@ -114,7 +114,7 @@ semlbci <- function(sem_out,
                       ciperc = ciperc,
                       ...)
       }
-    out <- do.call(rbind, out)
+    out <- do.call(rbind, out_raw)
     out_p <- ptable[, c("id", "lhs", "op", "rhs")]
     out_p$lbci_lb <- NA
     if (standardized) {
@@ -125,8 +125,47 @@ semlbci <- function(sem_out,
         out_p$est <- ptable[, c("est")]        
       }
     out_p$lbci_ub <- NA
+    
     out_p[i_selected, "lbci_lb"] <- out[, 1]
     out_p[i_selected, "lbci_ub"] <- out[, 2]
+
+    # Collect diagnostic info
+
+    lb_diag <- lapply(out_raw, attr, which = "lb_diag")
+    ub_diag <- lapply(out_raw, attr, which = "ub_diag")
+    lb_time <- sapply(out_raw, attr, which = "lb_time")
+    ub_time <- sapply(out_raw, attr, which = "ub_time")
+    ci_method <- sapply(out_raw, attr, which = "method")
+    p_names <- mapply(paste0, out_p[pars, "lhs"],
+                              out_p[pars, "op" ],
+                              out_p[pars, "rhs"],
+                      USE.NAMES = FALSE)
+    names(lb_diag) <- p_names
+    names(ub_diag) <- p_names
+    names(lb_time) <- p_names
+    names(ub_time) <- p_names
+    names(ci_method) <- p_names
+    
+    attr(out_p, "lb_diag") <- lb_diag
+    attr(out_p, "ub_diag") <- ub_diag
+    attr(out_p, "lb_time") <- lb_time
+    attr(out_p, "ub_time") <- ub_time
+    attr(out_p, "ci_method") <- ci_method
+
+    # Append diagnostic info
+
+    out_p[pars, "status_lb"] <- sapply(lb_diag, function(x) x$status)
+    out_p[pars, "status_ub"] <- sapply(ub_diag, function(x) x$status)
+    out_p[pars, "ci_org_lb"] <- sapply(lb_diag, function(x) x$ci_org_limit)
+    out_p[pars, "ci_org_ub"] <- sapply(ub_diag, function(x) x$ci_org_limit)
+    out_p[pars, "ratio_lb"] <- sapply(lb_diag, function(x) x$ci_limit_ratio)
+    out_p[pars, "ratio_ub"] <- sapply(ub_diag, function(x) x$ci_limit_ratio)
+    out_p[pars, "post_check_lb"] <- sapply(lb_diag, function(x) x$fit_post_check)
+    out_p[pars, "post_check_ub"] <- sapply(ub_diag, function(x) x$fit_post_check)
+    out_p[pars, "time_lb"] <- as.numeric(lb_time)
+    out_p[pars, "time_ub"] <- as.numeric(ub_time)
+    out_p[pars, "method"] <- ci_method
+
     class(out_p) <- c("semlbci", class(out_p))
     out_p
   }
\ No newline at end of file
diff --git a/R/set_constraint_nm.R b/R/set_constraint_nm.R
index a1d9d074..f6f46d71 100644
--- a/R/set_constraint_nm.R
+++ b/R/set_constraint_nm.R
@@ -46,36 +46,6 @@ set_constraint_nm <- function(i, sem_out, ciperc = .95, envir = NULL, get_fit_fr
     # Check if there are any equality constraints
     # if (sem_out@Model@eq.constraints) {
     if (FALSE) {
-        # Disabled. Not necessary for the Neale-Miller-1997 approach
-        # This part is to be deleted. It is copied from the code for the Wu-Nealed 2012 approach.
-        fn_constraint <- function(param, sem_out = NULL, debug = FALSE, lav_warn = FALSE) {
-            if (debug) {
-                cat(ls())
-                cat(ls(globalenv()))
-                }
-            start0 <- lavaan::parameterTable(sem_out)
-            start0[p_free, "est"] <- param
-            eq_out <- sem_out@Model@ceq.function(param)
-            eq_jac <- sem_out@Model@con.jac
-            if (lav_warn) {
-                    fit2 <- lavaan::update(sem_out, start = start0, do.fit = FALSE)
-                } else {
-                    suppressWarnings(fit2 <- lavaan::update(sem_out, start = start0, do.fit = FALSE))                    
-                }
-            if (lav_warn) {
-                    fit2_gradient <- rbind(lavaan::lavTech(fit2, "gradient"))
-                    fit2_jacobian <- rbind(eq_jac, lavaan::lavTech(fit2, "gradient"))
-                } else {
-                    suppressWarnings(fit2_gradient <- rbind(lavaan::lavTech(fit2, "gradient")))
-                    suppressWarnings(fit2_jacobian <- rbind(eq_jac, lavaan::lavTech(fit2, "gradient")))
-                }
-            list(
-                  objective = lavaan::lavTech(fit2, "optim")$fx,
-                  gradient = fit2_gradient,
-                  constraints = rbind(t(t(eq_out)), lavaan::lavTech(fit2, "optim")$fx - target),
-                  jacobian = fit2_jacobian,
-                  parameterTable = lavaan::parameterTable(fit2))
-          }
       } else {
         # if (length(i) == 1) {
         fn_constraint <- function(p_f, sem_out = NULL, debug = FALSE, lav_warn = FALSE) {
diff --git a/R/syntax_to_i.R b/R/syntax_to_i.R
index 6c908373..bea447d4 100644
--- a/R/syntax_to_i.R
+++ b/R/syntax_to_i.R
@@ -36,6 +36,7 @@ syntax_to_i <- function(syntax,
         l_model$req <- TRUE
         p_out <- merge(ptable, l_model[, c("lhs", "op", "rhs", "req")], 
                       by = c("lhs", "op", "rhs"), all.x = TRUE, sort = FALSE)
+        p_out <- p_out[match(ptable$id, p_out$id), ]
         i_par <- which(p_out$req)
       } else {
         i_par <- NULL
@@ -51,6 +52,5 @@ syntax_to_i <- function(syntax,
       } else {
         i_def <- NULL
       }
-    # which(p_out$req)
     c(i_par, i_def)
   }
\ No newline at end of file
diff --git a/R/update_model.R b/R/update_model.R
index a5218aca..2955621a 100644
--- a/R/update_model.R
+++ b/R/update_model.R
@@ -1,6 +1,6 @@
 #'@title Update a lavModel objec using user-supplied estimates
 #'
-#'@description Update a \code{lavaan model using user-supplied estimates
+#'@description Update a \code{lavaan} model using user-supplied estimates
 #'
 #'@details 
 #' 
diff --git a/data-raw/create_test_data.R b/data-raw/create_test_data.R
index b677a8f1..6a416763 100644
--- a/data-raw/create_test_data.R
+++ b/data-raw/create_test_data.R
@@ -52,4 +52,140 @@ head(dat)
 
 simple_med <- dat
 
-usethis::use_data(simple_med, overwrite = TRUE)
\ No newline at end of file
+usethis::use_data(simple_med, overwrite = TRUE)
+
+# # Data with a parameter estimate near its attainable boundary
+
+# library(MASS)
+
+# set.seed(3842093)
+# n <- 50
+# gamma <- matrix(c(.9995,
+#                   .5), 2, 1, byrow = TRUE)
+# beta  <- matrix(c(0, 0,
+#                   .5, 0), 2, 2, byrow = TRUE)
+# theta <- matrix(1, 1, 1)
+# # psi   <- diag(c(6, 6)^2)
+# xi    <- mvrnorm(n, 0, theta)
+# sigmay <-  solve(diag(2) - beta) %*% gamma %*% theta %*% t(gamma) %*% t(solve(diag(2) - beta))
+# sigmay
+# psi <- diag(1 - diag(sigmay))
+# eta   <- solve(diag(2) - beta) %*% (gamma %*% t(xi) + t(mvrnorm(n, c(0, 0), psi)))
+# eta   <- t(eta)
+# dat   <- cbind(xi, eta)
+# dat   <- as.data.frame(dat)
+# colnames(dat) <- c("x", "m", "y")
+# head(dat)
+# coef(lm.beta::lm.beta(lm(m ~ x, dat)))
+# coef(lm.beta::lm.beta(lm(y ~ m + x, dat)))
+# library(lavaan)
+# mod <- 
+# "
+# m ~ x
+# y ~ m + x
+# "
+# fit <- sem(mod, dat, fixed.x = FALSE)
+# # coef(fit)
+# dplyr::filter(as.data.frame(parameterEstimates(fit)), 
+#               lhs == "y", rhs == "y")
+# dplyr::filter(as.data.frame(standardizedSolution(fit)), 
+#               lhs == "m", rhs == "x")
+
+# Data based on a regression model with one correlation cloe to one.
+
+set.seed(8134704)
+n <- 100
+k <- 4
+xs <- mvrnorm(n, rep(0, k), diag(k))
+rxx <- .999
+xs2 <- rxx*xs[, k] + rnorm(n, 0, sqrt(1 - rxx))
+xs <- cbind(xs, xs2)
+colnames(xs) <- paste0("x", seq_len(k + 1))
+cor(xs)
+v_y <- .10
+bx <- sqrt((1 - v_y) / ((k + 1) + 2 * rxx))
+y <- xs %*% t(t(rep(bx, k + 1))) + rnorm(n, 0, sqrt(v_y))
+lm_out <- lm(y ~ xs)
+coef(lm.beta::lm.beta(lm_out))
+mod <- paste0("y ~ ", paste(colnames(xs), collapse = " + "))
+dat <- data.frame(xs, y)
+head(dat)
+library(lavaan)
+fit <- sem(mod, dat, fixed.x = FALSE)
+# coef(fit)
+dplyr::filter(as.data.frame(parameterEstimates(fit)), 
+              lhs == "y", rhs == "y")
+dplyr::filter(as.data.frame(standardizedSolution(fit)), 
+              lhs == paste0("x", k),
+              rhs == paste0("x", k + 1))
+mod2 <- 
+paste(mod, "\n",
+"
+x4 ~~ r45*x5
+r45 < 1
+")
+fit2 <- sem(mod2, dat, fixed.x = FALSE)
+dplyr::filter(as.data.frame(parameterEstimates(fit2)), 
+              lhs == "y", rhs == "y")
+dplyr::filter(as.data.frame(standardizedSolution(fit2)), 
+              lhs == paste0("x", k),
+              rhs == paste0("x", k + 1))
+
+reg_cor_near_one <- dat
+
+usethis::use_data(reg_cor_near_one, overwrite = TRUE)
+
+
+# Data based on a two factor model, with one error variance close to zero
+
+library(MASS)
+
+set.seed(5325235)
+n <- 120
+lambda <- matrix(c(  1,   0,
+                   1.2,   0,
+                   2.0,   0,
+                     0,   1,
+                     0, 0.8,
+                     0, 1.2), 6, 2, byrow = TRUE)
+f_var <- c(2, 3)
+f_cor <- .60
+phi <- matrix(c(1, f_cor,
+                f_cor, 1), 2, 2, byrow = TRUE)
+phi <- diag(sqrt(f_var)) %*% phi %*% diag(sqrt(f_var))
+psi <- diag(c(3, 2, .10, 4, 3, 2))
+dat_f <- mvrnorm(n, c(0, 0), phi)
+dat_e <- mvrnorm(n, rep(0, 6), psi)
+dat   <- lambda %*% t(dat_f) + t(dat_e)
+dat   <- t(dat)
+dat   <- as.data.frame(dat)
+colnames(dat) <- paste0("x", 1:6)
+head(dat)
+mod <- 
+"
+f1 =~ x1 + x2 + x3
+f2 =~ x4 + x5 + x6
+"
+fit <- cfa(mod, dat)
+dplyr::filter(as.data.frame(parameterEstimates(fit)), 
+              lhs == "x3", rhs == "x3")
+dplyr::filter(parameterEstimates(fit), 
+              lhs == rhs)
+
+mod2 <- 
+"
+f1 =~ x1 + x2 + x3
+f2 =~ x4 + x5 + x6
+x3 ~~ e3*x3
+e3 > 0
+"
+fit2 <- cfa(mod2, dat)
+dplyr::filter(as.data.frame(parameterEstimates(fit2)), 
+              lhs == "x3", rhs == "x3")
+dplyr::filter(parameterEstimates(fit2), 
+              lhs == rhs)
+
+
+cfa_evar_near_zero <- dat
+
+usethis::use_data(cfa_evar_near_zero, overwrite = TRUE)
diff --git a/data/cfa_evar_near_zero.rda b/data/cfa_evar_near_zero.rda
new file mode 100644
index 00000000..ea611b27
Binary files /dev/null and b/data/cfa_evar_near_zero.rda differ
diff --git a/data/reg_cor_near_one.rda b/data/reg_cor_near_one.rda
new file mode 100644
index 00000000..d5d203e8
Binary files /dev/null and b/data/reg_cor_near_one.rda differ
diff --git a/man/cfa_evar_near_zero.Rd b/man/cfa_evar_near_zero.Rd
new file mode 100644
index 00000000..22880010
--- /dev/null
+++ b/man/cfa_evar_near_zero.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{cfa_evar_near_zero}
+\alias{cfa_evar_near_zero}
+\title{Generated from a two-factor model, with one error variacne close to zero.}
+\format{
+An object of class \code{data.frame} with 120 rows and 6 columns.
+}
+\usage{
+cfa_evar_near_zero
+}
+\description{
+Generated from a two-factor model, with one error variacne close to zero.
+}
+\keyword{datasets}
diff --git a/man/cfa_two_factors.Rd b/man/cfa_two_factors.Rd
index a72a08ed..eb37af8d 100644
--- a/man/cfa_two_factors.Rd
+++ b/man/cfa_two_factors.Rd
@@ -3,7 +3,7 @@
 \docType{data}
 \name{cfa_two_factors}
 \alias{cfa_two_factors}
-\title{Sample dataset for testing}
+\title{Generated from a two-factor model with six variables, n 500}
 \format{
 An object of class \code{data.frame} with 500 rows and 6 columns.
 }
diff --git a/man/ci_bound_i.Rd b/man/ci_bound_i.Rd
index 12992f18..385ce97d 100644
--- a/man/ci_bound_i.Rd
+++ b/man/ci_bound_i.Rd
@@ -15,8 +15,9 @@ ci_bound_i(
   lb_var = -Inf,
   wald_ci_start = TRUE,
   standardized = FALSE,
-  opts = list(algorithm = "NLOPT_LD_SLSQP", xtol_rel = 1e-10, maxeval = 1000,
-    print_level = 0),
+  opts = list(),
+  ci_limit_ratio_tol = 1.5,
+  verbose = FALSE,
   ...
 )
 }
@@ -47,6 +48,14 @@ be used as the starting value for the target parameter if it is not a userd-defi
 
 \item{opts}{Options to be passed to \code{nloptr}}
 
+\item{ci_limit_ratio_tol}{The tolerance for the ratio of a to b, where a is the distance between an LBCI limit and the
+point estimate, and the b is the distance between the original \code{\link{lavaan}} CI limit (by default the Wald CI)
+and the point estimate. If the ratio is larger than this value, a warning is set in the status code.
+Default is 1.5.}
+
+\item{verbose}{If \code{\link{TRUE}}, will store more diagnostic information in the attribute \code{\link{diag}}.GlobalEnv
+Default is \code{\link{FALSE}}.}
+
 \item{...}{Optional arguments. Not used.}
 }
 \value{
diff --git a/man/ci_bound_nm_i.Rd b/man/ci_bound_nm_i.Rd
new file mode 100644
index 00000000..c4bbf265
--- /dev/null
+++ b/man/ci_bound_nm_i.Rd
@@ -0,0 +1,92 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ci_bound_nm_i.R
+\name{ci_bound_nm_i}
+\alias{ci_bound_nm_i}
+\title{Find the lower or upper bound for one parameter by the Neale-Miller-1997 approach}
+\usage{
+ci_bound_nm_i(
+  i = NULL,
+  npar = NULL,
+  sem_out = NULL,
+  f_constr = NULL,
+  which = NULL,
+  history = FALSE,
+  perturbation_factor = 0.9,
+  lb_var = -Inf,
+  wald_ci_start = TRUE,
+  standardized = FALSE,
+  opts = list(),
+  test_generic = TRUE,
+  ciperc = 0.95,
+  ci_limit_ratio_tol = 1.5,
+  verbose = FALSE,
+  ...
+)
+}
+\arguments{
+\item{i}{The position of the target parameters as in the parameter table of lavaan.}
+
+\item{npar}{The number of free parameters, including those constrained to be equal.}
+
+\item{sem_out}{The fit object.}
+
+\item{f_constr}{The constraint function generated by \code{set_constraint}.}
+
+\item{which}{Whether the lower bound or the upper bound is to be found. Must be "lbound" or "ubound".}
+
+\item{history}{If \code{TRUE}, return the full optimization output, stored in the \code{\link{diag}} attribute.
+Default is \code{FALSE}. This argument is no longer used. \code{\link{verbose}} will be used instead.}
+
+\item{perturbation_factor}{A number multiplied to the parameter estimates in the fit object.
+Using the parameter estimates as starting values may lead to errors
+in the first few iterations. Default is .90.}
+
+\item{lb_var}{The lower bound for free parameters that are variances. Defautl is -Inf. This is not an admissible
+value but seems to be necessary when we need to find the LBCI for user-defined parameter.}
+
+\item{wald_ci_start}{If TRUE and there are no equality constraints in the model, Wald confidence limit will
+be used as the starting value for the target parameter if it is not a userd-defined paramter.}
+
+\item{standardized}{If TRUE, the LBCI is for the standardized estimate.}
+
+\item{opts}{Options to be passed to \code{nloptr}}
+
+\item{test_generic}{If \code{\link{TRUE}}, will use the work-in-progress generic functions. Default is \code{\link{TRUE}}.
+The generic functions are now the default and other functions should not be used.
+Therefore, this argument will be ignored.}
+
+\item{ciperc}{The proportion of coverage for the confidence interval. Default
+is .95.}
+
+\item{ci_limit_ratio_tol}{The tolerance for the ratio of a to b, where a is the distance between an LBCI limit and the
+point estimate, and the b is the distance between the original \code{\link{lavaan}} CI limit (by default the Wald CI)
+and the point estimate. If the ratio is larger than this value, a warning is set in the status code.
+Default is 1.5.}
+
+\item{verbose}{If \code{\link{TRUE}}, will store more diagnostic information in the attribute \code{\link{diag}}.GlobalEnv
+Default is \code{\link{FALSE}}.}
+
+\item{...}{Optional arguments. Not used.}
+}
+\value{
+The requested bound.
+Can return the optimization history as an attribute.
+}
+\description{
+Find the lower or upper bound for one parameter by the Neale-Miller-1997 approach.
+}
+\details{
+Currently supports \code{lavaan} output only.
+}
+\examples{
+library(lavaan)
+data(cfa_two_factors)
+mod <-
+"
+f1 =~ x1 + x2 + a*x3
+f2 =~ x4 + a*x5 + equal('f1=~x2')*x6
+f1 ~~ 0*f2
+asq := a^2
+"
+fit <- sem(mod, cfa_two_factors)
+}
diff --git a/man/ci_i.Rd b/man/ci_i.Rd
index a7412c95..4d7426ee 100644
--- a/man/ci_i.Rd
+++ b/man/ci_i.Rd
@@ -4,11 +4,14 @@
 \alias{ci_i}
 \title{Find the lower bound and upper bound for one parameter}
 \usage{
-ci_i(i, ...)
+ci_i(i, method = "wn", ...)
 }
 \arguments{
 \item{i}{The position of the target parameters.}
 
+\item{method}{The approach to be used. Can be "wn" (Wu-Neale-2012) or "nm"
+(Neale-Miller-1997). Default is "wn".}
+
 \item{...}{Arguments to be passed to \code{ci_bound_i}.}
 }
 \value{
diff --git a/man/find_dependent.Rd b/man/find_dependent.Rd
new file mode 100644
index 00000000..33adfac8
--- /dev/null
+++ b/man/find_dependent.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/find_dependent.R
+\name{find_dependent}
+\alias{find_dependent}
+\title{Find the free parameters on which a parameter depends on}
+\usage{
+find_dependent(i = NULL, sem_out = NULL, standardized = FALSE, signed = FALSE)
+}
+\arguments{
+\item{i}{The position of the target parameter as in the parameter table of lavaan.}
+
+\item{sem_out}{The SEM output. Currently \code{lavaan} output only.}
+
+\item{standardized}{If TRUE, the LBCI is for the standardized estimate.}
+
+\item{signed}{If TRUE, return a vector of 1 or -1 to indicate the direction of the dependence.
+Default is \code{\link{FALSE}}.}
+}
+\value{
+A numeric vector of the positions of the free parameters in the \code{lavaan} parameter table.
+}
+\description{
+Find the free parameters on which a parameter depends on
+}
+\details{
+Currently supports \code{lavaan} output only.
+}
+\examples{
+library(lavaan)
+data(cfa_two_factors)
+mod <- 
+"
+f1 =~ x1 + x2 + a*x3
+f2 =~ x4 + a*x5 + equal('f1=~x2')*x6
+f1 ~~ 0*f2
+asq := a^2
+"
+fit <- sem(mod, cfa_two_factors)
+}
diff --git a/man/reg_cor_near_one.Rd b/man/reg_cor_near_one.Rd
new file mode 100644
index 00000000..cabfac87
--- /dev/null
+++ b/man/reg_cor_near_one.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{reg_cor_near_one}
+\alias{reg_cor_near_one}
+\title{Generated from a regression model, with one correlation close to one.}
+\format{
+An object of class \code{data.frame} with 100 rows and 6 columns.
+}
+\usage{
+reg_cor_near_one
+}
+\description{
+Generated from a regression model, with one correlation close to one.
+}
+\keyword{datasets}
diff --git a/man/semlbci.Rd b/man/semlbci.Rd
index d37aeb9c..1bb0e296 100644
--- a/man/semlbci.Rd
+++ b/man/semlbci.Rd
@@ -5,7 +5,16 @@
 \alias{semlbci}
 \title{semlbci: Likelihood Based Confidence Interval for SEM}
 \usage{
-semlbci(sem_out, pars = NULL, ciperc = 0.95, ..., parallel = FALSE, ncpu = 2)
+semlbci(
+  sem_out,
+  pars = NULL,
+  ciperc = 0.95,
+  standardized = FALSE,
+  method = "nm",
+  ...,
+  parallel = FALSE,
+  ncpu = 2
+)
 }
 \arguments{
 \item{sem_out}{The SEM output. Currently \code{lavaan} output only.}
@@ -19,6 +28,11 @@ parameter table.}
 \item{ciperc}{The proportion of coverage for the confidence interval. Default
 is .95.}
 
+\item{standardized}{If TRUE, the LBCI is for the standardized estimate.}
+
+\item{method}{The approach to be used. Can be "wn" (Wu-Neale-2012) or "nm"
+(Neale-Miller-1997). Default is "nm".}
+
 \item{...}{Arguments to be passed to \code{ci_bound_i}.}
 
 \item{parallel}{If \code{TRUE}, will use parallel. Currently disabled.
diff --git a/man/set_constraint_nm.Rd b/man/set_constraint_nm.Rd
new file mode 100644
index 00000000..dddd88b5
--- /dev/null
+++ b/man/set_constraint_nm.Rd
@@ -0,0 +1,51 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/set_constraint_nm.R
+\name{set_constraint_nm}
+\alias{set_constraint_nm}
+\title{Set the constraint for finding the LBCI by the Neale-Miller-1997 approach}
+\usage{
+set_constraint_nm(
+  i,
+  sem_out,
+  ciperc = 0.95,
+  envir = NULL,
+  get_fit_from_envir = FALSE
+)
+}
+\arguments{
+\item{i}{The position of the target parameters as in the parameter table of lavaan.}
+
+\item{sem_out}{The SEM output. Currently \code{lavaan} output only.}
+
+\item{ciperc}{The proportion of coverage for the confidence interval. Default
+is .95.}
+
+\item{envir}{The enviroment to stroe the state. Default is NULL}
+}
+\value{
+A constraint function for nloptr.
+}
+\description{
+Set the constraint for finding the LBCI by the Neale-Miller-1997 approach
+}
+\details{
+Currently supports \code{lavaan} output only.
+}
+\examples{
+library(lavaan)
+data(cfa_two_factors)
+mod <-
+"
+f1 =~ x1 + x2 + a*x3
+f2 =~ x4 + a*x5 + equal('f1=~x2')*x6
+f1 ~~ 0*f2
+asq := a^2
+"
+fit <- sem(mod, cfa_two_factors)
+#fn1 <- set_constraint(sem_out)
+
+#fn1(coef(sem_out))
+#fn1(runif(length(coef(sem_out))), .9, 1.1) * coef(sem_out))
+#fn1(coef(sem_out) + 2)
+#fn1(coef(sem_out) - .12)
+}
diff --git a/man/set_start.Rd b/man/set_start.Rd
index 4a95b9fd..f0037683 100644
--- a/man/set_start.Rd
+++ b/man/set_start.Rd
@@ -4,7 +4,7 @@
 \alias{set_start}
 \title{Set the starting values for optimization}
 \usage{
-set_start(i = NULL, sem_out = NULL, which = NULL)
+set_start(i = NULL, sem_out = NULL, which = NULL, standardized = FALSE)
 }
 \arguments{
 \item{i}{The position of the target parameters as in the parameter table of lavaan.}
@@ -12,6 +12,8 @@ set_start(i = NULL, sem_out = NULL, which = NULL)
 \item{sem_out}{The SEM output. Currently \code{lavaan} output only.}
 
 \item{which}{Whether the lower bound or the upper bound is to be found. Must be "lbound" or "ubound".}
+
+\item{standardized}{If TRUE, the LBCI is for the standardized estimate.}
 }
 \value{
 A lavaan parameter table, with parameters estimated with target fixed to its
diff --git a/man/update_model.Rd b/man/update_model.Rd
index 67464009..faa0b044 100644
--- a/man/update_model.Rd
+++ b/man/update_model.Rd
@@ -15,7 +15,7 @@ update_model(model, est)
 An updated lavModel object.
 }
 \description{
-
+Update a \code{lavaan} model using user-supplied estimates
 }
 \details{
 Modified from lav_model_set_parameters in lavaan.
diff --git a/tests/testthat/test-ci_bound_i_diagnostic_info.R b/tests/testthat/test-ci_bound_i_diagnostic_info.R
new file mode 100644
index 00000000..93236b8d
--- /dev/null
+++ b/tests/testthat/test-ci_bound_i_diagnostic_info.R
@@ -0,0 +1,62 @@
+library(testthat)
+library(semlbci)
+
+context("Check ci_bound_i: Check diagnostic info")
+
+data(simple_med)
+dat <- simple_med
+mod <- 
+"
+m ~ x
+y ~ m
+"
+fit_med <- lavaan::sem(mod, simple_med, fixed.x = FALSE)
+
+fn_constr0 <- set_constraint(fit_med)
+
+# opts0 <- list(print_level = 3)
+opts0 <- list()
+opts0 <- list(ftol_abs = 1e-7,
+              ftol_rel = 1e-7,
+              xtol_abs = 1e-7,
+              xtol_rel = 1e-7,
+              tol_constraints_eq = 1e-7
+              )
+system.time(out1l <- ci_bound_i(1, 5, sem_out = fit_med, f_constr = fn_constr0, which = "lbound", opts = opts0, verbose = TRUE))
+system.time(out1u <- ci_bound_i(1, 5, sem_out = fit_med, f_constr = fn_constr0, which = "ubound", opts = opts0, verbose = TRUE))
+# system.time(out1ls <- ci_bound_i(2, 5, sem_out = fit_med, f_constr = fn_constr0, which = "lbound", opts = opts0, standardized = TRUE, verbose = TRUE))
+# system.time(out1us <- ci_bound_i(2, 5, sem_out = fit_med, f_constr = fn_constr0, which = "ubound", opts = opts0, standardized = TRUE, verbose = TRUE))
+
+out1l_diag <- attr(out1l, "diag")
+# out1ls_diag <- attr(out1ls, "diag")
+
+test_that("Check diagnostic object: status", {
+    expect_equivalent(
+        out1l_diag$status, 0
+      )
+  })
+
+test_that("Check diagnostic object: est_org", {
+    expect_true(
+        is.numeric(out1l_diag$est_org) & length(out1l_diag$est_org) == 1
+      )
+  })
+
+test_that("Check diagnostic object: ci_org_limit", {
+    expect_true(
+        is.numeric(out1l_diag$ci_org_limit) & length(out1l_diag$ci_org_limit) == 1
+      )
+  })
+
+test_that("Check diagnostic object: ci_limit_ratio", {
+    expect_equivalent(
+        out1l_diag$ci_limit_ratio,
+        abs((out1l - out1l_diag$est_org) / (out1l_diag$ci_org_limit - out1l_diag$est_org))
+      )
+  })
+
+test_that("Check diagnostic object: fit_post_check", {
+    expect_true(
+        isTRUE(out1l_diag$fit_post_check) | isFALSE(out1l_diag$fit_post_check)
+      )
+  })
\ No newline at end of file
diff --git a/tests/testthat/test-ci_bound_nm_i_generic_eq_constraints.R b/tests/testthat/test-ci_bound_nm_i_generic_eq_constraints.R
index d1b86d48..cbb7ef63 100644
--- a/tests/testthat/test-ci_bound_nm_i_generic_eq_constraints.R
+++ b/tests/testthat/test-ci_bound_nm_i_generic_eq_constraints.R
@@ -29,8 +29,8 @@ system.time(out2u <- ci_bound_nm_i(2, 5, sem_out = fit_med, which = "ubound", op
 system.time(out3l <- ci_bound_nm_i(6, 5, sem_out = fit_med, which = "lbound", opts = opts0, test_generic = TRUE))
 system.time(out3u <- ci_bound_nm_i(6, 5, sem_out = fit_med, which = "ubound", opts = opts0, test_generic = TRUE))
 
-system.time(out1sl <- ci_bound_nm_i(3, 5, sem_out = fit_med, which = "lbound", opts = opts0, test_generic = TRUE, history = TRUE, standardized = TRUE))
-system.time(out1su <- ci_bound_nm_i(3, 5, sem_out = fit_med, which = "ubound", opts = opts0, test_generic = TRUE, history = TRUE, standardized = TRUE))
+# system.time(out1sl <- ci_bound_nm_i(3, 5, sem_out = fit_med, which = "lbound", opts = opts0, test_generic = TRUE, history = TRUE, standardized = TRUE))
+# system.time(out1su <- ci_bound_nm_i(3, 5, sem_out = fit_med, which = "ubound", opts = opts0, test_generic = TRUE, history = TRUE, standardized = TRUE))
 
 # (ci_semlbci <- c(out1l, out2l, out1u, out2u)) -
 #   unlist(ci_OpenMx[c("a", "b"), c("lbound", "ubound")])
diff --git a/tests/testthat/test-ci_i_diagnostic_info.R b/tests/testthat/test-ci_i_diagnostic_info.R
new file mode 100644
index 00000000..7c484d7b
--- /dev/null
+++ b/tests/testthat/test-ci_i_diagnostic_info.R
@@ -0,0 +1,57 @@
+library(testthat)
+library(semlbci)
+
+context("Check ci_i: Check diagnostic info")
+
+data(simple_med)
+dat <- simple_med
+mod <- 
+"
+m ~ x
+y ~ m
+"
+fit_med <- lavaan::sem(mod, simple_med, fixed.x = FALSE)
+
+fn_constr0 <- set_constraint(fit_med)
+
+# opts0 <- list(print_level = 3)
+opts0 <- list()
+opts0 <- list(ftol_abs = 1e-7,
+              ftol_rel = 1e-7,
+              xtol_abs = 1e-7,
+              xtol_rel = 1e-7,
+              tol_constraints_eq = 1e-7
+              )
+system.time(out1 <- ci_i(1, npar = 5, sem_out = fit_med, method = "wn", f_constr = fn_constr0, opts = opts0, verbose = TRUE))
+system.time(out1_nm <- ci_i(1, 5, sem_out = fit_med, method = "nm", opts = opts0, verbose = TRUE))
+# system.time(out1ls <- ci_bound_i(2, 5, sem_out = fit_med, f_constr = fn_constr0, which = "lbound", opts = opts0, standardized = TRUE, verbose = TRUE))
+# system.time(out1us <- ci_bound_i(2, 5, sem_out = fit_med, f_constr = fn_constr0, which = "ubound", opts = opts0, standardized = TRUE, verbose = TRUE))
+
+out1_lb_diag <- attr(out1, "lb_diag")
+out1_ub_diag <- attr(out1, "ub_diag")
+out1_nm_lb_diag <- attr(out1_nm, "lb_diag")
+out1_nm_ub_diag <- attr(out1_nm, "ub_diag")
+
+test_that("Check method", {
+    expect_true(
+        attr(out1, "method") == "wn"
+      )
+  })
+
+test_that("Check method", {
+    expect_true(
+        attr(out1_nm, "method") == "nm"
+      )
+  })
+
+test_that("Check existence of diagnostic objects", {
+    expect_true(
+        !is.null(out1_lb_diag) &
+        !is.null(out1_ub_diag) &
+        !is.null(out1_nm_lb_diag) &
+        !is.null(out1_nm_ub_diag) 
+      )
+  })
+
+
+
diff --git a/tests/testthat/test-extra_arguments.R b/tests/testthat/test-extra_arguments.R
index 5093abc9..f94e0527 100644
--- a/tests/testthat/test-extra_arguments.R
+++ b/tests/testthat/test-extra_arguments.R
@@ -37,12 +37,12 @@ fit_med_OpenMx <- mxRun(mod_mx, silent = TRUE, intervals = TRUE)
 ci_OpenMx <- summary(fit_med_OpenMx)$CI
 ci_OpenMx_target <- unlist(ci_OpenMx[c("a"), c("lbound", "ubound")])
 
-lbci_debug_TRUE <- semlbci(fit_med, pars = c(1), debug = TRUE)
-lbci_debug_FALSE <- semlbci(fit_med, pars = c(1), debug = FALSE)
-lbci_lav_warn_TRUE <- semlbci(fit_med, pars = c(1), lav_warn = TRUE)
-lbci_lav_warn_FALSE <- semlbci(fit_med, pars = c(1), lav_warn = FALSE)
-lbci_history_TRUE <- semlbci(fit_med, pars = c(1), history = TRUE)
-lbci_history_FALSE <- semlbci(fit_med, pars = c(1), history = FALSE)
+lbci_debug_TRUE <- semlbci(fit_med, pars = c(1), debug = TRUE, method = "wn")
+lbci_debug_FALSE <- semlbci(fit_med, pars = c(1), debug = FALSE, method = "wn")
+lbci_lav_warn_TRUE <- semlbci(fit_med, pars = c(1), lav_warn = TRUE, method = "wn")
+lbci_lav_warn_FALSE <- semlbci(fit_med, pars = c(1), lav_warn = FALSE, method = "wn")
+lbci_history_TRUE <- semlbci(fit_med, pars = c(1), history = TRUE, method = "wn")
+lbci_history_FALSE <- semlbci(fit_med, pars = c(1), history = FALSE, method = "wn")
 
 test_that("Can run with debug = TRUE", {
     expect_equivalent(
diff --git a/tests/testthat/test-semlbci.R b/tests/testthat/test-semlbci.R
index 2b0a0aad..8c098725 100644
--- a/tests/testthat/test-semlbci.R
+++ b/tests/testthat/test-semlbci.R
@@ -39,7 +39,7 @@ mod_mx <- mxModel("Mediation", type = "RAM",
 fit_med_OpenMx <- mxRun(mod_mx, silent = TRUE, intervals = TRUE)
 ci_OpenMx <- summary(fit_med_OpenMx)$CI
 
-lbci_med <- semlbci(fit_med)
+lbci_med <- semlbci(fit_med, method = "wn")
 
 test_that("Equal to OpenMx LBCIs for free parameters", {
     expect_equivalent(
@@ -49,7 +49,7 @@ test_that("Equal to OpenMx LBCIs for free parameters", {
       )
   })
 
-lbci_med2 <- semlbci(fit_med, pars = c(1, 3))
+lbci_med2 <- semlbci(fit_med, pars = c(1, 3), method = "wn")
 
 test_that("Check whether only selectd parameters were processed", {
     expect_equivalent(
diff --git a/tests/testthat/test-semlbci_derived_parameter.R b/tests/testthat/test-semlbci_derived_parameter.R
index ba9045a6..c3f2c6d3 100644
--- a/tests/testthat/test-semlbci_derived_parameter.R
+++ b/tests/testthat/test-semlbci_derived_parameter.R
@@ -43,7 +43,7 @@ mod_mx <- mxModel("Mediation", type = "RAM",
 fit_med_OpenMx <- mxRun(mod_mx, silent = TRUE, intervals = TRUE)
 ci_OpenMx <- summary(fit_med_OpenMx)$CI
 
-lbci_med <- semlbci(fit_med, pars = 6) 
+lbci_med <- semlbci(fit_med, pars = 6, method = "wn") 
 
 test_that("Equal to OpenMx LBCI", {
     expect_equivalent(
diff --git a/tests/testthat/test-semlbci_diagnostic.R b/tests/testthat/test-semlbci_diagnostic.R
new file mode 100644
index 00000000..d99c8f23
--- /dev/null
+++ b/tests/testthat/test-semlbci_diagnostic.R
@@ -0,0 +1,70 @@
+library(testthat)
+library(semlbci)
+
+context("Check semlbci: Diagnostic info")
+
+data(simple_med)
+dat <- simple_med
+# mod <- 
+# "
+# m ~ x
+# y ~ m
+# "
+# fit_med <- lavaan::sem(mod, simple_med, fixed.x = FALSE)
+
+# library(OpenMx)
+# cov_dat <- cov(dat)
+# n <- nrow(dat)
+# manifest_vars <- c("x", "m", "y")
+# mod_mx <- mxModel("Mediation", type = "RAM", 
+#     manifestVars = manifest_vars,
+#     mxPath(from = "x", to = "m", arrows = 1, free = TRUE, values = 1,
+#                   labels = "a"),
+#     mxPath(from = "m", to = "y", arrows = 1, free = TRUE, values = 1,
+#                   labels = "b"),
+#     mxPath(from = "x", arrows = 2, free = TRUE, values = 1,
+#                   labels = "var_x"),
+#     mxPath(from = "m", arrows = 2, free = TRUE, values = 1,
+#                   labels = "evar_m"),
+#     mxPath(from = "y", arrows = 2, free = TRUE, values = 1,
+#                   labels = "evar_y"),
+#     mxAlgebra(a * b , name = "ab"),
+#     mxCI(reference = c("a", "b", "ab", "evar_m"), 
+#                        interval = .95, type = "both"),
+#     mxData(observed = cov_dat, type = "cov", numObs = n)
+#   )
+# fit_med_OpenMx <- mxRun(mod_mx, silent = TRUE, intervals = TRUE)
+# ci_OpenMx <- summary(fit_med_OpenMx)$CI
+
+# lbci_med <- semlbci(fit_med, method = "nm")
+# lbci_med_std <- semlbci(fit_med, pars = c("m ~  x",
+#                                           "m ~~ m",
+#                                           "y ~~ y"),
+#                         method = "nm", standardized = TRUE)
+# lbci_med_wn <- semlbci(fit_med, method = "wn")
+# lbci_med_wn_std <- semlbci(fit_med, pars = c("m ~  x",
+#                                              "y ~  m",
+#                                              "y ~~ y"),
+#                             method = "wn", standardized = TRUE)
+
+mod2 <- 
+"
+m ~ a*x
+y ~ b*m
+ab := a * b
+"
+fit_med2 <- lavaan::sem(mod2, simple_med, fixed.x = FALSE)
+
+# lbci_med2 <- semlbci(fit_med2, method = "nm")
+# lbci_med2 <- semlbci(fit_med2, pars = c("m ~x", "y ~ m", "ab :="), method = "nm")
+# lbci_med2 <- semlbci(fit_med2, pars = c("m ~x", "y ~ m", "ab :="), method = "nm", verbose = TRUE)
+# lbci_med2_std <- semlbci(fit_med2, pars = c("m ~ x",
+#                                             "y ~ m",
+#                                             "ab := "),
+#                         method = "nm", standardized = TRUE)
+# lbci_med2_wn <- semlbci(fit_med2, method = "wn")
+# lbci_med2_wn_std <- semlbci(fit_med2, pars = c("m ~  x",
+#                                                "y ~  m",
+#                                                "ab := "),
+#                             method = "wn", standardized = TRUE)
+
diff --git a/tests/testthat/test-semlbci_eq_constr.R b/tests/testthat/test-semlbci_eq_constr.R
index 62578c4a..f00fe606 100644
--- a/tests/testthat/test-semlbci_eq_constr.R
+++ b/tests/testthat/test-semlbci_eq_constr.R
@@ -37,7 +37,7 @@ mod_mx <- mxModel("Mediation", type = "RAM",
 fit_med_OpenMx <- mxRun(mod_mx, silent = TRUE, intervals = TRUE)
 ci_OpenMx <- summary(fit_med_OpenMx)$CI
 
-lbci_med <- semlbci(fit_med)
+lbci_med <- semlbci(fit_med, method = "wn")
 
 test_that("Equality constraints are observed", {
     expect_equivalent(
diff --git a/tests/testthat/test-semlbci_text_input.R b/tests/testthat/test-semlbci_text_input.R
index 3a731ba2..6626a61e 100644
--- a/tests/testthat/test-semlbci_text_input.R
+++ b/tests/testthat/test-semlbci_text_input.R
@@ -47,7 +47,7 @@ lbci_med <- semlbci(fit_med, pars = c("m ~ x",
                                       "y ~ m",
                                       "asq := 1",
                                       "ab  := 2",
-                                      "not in table")) 
+                                      "not in table"), method = "wn") 
 
 test_that("Equal to OpenMx LBCI", {
     expect_equivalent(