diff --git a/R/main_function_documentation.R b/R/main_function_documentation.R
index 77dd170..40669bb 100644
--- a/R/main_function_documentation.R
+++ b/R/main_function_documentation.R
@@ -188,8 +188,6 @@ NULL
 #'  \item{\code{output} -- output of the model with information on the estimated population mean and standard errors.}
 #'  \item{\code{SE} -- standard error of the estimator of the population mean, divided into errors from probability and non-probability samples.}
 #'  \item{\code{confidence_interval} -- confidence interval of population mean estimator}
-#'  \item{\code{coeff_selection} -- `data.frame` of the estimated coefficients of the propensity score (inverse probability weighting) model and their standard errors.}
-#'  \item{\code{coeff_outcome} -- `data.frame` of the estimated coefficients of the outcome (mass imputation) model and their standard errors.}
 #'  \item{\code{nonprob_size} -- size of non-probability sample}
 #'  \item{\code{prob_size} -- size of probability sample}
 #'  \item{\code{pop_size} -- estimated population size derived from estimated weights (non-probability sample) or known design weights (probability sample)}
diff --git a/R/nonprobDR.R b/R/nonprobDR.R
index 5e42c33..27e1d4d 100644
--- a/R/nonprobDR.R
+++ b/R/nonprobDR.R
@@ -39,6 +39,7 @@ nonprobDR <- function(selection,
     message("Bootstrap variance only, analytical version during implementation")
     control_inference$var_method <- "bootstrap"
   }
+  outcome_init <- outcome
   h <- control_selection$h
   maxit <- control_selection$maxit
   optim_method <- control_selection$optim_method
@@ -633,7 +634,7 @@ nonprobDR <- function(selection,
     } else {
       stop("Please, provide only one of svydesign object or pop_totals/pop_means.")
     }
-    ys$k <- as.numeric(y_nons) # TODO name to change
+    ys[[k]] <- as.numeric(y_nons)
 
     if (se) {
       if (var_method == "analytic") { # TODO add estimator variance with model containg pop_totals to internal_varDR function
@@ -777,6 +778,7 @@ nonprobDR <- function(selection,
   names(OutcomeList) <- outcomes$f
   if (is.null(pop_size)) pop_size <- N_nons
   names(pop_size) <- "pop_size"
+  names(ys) <- all.vars(outcome_init[[2]])
 
   SelectionList <- list(
     coefficients = selection_model$theta_hat,
@@ -801,7 +803,7 @@ nonprobDR <- function(selection,
   structure(
     list(
       X = if (isTRUE(x)) X else NULL,
-      y = if (isTRUE(y)) as.numeric(y) else NULL,
+      y = if (isTRUE(y)) ys else NULL,
       prob = prop_scores,
       weights = as.vector(weights_nons),
       control = list(
diff --git a/R/nonprobIPW.R b/R/nonprobIPW.R
index 59b5866..e2a9c01 100644
--- a/R/nonprobIPW.R
+++ b/R/nonprobIPW.R
@@ -44,7 +44,7 @@ nonprobIPW <- function(selection,
   if (!(target[3] == "NULL()")) stop("ill-defined formula for the target")
   # formula for outcome variable if target defined
   dependents <- paste(selection, collapse = " ")
-  outcome <- stats::as.formula(paste(target[2], dependents))
+  outcome <- outcome_init <- stats::as.formula(paste(target[2], dependents))
   outcomes <- ff(outcome)
   output <- list()
   ys <- list()
@@ -365,7 +365,7 @@ nonprobIPW <- function(selection,
         flag = FALSE
       )$y_nons
     }
-    ys$k <- as.numeric(y_nons) # TODO name to change
+    ys[[k]] <- as.numeric(y_nons)
     mu_hat <- mu_hatIPW(
       y = y_nons,
       weights = weights,
@@ -503,6 +503,7 @@ nonprobIPW <- function(selection,
   rownames(output) <- rownames(confidence_interval) <- rownames(SE_values) <- outcomes$f
   if (is.null(pop_size)) pop_size <- N # estimated pop_size
   names(pop_size) <- "pop_size"
+  names(ys) <- all.vars(outcome_init[[2]])
 
   SelectionList <- list(
     coefficients = selection_model$theta_hat,
diff --git a/R/nonprobMI.R b/R/nonprobMI.R
index b137aa3..4926b48 100644
--- a/R/nonprobMI.R
+++ b/R/nonprobMI.R
@@ -32,6 +32,7 @@ nonprobMI <- function(outcome,
                       se,
                       ...) {
   var_selection <- control_inference$vars_selection
+  outcome_init <- outcome
   outcomes <- ff(outcome)
   output <- list()
   ys <- list()
@@ -210,7 +211,7 @@ nonprobMI <- function(outcome,
     } else {
       stop("Please, provide svydesign object or pop_totals/pop_means.")
     }
-    ys$k <- as.numeric(y_nons) # TODO name to change
+    ys[[k]] <- as.numeric(y_nons)
     if (se) {
       # design based variance estimation based on approximations of the second-order inclusion probabilities
       if (control_inference$var_method == "analytic") { # consider move variance implementation to internals
@@ -317,11 +318,12 @@ nonprobMI <- function(outcome,
   rownames(output) <- rownames(confidence_interval) <- rownames(SE_values) <- outcomes$f
   names(OutcomeList) <- outcomes$f
   names(pop_size) <- "pop_size"
+  names(ys) <- all.vars(outcome_init[[2]])
 
   structure(
     list(
       X = if (isTRUE(x)) X else NULL,
-      y = if (isTRUE(y)) as.numeric(y) else NULL,
+      y = if (isTRUE(y)) ys else NULL,
       control = list(
         control_outcome = control_outcome,
         control_inference = control_inference
diff --git a/man/nonprob.Rd b/man/nonprob.Rd
index 56c374a..b774df7 100644
--- a/man/nonprob.Rd
+++ b/man/nonprob.Rd
@@ -97,8 +97,6 @@ with type \code{list} containing:\cr
 \item{\code{output} -- output of the model with information on the estimated population mean and standard errors.}
 \item{\code{SE} -- standard error of the estimator of the population mean, divided into errors from probability and non-probability samples.}
 \item{\code{confidence_interval} -- confidence interval of population mean estimator}
-\item{\code{coeff_selection} -- \code{data.frame} of the estimated coefficients of the propensity score (inverse probability weighting) model and their standard errors.}
-\item{\code{coeff_outcome} -- \code{data.frame} of the estimated coefficients of the outcome (mass imputation) model and their standard errors.}
 \item{\code{nonprob_size} -- size of non-probability sample}
 \item{\code{prob_size} -- size of probability sample}
 \item{\code{pop_size} -- estimated population size derived from estimated weights (non-probability sample) or known design weights (probability sample)}