diff --git a/extreme_fit/estimator/abstract_estimator.py b/extreme_fit/estimator/abstract_estimator.py
index c8e93cad38e6a333e301105a652b9dfb7a053577..39a695afe6b7d77ffcb6f52791547582d85f644f 100644
--- a/extreme_fit/estimator/abstract_estimator.py
+++ b/extreme_fit/estimator/abstract_estimator.py
@@ -11,7 +11,7 @@ from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 
 class AbstractEstimator(object):
 
-    def __init__(self, dataset: AbstractDataset):
+    def __init__(self, dataset: AbstractDataset, **kwargs):
         self.dataset = dataset  # type: AbstractDataset
         self._result_from_fit = None  # type: Union[None, AbstractResultFromModelFit]
 
diff --git a/extreme_fit/estimator/quantile_estimator/abstract_quantile_estimator.py b/extreme_fit/estimator/quantile_estimator/abstract_quantile_estimator.py
index d43aae929945e9ec8b9f27b4a4c9cf657017b57c..5f0966b57abab724449247424b0424ecfda5775a 100644
--- a/extreme_fit/estimator/quantile_estimator/abstract_quantile_estimator.py
+++ b/extreme_fit/estimator/quantile_estimator/abstract_quantile_estimator.py
@@ -3,14 +3,17 @@ from cached_property import cached_property
 from extreme_fit.estimator.abstract_estimator import AbstractEstimator
 from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
 from extreme_fit.function.abstract_quantile_function import AbstractQuantileFunction
+from extreme_fit.function.margin_function.abstract_margin_function import AbstractMarginFunction
 from extreme_fit.model.margin_model.linear_margin_model.linear_margin_model import LinearMarginModel
+from extreme_fit.model.quantile_model.quantile_regression_model import AbstractQuantileRegressionModel
+from extreme_fit.model.result_from_model_fit.abstract_result_from_model_fit import AbstractResultFromModelFit
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 
 
 class AbstractQuantileEstimator(AbstractEstimator):
 
-    def __init__(self, quantile: float, **kwargs):
-        super().__init__(**kwargs)
+    def __init__(self, dataset: AbstractDataset, quantile: float, **kwargs):
+        super().__init__(dataset, **kwargs)
         assert 0 < quantile < 1
         self.quantile = quantile
 
@@ -26,5 +29,19 @@ class QuantileEstimatorFromMargin(AbstractQuantileEstimator, LinearMarginEstimat
 
     @cached_property
     def quantile_function_from_fit(self) -> AbstractQuantileFunction:
-        linear_margin_function = super().function_from_fit
+        linear_margin_function = super().function_from_fit  # type: AbstractMarginFunction
         return AbstractQuantileFunction(linear_margin_function, self.quantile)
+
+
+class QuantileRegressionEstimator(AbstractQuantileEstimator):
+
+    def __init__(self, dataset: AbstractDataset, quantile: float, quantile_regression_model_class: type, **kwargs):
+        super().__init__(dataset, quantile, **kwargs)
+        self.quantile_regression_model = quantile_regression_model_class(dataset, quantile) # type: AbstractQuantileRegressionModel
+
+    def _fit(self) -> AbstractResultFromModelFit:
+        return self.quantile_regression_model.fit()
+
+    @cached_property
+    def quantile_function_from_fit(self) -> AbstractQuantileFunction:
+        return self.result_from_model_fit.quantile_function
diff --git a/extreme_fit/function/abstract_quantile_function.py b/extreme_fit/function/abstract_quantile_function.py
index 4c06c95ec5f30bbc8ab287cbc517a73674faaf84..62596039e775d38726f1cf04600857ae3c7afbe6 100644
--- a/extreme_fit/function/abstract_quantile_function.py
+++ b/extreme_fit/function/abstract_quantile_function.py
@@ -12,4 +12,9 @@ class AbstractQuantileFunction(AbstractFunction):
 
     def get_quantile(self, coordinate: np.ndarray) -> float:
         gev_params = self.margin_function.get_gev_params(coordinate)
-        return gev_params.quantile(self.quantile)
\ No newline at end of file
+        return gev_params.quantile(self.quantile)
+
+    def visualize(self):
+        pass
+        # for coordine
+        # self.margin_function.
\ No newline at end of file
diff --git a/extreme_fit/model/quantile_model/quantile_regression_model.py b/extreme_fit/model/quantile_model/quantile_regression_model.py
new file mode 100644
index 0000000000000000000000000000000000000000..2652d5bbc843e30507fd7f4e052c0f0fc9e0dc80
--- /dev/null
+++ b/extreme_fit/model/quantile_model/quantile_regression_model.py
@@ -0,0 +1,57 @@
+import numpy as np
+from rpy2 import robjects
+
+from extreme_fit.model.abstract_model import AbstractModel
+from extreme_fit.model.result_from_model_fit.result_from_quantilreg import ResultFromQuantreg
+from extreme_fit.model.utils import r, safe_run_r_estimator, get_coord_df
+from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+
+
+class AbstractQuantileRegressionModel(AbstractModel):
+
+    def __init__(self, dataset: AbstractDataset, quantile: float):
+        self.dataset = dataset
+        self.quantile = quantile
+
+    @property
+    def data(self):
+        return get_coord_df(self.dataset.df_dataset)
+
+    @property
+    def first_column_of_observation(self):
+        return self.data.colnames[1]
+        # print(self.dataset.df_dataset.columns)
+        # return str(self.dataset.df_dataset.columns[0])
+
+    def fit(self):
+        print(self.data)
+        parameters = {
+            'tau': self.quantile,
+            'data': self.data,
+            'formula': self.formula
+
+        }
+        res = safe_run_r_estimator(r.rq, **parameters)
+        return ResultFromQuantreg(res)
+
+    @property
+    def formula_str(self):
+        raise NotImplementedError
+
+    @property
+    def formula(self):
+        return robjects.Formula(self.first_column_of_observation + '~ ' + self.formula_str)
+
+
+class ConstantQuantileRegressionModel(AbstractQuantileRegressionModel):
+
+    @property
+    def formula_str(self):
+        return '1'
+
+
+class AllCoordinatesQuantileRegressionModel(AbstractQuantileRegressionModel):
+
+    @property
+    def formula_str(self):
+        return  '+'.join(self.dataset.coordinates.coordinates_names[-1:])
diff --git a/extreme_fit/model/result_from_model_fit/result_from_quantilreg.py b/extreme_fit/model/result_from_model_fit/result_from_quantilreg.py
new file mode 100644
index 0000000000000000000000000000000000000000..2be5c564daeb1c67c20f60bbae58a67891a5eb86
--- /dev/null
+++ b/extreme_fit/model/result_from_model_fit/result_from_quantilreg.py
@@ -0,0 +1,14 @@
+from cached_property import cached_property
+
+from extreme_fit.model.result_from_model_fit.abstract_result_from_model_fit import AbstractResultFromModelFit
+
+
+class ResultFromQuantreg(AbstractResultFromModelFit):
+
+    @property
+    def coefficients(self):
+        return self.name_to_value['coefficients']
+
+    @cached_property
+    def quantile_function(self):
+        print(self.coefficients)
\ No newline at end of file
diff --git a/extreme_fit/model/utils.py b/extreme_fit/model/utils.py
index 658c2fcbade7865c3538945c60baac81638f665d..1582cf336f615ef09692bb4cc3ac06b8c305bfc1 100644
--- a/extreme_fit/model/utils.py
+++ b/extreme_fit/model/utils.py
@@ -23,6 +23,7 @@ numpy2ri.activate()
 pandas2ri.activate()
 r.library('SpatialExtremes')
 r.library('data.table')
+r.library('quantreg')
 # Desactivate temporarily warnings
 default_filters = warnings.filters.copy()
 warnings.filterwarnings("ignore")
@@ -125,12 +126,11 @@ def get_coord(df_coordinates: pd.DataFrame):
 def get_coord_df(df_coordinates: pd.DataFrame):
     coord = pandas2ri.py2ri_pandasdataframe(df_coordinates)
     # coord = r.transpose(coord)
-    colname = df_coordinates.columns[0]
+    colname = df_coordinates.columns
     coord.colnames = r.c(colname)
     coord = r('data.frame')(coord, stringsAsFactors=True)
     return coord
 
-
 def get_null():
     as_null = r['as.null']
     return as_null(1.0)
diff --git a/spatio_temporal_dataset/spatio_temporal_observations/abstract_spatio_temporal_observations.py b/spatio_temporal_dataset/spatio_temporal_observations/abstract_spatio_temporal_observations.py
index fe599b2c50a3b11e965a1a216b4aa24f1c253327..0b009da9e4aa7b152362e0e46399870d115a7354 100644
--- a/spatio_temporal_dataset/spatio_temporal_observations/abstract_spatio_temporal_observations.py
+++ b/spatio_temporal_dataset/spatio_temporal_observations/abstract_spatio_temporal_observations.py
@@ -53,7 +53,7 @@ class AbstractSpatioTemporalObservations(object):
                            (self.df_maxima_frech, self.OBSERVATIONS_FRECH)]:
             if df is not None:
                 df_maxima = df.copy()
-                df_maxima.columns = [str(c) + ' ' + suffix for c in df_maxima.columns]
+                df_maxima.columns = [str(c) for c in df_maxima.columns]
                 df_maxima_list.append(df_maxima)
         return pd.concat(df_maxima_list, axis=1)
 
diff --git a/test/test_extreme_fit/test_estimator/test_quantile_estimator.py b/test/test_extreme_fit/test_estimator/test_quantile_estimator.py
index bd517938b1dec9cb8c69ad594b9407b084942552..0a304b2750be3d1c1110ab154a41997da58c21e5 100644
--- a/test/test_extreme_fit/test_estimator/test_quantile_estimator.py
+++ b/test/test_extreme_fit/test_estimator/test_quantile_estimator.py
@@ -1,17 +1,18 @@
 import unittest
 
-from extreme_fit.estimator.quantile_estimator.abstract_quantile_estimator import QuantileEstimatorFromMargin
+from extreme_fit.estimator.quantile_estimator.abstract_quantile_estimator import QuantileEstimatorFromMargin, \
+    QuantileRegressionEstimator
 from spatio_temporal_dataset.dataset.simulation_dataset import MarginDataset
 from test.test_utils import load_test_1D_and_2D_spatial_coordinates, load_test_spatiotemporal_coordinates, \
-    load_smooth_margin_models
+    load_smooth_margin_models, load_smooth_quantile_model_classes
 
 
-class TestSmoothMarginEstimator(unittest.TestCase):
+class TestQuantileEstimator(unittest.TestCase):
     DISPLAY = False
 
     def test_smooth_margin_estimator_spatial(self):
-        self.nb_points = 2
-        self.nb_obs = 2
+        self.nb_points = 20
+        self.nb_obs = 1
         self.coordinates = load_test_1D_and_2D_spatial_coordinates(nb_points=self.nb_points)
 
     def test_smooth_margin_estimator_spatio_temporal(self):
@@ -21,21 +22,28 @@ class TestSmoothMarginEstimator(unittest.TestCase):
         self.coordinates = load_test_spatiotemporal_coordinates(nb_steps=self.nb_steps, nb_points=self.nb_points)
 
     def tearDown(self) -> None:
-        quantile = 0.98
+        quantile = 0.5
         for coordinates in self.coordinates:
             constant_margin_model = load_smooth_margin_models(coordinates=coordinates)[0]
             dataset = MarginDataset.from_sampling(nb_obs=self.nb_obs,
                                                   margin_model=constant_margin_model,
                                                   coordinates=coordinates)
-            quantile_estimators = [QuantileEstimatorFromMargin(dataset, constant_margin_model, quantile)]
-            help(QuantileEstimatorFromMargin)
-
+            print(dataset)
+            # Load quantile estimators
+            quantile_estimators = [
+                QuantileEstimatorFromMargin(dataset, constant_margin_model, quantile),
+            ]
+            for quantile_model_class in load_smooth_quantile_model_classes()[:1]:
+                quantile_estimator = QuantileRegressionEstimator(dataset, quantile, quantile_model_class)
+                quantile_estimators.append(quantile_estimator)
+
+            # Fit quantile estimators
             for quantile_estimator in quantile_estimators:
                 quantile_estimator.fit()
-                print(quantile_estimator.function_from_fit)
+                print(quantile_estimator.quantile_function_from_fit)
 
-        # self.assertTrue(True)
+        self.assertTrue(True)
 
 
 if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
+    unittest.main()
diff --git a/test/test_utils.py b/test/test_utils.py
index d5fe754c500eac5abe2bbd6e7372eb2e618dc6b3..bc2bf761591cb65a65c758c5a646395ce8a711d2 100644
--- a/test/test_utils.py
+++ b/test/test_utils.py
@@ -6,7 +6,8 @@ from experiment.meteo_france_data.scm_models_data.crocus.crocus import Crocus, C
 from extreme_fit.estimator.full_estimator.abstract_full_estimator import SmoothMarginalsThenUnitaryMsp, \
     FullEstimatorInASingleStepWithSmoothMargin
 from extreme_fit.estimator.max_stable_estimator.abstract_max_stable_estimator import MaxStableEstimator
-from extreme_fit.model.margin_model.linear_margin_model.linear_margin_model import LinearAllParametersAllDimsMarginModel, \
+from extreme_fit.model.margin_model.linear_margin_model.linear_margin_model import \
+    LinearAllParametersAllDimsMarginModel, \
     ConstantMarginModel
 from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import \
     NonStationaryLocationTemporalModel, NonStationaryScaleTemporalModel, NonStationaryShapeTemporalModel
@@ -16,6 +17,8 @@ from extreme_fit.model.max_stable_model.max_stable_models import Smith, BrownRes
     Geometric, ExtremalT, ISchlather
 from experiment.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall, Safran, SafranRainfall, \
     SafranTemperature, SafranPrecipitation
+from extreme_fit.model.quantile_model.quantile_regression_model import ConstantQuantileRegressionModel, \
+    AllCoordinatesQuantileRegressionModel
 from spatio_temporal_dataset.coordinates.spatial_coordinates.abstract_spatial_coordinates import \
     AbstractSpatialCoordinates
 from spatio_temporal_dataset.coordinates.spatial_coordinates.alps_station_3D_coordinates import \
@@ -40,7 +43,9 @@ TEST_3D_SPATIAL_COORDINATES = [AlpsStation3DCoordinatesWithAnisotropy]
 TEST_TEMPORAL_COORDINATES = [ConsecutiveTemporalCoordinates]
 TEST_SPATIO_TEMPORAL_COORDINATES = [UniformSpatioTemporalCoordinates, LinSpaceSpatial2DSpatioTemporalCoordinates]
 TEST_MARGIN_TYPES = [ConstantMarginModel, LinearAllParametersAllDimsMarginModel][:]
-TEST_NON_STATIONARY_TEMPORAL_MARGIN_TYPES = [NonStationaryLocationTemporalModel, NonStationaryScaleTemporalModel, NonStationaryShapeTemporalModel]
+TEST_QUANTILES_TYPES = [ConstantQuantileRegressionModel, AllCoordinatesQuantileRegressionModel][:]
+TEST_NON_STATIONARY_TEMPORAL_MARGIN_TYPES = [NonStationaryLocationTemporalModel, NonStationaryScaleTemporalModel,
+                                             NonStationaryShapeTemporalModel]
 TEST_MAX_STABLE_ESTIMATOR = [MaxStableEstimator]
 TEST_FULL_ESTIMATORS = [SmoothMarginalsThenUnitaryMsp, FullEstimatorInASingleStepWithSmoothMargin][:]
 
@@ -62,6 +67,10 @@ def load_smooth_margin_models(coordinates):
     return [margin_class(coordinates=coordinates) for margin_class in TEST_MARGIN_TYPES]
 
 
+def load_smooth_quantile_model_classes():
+    return [quantile_reg_class for quantile_reg_class in TEST_QUANTILES_TYPES]
+
+
 def load_test_max_stable_models(default_covariance_function=None):
     # Load all max stable model
     max_stable_models = []
@@ -84,14 +93,16 @@ def load_test_spatial_coordinates(nb_points, coordinate_types, train_split_ratio
             for coordinate_class in coordinate_types]
 
 
-def load_test_1D_and_2D_spatial_coordinates(nb_points, train_split_ratio=None, transformation_class=None) -> List[AbstractSpatialCoordinates]:
+def load_test_1D_and_2D_spatial_coordinates(nb_points, train_split_ratio=None, transformation_class=None) -> List[
+    AbstractSpatialCoordinates]:
     return load_test_spatial_coordinates(nb_points, TEST_1D_AND_2D_SPATIAL_COORDINATES,
                                          train_split_ratio=train_split_ratio,
                                          transformation_class=transformation_class)
 
 
 def load_test_3D_spatial_coordinates(nb_points, transformation_class=None) -> List[AbstractSpatialCoordinates]:
-    return load_test_spatial_coordinates(nb_points, TEST_3D_SPATIAL_COORDINATES, transformation_class=transformation_class)
+    return load_test_spatial_coordinates(nb_points, TEST_3D_SPATIAL_COORDINATES,
+                                         transformation_class=transformation_class)
 
 
 def load_test_temporal_coordinates(nb_steps, train_split_ratio=None, transformation_class=None):