def linear_model_fit(data_images, mask_images, design_matrix, vector):
    """
    Helper function for group data analysis using arbitrary design matrix
    """

    # Prepare arrays
    data, vardata, xyz, mask = prepare_arrays(data_images, None, mask_images)

    # Create glm instance
    G = glm(data, design_matrix)

    # Compute requested contrast
    c = G.contrast(vector)

    # Compute z-map image
    zmap = np.zeros(data_images[0].shape).squeeze()
    zmap[list(xyz)] = c.zscore()
    zimg = Image(zmap, get_affine(data_images[0]))

    return zimg