1
$\begingroup$

I have a set of 'activity' values for some enzyme assays I have been doing, that come out of some analysis I've been doing. The problem is, the data is fairly crap, and there aren't many points, but its for project for my MSc so I'm stuck with it. I appreciate that trying to regress 4 datapoints, and especially extrapolate from them is a statistical no-no, but it's what the PhD student who I'm emulating did (albeit with the help of our Uni's Biostats department - a resource I can't really utilise at this point) - and so I would like the results to be directly comparable.

DPConc DPActivity
0      100.000000
83     67.709971
166    6.296231
416    16.546593

I need to extrapolate from this an IC50 value, which is the concentration (x axis) value, for which 50% Activity is seen.

I want (have) to do this by fitting some kind of curve/spline to the points, but I'm at a bit of a loss. The best I've come up with so far is this:

require(splines)
library(splines)

plot(DPConc, DPActivity)
splineDP <- smooth.spline(DPConc, DPActivity, spar=0.45)
lines(splineDP)

Which for various values of spar gives me this graph (excuse its roughness - I'm really at the end of my tether!)

Dipeptide Activity vs. Concentration

It looks as though some kind of exponential decay is going on, but I'm not enough of a mathematician to be able to concoct the required function. I'd rather derive the function from the data, than make the data fit the function if possible (not least because I have other data that doesn't follow the same trend, though I suspect thats experimental error).

I've tried a couple of different approaches from around the web but my deadline is imminent and I really need the simplest solution. If the solution could be implemented with ggplot, that would be a bonus as my other plots so far have been made with it, the above was just a quick mockup.

$\endgroup$
  • 1
    $\begingroup$ This sounds less like a programing question and more like you need help on choosing an appropriate statistical method for small sample data. This type of question would probably be better on Cross Validated $\endgroup$ – MrFlick Sep 10 '14 at 17:58
  • $\begingroup$ A nuance here is that what you probably want is a (physics/chemistry/biology-based) model for the reaction or process, and that your curve fitting should be based on fitting model coefficients, not just randomly fitting using loess or splines. Once you have that model, you can find the 50% point easily enough. $\endgroup$ – Andy Clifton Sep 10 '14 at 18:16
  • $\begingroup$ I've seen similar questions posted here (SO) where people delved in to the statistics of it a bit. I suppose really the programming element I'm asking for is whether anyone here experienced in R knows of other spline-fitting tools that would spit out a polynomial/function from which I can calculate IC50. The 'scientific logic' of mass action (which is essentially what this is) would point to a linear or exponential relationship, and it's definitely non-linear. $\endgroup$ – Joe Healey Sep 10 '14 at 19:25
  • $\begingroup$ So, the question is then how we fit an exponential model to your data. You might want to look at e.g. stats.stackexchange.com/questions/11947/… $\endgroup$ – Andy Clifton Sep 10 '14 at 19:28
  • $\begingroup$ Thank you for that, it looks to be the right direction. I'd spent so long staring at the code that I think I couldn't really see the wood for the trees!! Do you know how/have any resources for how one would achieve the equivalent of that link in ggplot? $\endgroup$ – Joe Healey Sep 10 '14 at 19:53
6
$\begingroup$

According to the OP,

The 'scientific logic' of mass action (which is essentially what this is) would point to a linear or exponential relationship, and it's definitely non-linear

We can start by using ggplot to test out the best model. We'll do it purely visually, by plotting the two different options.

First up, we need to get the data in to a data frame. We'll use the usual method:

DP = data.frame(Activity = c(100,67.7,6.29,16.55),
                Conc = c(0, 83, 166, 416))

Now, we'll plot a linear fit to the data using ggplot's geom_smooth() function:

require(ggplot2)

p.dp.1 <- ggplot(data = DP,
               aes(x = Conc,
                   y = Activity)) + 
  ylab("Activity") + 
  xlab("Concentration") +
  geom_point() + 
  theme_bw(base_size = 8) +
  geom_smooth(method = "lm", 
              aes(colour = "Linear")) +
  scale_color_manual(name = "Fits",
                     breaks = c("Linear"),
                     values = c("blue"))

Which gives us this fit:

enter image description here

We can see the standard error (the grey-shaded area) is pretty large.

Because the OP suggested that this might be an exponential relationship, we'll now try adding a fit using an exponential. An exponential curve can be linearized by taking logs of both sides, and then doing a linear fit to the data, which would be very simple with ggplot. However, we have a problem; log(0) is -Inf, so we can't simply take the logs of both sides and do a linear fit. Instead, we have to use glm() to do the fit, and pass it through geom_smooth(). We can also take advantage of ggplot's ability to build up plots bit-by-bit, and just tack it on to the previous plot:

p.dp.2 <- p.dp.1 +
  geom_smooth(method = "glm",
              family = gaussian(link="log"), 
              aes(colour = "Exponential")) +
  scale_color_manual(name = "Fits",
                     breaks = c("Linear","Exponential"),
                     values = c("red","blue"))

Which gives us this:

enter image description here

Where we can see that the standard error (the grey bound) is much narrower than the linear fit. So, we have some justification of using the exponential, although it would still be nice to have a chemical model of activity versus concentration to base this on.

We now need the concentration at activity = 50. We can plot this using geom_abline():

p.dp.3 <- p.dp.2 +
  geom_abline(intercept = 50, slope =0,
              aes(colour = "Activity = 50")) +
  scale_color_manual(name = "Fits",
                     breaks = c("Linear","Exponential","Activity = 50"),
                     values = c("black","red","blue"))

which gives us this:

enter image description here

Our best estimate for 50% activity is the intersection of the black line of activity = 50 and the best-fit line using the exponential model. We can see that this point is slightly to the right of the data point at Conc = 83, and to the left of the line of concentration = 100.

To get the actual intersect, we create the model for activity as a function of concentration outside of ggplot, and approximate the fit using approx(). The model is set up like this:

dp.model <- glm(formula = Activity ~ Conc,
                family=gaussian(link="log"),
                data = DP)

Following this example we can get the concentration when activity is 50:

conc.est <- approx(x = dp.model$fitted.values,
      y = dp.model$data$Conc, xout=50)$y

and then the concentration is...

> conc.est
89.29575

So, using an exponential model for activity as a function of concentration, we can estimate that concentration = 89.3 when activity = 50. We can see from the figures that the standard error associated with this estimate is huge, though. Unfortunately I can't figure out how to get the error bands for this estimate.

$\endgroup$
  • 1
    $\begingroup$ According to this plot, the IC50 is somewhere between 50 and 400+, not 80-90. $\endgroup$ – whuber Sep 10 '14 at 21:07
  • $\begingroup$ By 80-90 I meant the approximate range of the intersect (I don't have the exact value). 50-400 would be correct if you believe the prediction bounds, but the bounds on this method are described elsewhere as being rough. I'm trying to find that link... $\endgroup$ – Andy Clifton Sep 10 '14 at 21:10
  • 2
    $\begingroup$ Right: "rough" in this instance includes (among other things) the fact the bounds do not evaluate lack of fit. Thus, to be realistic, one would need to widen these (fiducial) bounds even further. But at least they begin to indicate just how much uncertainty there is in the IC50 estimate in the absence of any information about the likely form of the response function. $\endgroup$ – whuber Sep 10 '14 at 21:13
  • $\begingroup$ Thank you Andy, that is perfect. One final question then in which case, is there any way you know of to get the exact IC50 value utilising the method you described here? @jarfa says below that stat_smooth doesn't have an equation quite like I'm getting at? $\endgroup$ – Joe Healey Sep 10 '14 at 21:23
  • 1
    $\begingroup$ I like this answer a lot better than the interpolating curve, but don't forget that you really need to believe that exponential model for this to be right. And that old saying that you can't make a silk purse from a sow's ear - in this case, the sow's ear has more going for it. $\endgroup$ – rvl Sep 10 '14 at 23:08
3
$\begingroup$

EDIT: I wrote this answer for the StackOverflow post, and was only attempting to answer the coding part of his question. Obviously a LOESS line isn't ideal for his needs.

You can use ggplot2 to give you smoothing line.

libary(ggplot2)
d = data.frame(DPConc= c(0, 83, 166, 416),
           DPActivity=c(100, 67.71, 6.3, 16.55))
ggplot(d, aes(x=DPConc, y=DPActivity)) + geom_point(size=3) + geom_smooth() + theme_minimal()

By default that uses a LOESS function to smooth the data points, but here you'll find information about other functions.

As others have said, go to CrossValidated for info on smoothing functions (although they may not be willing to help you interpolate 4 data points.)

enter image description here

$\endgroup$
  • $\begingroup$ The coding is good, but the physics is missing! For example, is there any reason for the locally very high gradient at ~90 conc? I know this won't impact the 50%$ level with this data set, but imagine if they were trying to get the value near the data points. $\endgroup$ – Andy Clifton Sep 10 '14 at 18:21
  • $\begingroup$ Agreed, I assumed that OP can figure out which function is best for his purposes. $\endgroup$ – jarfa Sep 10 '14 at 19:15
  • $\begingroup$ This looks to be the sort of thing I'd be after. Is there any way to 1) reduce the degree of fitting, since this is probably over-fit and 2) to pull out the equation of the line from geom_smooth? $\endgroup$ – Joe Healey Sep 10 '14 at 19:26
  • $\begingroup$ 1. Look at the ggplot help pages for stat_smooth and geom_s mooth. 2. A LOESS line doesn't really have an equation like you're thinking of. Try a regression with polynomial features (i.e. y ~ x1 + I(x1^2) + ...) $\endgroup$ – jarfa Sep 10 '14 at 19:27
  • 4
    $\begingroup$ But again, I should remind you that you should have a theory for this process if you're going to choose a given interpolation method. $\endgroup$ – jarfa Sep 10 '14 at 20:48

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.