improved parameter name handling if constraints are given

Author: davidsebfischer

batchglm/data.py

@@ -21,7 +21,7 @@ def design_matrix(
         as_categorical: Union[bool, list] = True,
         dmat: Union[pd.DataFrame, None] = None,
         return_type: str = "patsy",
-) -> Union[patsy.design_info.DesignMatrix, pd.DataFrame]:
+) -> Tuple[Union[patsy.design_info.DesignMatrix, pd.DataFrame], List[str]]:
     """
     Create a design matrix from some sample description.
 
@@ -62,6 +62,7 @@ def design_matrix(
                     sample_description[col] = sample_description[col].astype("category")
 
         dmat = patsy.dmatrix(formula, sample_description)
+        coef_names = dmat.design_info.column_names
 
         if return_type == "dataframe":
             df = pd.DataFrame(dmat, columns=dmat.design_info.column_names)
@@ -70,12 +71,12 @@ def design_matrix(
 
             return df
         elif return_type == "patsy":
-            return dmat
+            return dmat, coef_names
         else:
             raise ValueError("return type %s not recognized" % return_type)
     else:
         if return_type == "dataframe":
-            return dmat
+            return dmat, dmat.columns
         elif return_type == "patsy":
             raise ValueError("return type 'patsy' not supported for input (dmat is not None)")
         else:
@@ -134,7 +135,7 @@ def preview_coef_names(
 
 
 def constraint_system_from_star(
-        dmat: Union[None, np.ndarray] = None,
+        dmat: Union[None, patsy.design_info.DesignMatrix] = None,
         sample_description: Union[None, pd.DataFrame] = None,
         formula: Union[None, str] = None,
         as_categorical: Union[bool, list] = True,
@@ -202,7 +203,7 @@ def constraint_system_from_star(
         raise ValueError("supply either sample_description or dmat")
 
     if dmat is None and not isinstance(constraints, dict):
-       dmat = design_matrix(
+       dmat, coef_names = design_matrix(
             sample_description=sample_description,
             formula=formula,
             as_categorical=as_categorical,
@@ -213,39 +214,49 @@ def constraint_system_from_star(
         raise ValueError("dmat was supplied even though constraints were given as dict")
 
     if isinstance(constraints, dict):
-        dmat, cmat = constraint_matrix_from_dict(
+        dmat, coef_names, cmat, term_names = constraint_matrix_from_dict(
             sample_description=sample_description,
             formula=formula,
             as_categorical=as_categorical,
             constraints=constraints,
-            return_type="dataframe"
+            return_type="patsy"
         )
     elif isinstance(constraints, tuple) or isinstance(constraints, list):
-        cmat = constraint_matrix_from_string(
+        cmat, coef_names = constraint_matrix_from_string(
             dmat=dmat,
+            coef_names=dmat.design_info.column_names,
             constraints=constraints
         )
+        term_names = None  # not supported yet.
     elif isinstance(constraints, np.ndarray):
         cmat = constraints
+        term_names = None
     elif constraints is None:
         cmat = None
+        term_names = None
     else:
         raise ValueError("constraint format %s not recognized" % type(constraints))
 
-    return dmat, cmat
+    return dmat, coef_names, cmat, term_names
 
 
 def constraint_matrix_from_dict(
         sample_description: pd.DataFrame,
         formula: str,
         as_categorical: Union[bool, list] = True,
         constraints: dict = {},
-        return_type: str = "dataframe"
+        return_type: str = "patsy"
 ) -> Tuple:
     """
     Create a design matrix from some sample description and a constraint matrix
     based on factor encoding of constrained parameter sets.
 
+    Note that we build a dataframe instead of a pasty.DesignMatrix here if constraints are used.
+    This is done because we were not able to build a patsy.DesignMatrix of the constrained form
+    required in this context. In those cases in which the return type cannot be patsy, we encourage the
+    use of the returned term_names to perform term-wise slicing which is not supported by other
+    design matrix return types.
+
     :param sample_description: pandas.DataFrame of length "num_observations" containing explanatory variables as columns
     :param formula: model formula as string, describing the relations of the explanatory variables.
 
@@ -270,7 +281,9 @@ def constraint_matrix_from_dict(
 
         Can only group by non-constrained effects right now, use constraint_matrix_from_string
         for other cases.
-    :return: a model design matrix
+    :return:
+        - model design matrix
+        - term_names to allow slicing by factor if return type cannot be patsy.DesignMatrix
     """
     assert len(constraints) > 0, "supply constraints"
     sample_description: pd.DataFrame = sample_description.copy()
@@ -287,10 +300,11 @@ def constraint_matrix_from_dict(
     # absorption of the first level of each factor for each constrained factor onto the
     # core matrix.
     formula_unconstrained = formula.split("+")
-    formula_unconstrained = [x for x in formula_unconstrained if x not in constraints.keys()]
+    formula_unconstrained = [x for x in formula_unconstrained if x.strip(" ") not in constraints.keys()]
     formula_unconstrained = "+".join(formula_unconstrained)
     dmat = patsy.dmatrix(formula_unconstrained, sample_description)
     coef_names = dmat.design_info.column_names
+    term_names = dmat.design_info.term_names
 
     constraints_ls = string_constraints_from_dict(
         sample_description=sample_description,
@@ -301,6 +315,7 @@ def constraint_matrix_from_dict(
         dmat_constrained_temp = patsy.highlevel.dmatrix("0+" + x, sample_description)
         dmat = np.hstack([dmat, dmat_constrained_temp])
         coef_names.extend(dmat_constrained_temp.design_info.column_names)
+        term_names.extend(dmat_constrained_temp.design_info.term_names)
 
     # Build constraint matrix.
     constraints_ar = constraint_matrix_from_string(
@@ -312,8 +327,7 @@ def constraint_matrix_from_dict(
     # Format return type
     if return_type == "dataframe":
         dmat = pd.DataFrame(dmat, columns=coef_names)
-
-    return dmat, constraints_ar
+    return dmat, coef_names, constraints_ar, term_names
 
 
 def string_constraints_from_dict(
@@ -407,7 +421,7 @@ def constraint_matrix_from_string(
             constraint_mat[i, idx_unconstr_i] = 1
 
     # Test unconstrained subset design matrix for being full rank before returning constraints:
-    dmat_var =dmat[:, idx_unconstr]
+    dmat_var = dmat[:, idx_unconstr]
     if np.linalg.matrix_rank(dmat_var) != np.linalg.matrix_rank(dmat_var.T):
         logging.getLogger("batchglm").error("constrained design matrix is not full rank")
 

batchglm/models/base_glm/input.py

@@ -88,7 +88,7 @@ def __init__(
             design_matrix=design_loc,
             param_names=design_loc_names
         )
-        design_scale, design_scale_names =  parse_design(
+        design_scale, design_scale_names = parse_design(
             design_matrix=design_scale,
             param_names=design_scale_names
         )

batchglm/models/base_glm/simulator.py

@@ -3,7 +3,7 @@
 import numpy as np
 import pandas
 import patsy
-from typing import Union
+from typing import Union, Tuple
 
 from .model import _ModelGLM
 from .external import _SimulatorBase
@@ -14,7 +14,7 @@ def generate_sample_description(
         num_conditions: int = 2,
         num_batches: int = 4,
         shuffle_assignments=False
-) -> np.ndarray:
+) -> Tuple[patsy.DesignMatrix, pandas.DataFrame]:
     """ Build a sample description.
 
     :param num_observations: Number of observations to simulate.
@@ -45,7 +45,7 @@ def generate_sample_description(
             observations=np.random.permutation(sample_description.observations.values)
         )
 
-    return np.asarray(patsy.dmatrix("~1+condition+batch", sample_description)), sample_description
+    return patsy.dmatrix("~1+condition+batch", sample_description), sample_description
 
 
 class _SimulatorGLM(_SimulatorBase, metaclass=abc.ABCMeta):
@@ -73,6 +73,7 @@ def __init__(
         self.sample_description = None
         self.sim_a_var = None
         self.sim_b_var = None
+        self._size_factors = None
 
     def generate_sample_description(
             self,
@@ -139,6 +140,10 @@ def _generate_params(
             rand_fn_scale((self.sim_design_scale.shape[1], self.nfeatures))
         ], axis=0)
 
+    @property
+    def size_factors(self):
+        return self._size_factors
+
     @property
     def a_var(self):
         return self.sim_a_var

batchglm/models/base_glm/utils.py

@@ -37,8 +37,8 @@ def parse_design(
         assert False
 
     if param_names is not None:
-        if params is not None:
-            assert len(param_names) == len(params)
+        if params is None:
+            assert len(param_names) == dmat.shape[1]
             params = param_names
 
     return dmat, params
@@ -62,8 +62,11 @@ def parse_constraints(
         constraints = np.identity(n=dmat.shape[1])
         constraint_params = dmat_par_names
     else:
-        # Cannot use given parameter names if constraint matrix is not identity: Make up new ones.
-        constraint_params = ["var_"+str(x) for x in range(constraints.shape[1])]
+        # Cannot use all parameter names if constraint matrix is not identity: Make up new ones.
+        # Use variable names that can be mapped (unconstrained).
+        constraint_params = ["var_"+str(i) if np.sum(constraints[:, i] != 0) > 1
+                             else dmat_par_names[np.where(constraints[:, i] != 0)[0][0]]
+                             for i in range(constraints.shape[1])]
         assert constraints.shape[0] == dmat.shape[1], "constraint dimension mismatch"
 
     if constraint_par_names is not None:

batchglm/models/glm_nb/simulator.py

@@ -15,6 +15,10 @@ def __init__(
             num_observations=1000,
             num_features=100
     ):
+        Model.__init__(
+            self=self,
+            input_data=None
+        )
         _SimulatorGLM.__init__(
             self=self,
             model=None,
@@ -53,3 +57,4 @@ def generate_data(self):
             design_loc_names=None,
             design_scale_names=None
         )
+