fix: setting treatment to 'None' did not exclude it from the (g)lmm a…

…nalysis (close #224)

CHANGELOG

@@ -6,8 +6,11 @@
      ambiguities (#223)
    - Allow arbitrary number of repetitions for mixed model analysis
      with a spin control (#222)
+   - Setting the treatment to "None" for (g)lmm analysis did not
+     exclude it from the analysis, but introduced it as a different
+     treatment (#224)
  - Automatically identify treatment and repetition for mixed model
-   analysis
+   analysis (#222)
  - Automated release to PyPI with appveyor
  - Begin migration of docs to sphinx/readthedocs.org
 0.8.6

shapeout/lin_mix_mod.py

@@ -300,6 +300,24 @@ def linmixmod(xs, treatment, timeunit, model='lmm', RCMD=cran.rcmd):
               "Please select more treatment repetitions."
         raise ValueError(msg)
 
+    # Check that names are valid
+    for trt in treatment:
+        if trt not in ["None",
+                       "Control",
+                       "Reservoir Control",
+                       "Treatment",
+                       "Reservoir Treatment"]:
+            raise ValueError("Unknown treatment: '{}'".format(trt))
+
+    # Remove "None"s
+    treatment = np.array(treatment)
+    timeunit = np.array(timeunit)
+    xs = np.array(xs)
+    invalid = (treatment == "None")
+    treatment = list(treatment[~invalid])
+    timeunit = list(timeunit[~invalid])
+    xs = [xi for ii, xi in enumerate(xs) if ~invalid[ii]]
+
     ######################Differential Deformation#############################
     # If the user selected 'Control-Reservoir' and/or 'Treatment-Reservoir'
     Median_DiffDef = []
@@ -376,7 +394,7 @@ def linmixmod(xs, treatment, timeunit, model='lmm', RCMD=cran.rcmd):
         treatment = np.array(Treatment)
         timeunit = np.array(TimeUnit)
 
-    else:  # If there is no 'Reservoir Channel' selected dont apply bootstrapping
+    else:  # If there is no 'Reservoir Channel' selected don't apply bootstrapping
         if model == 'glmm':
             Head_string = "GENERALIZED LINEAR MIXED MODEL: \n" +\
                 "---Results are in log space (loglink was used)--- \n"

tests/test_linmixmod.py

@@ -16,15 +16,16 @@
 
 
 def test_linmixmod():
-    treatment = ['Control', 'Drug', 'Control', 'Drug']
+    treatment = ['Control', 'Treatment', 'Control', 'Treatment']
     timeunit = [1, 1, 2, 2]
     xs = [
         [100, 99, 80, 120, 140, 150, 100, 100, 110, 111, 140, 145],
-        [115, 110, 90, 110, 145, 155, 110, 120, 115, 120, 120, 150, 100, 90, 100],
+        [115, 110, 90, 110, 145, 155, 110, 120, 115, 120, 120, 150,
+         100, 90, 100],
         [150, 150, 130, 170, 190, 250, 150, 150, 160,
-            161, 180, 195, 130, 120, 125, 130, 125],
+         161, 180, 195, 130, 120, 125, 130, 125],
         [155, 155, 135, 175, 195, 255, 155, 155, 165, 165,
-            185, 200, 135, 125, 130, 135, 140, 150, 135, 140]
+         185, 200, 135, 125, 130, 135, 140, 150, 135, 140]
     ]
 
     res = linmixmod(xs=xs, treatment=treatment, timeunit=timeunit)

tests/test_linmixmod2.py

@@ -93,6 +93,22 @@ def test_linmixmod_4():
         [Result_4["p-Value (Likelihood Ratio Test)"]], [0.00365675950677214])
 
 
+def test_linmixmod_5():
+    # 5.: Add None values and get same result as in 4.
+    treatment5 = ['Treatment', 'Control', 'Treatment', 'Control',
+                  'None', 'Treatment', 'Control', 'Treatment', 'Control']
+    timeunit5 = [1, 1, 2, 2, 1, 3, 3, 4, 4]
+    xs2 = copy.deepcopy(xs)
+    xs2.insert(4, xs[0])
+    Result_5 = linmixmod(xs=xs2, treatment=treatment5,
+                         timeunit=timeunit5, model='glmm')
+    assert np.allclose([Result_5["Estimate"]], [2.71362344639])
+    assert 'BOOTSTAP-DISTRIBUTIONS' not in Result_5['Full Summary']
+    assert 'GENERALIZED' in Result_5['Full Summary']
+    assert np.allclose(
+        [Result_5["p-Value (Likelihood Ratio Test)"]], [0.00365675950677214])
+
+
 def test_diffdef():
     # Larger values (Channel1)
     y = np.array([100, 99, 80, 120, 140, 150, 100, 100, 110, 111, 140, 145])