stat_smooth: parameter specification for formula

701 views
Skip to first unread message

Keren Raiter

unread,
Jun 29, 2015, 6:39:41 PM6/29/15
to ggp...@googlegroups.com
Hi,

I'd like to add a smoother to my plot that is based on a hyperbolic function, y = ax/(b + x).

I'm happy to either let ggplot calculate the parameter values of 'a' and 'b', or to input them myself (I've calculated them using a nls fit function). 
I've tried 
p + stat_smooth(method="lm", formula= y~ a*x/(b + x), aes(x=mprojn, y=os$ft.fit)),
p + stat_smooth(method="lm", formula= y~ A1*x/(B1 + x), aes(x=mprojn, y=os$ft.fit)) (where I have A1 and B1 defined), 

and a few other combinations, all of which returned errors. 

I see in other examples of defining the formula of smoothers, such as:
# p + stat_smooth(method = "gam", formula = y ~ s(x), size = 1)
# p + stat_smooth(method = "lm", formula = y ~ x + I(x^2), size = 1)

That there do seem to be some letters that are accepted as parameters for certain formulas; how can I figure out what letters or otherwise how to tell ggplot what I want the formula of the smoother to be?

Thanks,
Keren

Dennis Murphy

unread,
Jun 30, 2015, 4:09:50 AM6/30/15
to Keren Raiter, ggplot2
Hi:

ggplot2 will neither fit nonlinear models for you on the fly, nor are
you allowed to pass unspecified parameters into a formula and expect
ggplot2 to estimate them. The formula specification in ggplot2 is
limited to specific types of models of one variable, although that is
not clear from the on-line documentation. In this case, you have to
fit the function outside of ggplot2, create a data frame with the
inputs necessary for the graph and then pass that data frame into
ggplot().

The model you cited is one parameterization of the Michaelis-Menten
model. Since you didn't provide a reproducible example, we can supply
one with the Puromycin data, one of the data sets built into R.

If all you want is the smoother, here is one way:

# Select the treatment group from the Puromycin data
PurTrt <- subset(Puromycin, state == "treated", select = -3)

# Fit a self-starting Michaelis-Menten model
mmod <- nls(rate ~ SSmicmen(conc, Vm, K), data = PurTrt)

# Values of x at which to predict
xpred <- seq(0.01, 1.1, by = 0.01)

mmpred <- data.frame(x = xpred,
y = predict(mmod, newdata = data.frame(conc = xpred)))

ggplot(PurTrt, aes(x = conc, y = rate)) +
geom_point() +
geom_line(data = mmpred, aes(x = x, y = y),
size = 1, color = "darkorange")


Since you didn't specify se = FALSE in your geom_smooth() call, I
assume you want 95% confidence bands around the fitted curve. To do
that requires more work using the bootstrap method. (In general, there
are few analytic formulas to compute confidence bands around fitted
models in nonlinear regression, so the bootstrap is commonly applied
for this purpose.)

Given the PurTrt data and the model mmod from above, this is a bit
tedious, but it gets the job done. It helps to install the nlstools
package for its nlsBoot() function.

# install.packages("nlstools")
library(nlstools)

# Fit 999 bootstrap iterations of the MM model where the
# mean-centered residuals are bootstrapped:
bootcoef <- nlsBoot(mmod)$coefboot

# Concentrations at which we want bootstrap predictions
xpred <- seq(0.01, 1.1, by = 0.01)

library(plyr)

# Function to compute predicted values from x
# and estimated coefficients
MMpred <- function(coef, x) coef[1] * x /(coef[2] + x)

# Apply above function to each row of bootcoef, yielding a
# matrix of bootstrap predicted values - each column corresponds
# to an x-value from xpred
bootpreds <- aaply(bootcoef, 1, function(b) MMpred(b, xpred))

# Find the quantile-based 95% CIs at each x in xpred
bootci <- aaply(bootpreds, 2, function(x) quantile(x, c(0.025, 0.975)))

# Put together xpred, the predictions from the original model for xpred
# and the bootstrap CIs at each value of xpred
bootlims <- data.frame(x = xpred,
yhat = predict(mmod, newdata = data.frame(conc = xpred)),
lcl95 = bootci[, 1], ucl95 = bootci[, 2])

# Finally....
ggplot(bootlims, aes(x = x, y= yhat)) +
theme_bw() +
geom_ribbon(aes(ymin = lcl95, ymax = ucl95), fill = "grey90") +
geom_point(data = PurTrt, aes(x = conc, y = rate), size = 3) +
geom_line(size = 1, color = "darkorange")

If you want a different type of bootstrap CI, then I recommend using
the boot package for this purpose.

HTH,
Dennis
> --
> --
> You received this message because you are subscribed to the ggplot2 mailing
> list.
> Please provide a reproducible example:
> https://github.com/hadley/devtools/wiki/Reproducibility
>
> To post: email ggp...@googlegroups.com
> To unsubscribe: email ggplot2+u...@googlegroups.com
> More options: http://groups.google.com/group/ggplot2
>
> ---
> You received this message because you are subscribed to the Google Groups
> "ggplot2" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to ggplot2+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

Tom

unread,
Jun 30, 2015, 6:36:27 PM6/30/15
to ggp...@googlegroups.com
I was about to ask a very similar question. Dennis, your answer helped immensely for my needs; thank you!

- Tom
Reply all
Reply to author
Forward
0 new messages