OnofriAndreaPG / aomisc

Utilities for stat courses
21 stars 7 forks source link

Calculate second derivative with DRC.PowerCurve #4

Open epstewart111 opened 3 years ago

epstewart111 commented 3 years ago

DRC.powerCurve() provides an estimate for the first derivative (i.e., slope), however, I was wondering if it's possible to get the second derivative? According to R documentation (http://ftp.uni-bayreuth.de/math/statlib/R/CRAN/doc/packages/drc.pdf) you can specify which parameters to estimate in the 'fixed' statement

LL.4(fixed = c(NA, NA, NA, NA), names = c("b", "c", "d", "e"), ...) However, when I try to do this in R with DRC.PowerCurve setting the parameter estimates to zero, I receive the error, "Not correct 'fixed' argument". I'm also unsure which initial correspond to the second derivative.

I can find little to no information on DRC.PowerCurve, let alone obtaining the second derivative.

OnofriAndreaPG commented 3 years ago

Hello, thanks for your message and sorry for delayed response. At the moment, you cannot put constraints on DRC.PowerCurve(), simply because they are not yet implemented (sorry, the aomisc package was mainly meant as a teaching tool; I decided to share it, but it is by no means a full blast R package...; that'll take some time...). As for the derivatives, try the following code as the starting point (you will need to update the aomisc package to its final version, though). Look forward to any comments...

Best regards

Andrea

# Use with drm and nls

rm(list=ls())
library(aomisc)
data(speciesArea)
model <- drm(numSpecies ~ Area, fct = DRC.powerCurve(),
             data = speciesArea)
model2 <- nls(numSpecies ~ NLS.powerCurve(Area, a, b),
             data = speciesArea)

# First derivative
D(expression(a * X^b), "X")
## a * (X^(b - 1) * b)

# Second derivative
D(D(expression(a * X^b), "X"), "X")
## a * (X^((b - 1) - 1) * (b - 1) * b)

funList <- list(~ a * (X^(b - 1) * b), ~ a * (X^((b - 1) - 1) * (b - 1) * b))
propdF <- data.frame(X = seq(1, 5, 1))
pred <- gnlht(model, funList, const = propdF, parameterNames = c("a", "b"))
pred
pred2 <- gnlht(model2, funList, const = propdF, parameterNames = c("a", "b"))
pred2