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

## ----proposed-data------------------------------------------------------------
library(sommer)
library(Matrix)
data(DT_example, package="enhancer")
DT <- DT_example
DT$Env <- droplevels(factor(DT$Env))
DT$Name <- droplevels(factor(DT$Name))

weight_center <- mean(DT$Weight)
DT$x <- DT$Weight - weight_center

yield_target <- 12
DT$HighYield <- as.integer(DT$Yield >= yield_target)

data.frame(
  observations=nrow(DT),
  individuals=nlevels(DT$Name),
  environments=nlevels(DT$Env),
  yield_weight_correlation=cor(DT$Yield, DT$Weight)
)
with(DT, table(Env, HighYield))

## ----proposed-relationship----------------------------------------------------
A <- A_example
name_levels <- levels(DT$Name)
stopifnot(
  !is.null(rownames(A)),
  !is.null(colnames(A)),
  all(name_levels %in% rownames(A)),
  all(name_levels %in% colnames(A))
)
A <- A[name_levels, name_levels, drop=FALSE]
stopifnot(isSymmetric(A), all(is.finite(A)))
invisible(chol(A))

Ai <- solve(A)
Ai <- as(as(as(Ai, "dMatrix"), "generalMatrix"), "CsparseMatrix")
attr(Ai, "inverse") <- TRUE

data.frame(
  min_diagonal=min(diag(A)),
  max_diagonal=max(diag(A)),
  max_absolute_off_diagonal=max(abs(A[row(A) != col(A)]))
)

## ----proposed-dg-declaration, eval=FALSE--------------------------------------
# ~ vsm(dsm(Env), ism(Name), Gu=Ai)

## ----proposed-methods, eval=FALSE---------------------------------------------
# library(sommer)
# 
# exact_ml <- mmes_method(objective="ml", integration="exact")
# laplace_ml <- mmes_method(objective="ml", integration="laplace")
# 
# fit_control <- mmes_control(
#   optimizer=list(maxit=1000L),
#   mode=list(maxit=100L, tol=1e-8)
# )

## ----proposed-linear, eval=FALSE----------------------------------------------
# fit_lm <- mmes(
#   Yield ~ Env + x,
#   family=gaussian(link="identity"),
#   data=DT,
#   method=exact_ml,
#   control=fit_control
# )

## ----proposed-linear-mixed, eval=FALSE----------------------------------------
# fit_lmm <- mmes(
#   Yield ~ Env + x,
#   random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
#   rcov=~ units,
#   family=gaussian(link="identity"),
#   data=DT,
#   method=exact_ml,
#   control=fit_control
# )

## ----proposed-generalized-linear, eval=FALSE----------------------------------
# fit_glm <- mmes(
#   HighYield ~ Env + x,
#   family=binomial(link="logit"),
#   data=DT,
#   method=exact_ml,
#   control=fit_control
# )

## ----proposed-generalized-linear-mixed, eval=FALSE----------------------------
# fit_glmm <- mmes(
#   HighYield ~ Env + x,
#   random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
#   family=binomial(link="logit"),
#   data=DT,
#   method=laplace_ml,
#   control=fit_control
# )

## ----proposed-nonlinear, eval=FALSE-------------------------------------------
# fit_nlm <- mmes(
#   Yield ~ Env + a * tanh(k * x),
#   parameters=list(
#     a=parameter(fixed=~ 1, random=NULL, link="identity"),
#     k=parameter(fixed=~ 1, random=NULL, link="log")
#   ),
#   family=gaussian(link="identity"),
#   data=DT,
#   start=c(a=3, k=1),
#   method=exact_ml,
#   control=fit_control
# )

## ----proposed-nonlinear-mixed-additive, eval=FALSE----------------------------
# fit_nlmm_additive <- mmes(
#   Yield ~ Env + a * tanh(k * x),
#   parameters=list(
#     a=parameter(fixed=~ 1, link="identity"),
#     k=parameter(fixed=~ 1, link="log")
#   ),
#   random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
#   rcov=~ units,
#   family=gaussian(link="identity"),
#   data=DT,
#   start=c(a=3, k=1),
#   method=exact_ml,
#   control=fit_control
# )

## ----proposed-nonlinear-mixed-rate, eval=FALSE--------------------------------
# fit_nlmm_rate <- mmes(
#   Yield ~ Env + a * tanh(k * x),
#   parameters=list(
#     a=parameter(fixed=~ 1, link="identity"),
#     k=parameter(
#       fixed=~ 1,
#       random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
#       link="log"
#     )
#   ),
#   rcov=~ units,
#   family=gaussian(link="identity"),
#   data=DT,
#   start=c(a=3, k=1),
#   method=laplace_ml,
#   control=fit_control
# )

## ----proposed-nonlinear-mixed-both, eval=FALSE--------------------------------
# fit_nlmm_both <- mmes(
#   Yield ~ Env + a * tanh(k * x),
#   parameters=list(
#     a=parameter(fixed=~ 1, link="identity"),
#     k=parameter(
#       fixed=~ 1,
#       random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
#       link="log"
#     )
#   ),
#   random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
#   rcov=~ units,
#   family=gaussian(link="identity"),
#   data=DT,
#   start=c(a=3, k=1),
#   method=laplace_ml,
#   control=fit_control
# )

## ----proposed-nonlinear-environment-rate, eval=FALSE--------------------------
# k=parameter(
#   fixed=~ Env,
#   random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
#   link="log"
# )

## ----proposed-generalized-nonlinear, eval=FALSE-------------------------------
# fit_gnlm <- mmes(
#   HighYield ~ Env + a * tanh(k * x),
#   parameters=list(
#     a=parameter(fixed=~ 1, link="identity"),
#     k=parameter(fixed=~ 1, link="log")
#   ),
#   family=binomial(link="logit"),
#   data=DT,
#   start=c(a=3, k=1),
#   method=exact_ml,
#   control=fit_control
# )

## ----proposed-generalized-nonlinear-mixed, eval=FALSE-------------------------
# fit_gnlmm <- mmes(
#   HighYield ~ Env + a * tanh(k * x),
#   parameters=list(
#     a=parameter(fixed=~ 1, link="identity"),
#     k=parameter(
#       fixed=~ 1,
#       random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
#       link="log"
#     )
#   ),
#   family=binomial(link="logit"),
#   data=DT,
#   start=c(a=3, k=1),
#   method=laplace_ml,
#   control=fit_control
# )

## ----proposed-gamma-nonlinear-mixed, eval=FALSE-------------------------------
# fit_gamma_nlmm <- mmes(
#   Yield ~ Env + a * tanh(k * x),
#   parameters=list(
#     a=parameter(fixed=~ 1, link="identity"),
#     k=parameter(fixed=~ 1, link="log")
#   ),
#   random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
#   family=Gamma(link="log"),
#   data=DT,
#   start=c(a=3, k=1),
#   method=laplace_ml,
#   control=fit_control
# )

## ----proposed-prediction, eval=FALSE------------------------------------------
# predict(fit_gnlmm, newdata=DT, type="response", target="conditional")
# predict(fit_gnlmm, newdata=DT, type="response", target="zero")
# predict(fit_gnlmm, newdata=DT, type="response", target="population")