LOESS and MA normalization in R? Attempting to do loess on two variables x and y in R using MA normalization(see MA-plot; see also Bland-Altman or Tukey mean-difference plot) like this:
> x = rnorm(100) + 5
> y = x + 0.6 + rnorm(100)*0.8
> m = log(x/y)
> a = 0.5*log(x*y)

I want to normalize x and y in such a way that the average m is 0, as in standard MA normalization, and then back-calculate the correct x and y values. First running loess on MA:
> l = loess(m ~ a)

What is the way to get corrected m values then? Is this correct?
> mc <- predict(l, a)
# original MA plot
> plot(a,m)
# corrected MA plot
> plot(a,m-mc)

not clear to me what predict actually does in the case of loess objects and how it's different from using l$residuals in the object l returned by loess - can someone explain?
finally, how can I back calculate new x and y values based on this correction?
Attempt to back-calculate x and y: the corrected m is, 
new_m = m - mc

so x,y can be derived from the definition of m:
m = log(x/y) = log(x) - log(y)

therefore,
x = exp(new_m + log(y))
y = exp(-1*(new_m - log(x)))

but this is wrong; there's a missing scaling factor of a half and I'd like to see a derivation of where that comes from. Probably from the definition of A but I don't see why I can't get the same be rewriting x,y in terms of A
edit if someone can explain where the half correction factor comes in when back calculating x -- i.e. why the formula x = exp(new_m + log(y)) is wrong -- i'd appreciate it.  
 A: first of all try and be consistent with logs (choose either log or log2 but not both with your m and a calculations)
The steps that you have outlined do look correct. predict() will take a model (generated by loess(m ~ a)) and give you the 'corrected' m (mc) for every value a. You can view the loess line by doing
#set.seed(12345)
x = rnorm(100) + 5
y = x + 0.6 + rnorm(100)*0.8

m = log(x/y)
a = 0.5*log(x*y)

l = loess(m ~ a) 
mc <- predict(l, a)

plot(a, m, ylim=c(-0.5, 0.5)); 
#make sure things are ordered for the lines() plot
lines(a[order(m)], mc[order(m)]) #this is the loess line through the points m
dev.new()

#if you want your m centered around loess line
plot(a, m-mc, ylim=c(-0.5, 0.5)); dev.new()

#rescaled values
x2 = exp(log(x) - mc/2)
y2 = exp(log(y) + mc/2)
m2 = log(x2/y2)
a2 = 0.5*log(x2*y2)

plot(a2, m2, ylim=c(-0.5, 0.5)) #same as second plot

In this case, l$residual and m-mc are the same values because you are interested in predicting everything that you fit. However, you could do mc2 = predict(l, a2) where a2 might be a superset or something else you want to fit a similar transform to, m2-mc2 will be different from l$residual. This is useful if you are only interested in using a subset of your data to fit to the model but interested in adjusting all of your data.
Also see this loess guide here
As for the back calculation, I think you just adjust X and Y by mc
x2 = exp(log(x) - mc/2)
y2 = exp(log(y) + mc/2)

If you are working on gene expression, the affy package has all of this implemented under normalize.loess()
A: yingw's answer is good and mostly correct - there is a minor mistake in the code, though. To correctly plot the loess curve, line 13 should read like the following (ordered by a rather than m):
lines(a[order(a)], mc[order(a)]) #this is the loess line through the points m


Sorry I didn't post this as a comment - not enought reputation, yet...
