diff --git a/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R b/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R
index 74d853e86abf6ecc88efac394ea4def6524131d9..92317a48e2225d56f8d84aeccc5a92de2b62cc66 100644
--- a/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R
+++ b/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R
@@ -27,7 +27,7 @@ coord = data.frame(coord, stringsAsFactors = TRUE)
 # res = fevd_fixed(x_gev, data=coord, method='MLE', verbose=TRUE, use.phi=FALSE)
 # res = fevd_fixed(x_gev, data=coord, location.fun= ~T, scale.fun= ~T, method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
 # res = fevd_fixed(x_gev, data=coord, location.fun= ~sin(X) + cos(T), method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
-res = fevd_fixed(x_gev, data=coord, location.fun= ~poly(X, 1, raw = TRUE) + 1 , method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
+res = fevd_fixed(x_gev, data=coord, location.fun= ~poly(X, 1, raw = TRUE) + poly(T, 2, raw = TRUE) , method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
 print(res)
 
 # Some display for the results
diff --git a/extreme_fit/function/margin_function/linear_margin_function.py b/extreme_fit/function/margin_function/linear_margin_function.py
index 3ef5baaecd8e6a0f6c15c6b3ee1e4566eeb258f6..d9b6b2733c6f62b92557b521f6e152677d3c8810 100644
--- a/extreme_fit/function/margin_function/linear_margin_function.py
+++ b/extreme_fit/function/margin_function/linear_margin_function.py
@@ -66,22 +66,26 @@ class LinearMarginFunction(ParametricMarginFunction):
 
     @property
     def form_dict(self) -> Dict[str, str]:
+        coordinate_name_to_dim = self.coefficient_name_to_dim(self.coordinates)
         form_dict = {}
         for param_name in self.params_class.PARAM_NAMES:
             linear_dims = self.param_name_to_dims.get(param_name, [])
             # Load spatial form_dict (only if we have some spatial coordinates)
             if self.coordinates.has_spatial_coordinates:
                 spatial_names = [name for name in self.coordinates.spatial_coordinates_names
-                                 if self.coefficient_name_to_dim(self.coordinates)[name] in linear_dims]
-                spatial_form = self.param_name_to_coef[param_name].spatial_form_dict(spatial_names)
+                                 if coordinate_name_to_dim[name] in linear_dims]
+                spatial_dims = [coordinate_name_to_dim[name] for name in spatial_names]
+                spatial_form = self.param_name_to_coef[param_name].spatial_form_dict(spatial_names, spatial_dims)
                 form_dict.update(spatial_form)
             # Load temporal form dict (only if we have some temporal coordinates)
+
             if self.coordinates.has_temporal_coordinates:
                 temporal_names = [name for name in self.coordinates.temporal_coordinates_names
-                                  if self.coefficient_name_to_dim(self.coordinates)[name] in linear_dims]
-                temporal_form = self.param_name_to_coef[param_name].temporal_form_dict(temporal_names)
+                                  if coordinate_name_to_dim[name] in linear_dims]
+                temporal_dims = [coordinate_name_to_dim[name] for name in temporal_names]
+                temporal_form = self.param_name_to_coef[param_name].temporal_form_dict(temporal_names, temporal_dims)
                 # Specifying a formula '~ 1' creates a bug in fitspatgev of SpatialExtreme R package
-                assert not any(['1' in formula for formula in temporal_form.values()])
+                assert not any(['~ 1' in formula for formula in temporal_form.values()])
                 form_dict.update(temporal_form)
         return form_dict
 
diff --git a/extreme_fit/function/param_function/abstract_coef.py b/extreme_fit/function/param_function/abstract_coef.py
index 28dfa2300bf816873c540c9a8df647919582f6e3..95df0f9d6422f02ae03b4821cf882c1a1ab72980 100644
--- a/extreme_fit/function/param_function/abstract_coef.py
+++ b/extreme_fit/function/param_function/abstract_coef.py
@@ -34,6 +34,6 @@ class AbstractCoef(object):
 
     """ Form dict """
 
-    def form_dict(self, names: List[str]) -> Dict[str, str]:
+    def form_dict(self, names: List[str], dims) -> Dict[str, str]:
         formula_str = '1' if not names else '+'.join(names)
         return {self.param_name + '.form': self.param_name + ' ~ ' + formula_str}
\ No newline at end of file
diff --git a/extreme_fit/function/param_function/linear_coef.py b/extreme_fit/function/param_function/linear_coef.py
index 5ea2d0cc609253137e804055a6f446d8f752ad86..f139dfe4596080c9996e5376cd25ac19b53a17c5 100644
--- a/extreme_fit/function/param_function/linear_coef.py
+++ b/extreme_fit/function/param_function/linear_coef.py
@@ -75,7 +75,7 @@ class LinearCoef(AbstractCoef):
             j += 1
         return coef_dict
 
-    def spatial_form_dict(self, coordinate_spatial_names: List[str]) -> Dict[str, str]:
+    def spatial_form_dict(self, coordinate_spatial_names: List[str], spatial_dims) -> Dict[str, str]:
         """
         Example of formula that could be specified:
         loc.form = loc ~ 1
@@ -85,9 +85,9 @@ class LinearCoef(AbstractCoef):
         :return:
         """
         assert all([name in AbstractCoordinates.COORDINATE_SPATIAL_NAMES for name in coordinate_spatial_names])
-        return self.form_dict(coordinate_spatial_names)
+        return self.form_dict(coordinate_spatial_names, spatial_dims)
 
-    def temporal_form_dict(self, coordinate_temporal_names: List[str]) -> Dict[str, str]:
+    def temporal_form_dict(self, coordinate_temporal_names: List[str], temporal_dims) -> Dict[str, str]:
         """
         Example of formula that could be specified:
         temp.form.loc = loc ~ coord_t
@@ -96,7 +96,7 @@ class LinearCoef(AbstractCoef):
         :return:
         """
         assert all([name in [AbstractCoordinates.COORDINATE_T] for name in coordinate_temporal_names])
-        k, v = self.form_dict(coordinate_temporal_names).popitem()
+        k, v = self.form_dict(coordinate_temporal_names, temporal_dims).popitem()
         k = 'temp.form.' + k.split('.')[0]
-        v = 'NULL' if '1' in v else v
+        v = 'NULL' if '~ 1' in v else v
         return {k: v}
diff --git a/extreme_fit/function/param_function/polynomial_coef.py b/extreme_fit/function/param_function/polynomial_coef.py
index e5784bd10811efd51ee3b451d8b11a118ffc6619..fa7c61984edb6c7cbd01c5cf08707d18a508416a 100644
--- a/extreme_fit/function/param_function/polynomial_coef.py
+++ b/extreme_fit/function/param_function/polynomial_coef.py
@@ -70,7 +70,7 @@ class PolynomialAllCoef(LinearCoef):
                 dim_to_polynomial_coef[dim] = PolynomialCoef(param_name=param_name, degree_to_coef=degree_to_coef)
         return cls(param_name=param_name, dim_to_polynomial_coef=dim_to_polynomial_coef, intercept=intercept)
 
-    def form_dict(self, coordinates_names: List[str]) -> Dict[str, str]:
+    def form_dict(self, coordinates_names: List[str], dims) -> Dict[str, str]:
         if len(coordinates_names) >= 2:
             raise NotImplementedError(
                 'Check how do we sum two polynomails without having two times an intercept parameter')
@@ -78,7 +78,7 @@ class PolynomialAllCoef(LinearCoef):
         if len(coordinates_names) == 0:
             formula_str = '1'
         else:
-            for dim, name in enumerate(coordinates_names):
+            for dim, name in zip(dims, coordinates_names):
                 polynomial_coef = self.dim_to_polynomial_coef[dim]
                 formula_list.append('poly({}, {}, raw = TRUE)'.format(name, polynomial_coef.max_degree))
             formula_str = ' '.join(formula_list)
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
index 4425b7ea2550dcd3c1a9e6a3890a863c8339ff7a..4346be81b686c411ce688a98454a09672ff9850d 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
@@ -19,6 +19,7 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
     @property
     def margin_coef_ordered_dict(self):
         values = self.name_to_value['results']
+        print(values)
         d = self.get_python_dictionary(values)
         if 'par' in d:
             values = {i: param for i, param in enumerate(np.array(d['par']))}
diff --git a/test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py b/test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py
index e5ac3c51719c589e21aec5f0e59bcf406cc939ba..c28f3503c598410b76a2d0b0769dfeb6c5f7fbf2 100644
--- a/test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py
+++ b/test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py
@@ -41,7 +41,6 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
         sample_fit = massif_fit.sample_id_to_sample_fit[0]
         model_fit = sample_fit.model_class_to_model_fit[model_class]  # type: TwoFoldModelFit
         estimator = model_fit.estimator_fold_1
-        print(estimator.dataset)
         return estimator
 
     # def test_location_spatio_temporal_linearity(self):