## -----------------------------------------------------------------------------
#| label: setup
#| include: false

required_pkgs <- c("glmmTMB", "lme4", "nlme")
have_all <- all(vapply(
  required_pkgs,
  requireNamespace,
  logical(1),
  quietly = TRUE
))
knitr::opts_chunk$set(eval = have_all)


## -----------------------------------------------------------------------------
#| label: simulated-data

library(seqwrap)
library(glmmTMB)


# Simulate data
dat <- seqwrap::simcounts(n_genes = 1000, 
                    n_samples = 50, 
                    beta_0 = 2, 
                    sigma_0 = 0.5, 
                    beta_1 = 0.5,, 
                    sigma_1 = 2, 
                    b_0 = 1,
                    clusters = 25,
                    overdispersion_min = c(1, 10))

# Calculate library sizes
dat$metadata$libsize <- log(colSums(dat$data[, -1]) )




## -----------------------------------------------------------------------------

swobject <- seqwrap_compose(
  data = dat$data, # A data frame with the first column as the target identifier
  metadata = dat$metadata, # Meta data bout samples, containing predictors
  samplename = "sample", # The name sample id column 
  modelfun = glmmTMB::glmmTMB, # The model fitting function
  # A list of argumnets supplied to the modelfunction
  arguments = list(
    formula = y ~ x + (1 | cluster),
    family = glmmTMB::nbinom1
  ),
  summary_fun = NULL,
  eval_fun = NULL,
  additional_vars = NULL
)



## -----------------------------------------------------------------------------

models <- seqwrap(swobject, 
                  return_models = TRUE,
                  subset = 1:10) 



## -----------------------------------------------------------------------------

modelsum <- seqwrap_summarise(models)



## -----------------------------------------------------------------------------
seqwrap::generic_summary


## -----------------------------------------------------------------------------
seqwrap::generic_evaluation


## -----------------------------------------------------------------------------

lib_sizes <- colSums(dat$data[,-1]) / 1e6

cpm <- sweep(dat$data[,-1], MARGIN = 2, 
             STATS = lib_sizes, 
             FUN = "/")

dat$data <- data.frame(cbind(dat$data[,1], cpm))



## -----------------------------------------------------------------------------

swobject <- seqwrap_compose(
  data = dat$data, # A data frame with the first column as the target identifier
  metadata = dat$metadata, # Meta data bout samples, containing predictors
  samplename = "sample", # The name sample id column 
  modelfun = lme4::lmer, # The model fitting function
  # A list of argumnets supplied to the modelfunction
  arguments = list(
    formula = y ~ x + (1 | cluster)
  ),
  summary_fun = NULL,
  eval_fun = NULL,
  additional_vars = NULL
)



## -----------------------------------------------------------------------------

models <- seqwrap(swobject, 
                  return_models = TRUE,
                  subset = 1:10) 



## -----------------------------------------------------------------------------

modelsum <- seqwrap_summarise(models)

modelsum$summaries



## -----------------------------------------------------------------------------

models_lme <- seqwrap(swobject, 
                  modelfun = nlme::lme,
                  
                  arguments = alist(fixed = y ~ x, 
                                    random = ~ 1|cluster), 
                  return_models = TRUE,
                  eval_fun = NULL,
                  subset = 1:10)

modelsum_lme <- seqwrap_summarise(models_lme)




## -----------------------------------------------------------------------------
lme_est <- data.frame(modelsum_lme$summaries[, c(1,4,5, 6,9)])
lmer_est <- data.frame(modelsum$summaries[, c(1,4,5, 6)])


## -----------------------------------------------------------------------------
#| echo: false
#| fig-cap: "Similar estimates between `lme4::lmer` and `nlme::lme` models"
#| label: fig-complme
plot(lme_est[lme_est$term == "x","estimate"], 
     lmer_est[lmer_est$term == "x","estimate"], 
     xlab = "nlme::lme estimate", 
     ylab = "lme4::lmer estimate")
abline(a = 0, b = 1)



## -----------------------------------------------------------------------------

models_gls <- seqwrap(swobject, 
                  modelfun = nlme::gls,
                  arguments = alist(
                  model = y ~ x,
                  correlation = nlme::corCompSymm(form = ~ 1|cluster) ),
                  additional_vars = "cluster",
                  return_models = TRUE,
                  eval_fun = NULL,
                  subset = 1:10)

modelsum_gls <- seqwrap_summarise(models_gls)



## -----------------------------------------------------------------------------


gls_est <- data.frame(modelsum_gls$summaries[, c(1,2,3,4,6)])




## -----------------------------------------------------------------------------
#| echo: false
#| fig-cap: "Similar estimates between `lme4::lmer` and `nlme::lme` models"
#| label: fig-compgls
plot(-log10(lme_est[lme_est$term == "x","p.value"]), 
     -log10(gls_est[gls_est$term == "x","p.value"]), 
     xlab = "nlme::lme p-value", 
     ylab = "nlme::gls p-value")
abline(a = 0, b = 1)