There really is no "right" answer, because it's a matter of how you view the problem. robustfit uses the residuals to calculate the weighting, but residual means the difference in the observed/predicted response variable. In your case, it looks like the points lie around a plane that is almost vertical ( normal(3) is small in all cases). This means that the residuals will be much greater when looking at z = f(x,y) than when looking at x = f(y,z) or y = f(x,z). Consequently, the weighting for a given point will be different in each formulation. In other words, different points are "outliers" in different formulations. To see this, try:
[p,s1] = robustfit([x y],z);
[p,s2] = robustfit([y z],x);
[p,s3] = robustfit([x z],y);
plotmatrix([s1.w,s2.w,s3.w])
This will plot the weightings against each other. I'm guessing you'll see reasonable correlation between s2 and s3, but not so great with s1.
This means that each formulation will give you a slightly different plane (because it's down-weighting different data points). I'd probably trust the last two, rather than the first.
