## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 5,
  fig.height = 4
)

## ----eval=FALSE---------------------------------------------------------------
# # Example code for fitting a Type 3 FR dynamical model:
# FR.formula = bf( NE | trials(N0) ~ Type3H_dyn(N0,P0,Time,b,h)/N0,
#                  b~1, h~1, nl=TRUE)
# FR.priors  = c(prior(exponential(1.0), nlpar="b", lb=0),
#                prior(exponential(1.0), nlpar="h", lb=0) )
# fit.1      = brm(FR.formula,
#                  family   = binomial(link="identity"),
#                  prior    = FR.priors,
#                  stanvars = stanvar(scode=Type3H_dyn_code, block="functions"),
#                  data     = df )

## ----warning=FALSE, message=FALSE---------------------------------------------
library(BayesFR)
library(brms)
library(ggplot2)

## -----------------------------------------------------------------------------
data(df_Michalko_and_Pekar_2017_AM_NAT)
df = subset(df_Michalko_and_Pekar_2017_AM_NAT, ID=="Figure 3c")
head(df)
ggplot(aes(N0,NE), data=df) +
  geom_jitter(width=0.1, height=0.1, alpha=0.6, size=2) +
  coord_cartesian(xlim=c(0,NA), ylim=c(0,NA))

## -----------------------------------------------------------------------------
FR.formula = bf( NE ~ a*N0/(1.0+a*h*N0)*Time,
                 a~1, h~1, 
                 nl = TRUE)

## -----------------------------------------------------------------------------
FR.formula = bf( NE ~ Type2H_fix(N0,1.0,7.0,a,h),
                 a~1, h~1, 
                 nl = TRUE)

## -----------------------------------------------------------------------------
FR.priors  = c(prior(exponential(1.0), nlpar="a", lb=0),
               prior(exponential(1.0), nlpar="h", lb=0))

## ----eval=FALSE---------------------------------------------------------------
# fit.1 = brm(FR.formula,
#             family   = poisson(link="identity"),
#             prior    = FR.priors,
#             data     = df,
#             # cores    = 4, # parallel computation of chains
#             stanvars = stanvar(scode = Type2H_fix_code, block = "functions")
# )
# expose_functions(fit.1, vectorize=TRUE)

## ----eval=FALSE---------------------------------------------------------------
# summary(fit.1)

## ----echo=FALSE---------------------------------------------------------------
cat("  Family: poisson 
   Links: mu = identity 
 Formula: NE ~ Type2H_fix(N0, 1, 7, a, h) 
          a ~ 1
          h ~ 1
    Data: df (Number of observations: 16) 
   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
          total post-warmup draws = 4000
 
 Regression Coefficients:
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
 a_Intercept     0.65      0.36     0.26     1.59 1.00     1035      936
 h_Intercept     0.70      0.18     0.37     1.06 1.00     1003      874
 
 Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
 and Tail_ESS are effective sample size measures, and Rhat is the potential
 scale reduction factor on split chains (at convergence, Rhat = 1).
")

## ----eval=FALSE---------------------------------------------------------------
# plot(fit.1)

## ----echo=FALSE, output=FALSE-------------------------------------------------
p1 = readRDS(
    system.file("extdata", "Tutorial_01_plot1.rds", package = "BayesFR")
)
p1[[1]]

## ----eval=FALSE---------------------------------------------------------------
# plot(conditional_effects(fit.1), points=TRUE)

## ----echo=FALSE, message=FALSE, warning=FALSE---------------------------------
p2 = readRDS(
    system.file("extdata", "Tutorial_01_plot2.rds", package = "BayesFR")
)
plot(p2, points=TRUE)

## ----eval=FALSE---------------------------------------------------------------
# plot(conditional_effects(fit.1,
#                          effects="N0",
#                          int_conditions=data.frame(N0=seq(0.01,12,length.out=100))),
#      points=TRUE, point_args=list(width=0.1, height=0.1, alpha=0.6, size=2))

## ----echo=FALSE, message=FALSE, warning=FALSE---------------------------------
p3 = readRDS(
    system.file("extdata", "Tutorial_01_plot3.rds", package = "BayesFR")
)
plot(p3, points=TRUE, point_args=list(width=0.1, height=0.1, alpha=0.6, size=2))

## -----------------------------------------------------------------------------
data(df_Hossie_and_Murray_2010_OECOLOGIA)
df = subset(df_Hossie_and_Murray_2010_OECOLOGIA, ID=="Figure 1e")
head(df)
ggplot(aes(N0,NE), data=df) +
  geom_jitter(width=0.1, height=0.0, alpha=0.6, size=2) +
  coord_cartesian(xlim=c(0,NA), ylim=c(0,NA))

## -----------------------------------------------------------------------------
FR.formula = bf( NE | trials(N0) ~ Type3H_dyn(N0,1.0,1.0,b,h)/N0,
                 b~1, h~1,
                 nl = TRUE)

## -----------------------------------------------------------------------------
FR.priors  = c(prior(exponential(1.0), nlpar="b", lb=0),
               prior(exponential(1.0), nlpar="h", lb=0))

## ----eval=FALSE---------------------------------------------------------------
# fit.1 = brm(FR.formula,
#             family   = binomial(link="identity"),
#             prior    = FR.priors,
#             data     = df,
#             # cores    = 4, # parallel computation of chains
#             stanvars = stanvar(scode=Type3H_dyn_code, block="functions")
# )
# expose_functions(fit.1, vectorize=TRUE)

## ----eval=FALSE---------------------------------------------------------------
# summary(fit.1)

## ----echo=FALSE---------------------------------------------------------------
cat("  Family: binomial 
   Links: mu = identity 
 Formula: NE | trials(N0) ~ Type3H_dyn(N0, 1, 1, b, h)/N0 
          b ~ 1
          h ~ 1
    Data: df (Number of observations: 29) 
   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
          total post-warmup draws = 4000
 
 Regression Coefficients:
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
 b_Intercept     0.26      0.06     0.16     0.41 1.00     1317     1495
 h_Intercept     0.05      0.00     0.04     0.05 1.00     1207     1572
 
 Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
 and Tail_ESS are effective sample size measures, and Rhat is the potential
 scale reduction factor on split chains (at convergence, Rhat = 1).
")

## ----eval=FALSE---------------------------------------------------------------
# plot(fit.1)

## ----echo=FALSE, output=FALSE-------------------------------------------------
p1 = readRDS(
    system.file("extdata", "Tutorial_02_plot1.rds", package = "BayesFR")
)
p1[[1]]

## ----eval=FALSE---------------------------------------------------------------
# plot(conditional_effects(fit.1), points=TRUE)

## ----echo=FALSE, message=FALSE, warning=FALSE---------------------------------
p2 = readRDS(
    system.file("extdata", "Tutorial_02_plot2.rds", package = "BayesFR")
)
plot(p2, points=TRUE)

## ----eval=FALSE---------------------------------------------------------------
# plot(conditional_effects(fit.1,
#                          effects="N0",
#                          int_conditions=data.frame(N0=1:60)),
#      points=TRUE, point_args=list(width=0.1, height=0.1, alpha=0.6, size=2))

## ----echo=FALSE, message=FALSE, warning=FALSE---------------------------------
p3 = readRDS(
    system.file("extdata", "Tutorial_02_plot3.rds", package = "BayesFR")
)
plot(p3, points=TRUE, point_args=list(width=0.1, height=0.1, alpha=0.6, size=2))