resde-vignette.R
In resde: Estimation in Reducible Stochastic Differential Equations

## ----setup, include=FALSE, cache=FALSE---------------------------
library(knitr)
opts_chunk$set(fig.align='center', # fig.show='hold',
dev='pdf', out.width='.45\\textwidth') # , highlight=FALSE)
old_options <- options(width=67)

## ----sdemodel, eval=FALSE----------------------------------------
#  sdemodel(phi=~x, phiprime=NULL, beta0=~beta0, beta1=~beta1, t0=0,
#           x0=0, mu0=0, mup=1, mum=1)

## ----------------------------------------------------------------
library(resde)
exmpl <- sdemodel(phi=~x^c, beta0=~b*a^c, beta1=~-b, mup=~sqrt(b))

## ----------------------------------------------------------------
lob301 <- Loblolly[Loblolly$Seed == 301, ]

## ----------------------------------------------------------------
m <- sdemodel(phi=~x^c, beta0=~b*a^c, beta1=~-b)  # else, defaults

## ----------------------------------------------------------------
f <- sdefit(m, x="height", t="age", data=lob301,
            start=c(a=60, b=0.1, c=1))
f

## ----------------------------------------------------------------
m <- sdemodel(phi=~x^c, beta0=~b*a^c, beta1=~-b, mum=0)
sdefit(m, x="height", t="age", data=lob301,
       start=c(a=60, b=0.1, c=1))

## ----------------------------------------------------------------
summary(f$fit)

## ----------------------------------------------------------------
m <- sdemodel(~log(abs(a^c - x^c)), beta0=~-b, beta1=0, mup=~b)
sdefit(m, x="height", t="age", data=lob301,
       start=c(a=70, b=0.1, c=1))

## ----------------------------------------------------------------
plot(c(0,height) ~ c(0, age), data=lob301)
curve(77.10687 * (1 - exp(-0.08405 * x))^(1/0.54946), add=T,
      col="red")
curve(72.5459 * (1 - exp(-0.0967 * x))^(1/0.5024), add=T,
      col="blue")

## ----------------------------------------------------------------
m <- sdemodel(phi=~bc(x/a, c), beta0=0, beta1=~-b,
              mup=~sqrt(abs(b)), mum=0)

## ----------------------------------------------------------------
alocal <- sdefit(m, x="height", t="age", unit="Seed",
          data=Loblolly, global=c(b=0.1, c=0.5), local=c(a=70))
alocal

## ----------------------------------------------------------------
(blocal <- sdefit(m, x="height", t="age", unit="Seed",
           data=Loblolly, global=c(a=70, c=0.5), local=c(b=0.1)))

## ----------------------------------------------------------------
(ablocal <- sdefit(m, x="height", t="age", unit="Seed",
            data=Loblolly, global=c(c=0.5), local=c(a=70, b=0.1)))

## ----------------------------------------------------------------
(blocal_mx <- sdefit(m, x="height", t="age", unit="Seed",
              data=Loblolly, global=c(a=70, c=0.5), local=c(b=0.1),
              method="nlme", control = nlme::nlmeControl(pnlsTol =
              0.01)))

## ----------------------------------------------------------------
coef(blocal_mx$fit)

## ----------------------------------------------------------------
m <- sdemodel(phi=~bc(x/a, c), beta0=0, beta1=~-b,
              mup=~sqrt(abs(b)))
sdefit(m, x="height", t="age", unit="Seed", data=Loblolly, 
       global=c(a=70, c=0.5), local=c(b=0.1))

## ----------------------------------------------------------------
lgLik <- eta <- seq(from=0, to=1, by=0.05)
for(i in seq_along(eta)){
  lgLik[i] <- (sdefit(m, x="height", t="age", unit="Seed",
               data=Loblolly, global=c(a=73, c=0.49),
               local=c(b=0.095), known=c(eta=eta[i]))$more
              )["logLik"]
}   
plot(lgLik ~ eta)

## ----------------------------------------------------------------
curve(unitran(x, alpha=-1, beta=0))

## ----------------------------------------------------------------
m <- sdemodel(phi=~unitran(x/a, alpha=alpha, beta=beta,
  reverse="no"), phiprime=~unitran_prime(x/a, alpha=alpha,
  beta=beta, reverse="no")/a, beta0=~b, beta1=0, mum=0)
(f <- sdefit(m, x="height", t="age", data=lob301, start=c(a=70,
  b=0.1, alpha=0.5, beta=0)))

## ----cleanup, include=FALSE, cache=FALSE--------------------------------------
options(old_options)