## ----include = FALSE----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(
  collapse = FALSE, 
  comment = '#>'
)

old = options(width = 200)

## ----warning=F, message=F-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library(dplyr)
library(tidyr)
library(glmmTMB)
library(powerutilities)
library(emmeans)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
ex1_trt_dict = expand.grid(Insecticide = c('Water', 'Spinosad', 'Sesame Oil'), 
                           Application = c('Curative', 'Prophylactic'))
ex1_trt_dict

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
ex1_trt_dict$Mortality = c(.05, 0.20, 0.15, 
                           .05, 0.35, 0.25) 

ex1a = ex1_trt_dict |> 
  merge(data.frame(Block = factor(1:4))) |> 
  mutate(Ttl_Thrips = 10)

ex1b = ex1_trt_dict |> 
  merge(data.frame(Block = factor(1:8))) |> 
  mutate(Ttl_Thrips = 5)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
ex1a_mod = set_glmm(Mortality ~ Insecticide*Application + (1|Block), 
                    data = ex1a, re_terms = c(0.3), REML = F,
                    family = binomial(link = 'logit'), 
                    weights = ex1a$Ttl_Thrips)

VarCorr(ex1a_mod)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# A warning is produced about non-integer number of successes, but we 
# can safely ignore this.
ex1b_mod = set_glmm(Mortality ~ Insecticide*Application + (1|Block), 
                  data = ex1b, family = binomial(link = 'logit'), 
                  weights = ex1b$Ttl_Thrips, REML = F,
                  re_terms = 0.3)

VarCorr(ex1b_mod)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
emmeans(ex1a_mod, ~ Insecticide:Application, type = 'response')

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
power_ftest(ex1a_mod, ddf = 'asymptotic')
power_ftest(ex1b_mod, ddf = 'asymptotic')

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
power_ftest(ex1a_mod, ddf = 'df.residual')
power_ftest(ex1b_mod, ddf = 'df.residual')

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
(ex1a_emm = emmeans(ex1a_mod, ~Insecticide:Application, type = 'response'))

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
(ex1b_emm = emmeans(ex1b_mod, ~Insecticide:Application, type = 'response'))

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
ex1_contr_list = list('Sesame Curative - Control' = c(-1, 0, 1, 0, 0, 0), 
                      'Sesame Prophylactic - Curative' = c(0, 0, -1, 0, 0, 1))

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
power_contrast(ex1a_emm, ex1_contr_list)

## ----include = FALSE----------------------------------------------------------
## Notes for future development of the vignette

# Example 2: A Randomized Complete Block Design and Gamma-Distributed Data
# Forthcoming

options(old)