diff --git a/NEWS.md b/NEWS.md
index c971fc7..8e36e2b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,7 +1,8 @@
 # (development version)
 
+* Help pages were revised for clarity, consistency, and improved user guidance [#21](https://github.com/PPaccioretti/paar/issues/21)
 * Now Summary index is calculated correctly. 
-* PCA results can be obtained within kmspc function
+* PCA results can be obtained within kmspc function.
 * Factor recoding.
 
 # paar 1.0.1
diff --git a/R/compare_zone.R b/R/compare_zone.R
index 0ee8245..f124189 100644
--- a/R/compare_zone.R
+++ b/R/compare_zone.R
@@ -1,23 +1,67 @@
-#' Compare spatial zone means
+#' Compare means between spatial zones
 #'
-#' @param data \code{sf} object with zones
-#' @param variable \code{character} or \code{sf} object to use for mean comparison
-#' @param zonesCol \code{character} colname from data were zone are specified
-#' @param alpha \code{numeric} Significance level to use for comparison
-#' @param returnLSD \code{logical} when LSD calculates with spatial variance should be returned
-#' @param join function to use for st_join if variable is \code{sf} object
-#' @param grid_dim \code{numeric} grid dimentins to estimate spatial variance
+#' @description
+#' Compares variable means across spatial zones using a spatially-adjusted
+#' least significant difference (LSD) approach based on kriging variance.
+#'
+#' The function accounts for spatial variability by estimating semivariograms
+#' and deriving a spatial variance component, which is then used to assess
+#' differences between zone means.
+#'
+#' @param data an \code{sf} object containing the spatial zones
+#'
+#' @param variable either:
+#' \itemize{
+#'   \item a \code{character} vector with column names in \code{data}, or
+#'   \item an \code{sf} object with external variables to be compared. In this
+#'   case, values are spatially joined to \code{data}.
+#' }
+#'
+#' @param zonesCol \code{character}. Column name in \code{data} defining zones
+#'
+#' @param alpha \code{numeric}. Significance level for mean comparison
+#'
+#' @param join function used in \code{sf::st_join} when \code{variable} is an
+#'   external \code{sf} object (default: \code{sf::st_nearest_feature})
+#'
+#' @param returnLSD \code{logical}. If \code{TRUE}, returns the LSD value used
+#'   for comparisons
+#'
+#' @param grid_dim \code{numeric}. Grid resolution used to estimate spatial
+#'   variance when interpolating external variables. If missing, it is
+#'   automatically determined.
+#'
+#' @details
+#' When \code{variable} is an external \code{sf} object, values are interpolated
+#' using ordinary kriging before comparison. Otherwise, cross-validation of the
+#' variogram model is used to estimate spatial variance.
+#'
+#' Pairwise comparisons between zones are evaluated using a spatially-adjusted
+#' LSD criterion:
+#'
+#' \deqn{LSD = z_{1-\alpha/2} \times \sigma_{spatial}}
+#'
+#' where \eqn{\sigma_{spatial}} is derived from kriging variance.
+#'
+#' Results are presented using compact letter displays to indicate groups
+#' of zones that are not significantly different.
+#'
+#' @return
+#' A list with:
+#' \describe{
+#'   \item{differences}{list of data frames with mean comparisons per variable}
+#'   \item{descriptive_stat}{data frame with descriptive statistics and spatial variance}
+#' }
 #'
-#' @return \code{list} with differences and descriptive_stat
 #' @references
 #' Paccioretti, P., Córdoba, M., & Balzarini, M. (2020).
 #' FastMapping: Software to create field maps and identify management zones
-#' in precision agriculture. Computers and Electronics in Agriculture, 175,
-#' 105556 https://doi.org/10.1016/j.compag.2020.105556.
-#' @export
+#' in precision agriculture. \emph{Computers and Electronics in Agriculture},
+#' 175, 105556. \doi{10.1016/j.compag.2020.105556}
 #'
 #' @example inst/examples/compare_zone.R
 #'
+#' @export
 
 compare_zone <- function(
   data,
@@ -37,7 +81,8 @@ compare_zone <- function(
       zonesCol,
       descriptive,
       returnLSD = returnLSD,
-      join = join
+      join = join,
+      alpha = alpha
     )
 
   list(
@@ -56,7 +101,8 @@ validVarKrig <-
     numCluster,
     descStat,
     join = sf::st_nearest_feature,
-    returnLSD
+    returnLSD,
+    alpha = 0.05
   ) {
     MedCV <-
       paste0(round(descStat$MediaVar, 2), " (", round(descStat$CVVar, 1), ")")
@@ -84,7 +130,8 @@ validVarKrig <-
             numCluster = numCluster,
             descStat = descStat,
             join = join,
-            retDMS = returnLSD
+            retDMS = returnLSD,
+            alpha = alpha
           )
         })
       names(MeansComparissons) <- colnames(sf::st_drop_geometry(dataToCompare))
@@ -96,7 +143,8 @@ validVarKrig <-
           Clasif = clusterData,
           numCluster = numCluster,
           descStat = descStat,
-          retDMS = returnLSD
+          retDMS = returnLSD,
+          alpha = alpha
         )
 
       names(MeansComparissons) <- dataToCompare
diff --git a/R/depurate.R b/R/depurate.R
index 3504260..3620273 100644
--- a/R/depurate.R
+++ b/R/depurate.R
@@ -1,61 +1,95 @@
-#' Remove errors from spatial data
+#' Spatial data depuration (error removal)
 #'
-#' @description Data can be filtered by null, edge values, global
-#'  outliers and spatial outliers or local defective observations. Default
-#'  values are optimized for precision agricultural data.
+#' @description
+#' Filters spatial point data by removing erroneous observations based on
+#' geometric, statistical, and spatial criteria. The function implements a
+#' sequential depuration workflow commonly used in precision agriculture.
 #'
-#' @param x an \code{sf} points object
-#' @param y \code{character} with the name of the variable to use for
-#'   depuration/filtering process
-#' @param toremove \code{character} vector specifying the procedure to
-#'   implement for errors removal. Default 'edges', 'outlier', 'inlier'.
-#'   See Details.
-#' @param crs coordinate reference system: integer with the EPSG code,
-#'   or character with proj4string to convert coordinates if \code{x} has
-#'   longitude/latitude data
-#' @param buffer \code{numeric} distance in meters to be removed. Negative
-#'   values are recommended
-#' @param ylimitmax \code{numeric} of length 1 indicating the maximum limit
-#'   for the \code{y} variable. If \code{NA} \code{Inf} is assumed
-#' @param ylimitmin \code{numeric} of length 1 indicating the minimum limit
-#'   for the \code{y} variable. If \code{NA} \code{-Inf} is assumed
-#' @param sdout \code{numeric} values outside the interval
-#'  \eqn{mean ± sdout × sdout} values will be removed
-#' @param ldist \code{numeric} lower distance bound to identify neighbors
-#' @param udist \code{numeric} upper distance bound to identify neighbors
-#' @param criteria \code{character} with "LM" and/or "MP" for methods to
-#'    identify spatial outliers
-#' @param zero.policy default NULL, use global option value;
-#'   if FALSE stop with error for any empty neighbors sets,
-#'   if TRUE permit the weights list to be formed with zero-length
-#'   weights vectors
-#' @param poly_border \code{sf} object with one polygon or NULL. Can be
-#' the result of \code{concaveman::concaveman}
+#' @param x An \code{sf} object with POINT geometries.
+#' @param y A \code{character} string indicating the variable name used
+#'   for filtering. If missing and only one attribute column is present,
+#'   it is used by default.
+#' @param toremove A \code{character} vector specifying which procedures
+#'   to apply. Options are \code{"edges"}, \code{"outlier"}, and
+#'   \code{"inlier"}. The order of execution is fixed and cannot be modified.
+#' @param crs Coordinate reference system used when transforming
+#'   longitude/latitude data. Can be an EPSG code or proj4string.
+#' @param buffer A \code{numeric} value indicating the distance (in meters)
+#'   for edge removal. Negative values are recommended to shrink boundaries.
+#' @param ylimitmax Numeric upper bound for \code{y}. If \code{NA}, \code{Inf}
+#'   is used.
+#' @param ylimitmin Numeric lower bound for \code{y}. If \code{NA}, \code{-Inf}
+#'   is used.
+#' @param sdout Numeric multiplier for standard deviation used to detect
+#'   global outliers.
+#' @param ldist Numeric lower distance bound for neighborhood definition.
+#' @param udist Numeric upper distance bound for neighborhood definition.
+#' @param criteria Character vector specifying spatial outlier detection
+#'   methods: \code{"LM"} (Local Moran) and/or \code{"MP"} (Moran Plot).
+#' @param zero.policy Logical. If \code{TRUE}, allows empty neighbor sets;
+#'   if \code{FALSE}, stops with an error.
+#' @param poly_border Optional \code{sf} polygon defining field boundaries.
+#'   If \code{NULL}, a hull is computed automatically.
 #'
 #' @details
-#' Possible values for \code{toremove} are one or more elements of:
+#' The depuration process is applied in a fixed sequence:
+#'
+#' \enumerate{
+#'   \item Edge removal (\code{"edges"})
+#'   \item Global outlier removal (\code{"outlier"})
+#'   \item Spatial outlier removal (\code{"inlier"})
+#' }
+#'
+#' The \code{toremove} argument controls which of these steps are applied,
+#' but **does not modify the order of execution**.
+#'
+#' Available procedures are:
+#'
 #' \describe{
-#'   \item{edges}{All data points for a distance of \code{buffer} m from data
-#'   edges are deleted.}
-#'   \item{outlier}{Values that are outside the mean±\code{sdout} are removed}
-#'   \item{inlier}{Local Moran index of spatial autocorrelation is calculated
-#'   for each datum as a tool to identify inliers}
+#'   \item{edges}{
+#'   Removes points located within a specified \code{buffer} distance from
+#'   the field boundary. The boundary is computed using a concave hull
+#'   (\code{concaveman}) or a convex hull if the package is not available.
+#'   }
+#'
+#'   \item{outlier}{
+#'   Removes global outliers based on:
+#'   \itemize{
+#'     \item user-defined limits (\code{ylimitmin}, \code{ylimitmax})
+#'     \item statistical thresholds defined as
+#'     \eqn{mean \pm sdout \times sd}
+#'   }
+#'   }
+#'
+#'   \item{inlier}{
+#'   Identifies and removes spatial outliers using:
+#'   \itemize{
+#'     \item Local Moran's I statistic ("LM")
+#'     \item Moran scatterplot influence ("MP")
+#'   }
+#'   }
 #' }
 #'
-#' @references
-#' Vega, A., Córdoba, M., Castro-Franco, M. et al. Protocol for
-#' automating error removal from yield maps. Precision Agric 20, 1030–1044
-#' (2019). https://doi.org/10.1007/s11119-018-09632-8
+#' Default parameter values are tuned for precision agriculture datasets
+#' (e.g., yield maps).
+#'
 #'
-#' @return an object of class \code{paar} with two elements:
+#' @return
+#' An object of class \code{paar} (list) with:
 #' \describe{
-#'  \item{depurated_data}{\code{sf} object with the data after the removal
-#'  process}
-#'  \item{condition}{\code{character} vector with the condition of each
-#'  observation}
-#'  }
-#' @export
+#'   \item{depurated_data}{Filtered \code{sf} object}
+#'   \item{condition}{Character vector indicating the reason each observation
+#'   was removed (or \code{NA} if retained)}
+#' }
+#'
+#' @references
+#' Vega, A., Córdoba, M., Castro-Franco, M. et al. (2019).
+#' Protocol for automating error removal from yield maps.
+#' \emph{Precision Agriculture}, 20, 1030–1044.
+#' \doi{10.1007/s11119-018-09632-8}
+#'
 #' @example inst/examples/depurate.R
+#' @export
 
 depurate <- function(
   x,
diff --git a/R/fuzzy_k_means.R b/R/fuzzy_k_means.R
index 908e156..40b0a3f 100644
--- a/R/fuzzy_k_means.R
+++ b/R/fuzzy_k_means.R
@@ -1,15 +1,42 @@
-#' Fuzzy k-means clustering
+#' Fuzzy k-means clustering (non-spatial)
+#'
+#' @description
+#' Performs fuzzy k-means clustering on tabular data (non-spatial).
+#' This function is a lightweight wrapper around \code{e1071::cmeans},
+#' providing a vectorized workflow and clustering quality indices.
+#'
+#' It is primarily intended as a fallback method when spatial clustering
+#' (e.g., \code{kmspc}) cannot be applied, such as when only one variable
+#' is available.
 #'
-#' @description Performs a vectorized fuzzy k-means clustering, this procedure
-#' it is not spatial. The function is almost a wrapper of the function cmeans
-#' from the package e1071. Is intended to be used when `KM-sPC` procedure is not
-#' possible because data set has only 1 variable.
 #'
 #' @inheritParams kmspc
-#' @return a list with classification results and indices to select best number of
-#' clusters.
-#' @export
+#' @details
+#' Missing values are removed prior to clustering. Observations with missing
+#' values are reintroduced in the output with \code{NA} cluster assignments.
+#'
+#' Clustering is performed for each value in \code{number_cluster}, and
+#' several indices are returned to assist in selecting the optimal number
+#' of clusters:
+#' \itemize{
+#'   \item Xie-Beni index
+#'   \item Partition coefficient
+#'   \item Partition entropy
+#'   \item Summary index
+#' }
+#'
+#' @return
+#' A list with:
+#' \describe{
+#'   \item{cluster}{\code{data.frame} with cluster assignments for each
+#'   evaluated number of clusters}
+#'   \item{indices}{\code{data.frame} with clustering validity indices}
+#'   \item{summaryResults}{\code{data.frame} with clustering metrics}
+#' }
+#'
+#' @seealso \code{\link{kmspc}}
 #' @example inst/examples/fuzzy_k_means.R
+#' @export
 #'
 fuzzy_k_means <- function(
   data,
@@ -22,10 +49,13 @@ fuzzy_k_means <- function(
     myNumVars <-
       unlist(lapply(sf::st_drop_geometry(data), is.numeric))
     if (sum(myNumVars) == 0) {
-      stop('Non numeric variables were found in data')
+      stop('No numeric variables found in data')
     }
-    warning("The numeric Variable will be used to make clusters", call. = FALSE)
+    warning("Numeric variables will be used for clustering", call. = FALSE)
     variables <- names(sf::st_drop_geometry(data))[myNumVars]
+    message(
+      paste("Using variables:", paste(variables, collapse = ", "))
+    )
   }
 
   # if (!inherits(data, "sf") & (length(variables) == 1)) {
diff --git a/R/kmspc.R b/R/kmspc.R
index a264a4d..7a34bcc 100644
--- a/R/kmspc.R
+++ b/R/kmspc.R
@@ -1,29 +1,75 @@
-#' MULTISPATI-PCA clustering
+#' Spatial PCA-based fuzzy clustering (MULTISPATI-PCA)
+#'
+#' @description
+#' Performs clustering of spatial data using a combination of spatial
+#' Principal Component Analysis (PCA), and fuzzy k-means clustering.
+#'
+#' The workflow consists of:
+#' \enumerate{
+#'   \item Dimensionality reduction using spatial PCA
+#'   \item Selection of components based on explained spatial variance
+#'   \item Fuzzy clustering over selected components
+#' }
 #'
 #'
 #' @inheritParams depurate
-#' @param data sf object
-#' @param variables variables to use for clustering, if missing, all numeric
-#' variables will be used
-#' @param number_cluster \code{numeric} vector with number of final clusters
-#' @param explainedVariance \code{numeric} number in percentage of explained variance
-#' from PCA analysis to keep and make cluster process
-#' @param fuzzyness A number greater than 1 giving the degree of fuzzification.
-#' @param center a logical or numeric value, centring option
-#' if TRUE, centring by the mean
-#' if FALSE no centring
-#' if a numeric vector, its length must be equal to the number of
-#' columns of the data frame df and gives the decentring
-#' @param distance \code{character} Must be one of the following:
-#' If "euclidean", the mean square error, if "manhattan", the mean
-#' absolute error is computed. Abbreviations are also accepted.
-#' @param only_spca_results \code{logical}; should return both PCA and sPCA
-#' results (\code{FALSE}), or only sPCA results (\code{TRUE})? This can be a
-#' time consuming process if there are multiple variables.
-#' @param all_results \code{logical}; should return the results from the
-#' sPCA and PCA call?
-#' @return a list with classification results and indices to select best number of
-#' clusters.
+#' @param data an \code{sf} object with point geometries
+#'
+#' @param variables \code{character} vector with variable names used for clustering.
+#'   If missing, all numeric variables in \code{data} are used.
+#'
+#' @param number_cluster \code{numeric} vector indicating the number of clusters
+#'   to evaluate (e.g., \code{3:5})
+#'
+#' @param explainedVariance \code{numeric}. Percentage (0–100) of cumulative
+#'   explained spatial variance used to select spatial principal components.
+#'   Values between 0 and 1 are interpreted as proportions.
+#'
+#' @param fuzzyness \code{numeric} value greater than 1 controlling the degree of
+#'   fuzziness in clustering (see \code{e1071::cmeans})
+#'
+#' @param center centering option passed to PCA:
+#'   \describe{
+#'     \item{TRUE}{center variables by their mean}
+#'     \item{FALSE}{no centering}
+#'     \item{numeric}{custom centering vector}
+#'   }
+#'
+#' @param distance \code{character} distance metric for clustering.
+#'   One of \code{"euclidean"} or \code{"manhattan"} (abbreviations allowed)
+#'
+#' @param only_spca_results \code{logical}. If \code{TRUE}, only spatial PCA results
+#'   are returned. If \code{FALSE}, both PCA and spatial PCA summaries are included.
+#'
+#' @param all_results \code{logical}. If \code{TRUE}, full PCA and spatial PCA
+#'   objects are returned (can increase computation time and memory use).
+#'
+#' @param ldist,udist \code{numeric}. Lower and upper distance thresholds used
+#'   to define spatial neighbors.
+#'
+#' @details
+#' Spatial relationships are defined using distance-based neighbors
+#' (\code{spdep::dnearneigh}). These relationships are incorporated into the
+#' spatial PCA analysis to extract spatially structured components.
+#'
+#' Clustering is performed using fuzzy c-means over selected spatial components.
+#' Several indices are computed to help determine the optimal number of clusters:
+#' \itemize{
+#'   \item Xie-Beni index
+#'   \item Partition coefficient
+#'   \item Partition entropy
+#'   \item Summary index (normalized combination)
+#' }
+#'
+#' @return
+#' A list with the following elements:
+#' \describe{
+#'   \item{cluster}{\code{data.frame} with cluster assignments for each evaluated number of clusters}
+#'   \item{indices}{\code{data.frame} with clustering validity indices}
+#'   \item{summaryResults}{\code{data.frame} with clustering metrics (iterations, SSDW)}
+#'   \item{pca_results}{(optional) PCA and/or spatial PCA summaries depending on arguments}
+#' }
+#'
 #' @export
 #' @example inst/examples/kmspc.R
 #'
diff --git a/R/spatial_t_test.R b/R/spatial_t_test.R
index 9e251f1..5fb90ed 100644
--- a/R/spatial_t_test.R
+++ b/R/spatial_t_test.R
@@ -1,15 +1,38 @@
 #' Modified t test
 #'
 #'
-#' @description Performs a modified version of the t test to assess the
-#' correlation between spatial processes. See SpatialPack::modified.ttest for details.
+#' @description
+#' Performs a modified t-test to assess the correlation between variables
+#' while accounting for spatial autocorrelation. This implementation wraps
+#' \code{SpatialPack::modified.ttest}.
 #'
-#' @param data \code{sf} data to extract coordinates or two
-#' columns \code{matrix} or \code{data.frame} specifying coordinates.
+#' @details
+#' The function computes pairwise correlations between the specified variables
+#' and adjusts the significance test to account for spatial dependence using
+#' coordinates. If \code{data} is an \code{sf} object, coordinates are extracted
+#' automatically. Otherwise, coordinates must be provided as an object with two
+#' columns.
+#'
+#' @param data An \code{sf} object containing geometry and variables, or a
+#'   \code{matrix}/\code{data.frame} with two columns representing
+#'   spatial coordinates (e.g., X and Y).
+#'
+#' @param variables A \code{character} vector with the names of the variables
+#'   to be tested. If \code{data} is not an \code{sf} object, this should be a
+#'   matrix or data.frame of variables to test.
+#'
+#' @return A \code{data.frame} with the following columns:
+#' \describe{
+#'   \item{Var1}{Name of the first variable}
+#'   \item{Var2}{Name of the second variable}
+#'   \item{corr}{Estimated correlation coefficient}
+#'   \item{p.value}{P-value adjusted for spatial autocorrelation}
+#' }
+#'
+#' @seealso \code{\link[SpatialPack]{modified.ttest}}
+#'
+#' @example inst/examples/spatial_t_test.R
 #'
-#' @param variables \code{character} vector with column names to perform ttest
-#' @return a data.frame with the correlation and p-value for each pair
-#'  of variables
 #' @export
 
 spatial_t_test <- function(data, variables) {
diff --git a/man/spatial_t_test.Rd b/man/spatial_t_test.Rd
index 7b9ed6d..507e839 100644
--- a/man/spatial_t_test.Rd
+++ b/man/spatial_t_test.Rd
@@ -7,16 +7,54 @@
 spatial_t_test(data, variables)
 }
 \arguments{
-\item{data}{\code{sf} data to extract coordinates or two
-columns \code{matrix} or \code{data.frame} specifying coordinates.}
+\item{data}{An \code{sf} object containing geometry and variables, or a
+\code{matrix}/\code{data.frame} with two columns representing
+spatial coordinates (e.g., X and Y).}
 
-\item{variables}{\code{character} vector with column names to perform ttest}
+\item{variables}{A \code{character} vector with the names of the variables
+to be tested. If \code{data} is not an \code{sf} object, this should be a
+matrix or data.frame of variables to test.}
 }
 \value{
-a data.frame with the correlation and p-value for each pair
- of variables
+A \code{data.frame} with the following columns:
+\describe{
+  \item{Var1}{Name of the first variable}
+  \item{Var2}{Name of the second variable}
+  \item{corr}{Estimated correlation coefficient}
+  \item{p.value}{P-value adjusted for spatial autocorrelation}
+}
 }
 \description{
-Performs a modified version of the t test to assess the
-correlation between spatial processes. See SpatialPack::modified.ttest for details.
+Performs a modified t-test to assess the correlation between variables
+while accounting for spatial autocorrelation. This implementation wraps
+\code{SpatialPack::modified.ttest}.
+}
+\details{
+The function computes pairwise correlations between the specified variables
+and adjusts the significance test to account for spatial dependence using
+coordinates. If \code{data} is an \code{sf} object, coordinates are extracted
+automatically. Otherwise, coordinates must be provided as an object with two
+columns.
+}
+\examples{
+if (requireNamespace("SpatialPack", quietly = TRUE)) {
+  library(sf)
+  data(wheat, package = 'paar')
+
+  # Transform the data.frame into a sf object
+  wheat_sf <- st_as_sf(wheat, coords = c('x', 'y'), crs = 32720)
+
+  # Run spatial t test
+  t_test_results <-
+    spatial_t_test(
+      wheat_sf,
+      variables = c('CE30', 'CE90')
+    )
+
+  # Print the t_test_results
+  t_test_results
+}
+}
+\seealso{
+\code{\link[SpatialPack]{modified.ttest}}
 }