Estoy analizando datos de un diseño factorial completo 2X2 de la siguiente forma:

    aovN <-aov(y ~ x + z + x:z + Bloque , data = A)
    anova(aovN)
    summary(aovN)
    aovNb<-aov(y~Tratamiento+Bloque, data = A)
    anova(aovNb)
    summary(aovNb)

Dándome resultados de este tipo:

    Response: y
                 Df  Sum Sq Mean Sq F value  Pr(>F)  
    Tratamiento  8  492.67   61.58  0.5283 0.81849  
    x            2  198.22   99.11  0.8502 0.44574  
    z            2   29.56   14.78  0.1268 0.88181  
    Bloque       2  758.22  379.11  3.2522 0.06528 .
    x:z          4  264.89   66.22  0.5681 0.68948  
    Residuals    16 1865.11  116.57 
Luego hacia una prueba de Tukey del paquete agricolae.

El problema es que un grupo de estos datos no es paramétrico y quisiera hacer un analisis que me pudiera mostrar los efectos individuales y la interacción entre ellos y luego desglosarlos con una prueba post-hoc como la de Tukey en R.
¿Qué procedimiento podría usar?

Dejo una muestra de los datos:

https://www.dropbox.com/s/kz87vgtkyj0u2ho/Duda.xlsx?dl=0