ggplot exponential smooth with tuning parameter inside exp

user_13487249 0
Viewed 543

ggplot provides various "smoothing methods" or "formulas" that determine the form of the trend line. However it is unclear to me how the parameters of the formula are specified and how I can get the exponential formula to fit my data. In other words how to tell ggplot that it should fit the parameter inside the exp.

df <- data.frame(x = c(65,53,41,32,28,26,23,19))
df$y <- c(4,3,2,8,12,8,20,15)

   x  y
1 65  4
2 53  3
3 41  2
4 32  8
5 28 12
6 26  8
7 23 20
8 19 15
p <- ggplot(data = df, aes(x = x, y = y)) +
  geom_smooth(method = "glm", se=FALSE, color="black", formula = y ~ exp(x)) +
  geom_point()

p

Problematic fit:

enter image description here

However if the parameter inside the exponential is fit then the form of the trend line becomes reasonable:

p <- ggplot(data = df, aes(x = x, y = y)) +
  geom_smooth(method = "glm", se=FALSE, color="black", formula = y ~ exp(-0.09 * x)) +
  geom_point()

p

enter image description here

ggplot2 r
2 Answers

Here is an approach with method nls instead of glm.

You can pass additional parameters to nls with a list supplied in method.args =. Here we define starting values for the a and r coefficients to be fit from. 

library(ggplot2)
ggplot(data = df, aes(x = x, y = y)) +
  geom_smooth(method = "nls", se = FALSE,
              formula = y ~ a * exp(r * x),
              method.args = list(start = c(a = 10, r = -0.01)),
              color = "black") +
  geom_point()

enter image description here

As discussed in the comments, the best way to get the coefficients on the graph is by fitting the model outside the ggplot call. 

model.coeff <- coef(nls( y ~ a * exp(r * x), data = df, start = c(a = 50, r = -0.04)))

ggplot(data = df, aes(x = x, y = y)) +
  geom_smooth(method = "nls", se = FALSE,
              formula = y ~ a * exp(r * x),
              method.args = list(start = c(a = 50, r = -0.04)),
              color = "black") +
  geom_point() + 
  geom_text(x = 40, y = 15,
            label = as.expression(substitute(italic(y) == a %.% italic(e)^(r %.% x),
                                             list(a = format(unname(model.coeff["a"]),digits = 3),
                                                  r = format(unname(model.coeff["r"]),digits = 3)))),
            parse = TRUE)

enter image description here

user_130953260

Firstly, to pass additional parameters to the function passed to the method param of geom_smooth, you can pass a list of named parameters to method.args.

Secondly, the problem you're seeing is that glm is placing the coefficient in front of the whole term: y ~ coef * exp(x) instead of inside: y ~ exp(coef * x) like you want. You could use optimization to solve the latter outside of glm, but you can fit it into the GLM paradigm by a transformation: a log link. This works because it's like taking the equation you want to fit, y = exp(coef * x), and taking the log of both sides, so you're now fitting log(y) = coef * x, which is equivalent to what you want to fit and works with the GLM paradigm. (This ignores the intercept. It also ends up in transformed link units, but it's easy enough to convert back if you like.)

You can run this outside of ggplot to see what the models look like: 

df <- data.frame(
    x = c(65,53,41,32,28,26,23,19), 
    y <- c(4,3,2,8,12,8,20,15)
)

bad_model <- glm(y ~ exp(x), family = gaussian(link = 'identity'), data = df)
good_model <- glm(y ~ x, family = gaussian(link = 'log'), data = df)

# this is bad
summary(bad_model)
#> 
#> Call:
#> glm(formula = y ~ exp(x), family = gaussian(link = "identity"), 
#>     data = df)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -7.7143  -2.9643  -0.8571   3.0357  10.2857  
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)   
#> (Intercept)  9.714e+00  2.437e+00   3.986  0.00723 **
#> exp(x)      -3.372e-28  4.067e-28  -0.829  0.43881   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for gaussian family taken to be 41.57135)
#> 
#>     Null deviance: 278.00  on 7  degrees of freedom
#> Residual deviance: 249.43  on 6  degrees of freedom
#> AIC: 56.221
#> 
#> Number of Fisher Scoring iterations: 2

# this is better
summary(good_model)
#> 
#> Call:
#> glm(formula = y ~ x, family = gaussian(link = "log"), data = df)
#> 
#> Deviance Residuals: 
#>    Min      1Q  Median      3Q     Max  
#> -3.745  -2.600   0.046   1.812   6.080  
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  3.93579    0.51361   7.663 0.000258 ***
#> x           -0.05663    0.02054  -2.757 0.032997 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for gaussian family taken to be 12.6906)
#> 
#>     Null deviance: 278.000  on 7  degrees of freedom
#> Residual deviance:  76.143  on 6  degrees of freedom
#> AIC: 46.728
#> 
#> Number of Fisher Scoring iterations: 6

From here, you can reproduce what geom_smooth is going to do: make a sequence of x values across the domain and use the predictions as the y values for the line:

# new data is a sequence across the domain of the model
new_df <- data.frame(x = seq(min(df$x), max(df$x), length = 501))

# `type = 'response'` because we want values for y back in y units
new_df$bad_pred <- predict(bad_model, newdata = new_df, type = 'response')
new_df$good_pred <- predict(good_model, newdata = new_df, type = 'response')

library(tidyr)
library(ggplot2)

new_df %>% 
    # reshape to long form for ggplot
    gather(model, y, contains('pred')) %>%
    ggplot(aes(x, y)) + 
    geom_line(aes(color = model)) + 
    # plot original points on top
    geom_point(data = df)

glm predictions plot

Of course, it's a lot easier to let ggplot handle all that for you:

ggplot(df, aes(x, y)) + 
    geom_smooth(
        method = 'glm', 
        formula = y ~ x, 
        method.args = list(family = gaussian(link = 'log'))
    ) + 
    geom_point()

glm smooth with log link

user_44970500
Related