def __init__(self, data, design_matrix, mask=None, formula=None,
                 model=def_model, method=None, niter=def_niter):

        # Convert input data and design into sequences
        if not hasattr(data, '__iter__'):
            data = [data]
        if not hasattr(design_matrix, '__iter__'):
            design_matrix = [design_matrix]

        # configure spatial properties
        # the 'sampling' direction is assumed to be the last
        # TODO: check that all input images have the same shape and
        # that it's consistent with the mask
        nomask = mask is None
        if nomask:
            self.xyz = None
            self.axis = len(data[0].shape) - 1
        else:
            self.xyz = np.where(mask.get_data() > 0)
            self.axis = 1

        self.spatial_shape = data[0].shape[0: -1]
        self.affine = get_affine(data[0])

        self.glm = []
        for i in range(len(data)):
            if not isinstance(design_matrix[i], np.ndarray):
                raise ValueError('Invalid design matrix')
            if nomask:
                Y = data[i].get_data()
            else:
                Y = data[i].get_data()[self.xyz]
            X = design_matrix[i]

            self.glm.append(glm(Y, X, axis=self.axis,
                                formula=formula, model=model,
                                method=method, niter=niter))