diff --git a/NAMESPACE b/NAMESPACE
index eae8719b..4393f6ed 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -276,3 +276,4 @@ importFrom(nlme,getGroups)
 importFrom(rlang,.data)
 importFrom(stats,formula)
 importFrom(stats,model.frame)
+importFrom(stats,na.omit)
diff --git a/NEWS.md b/NEWS.md
index 0781cfc1..95252ead 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -87,6 +87,7 @@ when the issue is due to an excluded point (#310)
   `clr.pred.dn` (#433)
 * `PKNCA.set.summary(reset = TRUE)` warns that it may break the use of
   `summary()` (#477)
+* `pk.nca` output now includes a `PPANMETH` column describing the analysis methods used for each parameter regarding imputations, AUC and half.life calculations (#457)
 * Added new `tmin` parameter
 * New post-processing functions to normalize PKNCA result parameters based on any column in PKNCAconc data.frame (`normalize_by_col()`) or by using a custom normalization table (`normalize()`)
 * New excretion rate parameters: `ermax`  (Maximum excretion rate), `ertmax` (Midpoint time of maximum excretion rate) and `ertlst` (Time of last excretion rate measurement) (#433)
diff --git a/R/auc.R b/R/auc.R
index 19e4f685..3d7e6c6b 100644
--- a/R/auc.R
+++ b/R/auc.R
@@ -179,6 +179,8 @@ pk.calc.auxc <- function(conc, time, interval=c(0, Inf),
         fun_linear = fun_linear, fun_log = fun_log, fun_inf = fun_inf
       )
   }
+  # Add method details as an attribute
+  attr(ret, "method") <- paste0("AUC: ", method)
   ret
 }
 
diff --git a/R/aucint.R b/R/aucint.R
index d9978bdd..1acad839 100644
--- a/R/aucint.R
+++ b/R/aucint.R
@@ -175,6 +175,9 @@ pk.calc.auxcint <- function(conc, time,
       fun_log = fun_log,
       fun_inf = fun_inf
     )
+  # Add method details as an attribute
+  attr(ret, "method") <- paste0("AUC: ", method)
+
   ret
 }
 
diff --git a/R/class-PKNCAresults.R b/R/class-PKNCAresults.R
index e83f2dfb..cb63e4b8 100644
--- a/R/class-PKNCAresults.R
+++ b/R/class-PKNCAresults.R
@@ -89,7 +89,7 @@ as.data.frame.PKNCAresults <- function(x, ..., out_format = c('long', 'wide'), f
     # Since we moved the results into PPTESTCD and PPORRES regardless of what
     # they really are in the source data, remove the extra units and unit
     # conversion columns to allow spread to work.
-    ret <- ret[, setdiff(names(ret), c("PPSTRES", "PPSTRESU", "PPORRESU"))]
+    ret <- ret[, setdiff(names(ret), c("PPSTRES", "PPSTRESU", "PPORRESU", "PPANMETH"))]
     ret <- tidyr::spread(ret, key="PPTESTCD", value="PPORRES")
   }
   ret
diff --git a/R/half.life.R b/R/half.life.R
index 989775c1..79f50bd0 100644
--- a/R/half.life.R
+++ b/R/half.life.R
@@ -271,6 +271,7 @@ pk.calc.half.life <- function(conc, time, tmax, tlast,
     }
 
     if (manually.selected.points) {
+      attr(ret, "method") <- "Lambda Z: Manual selection"
       if (nrow(data) > 0) {
         fit <- fit_half_life(data=data, tlast=ret$tlast)
         ret[, ret_replacements] <- fit[, ret_replacements]
@@ -360,6 +361,7 @@ pk.calc.half.life <- function(conc, time, tmax, tlast,
     n_above_lloq <- sum(!dfK_all$mask_blq)
 
     if (manually.selected.points) {
+      attr(ret, "method") <- "Lambda Z: Manual selection"
       # Use data_tobit as-is (all non-NA points, no tmax filter applied again)
       if (nrow(data_tobit) > 0) {
         fit <- fit_half_life_tobit(
diff --git a/R/pk.calc.all.R b/R/pk.calc.all.R
index c5faddba..6021aa40 100644
--- a/R/pk.calc.all.R
+++ b/R/pk.calc.all.R
@@ -361,6 +361,7 @@ pk.nca.intervals <- function(data_conc, data_dose, data_intervals, sparse,
 #'
 #' @seealso [check.interval.specification()]
 #' @export
+#' @importFrom stats na.omit
 pk.nca.interval <- function(conc, time, volume, duration.conc,
                             dose, time.dose, duration.dose, route,
                             conc.group=NULL, time.group=NULL, volume.group=NULL, duration.conc.group=NULL,
@@ -389,6 +390,9 @@ pk.nca.interval <- function(conc, time, volume, duration.conc,
     }
     conc <- impute_data$conc
     time <- impute_data$time
+    tmp_imp_method <- paste0("Imputation: ", paste(na.omit(impute_method), collapse = ", "))
+  } else {
+    tmp_imp_method <- character()
   }
   # Prepare the return value using SDTM names
   ret <- data.frame(PPTESTCD=NA, PPORRES=NA)[-1,]
@@ -528,6 +532,10 @@ pk.nca.interval <- function(conc, time, volume, duration.conc,
         } else {
           NA_character_
         }
+      # The handling of the method column (PPANMETH)
+      tmp_method <- c(tmp_imp_method, attr(tmp_result, "method"))
+      attr(tmp_result, "method") <- NULL
+
       # If the function returns a data frame, save all the returned values,
       # otherwise, save the value returned.
       if (is.data.frame(tmp_result)) {
@@ -540,6 +548,7 @@ pk.nca.interval <- function(conc, time, volume, duration.conc,
         data.frame(
           PPTESTCD=tmp_testcd,
           PPORRES=tmp_result,
+          PPANMETH=paste(tmp_method, collapse=". "),
           exclude=exclude_reason,
           stringsAsFactors=FALSE
         )
diff --git a/R/sparse.R b/R/sparse.R
index 8157acee..dbddf474 100644
--- a/R/sparse.R
+++ b/R/sparse.R
@@ -302,7 +302,7 @@ sparse_to_dense_pk <- function(sparse_pk) {
 #' @family Sparse Methods
 #' @export
 pk.calc.sparse_auc <- function(conc, time, subject,
-                               method=NULL,
+                               method="linear",
                                auc.type="AUClast",
                                ...,
                                options=list()) {
@@ -314,15 +314,23 @@ pk.calc.sparse_auc <- function(conc, time, subject,
       conc=sparse_pk_attribute(sparse_pk_mean, "mean"),
       time=sparse_pk_attribute(sparse_pk_mean, "time"),
       auc.type=auc.type,
-      method="linear"
+      method=method
     )
+
   var_auc <- var_sparse_auc(sparse_pk_mean)
-  data.frame(
+  ret <- data.frame(
     sparse_auc=auc,
     # as.numeric() drops the "df" attribute
     sparse_auc_se=sqrt(as.numeric(var_auc)),
     sparse_auc_df=attr(var_auc, "df")
   )
+
+  # Add method details as an attribute
+  for (col in names(ret)) {
+    attr(ret[[col]], "method") <- c(paste0("AUC: ", method), "Sparse: arithmetic mean, <=50% BLQ")
+  }
+
+  ret
 }
 
 #' @describeIn pk.calc.sparse_auc Compute the AUClast for sparse PK
diff --git a/man/pk.calc.sparse_auc.Rd b/man/pk.calc.sparse_auc.Rd
index e2677cbb..3f3189da 100644
--- a/man/pk.calc.sparse_auc.Rd
+++ b/man/pk.calc.sparse_auc.Rd
@@ -9,7 +9,7 @@ pk.calc.sparse_auc(
   conc,
   time,
   subject,
-  method = NULL,
+  method = "linear",
   auc.type = "AUClast",
   ...,
   options = list()
diff --git a/tests/testthat/test-auc.R b/tests/testthat/test-auc.R
index a6572db3..887dbc2c 100644
--- a/tests/testthat/test-auc.R
+++ b/tests/testthat/test-auc.R
@@ -54,9 +54,9 @@ test_that("pk.calc.auxc", {
 
 test_that("pk.calc.auc: Linear AUC when the conc at the end of the interval is above LOQ", {
   # lambda.z is unused
-  tests <- list(AUCinf=as.numeric(NA),
-                AUClast=1.5,
-                AUCall=1.5)
+  tests <- list(AUCinf=structure(as.numeric(NA), method="AUC: linear"),
+                AUClast=structure(1.5, method="AUC: linear"),
+                AUCall=structure(1.5, method="AUC: linear"))
   for (t in names(tests)) {
     # Note: using this structure ensures that there will not be
     # excessive warnings during testing.
@@ -78,9 +78,9 @@ test_that("pk.calc.auc: Linear AUC when the conc at the end of the interval is a
 test_that("pk.calc.auc: Linear AUC when the conc at the end of the interval is BLQ, lambda.z missing", {
   # lambda.z is used to extrapolate to the end of the interval.
   # Since lambda.z is NA, the result is NA.
-  tests <- list(AUCinf=as.numeric(NA),
-                AUClast=0.5,
-                AUCall=1)
+  tests <- list(AUCinf=structure(as.numeric(NA), method="AUC: linear"),
+                AUClast=structure(0.5, method="AUC: linear"),
+                AUCall=structure(1, method="AUC: linear"))
   for (t in names(tests)) {
     # Note: using this structure ensures that there will not be
     # excessive warnings during testing.
@@ -101,9 +101,9 @@ test_that("pk.calc.auc: Linear AUC when the conc at the end of the interval is B
 
 test_that("pk.calc.auc: Linear AUC when the conc at the end of the interval is BLQ, lambda.z given", {
   # The same when lambda.z is given
-  tests <- list(AUCinf=1.5,
-                AUClast=0.5,
-                AUCall=1)
+  tests <- list(AUCinf=structure(1.5, method="AUC: linear"),
+                AUClast=structure(0.5, method="AUC: linear"),
+                AUCall=structure(1, method="AUC: linear"))
   for (t in names(tests)) {
     # Note: using this structure ensures that there will not be
     # excessive warnings during testing.
@@ -123,9 +123,9 @@ test_that("pk.calc.auc: Linear AUC when the conc at the end of the interval is B
 })
 
 test_that("pk.calc.auc: Linear AUC when when there are multiple BLQ values at the end, lambda.z given", {
-  tests <- list(AUCinf=1.5,
-                AUClast=0.5,
-                AUCall=1)
+  tests <- list(AUCinf=structure(1.5, method="AUC: linear"),
+                AUClast=structure(0.5, method="AUC: linear"),
+                AUCall=structure(1, method="AUC: linear"))
   for (t in names(tests)) {
     # Note: using this structure ensures that there will not be
     # excessive warnings during testing.
@@ -148,13 +148,13 @@ test_that("pk.calc.auc: Confirm that center BLQ points are dropped, kept, or imp
   # Do this with both "linear" and "lin up/log down"
   tests <- list(
     "linear"=list(
-      AUCinf=1+1+0.5+1.5+1.5+1,
-      AUClast=1+1+0.5+1.5+1.5,
-      AUCall=1+1+0.5+1.5+1.5),
+      AUCinf=structure(1+1+0.5+1.5+1.5+1, method="AUC: linear"),
+      AUClast=structure(1+1+0.5+1.5+1.5, method="AUC: linear"),
+      AUCall=structure(1+1+0.5+1.5+1.5, method="AUC: linear")),
     "lin up/log down"=list(
-      AUCinf=1+1+0.5+1.5+1/log(2)+1,
-      AUClast=1+1+0.5+1.5+1/log(2),
-      AUCall=1+1+0.5+1.5+1/log(2)))
+      AUCinf=structure(1+1+0.5+1.5+1/log(2)+1, method="AUC: lin up/log down"),
+      AUClast=structure(1+1+0.5+1.5+1/log(2), method="AUC: lin up/log down"),
+      AUCall=structure(1+1+0.5+1.5+1/log(2), method="AUC: lin up/log down")))
   for (t in names(tests)) {
     for (n in names(tests[[t]])) {
       # Note: using this structure ensures that there will not be
@@ -180,13 +180,13 @@ test_that("pk.calc.auc: Confirm BLQ in the middle or end are calculated correctl
   # AUCall looks different when there are BLQs at the end
   tests <- list(
     "linear"=list(
-      AUCinf=1+1+0.5+1.5+1.5+1,
-      AUClast=1+1+0.5+1.5+1.5,
-      AUCall=1+1+0.5+1.5+1.5+0.5),
+      AUCinf=structure(1+1+0.5+1.5+1.5+1, method="AUC: linear"),
+      AUClast=structure(1+1+0.5+1.5+1.5, method="AUC: linear"),
+      AUCall=structure(1+1+0.5+1.5+1.5+0.5, method="AUC: linear")),
     "lin up/log down"=list(
-      AUCinf=1+1+0.5+1.5+1/log(2)+1,
-      AUClast=1+1+0.5+1.5+1/log(2),
-      AUCall=1+1+0.5+1.5+1/log(2)+0.5))
+      AUCinf=structure(1+1+0.5+1.5+1/log(2)+1, method="AUC: lin up/log down"),
+      AUClast=structure(1+1+0.5+1.5+1/log(2), method="AUC: lin up/log down"),
+      AUCall=structure(1+1+0.5+1.5+1/log(2)+0.5, method="AUC: lin up/log down")))
   for (t in names(tests)) {
     for (n in names(tests[[t]])) {
       # Note: using this structure ensures that there will not be
@@ -212,13 +212,13 @@ test_that("pk.calc.auc: Confirm BLQ in the middle or end are calculated correctl
   # starting times differing, so not tested here.)
   tests <- list(
     "linear"=list(
-      AUCinf=1+3+1.5+1.5+1,
-      AUClast=1+3+1.5+1.5,
-      AUCall=1+3+1.5+1.5+0.5),
+      AUCinf=structure(1+3+1.5+1.5+1, method="AUC: linear"),
+      AUClast=structure(1+3+1.5+1.5, method="AUC: linear"),
+      AUCall=structure(1+3+1.5+1.5+0.5, method="AUC: linear")),
     "lin up/log down"=list(
-      AUCinf=1+2/log(2)+1.5+1/log(2)+1,
-      AUClast=1+2/log(2)+1.5+1/log(2),
-      AUCall=1+2/log(2)+1.5+1/log(2)+0.5))
+      AUCinf=structure(1+2/log(2)+1.5+1/log(2)+1, method="AUC: lin up/log down"),
+      AUClast=structure(1+2/log(2)+1.5+1/log(2), method="AUC: lin up/log down"),
+      AUCall=structure(1+2/log(2)+1.5+1/log(2)+0.5, method="AUC: lin up/log down")))
   for (t in names(tests)) {
     for (n in names(tests[[t]])) {
       # Note: using this structure ensures that there will not be
@@ -248,13 +248,13 @@ test_that("pk.calc.auc: Confirm BLQ in the middle or end are calculated correctl
 test_that("pk.calc.auc: When AUCinf is requested with NA for lambda.z, the result is NA", {
   tests <- list(
     "linear"=list(
-      AUCinf=as.numeric(NA),
-      AUClast=1+3+1.5+1.5,
-      AUCall=1+3+1.5+1.5+0.5),
+      AUCinf=structure(as.numeric(NA), method="AUC: linear"),
+      AUClast=structure(1+3+1.5+1.5, method="AUC: linear"),
+      AUCall=structure(1+3+1.5+1.5+0.5, method="AUC: linear")),
     "lin up/log down"=list(
-      AUCinf=as.numeric(NA),
-      AUClast=1+2/log(2)+1.5+1/log(2),
-      AUCall=1+2/log(2)+1.5+1/log(2)+0.5))
+      AUCinf=structure(as.numeric(NA), method="AUC: lin up/log down"),
+      AUClast=structure(1+2/log(2)+1.5+1/log(2), method="AUC: lin up/log down"),
+      AUCall=structure(1+2/log(2)+1.5+1/log(2)+0.5, method="AUC: lin up/log down")))
   for (t in names(tests)) {
     for (n in names(tests[[t]])) {
       # Note: using this structure ensures that there will not be
@@ -288,13 +288,13 @@ test_that("pk.calc.auc: Test NA at the end", {
   # Test NA at the end
   tests <- list(
     "linear"=list(
-      AUCinf=as.numeric(NA),
-      AUClast=1+3+1.5+1.5,
-      AUCall=1+3+1.5+1.5+1),
+      AUCinf=structure(as.numeric(NA), method="AUC: linear"),
+      AUClast=structure(1+3+1.5+1.5, method="AUC: linear"),
+      AUCall=structure(1+3+1.5+1.5+1, method="AUC: linear")),
     "lin up/log down"=list(
-      AUCinf=as.numeric(NA),
-      AUClast=1+2/log(2)+1.5+1/log(2),
-      AUCall=1+2/log(2)+1.5+1/log(2)+1))
+      AUCinf=structure(as.numeric(NA), method="AUC: lin up/log down"),
+      AUClast=structure(1+2/log(2)+1.5+1/log(2), method="AUC: lin up/log down"),
+      AUCall=structure(1+2/log(2)+1.5+1/log(2)+1, method="AUC: lin up/log down")))
   for (t in names(tests))
     for (n in names(tests[[t]])) {
       # Note: using this structure ensures that there will not be
@@ -323,13 +323,13 @@ test_that("pk.calc.auc: Test NA at the end", {
 test_that("pk.calc.auc: interpolation of times within the time interval", {
   tests <- list(
     "linear"=list(
-      AUCinf=1+3+1.5+1.5+1,
-      AUClast=1+3+1.5+1.5,
-      AUCall=1+3+1.5+1.5+0.75),
+      AUCinf=structure(1+3+1.5+1.5+1, method="AUC: linear"),
+      AUClast=structure(1+3+1.5+1.5, method="AUC: linear"),
+      AUCall=structure(1+3+1.5+1.5+0.75, method="AUC: linear")),
     "lin up/log down"=list(
-      AUCinf=1+2/log(2)+1.5+1/log(2)+1,
-      AUClast=1+2/log(2)+1.5+1/log(2),
-      AUCall=1+2/log(2)+1.5+1/log(2)+0.5/log(2)))
+      AUCinf=structure(1+2/log(2)+1.5+1/log(2)+1, method="AUC: lin up/log down"),
+      AUClast=structure(1+2/log(2)+1.5+1/log(2), method="AUC: lin up/log down"),
+      AUCall=structure(1+2/log(2)+1.5+1/log(2)+0.5/log(2), method="AUC: lin up/log down")))
   for (t in names(tests)) {
     for (n in names(tests[[t]])) {
       # Note: using this structure ensures that there will not be
@@ -394,11 +394,9 @@ test_that("pk.calc.auc: warning with beginning of interval before the beginning
                           method=t),
         class="pknca_warn_auc_before_first"
       )
-      expect_equal(
-        v1,
-        tests[[t]][[n]],
-        info=paste(t, n)
-      )
+      expect_equal(v1,
+                   tests[[t]][[n]],
+                   info=paste(t, n))
     }
   }
 })
@@ -544,7 +542,7 @@ test_that("pk.calc.auc.inf.obs returns NA when lambda.z is NA", {
       interval = c(0, Inf),
       method = "linear"
     ),
-    NA_real_
+    structure(NA_real_, method = "AUC: linear")
   )
 })
 
@@ -577,14 +575,14 @@ test_that("pk.calc.aumc", {
       time=0:3,
       interval=c(0, 3),
       method="linear"),
-    3.75)
+    structure(3.75, method="AUC: linear"))
   expect_equal(
     pk.calc.aumc(
       conc=c(0, 1, 1, 0.5),
       time=0:3,
       interval=c(0, 3),
       method="lin up/log down"),
-    2-0.5/log(0.5)+0.5/(log(0.5)^2))
+    structure(2-0.5/log(0.5)+0.5/(log(0.5)^2), method="AUC: lin up/log down"))
   expect_equal(
     pk.calc.aumc(
       conc=c(0, 1, 1, 0.5),
@@ -593,7 +591,7 @@ test_that("pk.calc.aumc", {
       auc.type="AUCinf",
       lambda.z=1,
       method="lin up/log down"),
-    2 - 0.5/log(0.5) + 0.5/(log(0.5)^2) + 1.5 + 0.5)
+    structure(2 - 0.5/log(0.5) + 0.5/(log(0.5)^2) + 1.5 + 0.5, method="AUC: lin up/log down"))
 })
 
 
@@ -742,3 +740,34 @@ test_that("AUC with a single concentration measured should return NA (fix #176)"
     )
   )
 })
+
+test_that("pk.calc.auc and wrappers: method attribute is set and propagated", {
+
+  auc_params <- c(
+    "auc", "auc.last", "auc.inf.obs", "auc.inf.pred", "auc.all",
+    "aumc.last", "aumc.inf.obs", "aumc.inf.pred", "aumc.all"
+  )
+  auc_methods <- c("linear", "lin up/log down", "lin-log")
+  auc_args <- list(
+    conc=c(0,1,1),
+    time=0:2,
+    interval=c(0,2),
+    lambda.z=1,
+    clast.pred = 1,
+    clast.obs = 1
+  )
+
+  for (param in auc_params) {
+    auc_fun <- get(paste0("pk.calc.", param))
+    args_fun <- auc_args[intersect(names(auc_args), names(formals(auc_fun)))]
+    for (method in auc_methods) {
+      args_fun$method <- method
+      v <- do.call(auc_fun, args_fun)
+      expect_equal(
+        attr(v, "method"),
+        paste0("AUC: ", method),
+        info=paste("pk.calc.param sets method attribute for", param, "with method", method)
+      )
+    }
+  }
+})
diff --git a/tests/testthat/test-aucint.R b/tests/testthat/test-aucint.R
index 545118a0..4d70d293 100644
--- a/tests/testthat/test-aucint.R
+++ b/tests/testthat/test-aucint.R
@@ -469,11 +469,12 @@ test_that("aucint uses log extrapolation regardless of the interpolation method
   # the second is more directly mathematical.
   expect_equal(
     aucinf_obs6_lin - aucinf_obs5_lin,
-    aucinf_obs6_log - aucinf_obs5_log
+    aucinf_obs6_log - aucinf_obs5_log,
+    ignore_attr = TRUE
   )
   expect_equal(
     aucinf_obs6_lin,
-    aucinf_obs5_lin + (6-5)*(clast-ctau_extrap)/log(clast/ctau_extrap)
+    structure(aucinf_obs5_lin + (6-5)*(clast-ctau_extrap)/log(clast/ctau_extrap), method="AUC: linear")
   )
 })
 
@@ -498,6 +499,34 @@ test_that("aucint.inf.pred returns NA when half-life is not estimable (#450)", {
   expect_equal(aucint_inf_pred, NA_real_)
 })
 
+test_that("pk.calc.aucint and wrappers: method attribute is set and propagated", {
+  aucint_params <- c("aucint", "aucint.last", "aucint.inf.obs", "aucint.inf.pred", "aucint.all")
+  auc_methods <- c("linear", "lin up/log down", "lin-log")
+  auc_args <- list(
+    conc = c(0,1,1),
+    time = 0:2,
+    interval = c(0,2),
+    lambda.z = 1,
+    clast.pred = 1,
+    clast.obs = 1,
+    start = 0,
+    end = 2
+  )
+
+  for (param in aucint_params) {
+    auc_fun <- get(paste0("pk.calc.", param))
+    args_fun <- auc_args[intersect(names(auc_args), names(formals(auc_fun)))]
+    for (method in auc_methods) {
+      args_fun$method <- method
+      v <- do.call(auc_fun, args_fun)
+      expect_equal(
+        attr(v, "method"),
+        paste0("AUC: ", method),
+        info = paste("pk.calc.param sets method attribute for", param, "with method", method)
+      )
+    }
+  }
+})
 # ============================================================================
 # AUMC Tests (parallel to AUC tests)
 # ============================================================================
@@ -705,4 +734,4 @@ test_that("Integration functions are passed correctly through wrapper", {
   # Both should be positive numeric values
   expect_true(is.numeric(auc_result) && auc_result > 0)
   expect_true(is.numeric(aumc_result) && aumc_result > 0)
-})
\ No newline at end of file
+})
diff --git a/tests/testthat/test-auciv.R b/tests/testthat/test-auciv.R
index ee96b1df..cb0a854f 100644
--- a/tests/testthat/test-auciv.R
+++ b/tests/testthat/test-auciv.R
@@ -13,14 +13,16 @@ test_that("pk.calc.auciv works correctly", {
   expect_equal(
     # No check is done to confirm that the auc argument matches the data
     pk.calc.auciv(conc = 0:5, time = 0:5, c0 = 1, auc = 2.75),
-    2.75 + 1 - 0.5
+    2.75 + 1 - 0.5,
+    ignore_attr = TRUE
   )
   
   # With check = FALSE
   expect_equal(
     # No verifications are made on the data
-    pk.calc.auciv(conc = 0:5, time = 0:5, c0 = 1, auc = 2.75, check = FALSE),
-    2.75 + 1 - 0.5
+    pk.calc.auciv(conc = 0:5, time = 0:5, c0 = 1, auc = 2.75, check=FALSE),
+    2.75 + 1 - 0.5,
+    ignore_attr = TRUE
   )
   
   # With NA c0
@@ -111,6 +113,30 @@ test_that("missing dose information does not cause NA time (#353)", {
   expect_s3_class(o_nca, "PKNCAresults")
 })
 
+test_that("pk.calc.auciv: method attribute is set and propagated", {
+
+  auc_params <- c("auciv")
+  auc_methods <- c("linear", "lin up/log down", "lin-log")
+  auc_args <- list(
+    conc=3:1,
+    time=0:2,
+    c0 = 1,
+    auc = 2
+  )
+
+  for (param in auc_params) {
+    auc_fun <- get(paste0("pk.calc.", param))
+    args_fun <- auc_args[intersect(names(auc_args), names(formals(auc_fun)))]
+    for (method in auc_methods) {
+      args_fun$method <- method
+      v <- do.call(auc_fun, args_fun)
+      expect_equal(
+        attr(v, "method"),
+        paste0("AUC: ", method),
+        info=paste("pk.calc.param sets method attribute for", param, "with method", method)
+      )
+    }
+  }
 # ============================================================================
 # Wrapper Function Tests (pk.calc.auxciv)
 # ============================================================================
diff --git a/tests/testthat/test-class-PKNCAresults.R b/tests/testthat/test-class-PKNCAresults.R
index 09faa6b1..f214ac85 100644
--- a/tests/testthat/test-class-PKNCAresults.R
+++ b/tests/testthat/test-class-PKNCAresults.R
@@ -55,6 +55,16 @@ test_that("PKNCAresults generation", {
                 24.00, 0.3148, 0.05689, 0.9000,
                 0.8944, -0.9487, 5.000, 24.00, 20.00, 0.3011,
                 12.18, 1.560, 19.56),
+      PPANMETH=c(
+        "AUC: lin up/log down",
+        rep("", 4),
+        rep("", 10),
+        "AUC: lin up/log down",
+        "AUC: lin up/log down",
+        rep("", 4),
+        rep("", 10),
+        "AUC: lin up/log down"
+      ),
       exclude=NA_character_
     )
   expect_equal(
@@ -79,7 +89,7 @@ test_that("PKNCAresults generation", {
   )
   expect_equal(
     as.data.frame(o_result, out_format="wide"),
-    tidyr::spread(verify.result, key="PPTESTCD", value="PPORRES"),
+    tidyr::spread(verify.result[names(verify.result) != "PPANMETH"], key="PPTESTCD", value="PPORRES"),
     tolerance=0.001,
     info="Conversion of PKNCAresults to a data.frame in wide format (specifying wide format)"
   )
diff --git a/tests/testthat/test-half.life.R b/tests/testthat/test-half.life.R
index da4bfa10..558d68fd 100644
--- a/tests/testthat/test-half.life.R
+++ b/tests/testthat/test-half.life.R
@@ -122,15 +122,15 @@ test_that("pk.calc.half.life", {
 })
 
 test_that("half-life manual point selection", {
-  expect_equal(
-      pk.calc.half.life(conc=c(3, 1, 0.5, 0.13, 0.12, 0.113),
-                        time=c(0, 1, 2, 3, 4, 5),
-                        manually.selected.points=TRUE,
-                        min.hl.points=3,
-                        allow.tmax.in.half.life=FALSE,
-                        check=FALSE)$half.life,
-      1.00653,
-      tolerance=0.0001,
+  expect_equal(
+      pk.calc.half.life(conc=c(3, 1, 0.5, 0.13, 0.12, 0.113),
+                        time=c(0, 1, 2, 3, 4, 5),
+                        manually.selected.points=TRUE,
+                        min.hl.points=3,
+                        allow.tmax.in.half.life=FALSE,
+                        check=FALSE)$half.life,
+      1.00653,
+      tolerance=0.0001,
       info="manual selection uses the given points as is")
   expect_true(
     pk.calc.half.life(conc=c(3, 1, 0.5, 0.13, 0.12, 0.113),
@@ -194,6 +194,7 @@ test_that("half-life manual point selection", {
       tlast = 3L
     )
   attr(excluded_result, "exclude") <- "Negative half-life estimated with manually-selected points"
+  attr(excluded_result, "method") <- "Lambda Z: Manual selection"
   expect_equal(
     pk.calc.half.life(conc = 2^(1:3), time = 1:3, manually.selected.points = TRUE),
     excluded_result
diff --git a/tests/testthat/test-pk.calc.all.R b/tests/testthat/test-pk.calc.all.R
index 7a4adb34..dd72dd8f 100644
--- a/tests/testthat/test-pk.calc.all.R
+++ b/tests/testthat/test-pk.calc.all.R
@@ -63,6 +63,16 @@ test_that("pk.nca", {
                 24.00, 0.3148, 0.05689, 0.9000, 0.8944, -0.952,
                 5.000, 24.00, 20.00, 0.3011, 12.18,
                 1.560, 19.56),
+      PPANMETH = c(
+        "AUC: lin up/log down",
+        rep("", 4),
+        rep("", 10),
+        "AUC: lin up/log down",
+        "AUC: lin up/log down",
+        rep("", 4),
+        rep("", 10),
+        "AUC: lin up/log down"
+      ),
       exclude=NA_character_
     )
   expect_equal(
@@ -244,6 +254,12 @@ test_that("Calculations when no dose info is given", {
       PPTESTCD=rep(c("auclast", "cmax", "cl.last"), 2),
       PPORRES=c(13.5417297156528, 0.999812606062292, NA,
                 14.0305397438242, 0.94097296083447, NA),
+      PPANMETH = c(
+        "AUC: lin up/log down",
+        rep("", 2),
+        "AUC: lin up/log down",
+        rep("", 2)
+      ),
       exclude=NA_character_
     )
   )
@@ -757,6 +773,94 @@ test_that("do not give rbind error when interval columns have attributes (#381)"
   )
 })
 
+test_that("pk.nca produces the PPANMETH column", {
+  # --- Setup shared concentration and dose data ---
+  tmpconc <- generate.conc(1, 1, 0:24)
+  tmpdose <- generate.dose(tmpconc)
+  myconc <- PKNCAconc(tmpconc, formula=conc~time|treatment+ID)
+  mydose <- PKNCAdose(tmpdose, formula=dose~time|treatment+ID)
+
+  # --- PPANMETH differentiates based on the AUC method used ---
+  mydata_linear <- PKNCAdata(myconc, mydose, intervals=data.frame(start=0, end=24, auclast=TRUE), options=list(auc.method="linear"))
+  mydata_linlog <- PKNCAdata(myconc, mydose, intervals=data.frame(start=0, end=24, auclast=TRUE), options=list(auc.method="lin up/log down"))
+  res_linear <- pk.nca(mydata_linear)
+  res_linlog <- pk.nca(mydata_linlog)
+  expect_true("PPANMETH" %in% names(res_linear$result))
+  expect_true("PPANMETH" %in% names(res_linlog$result))
+  expect_true(any(grepl("AUC: linear", res_linear$result$PPANMETH, fixed=TRUE)))
+  expect_true(any(grepl("AUC: lin up/log down", res_linlog$result$PPANMETH, fixed=TRUE)))
+
+  # --- PPANMETH distinguishes how the half.life was adjusted ---
+  tmpconc$include_hl <- tmpconc$time <= 22
+  tmpconc$exclude_hl <- tmpconc$time == 22
+  myconc_base <- PKNCAconc(tmpconc, formula=conc~time|treatment+ID)
+  myconc_incl <- PKNCAconc(tmpconc, formula=conc~time|treatment+ID, include_half.life="include_hl")
+  myconc_excl <- PKNCAconc(tmpconc, formula=conc~time|treatment+ID, exclude_half.life="exclude_hl")
+  mydata_base <- PKNCAdata(myconc_base, mydose, intervals=data.frame(start=0, end=24, lambda.z=TRUE))
+  mydata_incl <- PKNCAdata(myconc_incl, mydose, intervals=data.frame(start=0, end=24, lambda.z=TRUE))
+  mydata_excl <- PKNCAdata(myconc_excl, mydose, intervals=data.frame(start=0, end=24, lambda.z=TRUE))
+  res_base <- pk.nca(mydata_base)
+  res_incl <- pk.nca(mydata_incl)
+  res_excl <- pk.nca(mydata_excl)
+  expect_true("PPANMETH" %in% names(res_base$result))
+  expect_true("PPANMETH" %in% names(res_incl$result))
+  expect_true("PPANMETH" %in% names(res_excl$result))
+  expect_equal(
+    unique(res_base$result$PPANMETH[res_base$result$PPTESTCD %in% c("lambda.z", "half.life", "r.squared")]),
+    ""
+  )
+  expect_equal(
+    unique(res_incl$result$PPANMETH[res_incl$result$PPTESTCD %in% c("lambda.z", "half.life", "r.squared")]),
+    "Lambda Z: Manual selection"
+  )
+  expect_equal(
+    unique(res_excl$result$PPANMETH[res_excl$result$PPTESTCD %in% c("lambda.z", "half.life", "r.squared")]),
+    ""
+  )
+  expect_equal(
+    unique(res_base$result$PPANMETH[res_base$result$PPTESTCD %in% c("tmax", "cmax")]),
+    ""
+  )
+  expect_equal(
+    unique(res_incl$result$PPANMETH[res_incl$result$PPTESTCD %in% c("tmax", "cmax")]),
+    ""
+  )
+  expect_equal(
+    unique(res_excl$result$PPANMETH[res_excl$result$PPTESTCD %in% c("tmax", "cmax")]),
+    ""
+  )
+
+  # --- PPANMETH specifies if an imputation method was used in the interval ---
+  o_data <- PKNCAdata(myconc, mydose, intervals=data.frame(start=0, end=24, c0=TRUE))
+  o_data_impute <- PKNCAdata(myconc, mydose, intervals=data.frame(start=0, end=24, c0=TRUE), impute="start_conc0")
+  res <- pk.nca(o_data)
+  res_impute <- pk.nca(o_data_impute)
+  expect_equal(res$result$PPANMETH, "")
+  expect_true("PPANMETH" %in% names(res$result))
+  expect_equal(res$result$PPANMETH, "")
+  expect_equal(res_impute$result$PPANMETH, "Imputation: start_conc0")
+
+  # --- PPANMETH reports based on the parameter dependencies ---
+  mydata <- PKNCAdata(
+    myconc_incl, mydose,
+    intervals=data.frame(start=0, end=24, c0 = TRUE, half.life = TRUE, aucinf.pred=TRUE),
+    impute = "start_conc0"
+  )
+  res <- pk.nca(mydata)
+  expect_equal(
+    res$result$PPANMETH[res$result$PPTESTCD == "c0"],
+    "Imputation: start_conc0"
+  )
+  expect_equal(
+    res$result$PPANMETH[res$result$PPTESTCD == "half.life"],
+    "Imputation: start_conc0. Lambda Z: Manual selection"
+  )
+  expect_equal(
+    res$result$PPANMETH[res$result$PPTESTCD == "aucinf.pred"],
+    "Imputation: start_conc0. AUC: lin up/log down"
+  )
+})
+
 test_that("pk.nca can be run for each parameter independently (#473)", {
   
   # ── Dense data setup ──────────────────────────────────────────────────────
diff --git a/tests/testthat/test-pk.calc.simple.R b/tests/testthat/test-pk.calc.simple.R
index a4ae05cc..7074910f 100644
--- a/tests/testthat/test-pk.calc.simple.R
+++ b/tests/testthat/test-pk.calc.simple.R
@@ -489,6 +489,7 @@ test_that("pk.calc.aucabove", {
           pk.calc.aucabove(conc = d_conc$conc, time = d_conc$time, conc_above = 2),
           pk.calc.aucabove(conc = d_conc$conc, time = d_conc$time, conc_above = 3)
         ),
+      PPANMETH = c("", "", "AUC: lin up/log down", "AUC: lin up/log down"),
       exclude = NA_character_
     )
   )
diff --git a/tests/testthat/test-sparse.R b/tests/testthat/test-sparse.R
index efe3bfd3..6eaa2c5e 100644
--- a/tests/testthat/test-sparse.R
+++ b/tests/testthat/test-sparse.R
@@ -23,14 +23,14 @@ test_that("sparse_auc", {
     sparse_batch <- pk.calc.sparse_auc(conc = d_sparse$conc, time = d_sparse$time, subject = d_sparse$id),
     regexp = "Cannot yet calculate sparse degrees of freedom for multiple samples per subject"
   )
-  expect_equal(sparse_batch$sparse_auc, auclast)
-  expect_equal(sparse_batch$sparse_auc_se, auclast_se_batch)
-  expect_equal(sparse_batch$sparse_auc_df, NA_real_)
-  
-  sparse_serial <- pk.calc.sparse_auc(conc = d_sparse$conc, time = d_sparse$time, subject = seq_len(nrow(d_sparse)))
-  expect_equal(sparse_serial$sparse_auc, auclast)
+  expect_equal(sparse_batch$sparse_auc, structure(auclast, method=c("AUC: linear", "Sparse: arithmetic mean, <=50% BLQ")))
+  expect_equal(sparse_batch$sparse_auc_se, structure(auclast_se_batch, method=c("AUC: linear", "Sparse: arithmetic mean, <=50% BLQ")))
+  expect_equal(sparse_batch$sparse_auc_df, structure(NA_real_, method=c("AUC: linear", "Sparse: arithmetic mean, <=50% BLQ")))
+
+  sparse_serial <- pk.calc.sparse_auc(conc=d_sparse$conc, time=d_sparse$time, subject=seq_len(nrow(d_sparse)))
+  expect_equal(sparse_serial$sparse_auc, structure(auclast, method=c("AUC: linear", "Sparse: arithmetic mean, <=50% BLQ")))
   expect_equal(as.numeric(sparse_serial$sparse_auc_se), auclast_se_serial)
-  expect_equal(sparse_serial$sparse_auc_df, auclast_df_serial)
+  expect_equal(sparse_serial$sparse_auc_df, structure(auclast_df_serial, method=c("AUC: linear", "Sparse: arithmetic mean, <=50% BLQ")))
 })
 
 test_that("sparse_auclast expected errors", {
@@ -65,6 +65,22 @@ test_that("sparse_mean", {
   )
 })
 
+test_that("sparse_auc and sparse_auclast method attribute", {
+  d_sparse <-
+    data.frame(
+      id = c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 4L, 5L, 6L, 4L, 5L, 6L, 7L, 8L, 9L, 7L, 8L, 9L),
+      conc = c(0, 0, 0,  1.75, 2.2, 1.58, 4.63, 2.99, 1.52, 3.03, 1.98, 2.22, 3.34, 1.3, 1.22, 3.54, 2.84, 2.55, 0.3, 0.0421, 0.231),
+      time = c(0, 0, 0, 1, 1, 1, 6, 6, 6, 2, 2, 2, 10, 10, 10, 4, 4, 4, 24, 24, 24),
+      dose = c(100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100)
+    )
+  auc <- pk.calc.sparse_auc(conc=d_sparse$conc, time=d_sparse$time, subject=seq_len(nrow(d_sparse)))
+  expect_equal(attr(auc$sparse_auc, "method"),
+               c("AUC: linear", "Sparse: arithmetic mean, <=50% BLQ"))
+
+  auclast <- pk.calc.sparse_auclast(conc=d_sparse$conc, time=d_sparse$time, subject=seq_len(nrow(d_sparse)))
+  expect_equal(attr(auclast$sparse_auclast, "method"),
+               c("AUC: linear", "Sparse: arithmetic mean, <=50% BLQ"))
+})
 test_that("cov_holder clips covariance to Cauchy-Schwartz bound", {
   # Construct data where the Holder covariance formula exceeds sqrt(var1*var2).
   # Time 1: subjects 1 & 2, concentrations 0 & 10 → var = 50