Nonlinear subdivision schemes arise from, among other applications, nonlinear multiscale signal processing and shape preserving interpolation. For the univariate homogeneous subdivision operator $S:\ell(\bZ) \goto \ell(\bZ)$ we establish a set of commutation/recurrence relations which can be used to analyze the asymptotic decay rate of $\|\Delta^r S^j m\|_{\ell^\infty}$, $j=1,2,\ldots,$ the latter in turn determines the convergence and H\”older regularity of $S$. We apply these results to prove that the critical H\”older regularity exponent of a nonlinear subdivision scheme based on median-interpolation is equal to that of an approximating linear subdivision scheme, resolving a conjecture by Donoho and Yu. We also consider a family of nonlinear but affine invariant subdivision operators based on interpolation-imputation of $p$-mean (of which median corresponds to the special case $p=1$) as well as general continuous $M$-estimators. We propose a linearization principle which, when applied to $p$-mean subdivision operators, leads to a family of linear subdivision schemes. Numerical evidence indicates that in at least many cases the critical smoothness of a $p$-mean subdivision scheme is the same as that of the corresponding linear scheme. This suggests a more coherent view of the result obtained in this paper.