Introduction

knitr::opts_chunk$set(echo = TRUE)

The overarching aim of this package is to fit both a sigmoidal and a double-sigmoidal curve to time--intensity data and then determine which of the two models fit the data best. The sigmoidal model is assumed to start at zero and rise to an asymptotic value. The double-sigmoidal model first increases from 0 to a maximum intensity and the decreases to a final asymptotic intensity. After calculating the best candidate curves for both of models a decision algorithm determines which of the two candidates (if any) fits the data best. Finally, the package calculates key parameters that quantify the time--intensity data with respect to the best-fitting model.

For the sigmoidal model, key parameters describing the curve are:

For the double-sigmoidal model, key parameters describing the curve are:

Importantly, the parameters used to mathematically describe the fitted curves do not necessarily directly correspond to these intuitively meaningful parameters. Therefore, the sicegar package calculates the intuitively meaningful parameters from the estimated parameters. See the vignette on additional parameters for details.

Example fit on simulated input data

###*****************************
# INITIAL COMMANDS TO RESET THE SYSTEM
seedNo=14159
set.seed(seedNo)
###*****************************

###*****************************
require("sicegar")
require("dplyr")
require("ggplot2")
require("cowplot")
###*****************************

The input to the fitting function must be in the form of a data frame with two columns called time and intensity. We will here use double-sigmoidal data generated with some arbitrarily chosen parameters. We can use the function doublesigmoidalFitFormula() to generate a double-sigmoidal curve, to which we add random noise.

time <- seq(3, 24, 0.5)
noise_parameter <- 0.2
intensity_noise <- runif(n = length(time), min = 0, max = 1) * noise_parameter
intensity <- doublesigmoidalFitFormula(time,
                                    finalAsymptoteIntensityRatio = .3,
                                    maximum = 4,
                                    slope1Param = 1,
                                    midPoint1Param = 7,
                                    slope2Param = 1,
                                    midPointDistanceParam = 8)
intensity <- intensity+intensity_noise
dataInput <- data.frame(time, intensity)

ggplot(dataInput, aes(time, intensity)) + geom_point() + theme_bw()

We can now fit the two models to the data and determine which is the better fit. This is done with the function fitAndCategorize(). The three provided threshold parameters are used in the categorization process and depend on the units in which the data are measured. See the vignette on categorizing fits for details.

fitObj <- fitAndCategorize(dataInput,
                           threshold_minimum_for_intensity_maximum = 0.3,
                           threshold_intensity_range = 0.1,
                           threshold_t0_max_int = 0.05)

The two fitted curves can be visualized with the function figureModelCurves(), which returns a ggplot2 plot.

# Double-sigmoidal fit with parameter related lines
fig_a <- figureModelCurves(dataInput = fitObj$normalizedInput,
                                  sigmoidalFitVector = fitObj$sigmoidalModel,
                                  showParameterRelatedLines = TRUE)

fig_b <- figureModelCurves(dataInput = fitObj$normalizedInput,
                                  doubleSigmoidalFitVector = fitObj$doubleSigmoidalModel,
                                  showParameterRelatedLines = TRUE)

plot_grid(fig_a, fig_b, ncol = 2) # function from the cowplot package

Clearly the regular sigmoidal curve does not provide a good fit but the double-sigmoidal curve does. This information is available from the returned fit object:

fitObj$decisionProcess$decision # final decision

The fit object

The fit object returned by fitAndCategorize() contains all the information potentially of interest in the course of this type of an analysis. It consists of five distinct components:

names(fitObj)

Each component holds numerous parameters of interest:

Here, the contents of the summary vector is as follows:

str(fitObj$summaryVector)

All these parameters are defined in the vignette on additional parameters.



Try the sicegar package in your browser

Any scripts or data that you put into this service are public.

sicegar documentation built on Aug. 23, 2019, 5:05 p.m.