diff --git a/NEWS.md b/NEWS.md
index dcd301d..9b5a920 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -11,6 +11,7 @@
 - format_data_PPBstats.data_network: take into accuont when there are not lat and long columns
 - pmap in common_functions : use stamen instead of google on map, cf #90
 - implement spatial analysis, cf #20 
+- mean comparisons with spatial model, cf #124 
 
 ### add new functions
 - HA_to_LF, cf #25
diff --git a/R/mean_comparisons.check_model_spatial.R b/R/mean_comparisons.check_model_spatial.R
index 85f1124..0a20eda 100644
--- a/R/mean_comparisons.check_model_spatial.R
+++ b/R/mean_comparisons.check_model_spatial.R
@@ -39,6 +39,8 @@
 #' 
 #' @import agricolae
 #' 
+#' @importFrom SpATS predict.SpATS
+#' 
 mean_comparisons.check_model_spatial <- function(
   x, 
   alpha = 0.05,
@@ -51,10 +53,16 @@ mean_comparisons.check_model_spatial <- function(
   variable = info$variable
   summary_model = x$spatial$model$summary
   
-  # 2. Mean comparison on germplasm ----------
+  # 2. Get prediction of BLUPs
+  pre = predict.SpATS(out_spatial$model$model, which = "germplasm")
+  pre$germplasm = factor(pre$germplasm, levels = pre$germplasm[order(pre$predicted.values)])
+  pre$lower = pre$predicted.values-pre$standard.errors
+  pre$upper = pre$predicted.values+pre$standard.errors
+
+  # 3. Mean comparison on germplasm ----------
   lsd = agricolae::LSD.test(
-    y = data[,variable],
-    trt = data$germplasm,
+    y = pre$predicted.values,
+    trt = pre$germplasm,
     DFerror = x$spatial$model$df_residual, 
     MSerror = x$spatial$model$MSerror, 
     alpha = alpha, 
@@ -73,6 +81,7 @@ mean_comparisons.check_model_spatial <- function(
   # 3. return results
   out <- list(
     "info" = info,
+    "blups_prediction" = pre,
     "data_ggplot_LSDbarplot_germplasm" = data_ggplot_LSDbarplot_germplasm
   )
   
diff --git a/R/plot.mean_comparisons_model_spatial.R b/R/plot.mean_comparisons_model_spatial.R
index 05e498f..259c7c1 100644
--- a/R/plot.mean_comparisons_model_spatial.R
+++ b/R/plot.mean_comparisons_model_spatial.R
@@ -71,9 +71,15 @@ plot.mean_comparisons_model_spatial <- function(
     ggplot_LSDbarplot_germplasm = NULL
   }
   
-  # 6. return results ----------
+  # 4. Germplasm BLUPs with confidence intervalle
+  pre = x$blups_prediction
+  p_blup = ggplot(pre, aes(y = predicted.values, x = germplasm)) + geom_point() 
+  p_blup = p_blup + geom_errorbar(aes(ymin = lower, ymax = upper)) + coord_flip()
+  
+  # 5. return results ----------
   out = list(
-    "germplasm" = ggplot_LSDbarplot_germplasm
+    "germplasm_blup" = p_blup,
+    "germplasm_barplot" = ggplot_LSDbarplot_germplasm
   )
   
   return(out)
diff --git a/inst/bookdown/307-agro_spatial_analysis.Rmd b/inst/bookdown/307-agro_spatial_analysis.Rmd
index 94ae912..f6e67f6 100644
--- a/inst/bookdown/307-agro_spatial_analysis.Rmd
+++ b/inst/bookdown/307-agro_spatial_analysis.Rmd
@@ -151,16 +151,20 @@ out_mean_comparisons_spatial = mean_comparisons(out_check_spatial, p.adj = "bonf
 p_out_mean_comparisons_spatial = plot(out_mean_comparisons_spatial)
 ```
 
-`p_out_mean_comparisons_spatial` is a list of one elements with barplots :
+`p_out_mean_comparisons_spatial` is a list of two elements
 
+- `germplasm_blup` whih represent the BLUPs with their confidence intervalle
+```{r}
+p_out_mean_comparisons_spatial$germplasm_blup
+```
+
+- `germplasm_barplot`. 
 For each element of the list, there are as many graph as needed with `nb_parameters_per_plot` parameters per graph.
 Letters are displayed on each bar. Parameters that do not share the same letters are different regarding type I error (alpha) and alpha correction. 
 The error I (alpha) and the alpha correction are displayed in the title. 
 
-`germplasm` : mean comparison for germplasm
 ```{r}
-pg = p_out_mean_comparisons_spatial$germplasm
-names(pg)
+pg = p_out_mean_comparisons_spatial$germplasm_barplot
 pg$`23`
 ```
 
@@ -195,4 +199,3 @@ workflow_spatial = function(x, data){
 # names(out) = vec_variables
 ```
 
-