A slightly more satisfying procedure:
Produce a cartogram on a 'diffusion-based' method (like scapetoad).
Create a hexagonal-layer [Plugins > Qmarxan > Create Planning Grid] (use appropriate/ desired number of hexagons). edit: Depending on the purpose of the map, Qmarxans grid could be of dubious quality. The mmqgis plugin produces much better results (polygons align perfectly).
Join the the hexagonal-layer with the cartogram [Vector > Data Managment Tools > Join attributes by location]
Now the new hexagonal-layer has all the properties/attributes from the cartogram, and it is possible to access each subunit individually.
For better shape-preservance of subunits some manual tweeking is possible: Select the hexagon(s) and change the attribute(s) directly in the attribute-table. This step works best, if the layer is semi-transparent and the original diffusion-based cartogram is shown below.
Merging the hexagons of a subunit to one unit: 1. Select all hexagons of a subunit (if many: use attribute-table for selection). 2. 'toggle edit' > 'merge selected features'.
Mainly shortcoming b) remains - but a) is solved and also to some degree c) isn't a big problem anymore, b/c this method doesn't take that much time.
Edit: Step 4. (and colored afterwards) leads to this result (transparent hexagon-layer on top, 'diffusion-based' cartogram below):
The bright-green hexagons clearly take one row to much...
Any ideas how it is possible to achieve a better/ a custom fitting of the "join attributes"? - I am looking for a method to give the algorithm some kind of threshold (like: >50% of area).