4
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – Glen_b -Reinstate Monica Mar 20 '13 at 22:02
  • $\begingroup$ @Gleb_b: sorry for confusion, stands for this: en.wikipedia.org/wiki/MA_plot $\endgroup$ – user248237 Mar 20 '13 at 22:05
  • $\begingroup$ Thanks. Never heard of it before. Sounds a lot like a Tukey mean-difference plot (/sum-difference plot). $\endgroup$ – Glen_b -Reinstate Monica Mar 20 '13 at 22:11
  • 1
    $\begingroup$ Please do not cross post on Stackoverflow $\endgroup$ – mnel Mar 20 '13 at 22:36
  • 1
    $\begingroup$ @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. $\endgroup$ – cashoes Jan 13 '14 at 18:32
4
$\begingroup$

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()

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – user248237 Mar 20 '13 at 22:58
  • $\begingroup$ 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. $\endgroup$ – yingw Mar 20 '13 at 23:16
  • $\begingroup$ 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? $\endgroup$ – user248237 Mar 20 '13 at 23:27
  • $\begingroup$ 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. $\endgroup$ – yingw Mar 20 '13 at 23:53
  • $\begingroup$ 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. $\endgroup$ – user248237 Mar 21 '13 at 0:09
1
$\begingroup$

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

enter image description here

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.