diff --git a/NAMESPACE b/NAMESPACE
index 8193127..eb7eb03 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -6,6 +6,7 @@ export(broken_stick_stevens)
 export(fake_crustaceans)
 export(infl_pt)
 export(regrans_fun)
+export(somerton_fun)
 export(two_line_logistic)
 export(two_line_stevens)
 import(ggplot2)
diff --git a/R/broken_stick.R b/R/broken_stick.R
index dd7ca6c..a0e7346 100644
--- a/R/broken_stick.R
+++ b/R/broken_stick.R
@@ -61,11 +61,11 @@ broken_stick <- function(dat,
   }
 
   if (is.null(lower)) {
-    lower <- stats::quantile(dat$xvar, 0.2)
+    lower <- stats::quantile(dat[[xvar]], 0.2)
   }
 
   if (is.null(upper)) {
-    upper <- stats::quantile(dat$xvar, 0.8)
+    upper <- stats::quantile(dat[[xvar]], 0.8)
   }
 
   if ("chngpt" %in% method | "all" %in% method) {
diff --git a/vignettes/classification.Rmd b/vignettes/classification.Rmd
index 0c9686b..be35704 100644
--- a/vignettes/classification.Rmd
+++ b/vignettes/classification.Rmd
@@ -83,6 +83,17 @@ ggplot() +
 
 hclust, kmeans, hkmeans, PAM, OPTICS, DBSCAN, HDBSCAN, Somerton method, mclust
 
+## Somerton method
+
+```{r}
+out_df <- somerton_fun(fc, xvar = "x", yvar = "y")[[1]]
+mod <- glm(data = out_df,
+           pred_mat_num ~ x,
+           family = binomial(link = "logit"))
+unname(-coef(mod)[1] / coef(mod)[2])
+```
+
+
 # Other classification methods
 
 ## Inflection point discriminant method
@@ -103,8 +114,8 @@ ggplot() +
 
 The line can then be used to classify each point as immature and mature, enabling logistic regression to be applied to estimate SM50:
 ```{r}
-fc_infl_pt <- fc %>% mutate(pred_mat = as.factor(if_else(y/x > disc, 1, 0)),
-                    correct = if_else(pred_mat == mature, TRUE, FALSE))
+fc_infl_pt <- fc %>% mutate(pred_mat = as.factor(if_else(y / x > disc, 1, 0)),
+                            correct = if_else(pred_mat == mature, TRUE, FALSE))
 ```
 
 ```{r}
@@ -132,10 +143,12 @@ if (rlang::is_installed("patchwork")) {
 ```
 
 ```{r}
-infl_pt_mod <- glm(pred_mat ~ x, family=binomial(link="logit"), data=fc_infl_pt)
+infl_pt_mod <- glm(pred_mat ~ x,
+                   family = binomial(link = "logit"),
+                   data = fc_infl_pt)
 
 # SM50 = -A/B
-unname(-coef(infl_pt_mod)[1]/coef(infl_pt_mod)[2])
+unname(-coef(infl_pt_mod)[1] / coef(infl_pt_mod)[2])
 ```