kurtosis, skewness, start with ML

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: AMR
-Version: 0.2.0.9008
-Date: 2018-07-04
+Version: 0.2.0.9009
+Date: 2018-07-06
 Title: Antimicrobial Resistance Analysis
 Authors@R: c(
     person(
@@ -28,6 +28,7 @@ Depends:
     R (>= 3.0.0)
 Imports:
     backports,
+    broom,
     clipr,
     curl,
     dplyr (>= 0.7.0),

NAMESPACE

@@ -6,6 +6,11 @@ S3method(as.integer,mic)
 S3method(as.numeric,mic)
 S3method(barplot,mic)
 S3method(barplot,rsi)
+S3method(hist,frequency_tbl)
+S3method(kurtosis,data.frame)
+S3method(kurtosis,default)
+S3method(kurtosis,matrix)
+S3method(plot,frequency_tbl)
 S3method(plot,mic)
 S3method(plot,rsi)
 S3method(print,data.table)
@@ -14,6 +19,9 @@ S3method(print,mic)
 S3method(print,rsi)
 S3method(print,tbl)
 S3method(print,tbl_df)
+S3method(skewness,data.frame)
+S3method(skewness,default)
+S3method(skewness,matrix)
 S3method(summary,mic)
 S3method(summary,rsi)
 export("%like%")
@@ -43,6 +51,7 @@ export(interpretive_reading)
 export(is.mic)
 export(is.rsi)
 export(key_antibiotics)
+export(kurtosis)
 export(left_join_microorganisms)
 export(like)
 export(mo_property)
@@ -54,6 +63,7 @@ export(rsi)
 export(rsi_df)
 export(rsi_predict)
 export(semi_join_microorganisms)
+export(skewness)
 export(top_freq)
 exportMethods(as.data.frame.frequency_tbl)
 exportMethods(as.double.mic)
@@ -61,6 +71,12 @@ exportMethods(as.integer.mic)
 exportMethods(as.numeric.mic)
 exportMethods(barplot.mic)
 exportMethods(barplot.rsi)
+exportMethods(hist.frequency_tbl)
+exportMethods(kurtosis)
+exportMethods(kurtosis.data.frame)
+exportMethods(kurtosis.default)
+exportMethods(kurtosis.matrix)
+exportMethods(plot.frequency_tbl)
 exportMethods(plot.mic)
 exportMethods(plot.rsi)
 exportMethods(print.data.table)
@@ -69,8 +85,13 @@ exportMethods(print.mic)
 exportMethods(print.rsi)
 exportMethods(print.tbl)
 exportMethods(print.tbl_df)
+exportMethods(skewness)
+exportMethods(skewness.data.frame)
+exportMethods(skewness.default)
+exportMethods(skewness.matrix)
 exportMethods(summary.mic)
 exportMethods(summary.rsi)
+importFrom(broom,tidy)
 importFrom(clipr,read_clip_tbl)
 importFrom(clipr,write_clip)
 importFrom(curl,nslookup)
@@ -107,6 +128,7 @@ importFrom(dplyr,vars)
 importFrom(grDevices,boxplot.stats)
 importFrom(graphics,axis)
 importFrom(graphics,barplot)
+importFrom(graphics,hist)
 importFrom(graphics,plot)
 importFrom(graphics,text)
 importFrom(knitr,kable)
@@ -119,6 +141,7 @@ importFrom(rvest,html_nodes)
 importFrom(rvest,html_table)
 importFrom(stats,fivenum)
 importFrom(stats,mad)
+importFrom(stats,na.omit)
 importFrom(stats,pchisq)
 importFrom(stats,sd)
 importFrom(tibble,tibble)

NEWS.md

@@ -1,10 +1,14 @@
 # 0.2.0.90xx (development version)
 #### New
 * Support for Addins menu in RStudio to quickly insert `%in%` or `%like%` (and give them keyboard shortcuts), or to view the datasets that come with this package
-* Function `top_freq` function to get the top/below *n* items of frequency tables
-* Vignette about frequency tables
-* Header of frequency tables now also show MAD and IQR
-* Possibility to globally set the default for the amount of items to print in frequency tables (`freq` function), with `options(max.print.freq = n)`
+* For convience, descriptive statistical functions `kurtosis` and `skewness` that are lacking in base R - they are generic functions and have support for vectors, data.frames and matrices
+* New for frequency tables (function `freq`):
+  * A vignette to explain its usage
+  * Support for existing functions `hist` and `plot` to use a frequency table as input: `hist(freq(df$age))`
+  * Support for quasiquotation: `freq(mydata, mycolumn)` is the same as `mydata %>% freq(mycolumn)`
+  * Function `top_freq` function to return the top/below *n* items as vector
+  * Header of frequency tables now also show Mean Absolute Deviaton (MAD) and Interquartile Range (IQR)
+  * Possibility to globally set the default for the amount of items to print, with `options(max.print.freq = n)` where *n* is your preset value
 * Functions `clipboard_import` and `clipboard_export` as helper functions to quickly copy and paste from/to software like Excel and SPSS
 * Function `g.test` to perform the Χ<sup>2</sup> distributed [*G*-test](https://en.wikipedia.org/wiki/G-test)
 * Function `ratio` to transform a vector of values to a preset ratio (convenient to use with `g.test`). For example:
@@ -16,7 +20,6 @@ ratio(c(772, 1611, 737), ratio = "1:2:1")
 
 #### Changed
 * `%like%` now supports multiple patterns
-* Frequency tables (function `freq`) now supports quasiquotation: `freq(mydata, mycolumn)`, or `mydata %>% freq(mycolumn)`
 * Frequency tables are now actual `data.frame`s with altered console printing to make it look like a frequency table. Because of this, the parameter `toConsole` is not longer needed.
 * Small translational improvements to the `septic_patients` dataset
 * Small improvements to the `microorganisms` dataset, especially for *Salmonella*

R/freq.R

@@ -46,8 +46,8 @@
 #'
 #' For dates and times of any class, these additional values will be calculated with \code{na.rm = TRUE} and shown into the header:
 #' \itemize{
-#'   \item{Oldest, using \code{\link[base]{min}}}
-#'   \item{Newest, using \code{\link[base]{max}}, with difference between newest and oldest}
+#'   \item{Oldest, using \code{\link{min}}}
+#'   \item{Newest, using \code{\link{max}}, with difference between newest and oldest}
 #'   \item{Median, using \code{\link[stats]{median}}, with percentage since oldest}
 #' }
 #'
@@ -522,3 +522,42 @@ as.data.frame.frequency_tbl <- function(x, ...) {
   attr(x, 'opt') <- NULL
   as.data.frame.data.frame(x, ...)
 }
+
+#' @noRd
+#' @exportMethod hist.frequency_tbl
+#' @export
+#' @importFrom dplyr %>% pull
+#' @importFrom graphics hist
+hist.frequency_tbl <- function(x, ...) {
+
+  opt <- attr(x, 'opt')
+
+  if (!is.null(opt$vars)) {
+    title <- opt$vars
+  } else {
+    title <- ""
+  }
+
+  items <- x %>% pull(item)
+  counts <- x %>% pull(count)
+  vect <- rep(items, counts)
+  hist(vect, main = paste("Histogram of", title), xlab = title, ...)
+}
+
+#' @noRd
+#' @exportMethod plot.frequency_tbl
+#' @export
+#' @importFrom dplyr %>% pull
+plot.frequency_tbl <- function(x, y, ...) {
+  opt <- attr(x, 'opt')
+
+  if (!is.null(opt$vars)) {
+    title <- opt$vars
+  } else {
+    title <- ""
+  }
+
+  items <- x %>% pull(item)
+  counts <- x %>% pull(count)
+  plot(x = items, y = counts, ylab = "Count", xlab = title, ...)
+}

R/globals.R

@@ -35,12 +35,14 @@ globalVariables(c('abname',
                   'key_ab',
                   'key_ab_lag',
                   'key_ab_other',
+                  'labs',
                   'median',
                   'mic',
                   'microorganisms',
                   'mocode',
                   'molis',
                   'n',
+                  'na.omit',
                   'other_pat_or_mo',
                   'patient_id',
                   'quantile',

R/kurtosis.R

@@ -0,0 +1,40 @@
+#' Kurtosis of the sample
+#'
+#' @description Kurtosis is a measure of the "tailedness" of the probability distribution of a real-valued random variable.
+#'
+#' @param x a vector of values, a \code{matrix} or a \code{data frame}
+#' @param na.rm a logical value indicating whether \code{NA} values should be stripped before the computation proceeds.
+#' @exportMethod kurtosis
+#' @seealso \code{\link{skewness}}
+#' @rdname kurtosis
+#' @export
+kurtosis <- function(x, na.rm = FALSE) {
+  UseMethod("kurtosis")
+}
+
+#' @exportMethod kurtosis.default
+#' @rdname kurtosis
+#' @export
+kurtosis.default <- function (x, na.rm = FALSE) {
+  x <- as.vector(x)
+  if (na.rm == TRUE) {
+    x <- x[!is.na(x)]
+  }
+  n <- length(x)
+  n * base::sum((x - base::mean(x, na.rm = na.rm))^4, na.rm = na.rm) /
+    (base::sum((x - base::mean(x, na.rm = na.rm))^2, na.rm = na.rm)^2)
+}
+
+#' @exportMethod kurtosis.matrix
+#' @rdname kurtosis
+#' @export
+kurtosis.matrix <- function (x, na.rm = FALSE) {
+  base::apply(x, 2, kurtosis.default, na.rm = na.rm)
+}
+
+#' @exportMethod kurtosis.data.frame
+#' @rdname kurtosis
+#' @export
+kurtosis.data.frame <- function (x, na.rm = FALSE) {
+  base::sapply(x, kurtosis.default, na.rm = na.rm)
+}

R/like.R

@@ -18,7 +18,7 @@
 
 #' Pattern Matching
 #'
-#' Convenient wrapper around \code{\link[base]{grepl}} to match a pattern: \code{a \%like\% b}. It always returns a \code{logical} vector and is always case-insensitive. Also, \code{pattern} (\code{b}) can be as long as \code{x} (\code{a}) to compare items of each index in both vectors.
+#' Convenient wrapper around \code{\link[base]{grep}} to match a pattern: \code{a \%like\% b}. It always returns a \code{logical} vector and is always case-insensitive. Also, \code{pattern} (\code{b}) can be as long as \code{x} (\code{a}) to compare items of each index in both vectors.
 #' @inheritParams base::grepl
 #' @return A \code{logical} vector
 #' @name like

R/skewness.R

@@ -0,0 +1,40 @@
+#' Skewness of the sample
+#'
+#' @description Skewness is a measure of the asymmetry of the probability distribution of a real-valued random variable about its mean.
+#'
+#' When negative: the left tail is longer; the mass of the distribution is concentrated on the right of the figure. When positive: the right tail is longer; the mass of the distribution is concentrated on the left of the figure.
+#' @param x a vector of values, a \code{matrix} or a \code{data frame}
+#' @param na.rm a logical value indicating whether \code{NA} values should be stripped before the computation proceeds.
+#' @exportMethod skewness
+#' @seealso \code{\link{kurtosis}}
+#' @rdname skewness
+#' @export
+skewness <- function(x, na.rm = FALSE) {
+  UseMethod("skewness")
+}
+
+#' @exportMethod skewness.default
+#' @rdname skewness
+#' @export
+skewness.default <- function (x, na.rm = FALSE) {
+  x <- as.vector(x)
+  if (na.rm == TRUE) {
+    x <- x[!is.na(x)]
+  }
+  n <- length(x)
+  (base::sum((x - base::mean(x))^3) / n) / (base::sum((x - base::mean(x))^2) / n)^(3/2)
+}
+
+#' @exportMethod skewness.matrix
+#' @rdname skewness
+#' @export
+skewness.matrix <- function (x, na.rm = FALSE) {
+  base::apply(x, 2, skewness.default, na.rm = na.rm)
+}
+
+#' @exportMethod skewness.data.frame
+#' @rdname skewness
+#' @export
+skewness.data.frame <- function (x, na.rm = FALSE) {
+  base::sapply(x, skewness.default, na.rm = na.rm)
+}

R/trends.R

@@ -0,0 +1,123 @@
+#' Detect trends using Machine Learning
+#'
+#' Test text
+#' @param data a \code{data.frame}
+#' @param threshold_unique do not analyse more unique \code{threshold_unique} items per variable
+#' @param na.rm a logical value indicating whether \code{NA} values should be stripped before the computation proceeds.
+#' @param info print relevant combinations to console
+#' @return A \code{list} with class \code{"trends"}
+#' @importFrom stats na.omit
+#' @importFrom broom tidy
+# @export
+trends <- function(data, threshold_unique = 30, na.rm = TRUE, info = TRUE) {
+
+  cols <- colnames(data)
+  relevant <- list()
+  count <- 0
+  for (x in 1:length(cols)) {
+    for (y in 1:length(cols)) {
+      if (x == y) {
+        next
+      }
+      if (n_distinct(data[, x]) > threshold_unique | n_distinct(data[, y]) > threshold_unique) {
+        next
+      }
+      count <- count + 1
+      df <- data %>%
+        group_by_at(c(cols[x], cols[y])) %>%
+        summarise(n = n())
+      n <- df %>% pull(n)
+      # linear regression model
+      lin <- stats::lm(1:length(n) ~ n, na.action = ifelse(na.rm == TRUE, na.omit, NULL))
+
+      res <- list(
+        df = df,
+        x = cols[x],
+        y = cols[y],
+        m = base::mean(n, na.rm = na.rm),
+        sd = stats::sd(n, na.rm = na.rm),
+        cv = cv(n, na.rm = na.rm),
+        cqv = cqv(n, na.rm = na.rm),
+        kurtosis = kurtosis(n, na.rm = na.rm),
+        skewness = skewness(n, na.rm = na.rm),
+        lin.p = broom::tidy(lin)[2, 'p.value']
+        #binom.p <- broom::tidy(binom)[2, 'p.value']
+      )
+
+      include <- TRUE
+      # ML part
+      if (res$cv > 0.25) {
+        res$reason <- "cv > 0.25"
+      } else if (res$cqv > 0.75) {
+        res$reason <- "cqv > 0.75"
+      } else {
+        include <- FALSE
+      }
+
+      if (include == TRUE) {
+        relevant <- c(relevant, list(res))
+        if (info == TRUE) {
+          # minus one because the whole data will be added later
+          cat(paste0("[", length(relevant), "]"), "Relevant:", cols[x], "vs.", cols[y], "\n")
+        }
+      }
+
+    }
+  }
+
+  cat("Total of", count, "combinations analysed;", length(relevant), "seem relevant.\n")
+  class(relevant) <- 'trends'
+  relevant <- c(relevant, list(data = data))
+  relevant
+
+}
+
+# @exportMethod print.trends
+# @export
+#' @noRd
+print.trends <- function(x, ...) {
+  cat(length(x) - 1, "relevant trends, out of", length(x$data)^2, "\n")
+}
+
+# @exportMethod plot.trends
+# @export
+#' @noRd
+# plot.trends <- function(x, n = NULL, ...) {
+#   if (is.null(n)) {
+#     oask <- devAskNewPage(TRUE)
+#     on.exit(devAskNewPage(oask))
+#     n <- c(1:(length(x) - 1))
+#   } else {
+#     if (n > length(x) - 1) {
+#       stop('trend unavailable, max is ', length(x) - 1, call. = FALSE)
+#     }
+#     oask <- NULL
+#   }
+#   for (i in n) {
+#     data <- x[[i]]$df
+#     if (as.character(i) %like% '1$') {
+#       suffix <- "st"
+#     } else if (as.character(i) %like% '2$') {
+#       suffix <- "nd"
+#     } else if (as.character(i) %like% '3$') {
+#       suffix <- "rd"
+#     } else {
+#       suffix <- "th"
+#     }
+#     if (!is.null(oask)) {
+#       cat(paste("Coming up:", colnames(data)[1], "vs.", colnames(data)[2]), "\n")
+#     }
+#     print(
+#       ggplot(
+#         data,
+#         aes_string(x = colnames(data)[1],
+#                    y = colnames(data)[3],
+#                    group = colnames(data)[2],
+#                    fill = colnames(data)[2])) +
+#         geom_col(position = "dodge") +
+#         theme_minimal() +
+#         labs(title = paste(colnames(data)[1], "vs.", colnames(data)[2]),
+#              subtitle = paste0(i, suffix, " trend"))
+#     )
+#   }
+# }

README.md

@@ -32,7 +32,7 @@ With `AMR` you can also:
 * Get antimicrobial ATC properties from the WHO Collaborating Centre for Drug Statistics Methodology ([WHOCC](https://www.whocc.no/atc_ddd_methodology/who_collaborating_centre/)), to be able to:
   * Translate antibiotic codes (like *AMOX*), official names (like *amoxicillin*) and even trade names (like *Amoxil* or *Trimox*) to an [ATC code](https://www.whocc.no/atc_ddd_index/?code=J01CA04&showdescription=no) (like *J01CA04*) and vice versa with the `abname` function
   * Get the latest antibiotic properties like hierarchic groups and [defined daily dose](https://en.wikipedia.org/wiki/Defined_daily_dose) (DDD) with units and administration form from the WHOCC website with the `atc_property` function
-* Create frequency tables with the `freq` function
+* Conduct descriptive statistics: calculate kurtosis, skewness and create frequency tables
 
 And it contains:
 * A recent data set with ~2500 human pathogenic microorganisms, including family, genus, species, gram stain and aerobic/anaerobic
@@ -41,13 +41,17 @@ And it contains:
 
 With the `MDRO` function (abbreviation of Multi Drug Resistant Organisms), you can check your isolates for exceptional resistance with country-specific guidelines or EUCAST rules. Currently guidelines for Germany and the Netherlands are supported. Please suggest addition of your own country here: [https://github.com/msberends/AMR/issues/new](https://github.com/msberends/AMR/issues/new?title=New%20guideline%20for%20MDRO&body=%3C--%20Please%20add%20your%20country%20code,%20guideline%20name,%20version%20and%20source%20below%20and%20remove%20this%20line--%3E).
 
+#### Read all changes and new functions in [NEWS.md](NEWS.md).
+
 ## How to get it?
 This package is available on CRAN and also here on GitHub.
 
 ### From CRAN (recommended)
 [![CRAN_Badge](https://img.shields.io/cran/v/AMR.svg?label=CRAN&colorB=3679BC)](http://cran.r-project.org/package=AMR)
+
+Downloads via RStudio CRAN server (downloads by all other CRAN mirrors not measured):
 [![CRAN_Downloads](https://cranlogs.r-pkg.org/badges/grand-total/AMR)](http://cran.r-project.org/package=AMR)
-[![CRAN_Downloads](https://cranlogs.r-pkg.org/badges/AMR)](http://cran.r-project.org/package=AMR)
+[![CRAN_Downloads](https://cranlogs.r-pkg.org/badges/AMR)](https://cranlogs.r-pkg.org/downloads/daily/last-month/AMR)
   
 - <img src="http://www.rstudio.com/favicon.ico" alt="RStudio favicon" height="20px"> In [RStudio](http://www.rstudio.com) (recommended):
   - Click on `Tools` and then `Install Packages...`

man/freq.Rd

@@ -66,8 +66,8 @@ For numeric values of any class, these additional values will all be calculated
 
 For dates and times of any class, these additional values will be calculated with \code{na.rm = TRUE} and shown into the header:
 \itemize{
-  \item{Oldest, using \code{\link[base]{min}}}
-  \item{Newest, using \code{\link[base]{max}}, with difference between newest and oldest}
+  \item{Oldest, using \code{\link{min}}}
+  \item{Newest, using \code{\link{max}}, with difference between newest and oldest}
   \item{Median, using \code{\link[stats]{median}}, with percentage since oldest}
 }
 

man/kurtosis.Rd

@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/kurtosis.R
+\name{kurtosis}
+\alias{kurtosis}
+\alias{kurtosis.default}
+\alias{kurtosis.matrix}
+\alias{kurtosis.data.frame}
+\title{Kurtosis of the sample}
+\usage{
+kurtosis(x, na.rm = FALSE)
+
+\method{kurtosis}{default}(x, na.rm = FALSE)
+
+\method{kurtosis}{matrix}(x, na.rm = FALSE)
+
+\method{kurtosis}{data.frame}(x, na.rm = FALSE)
+}
+\arguments{
+\item{x}{a vector of values, a \code{matrix} or a \code{data frame}}
+
+\item{na.rm}{a logical value indicating whether \code{NA} values should be stripped before the computation proceeds.}
+}
+\description{
+Kurtosis is a measure of the "tailedness" of the probability distribution of a real-valued random variable.
+}
+\seealso{
+\code{\link{skewness}}
+}

man/like.Rd

@@ -29,7 +29,7 @@ x \%like\% pattern
 A \code{logical} vector
 }
 \description{
-Convenient wrapper around \code{\link[base]{grepl}} to match a pattern: \code{a \%like\% b}. It always returns a \code{logical} vector and is always case-insensitive. Also, \code{pattern} (\code{b}) can be as long as \code{x} (\code{a}) to compare items of each index in both vectors.
+Convenient wrapper around \code{\link[base]{grep}} to match a pattern: \code{a \%like\% b}. It always returns a \code{logical} vector and is always case-insensitive. Also, \code{pattern} (\code{b}) can be as long as \code{x} (\code{a}) to compare items of each index in both vectors.
 }
 \details{
 Using RStudio? This function can also be inserted from the Addins menu and can have its own Keyboard Shortcut like Ctrl+Shift+L or Cmd+Shift+L (see Tools > Modify Keyboard Shortcuts...).

man/skewness.Rd

@@ -0,0 +1,30 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/skewness.R
+\name{skewness}
+\alias{skewness}
+\alias{skewness.default}
+\alias{skewness.matrix}
+\alias{skewness.data.frame}
+\title{Skewness of the sample}
+\usage{
+skewness(x, na.rm = FALSE)
+
+\method{skewness}{default}(x, na.rm = FALSE)
+
+\method{skewness}{matrix}(x, na.rm = FALSE)
+
+\method{skewness}{data.frame}(x, na.rm = FALSE)
+}
+\arguments{
+\item{x}{a vector of values, a \code{matrix} or a \code{data frame}}
+
+\item{na.rm}{a logical value indicating whether \code{NA} values should be stripped before the computation proceeds.}
+}
+\description{
+Skewness is a measure of the asymmetry of the probability distribution of a real-valued random variable about its mean.
+
+When negative: the left tail is longer; the mass of the distribution is concentrated on the right of the figure. When positive: the right tail is longer; the mass of the distribution is concentrated on the left of the figure.
+}
+\seealso{
+\code{\link{kurtosis}}
+}

man/trends.Rd

@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/trends.R
+\name{trends}
+\alias{trends}
+\title{Detect trends using Machine Learning}
+\usage{
+trends(data, threshold_unique = 30, na.rm = TRUE, info = TRUE)
+}
+\arguments{
+\item{data}{a \code{data.frame}}
+
+\item{threshold_unique}{do not analyse more unique \code{threshold_unique} items per variable}
+
+\item{na.rm}{a logical value indicating whether \code{NA} values should be stripped before the computation proceeds.}
+
+\item{info}{print relevant combinations to console}
+}
+\value{
+A \code{list} with class \code{"trends"}
+}
+\description{
+Test text
+}

tests/testthat/test-freq.R

@@ -49,5 +49,9 @@ test_that("frequency table works", {
   # input must be freq tbl
   expect_error(septic_patients %>% top_freq(1))
 
+  # charts from plot and hist, should not raise errors
+  plot(freq(septic_patients, age))
+  hist(freq(septic_patients, age))
+
 })
 

tests/testthat/test-kurtosis.R

@@ -0,0 +1,13 @@
+context("kurtosis.R")
+
+test_that("kurtosis works", {
+  expect_equal(kurtosis(septic_patients$age),
+               6.423118,
+               tolerance = 0.00001)
+  expect_equal(unname(kurtosis(data.frame(septic_patients$age))),
+               6.423118,
+               tolerance = 0.00001)
+  expect_equal(kurtosis(matrix(septic_patients$age)),
+               6.423118,
+               tolerance = 0.00001)
+})

tests/testthat/test-skewness.R

@@ -0,0 +1,13 @@
+context("skewness.R")
+
+test_that("skewness works", {
+  expect_equal(skewness(septic_patients$age),
+               -1.637164,
+               tolerance = 0.00001)
+  expect_equal(unname(skewness(data.frame(septic_patients$age))),
+               -1.637164,
+               tolerance = 0.00001)
+  expect_equal(skewness(matrix(septic_patients$age)),
+               -1.637164,
+               tolerance = 0.00001)
+})