diff --git a/clinica/pipelines/statistics_surface/_model.py b/clinica/pipelines/statistics_surface/_model.py
index bb496c8dec..3b974b17b5 100644
--- a/clinica/pipelines/statistics_surface/_model.py
+++ b/clinica/pipelines/statistics_surface/_model.py
@@ -670,9 +670,12 @@ def _write_to_json(self, result: StatisticsResults):
         ----------
         results : Results to write to disk in JSON format.
         """
+        import os
         import json
 
-        out_json_file = str(self.output_file) + self.json_extension
+        out_json_file = Path(str(self.output_file) + self.json_extension)
+        if not os.path.exists(out_json_file.parents[0]):
+            os.makedirs(out_json_file.parents[0])
         cprint(
             msg=f"Writing results to JSON in {out_json_file}...",
             lvl="info",
diff --git a/test/unittests/pipelines/statistics_surface/test_model.py b/test/unittests/pipelines/statistics_surface/test_model.py
index 2fa0498809..efeea456bf 100644
--- a/test/unittests/pipelines/statistics_surface/test_model.py
+++ b/test/unittests/pipelines/statistics_surface/test_model.py
@@ -253,3 +253,76 @@ def test_group_glm_with_interaction_instantiation(df):
     ):
         model.filename_root("foo")
     assert model.filename_root("age * sex") == "interaction-age * sex_measure-feature_label_fwhm-2.0"
+
+
+def test_glm_factory(df):
+    from clinica.pipelines.statistics_surface._model import (
+        GLMFactory, GroupGLM, GroupGLMWithInteraction, CorrelationGLM
+    )
+    factory = GLMFactory("feature_label")
+    assert factory.feature_label == "feature_label"
+    parameters = {
+        "group_label": "group_label",
+        "sizeoffwhm": 2.0,
+    }
+    model = factory.create_model("correlation", "age", df, "age", **parameters)
+    assert isinstance(model, CorrelationGLM)
+    model = factory.create_model("group_comparison", "age", df, "sex", **parameters)
+    assert isinstance(model, GroupGLM)
+    model = factory.create_model("group_comparison", "age * sex", df, "age * sex", **parameters)
+    assert isinstance(model, GroupGLMWithInteraction)
+
+
+def test_p_value_results():
+    from clinica.pipelines.statistics_surface._model import PValueResults
+    pvalues = np.random.random((5, 10))
+    threshold = .3
+    mask = pvalues >= threshold
+    results = PValueResults(pvalues, mask, threshold)
+    d = results.to_dict(jsonable=False)
+    assert_array_equal(d["P"], pvalues)
+    assert_array_equal(d["mask"], mask)
+    assert d["thresh"] == threshold
+    d = results.to_dict(jsonable=True)
+    assert d["P"] == pvalues.tolist()
+    assert d["mask"] == mask.tolist()
+    assert d["thresh"] == threshold
+
+
+def test_statistics_results_serializer(tmp_path):
+    import json
+    from scipy.io import loadmat
+    from clinica.pipelines.statistics_surface._model import (
+        StatisticsResults, StatisticsResultsSerializer,
+        PValueResults, CorrectedPValueResults
+    )
+    dummy_input = np.empty([3, 6])
+    uncorrected = PValueResults(*[dummy_input]*2, 0.1)
+    corrected = CorrectedPValueResults(*[dummy_input]*3, 0.2)
+    results = StatisticsResults(*[dummy_input]*2, uncorrected, dummy_input, corrected)
+    serializer = StatisticsResultsSerializer(str(tmp_path / Path("out/dummy")))
+    assert serializer.output_file == str(tmp_path / Path("out/dummy"))
+    assert serializer.json_extension == "_results.json"
+    assert serializer.json_indent == 4
+    assert serializer.mat_extension == ".mat"
+    with pytest.raises(
+        NotImplementedError,
+        match="Serializing method foo is not implemented."
+    ):
+        serializer.save(results, "foo")
+    serializer.save(results, "json")
+    assert os.path.exists(tmp_path / Path("out/dummy_results.json"))
+    with open(tmp_path / Path("out/dummy_results.json"), "r") as fp:
+        serialized = json.load(fp)
+    assert results.to_dict(jsonable=True) == serialized
+    serializer.save(results, "mat")
+    names = ["coefficients", "TStatistics", "uncorrectedPValues", "correctedPValues", "FDR"]
+    keys = ["coef", "tvaluewithmask", "uncorrectedpvaluesstruct", "correctedpvaluesstruct", "FDR"]
+    for name, key in zip(names, keys):
+        assert os.path.exists(tmp_path / Path(f"out/dummy_{name}.mat"))
+        mat = loadmat(tmp_path / Path(f"out/dummy_{name}.mat"))
+        if key in ["uncorrectedpvaluesstruct", "correctedpvaluesstruct"]:
+            assert_array_equal(mat[key]['P'][0, 0], dummy_input)
+            assert_array_equal(mat[key]['mask'][0, 0], dummy_input)
+        else:
+            assert_array_equal(mat[key], dummy_input)