# 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. • Generally, if I see the abbreviation "MA", I would expect this is being discussed. Could you explain (in your question) what MA stands for the first time you mention it, please? Mar 20, 2013 at 22:02 • @Gleb_b: sorry for confusion, stands for this: en.wikipedia.org/wiki/MA_plot Mar 20, 2013 at 22:05 • Thanks. Never heard of it before. Sounds a lot like a Tukey mean-difference plot (/sum-difference plot). Mar 20, 2013 at 22:11 • Please do not cross post on Stackoverflow – mnel Mar 20, 2013 at 22:36 • @Glen_b: The MA-plot is an application of a Bland–Altman plot for visual representation of microarray gene expression data. The Bland-Altman plot is known as Tukey's mean-difference plot in other fields. Jan 13, 2014 at 18:32 ## 2 Answers 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() • thank you! final question: can you explain the difference between m+predict(l,a) and m+l$residuals? shouldn't the corrected m just be m plus its residuals? Mar 20, 2013 at 22:58
• I tried to explain this in the 3rd paragraph. The residual is the distance from the point to the model. When you run loess command with m and a, the l\$residual will be the distance between m and the model. When you are predicting a, what the function will give back is: Given a model l and these points a, how much would I adjust to get to the line from model. Maybe its easier to think of it this way, if your a was always integers and you wanted to predict 5.5, you aren't going to find 5.5 in your a but you can figure it out using the model. Mar 20, 2013 at 23:16
• ok so just to make it explicit could you say what the notation would be for expressing the corrected m value as a function of predict's output? is it: new_m = g(a) + error? Mar 20, 2013 at 23:27
• I think I am confusing myself and confusing you, the best thing to look at is to plot the lines(a, mc) in in my answer and I hope that will make things easier to see. Mar 20, 2013 at 23:53
• it all pretty much makes sense except getting the scaling back to x and y both corrected. there's some sort of 0.5 scaling factor that I am missing. I don't think the above definition of X_rescale and Y_rescale is correct. Mar 21, 2013 at 0:09

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...