We propose an intrinsic multi-factorial model for data on Riemannian manifolds as typically occur in the statistical analysis of shape. Due to the lack of a linear structure, linear models cannot be defined in general; to date only one-way MANOVA is available. For a general multi-factorial model, we assume that variation not explained by the model is concentrated near elements defining the effects. By determining the asymptotic distributions of respective sample covariances under parallel transport, we show that they can be compared by standard MANOVA. Often in applications, manifolds are only implicitly given as quotients where the bottom space parallel transport can be expressed through a differential equation. For Kendall's space of planar shapes, we provide an explicit solution. We illustrate our method by an intrinsic two-way MANOVA for a set of leaf shapes. While biologists can identify genotype effects by sight, we can detect height effects that are otherwise not identifiable.