R/util-basicandfitfunc.R

In DRomics: Dose Response for Omics

###########################################
# Basic functions defining models 
# for plot, fit, BMD calculation
###########################################

flin <- function(x, b, d)
  # x the dose
  # d the y value for a null dose
  # b the slope
{
  d  + b*x
}
## inverse lin (X for an y value)
invlin <- function(y, b, d)
{
  (y - d) /b
}

### Expo model and starting values
formExp3p <- stats::as.formula(signal ~ d + b * (exp(dose/e)  - 1) )
startvalExp3pnls.1 <- function(xm, ym, increase, Ushape)
  # inputs
  # - xm unique values of the dose (sorted by dose)
  # - ym means of the signal at each value of xm (sorted by dose)
  # - Ushape TRUE if U shape FALSE if umbrella shape
{
  # initial value of d
  d <- ym[1] 
  # initial value of absolute valeur of e
  eabs <- 0.1*max(xm)
  
  if ((increase & !Ushape) | (!increase & Ushape))
  {
    # initial value of c even if it does not appear as a parameter
    c <- ym[which.max(xm)] # mean of the response at the highest dose
    b <- d - c
    # initial value of e that corresponds to its value in the 2 parameter model
    e <- -eabs
  } else
  {
    # initial value of e that corresponds to its value in the 2 parameter model
    e <- eabs
    # initial value of b 
    reg <- stats::lm(ym ~ exp(xm / e))
    b <- stats::coef(reg)[2]
  }
  startval <- list(b = b, d = d, e = e)
}

startvalExp3pnls.2 <- function(xm, ym, increase, Ushape)
  # inputs
  # - xm unique values of the dose (sorted by dose)
  # - ym means of the signal at each value of xm (sorted by dose)
  # - Ushape TRUE if U shape FALSE if umbrella shape
{
  # initial value of d
  d <- ym[1] 
  # initial value of absolute valeur of e
  eabs <- max(xm)     # WHAT CHANGES IN THE SECOND TRIAL !!!!!!!!!!!
  
  if ((increase & !Ushape) | (!increase & Ushape))
  {
    # initial value of c even if it does not appear as a parameter
    c <- ym[which.max(xm)] # mean of the response at the highest dose
    b <- d - c
    # initial value of e that corresponds to its value in the 2 parameter model
    e <- -eabs
  } else
  {
    # initial value of e that corresponds to its value in the 2 parameter model
    e <- eabs
    # initial value of b 
    reg <- stats::lm(ym ~ exp(xm / e)) 
    b <- stats::coef(reg)[2]
  }
  startval <- list(b = b, d = d, e = e)
}


# for plot
fExpo <- function(x, b, d, e)
{
  d + b * (exp(x/e)  - 1)
}
# x the dose
# d the y value for a null dose
# b a shape parameter 
# e another shape parameter 
# if e < 0 the model tends to a non infinite limit (d - b) when x tends toward infinity
# the model is increasing if e*b > 0

## inverse Expo (X for an y value)
invExpo <- function(y, b, d, e)
{
  if ( ((e < 0) & (b < 0) & (y > d - b)) | ((e < 0) & (b > 0) & (y < d - b)) )
    {return(NaN)} else {return(e * log(1 + (y - d) / b))}
}

### Hill model and starting values
formHill <- stats::as.formula(signal ~ c + (d - c) / (1 + (dose/e)^b ) )
startvalHillnls2 <- function(x, y, xm, ym, increase) # requires the definition of increase from min and max values
  # inputs
  # - x values of the dose
  # - y values the corresponding signal
  # - xm unique values of the dose (sorted by dose)
  # - ym means of the signal at each value of xm (sorted by dose)
  # 
{
  maxi <- max(y, na.rm = TRUE)
  mini <- min(y, na.rm = TRUE)
  ampl <- maxi - mini
  # inflate maxi and mini so as all values are strictly inside the interval [mini; maxi]
  maxi <- maxi + 0.001 * ampl
  mini <- mini - 0.001 * ampl
  # initial value of c
  c <- ifelse(increase, maxi, mini) 
  # initial value of d
  d <-ifelse(increase, mini, maxi) 
  # initial value of e and b from regression
  yreg <- log((d - c) / (y[x!=0] - c) - 1)
  xreg <- log(x[x!=0])
  reg <- stats::lm(yreg ~ xreg)
  b <- reg$coefficients[2]
  e <- reg$coefficients[1] / (-b)
  startval <- list(b = b, c = c, d = d, e = e)
}


# for plot
fHill <- function(x, b, c, d, e)
{
  c + (d - c) / (1 + (x/e)^b)
}

# x the dose
# c the y value for a high dose
# d the y value for a null dose
# b a shape parameter > 0
# e the dose at the inflexion point 

## inverse Hill (X for an y value)
invHill <- function(y, b, c, d, e)
{
  if ( ((d < c) & (y > c)) | ((d > c) & (y < c)) )
    {return(NaN)} else
      {return(e * ((d - y) / (y - c))^(1/b))}
}

### Gaussian model 5 p and starting values
formGauss5p <- stats::as.formula(signal ~ f * exp(-0.5 * ((dose-e)/b)^2) + d + (c - d) * pnorm((dose-e)/b)) 
startvalGauss5pnls <- function(xm, ym, Ushape)
  # inputs
  # - xm unique values of the dose (sorted by dose)
  # - ym means of the signal at each value of xm (sorted by dose)
  # - Ushape TRUE if U shape FALSE if umbrella shape
{
  # initial value of d
  d <- ym[1]# mean of the response at the first dose
  # initial value of c
  c <- ym[which.max(xm)]# mean of the response at the highest dose
  # initial value of f
  yextremum <- ifelse(Ushape, min(ym), max(ym))
  f <-  yextremum - (c + d) / 2 # amplitude of the gaussian part
  # initial value of b (standard error) assuming (after having looked many curves) the
  # bell shape extends to the whole dose range so dose max = 4 sd = 4 * b
  # 
  b <- max(xm) / 4
  # initial value of e (dose corresponding to the maximal (or minimal) signal)
  xextremum <- stats::median(xm[which(ym == yextremum)]) # just in case there is more than one dose at which ym == yextremum
  e <- min(xextremum - (c - d)*b/(f*sqrt(2*pi)), 1e-6) 
  startval <- list(b = b, c = c, d = d, e = e, f = f)
}

## simplified version for cases where Gauss 5p does not converge (c = d)
formGauss4p <- stats::as.formula(signal ~ f * exp(-0.5 * ((dose-e)/b)^2) + d ) # without c and without probit part
startvalGauss4pnls <- function(xm, ym, Ushape)
  # inputs
  # - xm unique values of the dose (sorted by dose)
  # - ym means of the signal at each value of xm (sorted by dose)
  # - Ushape TRUE if U shape FALSE if umbrella shape
{
  # initial value of d
  d <- ym[which.max(xm)]# mean of the response at the highest dose (supposed identical to the other asymptote)
  # initial value of f
  yextremum <- ifelse(Ushape, min(ym), max(ym))
  f <-  yextremum - d # amplitude of the gaussian part
  # initial value of b (standard error) assuming (after having looked many curves) the
  # bell shape extends to the whole dose range so dose max = 4 sd = 4 * b
  # with 2sd between 0 and e -> b = e /2
  b <- max(xm) / 4
  # initial value of e (dose corresponding to the maximal (or minimal) signal)
  xextremum <- stats::median(xm[which(ym == yextremum)] )# just in case there is more than one dose with ym = yextremum
  e <- min(xextremum, 1e-6)
  startval <- list(b = b, d = d, e = e, f = f)
}

### probit model
formprobit <- stats::as.formula(signal ~ d + (c - d) * pnorm((dose-e)/b)) 



# for plot

fGauss5p <- function(x, b, c, d, e, f)
{
  f * exp(-0.5 * ((x-e)/b)^2) + # gaussian part
    d + (c - d) * stats::pnorm((x-e)/b) # probit part
}

fGauss5pBMR <- function(x, b, c, d, e, g, threshold)
{
  g * exp(-0.5 * ((x-e)/b)^2) + # gaussian part
    d + (c - d) * stats::pnorm((x-e)/b) - threshold # probit part
  
}

fGauss5pBMR_xinlog <- function(xinlog, b, c, d, e, g, threshold)
{
  x <- exp(xinlog)
  g * exp(-0.5 * ((x-e)/b)^2) + # gaussian part
    d + (c - d) * stats::pnorm((x-e)/b) - threshold # probit part
  
}

fprobit <- function(x, b, c, d, e) 
{
  d + (c - d) * stats::pnorm((x - e)/b)
}

## inverse probit (X for an y value)
invprobit <- function(y, b, c, d, e)
{
  if ( ((d < c) & (y > c)) | ((d > c) & (y < c)) )
    {return(NaN)} else
      {return(e + b *stats::qnorm((y - d) / (c - d)))}
}


fGauss5poutofrange <- function(fit, signalmin, signalmax)
# TRUE if the fit gives an extremum value out of the signal range 
# function that takes the result of nls as first argument
# to be used before the choice of the best model
{
  par <- stats::coef(fit)
  b.i <- par["b"]
  c.i <- par["c"]
  d.i <- par["d"]
  e.i <- par["e"]
  f.i <- par["f"]
  xextr.i <- e.i + (c.i - d.i)*b.i / (f.i*sqrt(2*pi)) 
  yextr.i <- fGauss5p(xextr.i, b = b.i, c = c.i, d = d.i, e = e.i, f = f.i)
  outofrange <- (yextr.i > signalmax) | (yextr.i < signalmin)
}

fGauss4poutofrange <- function(fit, signalmin, signalmax)
  # TRUE if the fit gives an extremum value out of the signal range 
  # function that takes the result of nls as first argument
  # to be used before the choice of the best model
{
  par <- stats::coef(fit)
  b.i <- par["b"]
  c.i <- par["d"]
  d.i <- par["d"]
  e.i <- par["e"]
  f.i <- par["f"]
  xextr.i <- e.i + (c.i - d.i)*b.i / (f.i*sqrt(2*pi)) 
  yextr.i <- fGauss5p(xextr.i, b = b.i, c = c.i, d = d.i, e = e.i, f = f.i)
  outofrange <- (yextr.i > signalmax) | (yextr.i < signalmin)
}


### Gaussian model 5 p and starting values in log scale
formLGauss5p <- stats::as.formula(signal ~ f * exp(-0.5 * (log(dose/e)/b)^2) + d + (c - d) * pnorm(log(dose/e)/b)) 
startvalLGauss5pnls <- function(xm, ym, Ushape)
  # inputs
  # - xm unique values of the dose (sorted by dose)
  # - ym means of the signal at each value of xm (sorted by dose)
  # - Ushape TRUE if U shape FALSE if umbrella shape
{
  # initial value of d
  d <- ym[1]# mean of the response at the first dose
  # initial value of c
  c <- ym[which.max(xm)]# mean of the response at the highest dose
  # initial value of f
  yextremum <- ifelse(Ushape, min(ym), max(ym))
  f <-  yextremum - (c + d) / 2 # amplitude of the gaussian part
  # initial value of b (standard error) assuming (after having looked many curves) the
  # bell shape extends to the whole dose range so dose max = 4 sd = 4 * b IN LOG SCALE !!!
  b <- (log(max(xm)) - log(min(xm[xm!=0]))) / 4
  # initial value of e (dose corresponding to the maximal (or minimal) signal)
  xextremum <- stats::median(xm[which(ym == yextremum)]) # median just in case there is more than one dose with ym = yextremum
  e <- xextremum * exp(- (c - d)*b/(f*sqrt(2*pi))) 
  startval <- list(b = b, c = c, d = d, e = e, f = f)
}


## simplified version for cases where Gauss 5p dose not converge (c = d)
formLGauss4p <- stats::as.formula(signal ~ f * exp(-0.5 * (log(dose/e)/b)^2) + d) # without c and without probit part
startvalLGauss4pnls <- function(xm, ym, Ushape)
  # inputs
  # - xm unique values of the dose (sorted by dose)
  # - ym means of the signal at each value of xm (sorted by dose)
  # - Ushape TRUE if U shape FALSE if umbrella shape
{
  # initial value of d
  d <- ym[1]# mean of the response at the smallest dose
  # initial value of f
  yextremum <- ifelse(Ushape, min(ym), max(ym))
  f <-  yextremum - d # amplitude of the gaussian part
  # initial value of b (standard error) assuming (after having looked many curves) the
  # bell shape extends to the whole dose range so dose max = 4 sd = 4 * b IN LOG SCALE !!!
  b <- (log(max(xm)) - log(min(xm[xm!=0]))) / 4
  # initial value of e (dose corresponding to the maximal (or minimal) signal)
  xextremum <- stats::median(xm[which(ym == yextremum)]) # just in case there is more than one dose with ym = yextremum
  e <- xextremum 
  startval <- list(b = b, d = d, e = e, f = f)
}


### log-probit model - used to simplify the LGP model
formLprobit <- stats::as.formula(signal ~ d + (c - d) * pnorm(log(dose/e)/b)) 

### log-probit starting values - no more used
# startvalLprobitnls1 <- function(xm, ym) # to suppress
#   # inputs
#   # - xm unique values of the dose (sorted by dose)
#   # - ym means of the signal at each value of xm (sorted by dose)
#   # - Ushape TRUE if U shape FALSE if umbrella shape
# {
#   # initial value of d
#   d <- ym[1]# mean of the response at the first dose
#   # initial value of c
#   c <- ym[which.max(xm)]# mean of the response at the highest dose
#   # initial value of b (standard error) assuming (after having looked many curves) the
#   # bell shape extends to the whole dose range so dose max = 4 sd = 4 * b IN LOG SCALE !!!
#   ldosemin <- log(min(xm[xm!=0]))
#   ldosemax <- log(max(xm))
#   b <- (ldosemax - ldosemin) / 4
#   # initial value of e (in the middle of the dose range in log)
#   e <- exp( ldosemin + (ldosemax - ldosemin) /2 )
#   startval <- list(b = b, c = c, d = d, e = e)
# }
# 
# startvalLprobitnls2 <- function(x, y, xm, ym, increase) # requires the definition of increase from min and max values
#   # inputs
#   # - x values of the dose
#   # - y values the corresponding signal
#   # - xm unique values of the dose (sorted by dose)
#   # - ym means of the signal at each value of xm (sorted by dose)
#   # 
# {
#   maxi <- max(y, na.rm = TRUE)
#   mini <- min(y, na.rm = TRUE)
#   ampl <- maxi - mini
#   # inflate maxi and mini so as all values are strictly inside the interval [mini; maxi]
#   maxi <- maxi + 0.001 * ampl
#   mini <- mini - 0.001 * ampl
#   # initial value of c
#   c <- ifelse(increase, maxi, mini) 
#   # initial value of d
#   d <-ifelse(increase, mini, maxi) 
#   # initial value of e and b from regression
#   Y <- (y[x!=0] - d) / (c - d)
#   yreg <- qnorm( Y )
#   xreg <- log(x[x!=0])
#   reg <- lm(yreg ~ xreg)
#   b <- 1/ reg$coefficients[2]
#   e <- reg$coefficients[1] * (-b)
#   startval <- list(b = b, c = c, d = d, e = e)
# }


# for plot

fLGauss5p <- function(x, b, c, d, e, f)
{
  f * exp(-0.5 * (log(x/e)/b)^2) + # gaussian part
    d + (c - d) * stats::pnorm(log(x/e)/b) # probit part
}

fLGauss5pBMR <- function(x, b, c, d, e, g, threshold)
{
  g * exp(-0.5 * (log(x/e)/b)^2) + # gaussian part
    d + (c - d) * stats::pnorm(log(x/e)/b) - threshold # probit part
  
}

fLGauss5pBMR_xinlog <- function(xinlog, b, c, d, e, g, threshold)
{
  g * exp(-0.5 * ((xinlog - log(e))/b)^2) + # gaussian part
    d + (c - d) * stats::pnorm((xinlog - log(e))/b) - threshold # probit part
  
}


fLprobit <- function(x, b, c, d, e)
{
  d + (c - d) * stats::pnorm(log(x/e)/b)
}

## inverse Lprobit (X for an y value)
invLprobit <- function(y, b, c, d, e)
{
  if ( ((d < c) & (y > c)) | ((d > c) & (y < c)) )
    {return(NaN)} else
      {return(e * exp(stats::qnorm((y - d) / (c - d)) *b))}
}


fLGauss5poutofrange <- function(fit, signalmin, signalmax)
  # TRUE if the fit gives an extremum value out of the signal range 
  # function that takes the result of nls as first argument
  # to be used before the choice of the best model
{
  par <- stats::coef(fit)
  b.i <- par["b"]
  c.i <- par["c"]
  d.i <- par["d"]
  e.i <- par["e"]
  f.i <- par["f"]
  xextr.i <- exp(log(e.i) + (c.i - d.i)*b.i/(f.i*sqrt(2*pi))) 
  yextr.i <- fLGauss5p(xextr.i, b = b.i, c = c.i, d = d.i, e = e.i, f = f.i)
  outofrange <- (yextr.i > signalmax) | (yextr.i < signalmin)
}

fLGauss4poutofrange <- function(fit, signalmin, signalmax)
  # TRUE if the fit gives an extremum value out of the signal range 
  # function that takes the result of nls as first argument
  # to be used before the choice of the best model
{
  par <- stats::coef(fit)
  b.i <- par["b"]
  c.i <- par["d"]
  d.i <- par["d"]
  e.i <- par["e"]
  f.i <- par["f"]
  xextr.i <- exp(log(e.i) + (c.i - d.i)*b.i/(f.i*sqrt(2*pi))) 
  yextr.i <- fLGauss5p(xextr.i, b = b.i, c = c.i, d = d.i, e = e.i, f = f.i)
  outofrange <- (yextr.i > signalmax) | (yextr.i < signalmin)
}