Statistics on LOESS smoothing I'm trying to make sense of a data set which contains thousands of independent measurements of intensity in the form of a scatterplot. These measurements are dependent on two major variables: treatment and genotype, so I represented them on a grid. And because I want to see the trend of my data, I add a LOESS regression. I'm not an R pro so it's just something very simple:
ggplot(data, aes(Distance, Response, color = Treatment))  + 
  geom_point(alpha = 1/8, show.legend = FALSE) + 
  stat_smooth(method = "loess", formula = "y ~ x", show.legend = FALSE) + 
  labs(x = "Distance", y = "Response") + 
  theme_classic () +
  facet_grid(Genotype ~ Treatment)


So, after adding the LOESS curve, it seems that the treatments have a different effect in each genotype. And it is much clearer to see it with the curve than with the scatterplot. Therefore, I would like to know the statistical significance of that effect, ideally comparing the curves instead of the data points. If I do for example an ANOVA, I see that both the effects of "Treatment" and "Genotype" are statistically significant. But I would like to see the effect of Treatment2, Treatment3 and Treatment4 in each genotype in comparison with Treatment1 (mock treatment). And also to examine the differences between genotypes under the same treatment. I do not know if I'm explaining it correctly but it would be ideally something like this.
The (huge) problem is that I have no idea if that is possible, and what to compare. I have been all the day reading documentation about LOESS regression, ggplot, ggpubr ... But nothing seems to do what I want. The only thing I can think of is to obtain the fitting values for the LOESS curves, integrate them and compare that value as a measure of how "different" they are, but before doing that I would love to know if there is a more or less direct way to obtain my desired graph. Any idea or advice??
Thank you very much in advance!
EDIT: Thanks for migrating this here. I'm sorry I didn't posted it on the right place since the begging.
 A: I'm sure there are better ways to do it but maybe this is an idea to get started.
Divide your "Distance" variable in bins, each containing a reasonable number of datapoints. At a glance, bins from -2 to +2 in 0.5 step could do [i.e. seq(-2, 2, by= 0.5)].
Your data table now should have columns: "response", "bin", "treatment", "genotype".
Then fit the ANOVA model with interaction between bin, treatment, genotype:
aov1 <- aov(data= dat, response ~ as.factor(bin) * treatment * genotype)
summary.lm(aov1)

This should pick up that bins around 0 in treatment 1, genotype 4 are different from the baseline.
You can then check for comparisons between bins, treatments and genotypes with:
TukeyHSD(aov1)


EDIT after whuber's comment: A simple improvement to the above solution may be to use the absolute distance from 0 for the binning, e.g. use abs(seq(-2, 2, by= 0.5)), since bins to the left and right of 0 are assumed to be equivalent with respect to the response. This will half the number of bins and increase power. Bins could either be treated as a nominal variable or as an ordinal variable to reflect that there is an increasing trend moving towards 0.
