Non-parametric for two-way ANOVA (3x3) My dependent variable is continuous, non-normal (skewed left according to Shapiro-Wilk test). I have two independent variables (treatment group by colour, food type). There are 3 levels within each independent variable. The number of observations for each independent variables are not equal.
I have looked up non-parametric tests such as Friedman's Test and Scheirer-Ray-Hare Test, neither of which seem suitable (due to unequal number of observations).
Are there alternative tests that anyone could suggest? I am using SAS.
 A: What question are you trying to answer?
If you want an overall test of anything going on, The null is that both main effects and the interaction are all 0, then you can replace all the data points with their ranks and just do a regular ANOVA to compare against an intercept/grand mean only model.  This is basically how many of the non-parametric tests work, using the ranks transforms the data to a uniform distribution (under the null) and you get a good approximation by treating it as normal (the Central Limit Theorem applies for the uniform for sample sizes above about 5 or 6).
For other questions you could use permutation tests.  If you want to test one of the main effects and the interaction together (but allow the other main effect to be non-zero) then you can permute the predictor being tested.  I you want to test the interaction while allowing both main effects to be non-zero then you can fit the reduced model of main-effects only and compute the fitted values and residuals, then randomly permute the residuals and add the permuted residuals back to the fitted values and fit the full anova model including the interaction.  Repeat this a bunch of times to get the null distribution for the size of the interaction effect to compare with the size of the interaction effect from the original data. 
There may be existing SAS code for doing things like this, I have seen some basic tutorials on using SAS for bootstrap and permutation tests (the quickest way seems to be using the data step to create all the datasets in one big table, then using by processing to do the analyses).  Personally I use R for this type of thing so can't be of more help in using SAS.

Edit
Here is an example using R code:  
> fit1 <- aov(breaks ~ wool*tension, data=warpbreaks)
> summary(fit1)
             Df Sum Sq Mean Sq F value   Pr(>F)    
wool          1    451   450.7   3.765 0.058213 .  
tension       2   2034  1017.1   8.498 0.000693 ***
wool:tension  2   1003   501.4   4.189 0.021044 *  
Residuals    48   5745   119.7                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> fit2 <- aov(breaks ~ wool + tension, data=warpbreaks)
> 
> tmpfun <- function() {
+   new.df <- data.frame(breaks = fitted(fit2) + sample(resid(fit2)),
+                        wool = warpbreaks$wool,
+                        tension = warpbreaks$tension)
+   fitnew <- aov(breaks ~ wool*tension, data=new.df)
+   fitnew2 <- update(fitnew, .~ wool + tension)
+   c(coef(fitnew), F=anova(fitnew2,fitnew)[2,5])
+ }
> 
> out <- replicate(10000, tmpfun())
> 
> # based on only the interaction coefficients
> mean(out[5,] >= coef(fit1)[5])
[1] 0.002
> mean(out[6,] >= coef(fit1)[6])
[1] 0.0796
> 
> # based on F statistic from full-reduced model
> mean(out[7,] >= anova(fit2,fit1)[2,5])
[1] 0.022

A: +1 to @Greg Snow. In keeping with his non-parametric use ranks strategy, you could use ordinal logistic regression.  This will allow you to fit a model with multiple factors and interactions between them.  The relevant SAS documentation is here.  There is a tutorial on how to do this in SAS at UCLA's excellent stats help website here, and another tutorial from Indiana University here.  
