A linear activation function (or none at all) should only be used when the relation between input and output is linear. Why doesn't the same rule apply for other activation functions? For example, why doesn't sigmoid only work when the relation between input and output is "of sigmoid shape"?
It is correct to say that a sigmoid activation function would only work well as a model if the desired output is close to the sigmoid function applied to the input. This is a trivial fact that applies to a single layer perceptron. This is true for for the single layer case for any activation function, also a trivial fact.
When the layer number is between one and infinity (two or more), the theory bifurcates. The identity function becomes a special case: Any number of layers that conform to a first degree polynomial, $ax + b$, can be replaced with a single identity layers. The alternative case, where there are multiple layers that functionally do not conform to a first degree polynomial, $ax + b$, cannot be replaced by a single layer of some equally simple function. The complexity increases geometrically, which is the entire point of multilayer perceptrons.
Under particular constraints, the multilayer perceptron can produce a wide variety of functional behaviors that do not resemble the activation functions of the layers.
For instance, a properly trained network using sigmoid activation functions, with sufficient layer depth and sufficient massive allocation of computing resources, could theoretically approximate the topography of the Himalayas.