#+TITLE: Overview of the package LMMstar
#+Author: Brice Ozenne
#+BEGIN_SRC R :exports none :results output :session *R* :cache no
options(width = 100)
if(system("whoami",intern=TRUE)=="bozenne"){  
  setwd("~/Documents/GitHub/LMMstar/inst/doc-software/")
}else if(system("whoami",intern=TRUE)=="unicph\\hpl802"){  
  setwd("c:/Users/hpl802/Documents/Github/LMMstar/inst/doc-software/")
}
library(ggplot2, quietly = TRUE, verbose = FALSE, warn.conflicts = FALSE)
library(emmeans, quietly = TRUE, verbose = FALSE, warn.conflicts = FALSE)
library(qqtest, quietly = TRUE, verbose = FALSE, warn.conflicts = FALSE)
library(multcomp, quietly = TRUE, verbose = FALSE, warn.conflicts = FALSE)
library(ggpubr, quietly = TRUE, verbose = FALSE, warn.conflicts = FALSE)
library(LMMstar, quietly = TRUE, verbose = FALSE)
#+END_SRC

#+RESULTS:

This vignette describes the main functionalities of the *LMMstar*
package. This package implements specific types of linear mixed models
mainly useful when having repeated observations over a discrete
variable (e.g. time, brain region, ...). Key assumptions are that at
the cluster level, observation are independent and identically
distributed and that the mean and variance are independent
(conditionally on covariates). In particular, in large samples the
residuals do not have to be normally distributed.

\bigskip

The *LMMstar* package contains four main functions:
- the function =lmm= is the main function of the package which fits
  linear mixed models. The user can interact with /lmm/ objects using:
    + =anova= to test linear combinations of coefficients (Wald test
      or Likelihood ratio tests). Can be combined via =rbind=.
    + =coef= to extract the estimates.
    + =dummy.coef= to extract the estimated (marginal) mean for each combination of categorical covariate.
    + =estimate= to test non-linear combinations of coefficients (Wald test via a first order delta method).
    + =levels= to extract the reference level for the mean structure.
      (i.e. what =(Intercept)= refers to in presence of categorical.
      covariates).
    + =getVarCov= to extract the modeled residual variance covariance matrix.
    + =logLik= to output the log-likelihood of the estimated model.
    + =model.tables= to extract table containing estimates with the corresponding uncertainty.
    + =plot= to obtain a diagnostic plots, partial residual plots, or a graphical display of the fitted values.
    + =predict= to compute the conditional mean for new observations.
    + =residuals= to extract the observed residuals of the fitted model.
    + =summary= to obtain a summary of the input, model fit, and estimated values.
- the =summarize= function to compute summary statistics stratified on a categorical variable (typically time).
- the =sampleRem= function to simulate longitudinal data.
- the =LMMstar.options= function enables the user to display the
  default values used in the *LMMstar* package. The function
  can also change the default values to better match the user needs.

\clearpage

Before going further we need to load the *LMMstar* package in the R
session:
#+BEGIN_SRC R  :results silent   :exports code  :session *R* :cache no
library(LMMstar)
#+END_SRC

To illustrate the functionalities of the package, we will use the
gastricbypass dataset:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
data(gastricbypassL, package = "LMMstar")
head(gastricbypassL)
#+END_SRC

#+RESULTS:
:   id visit            time weight glucagonAUC
: 1  1     1 3 months before  127.2     5032.50
: 2  2     1 3 months before  165.2    12142.50
: 3  3     1 3 months before  109.7    10321.35
: 4  4     1 3 months before  146.2     6693.00
: 5  5     1 3 months before  113.1     7090.50
: 6  6     1 3 months before  158.8    10386.00

See =?gastricbypassL= for a presentation of the database. We will use a shorter version of the time variable:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
gastricbypassL$time <- factor(gastricbypassL$time,
                              levels = c("3 months before", "1 week before",
                                         "1 week after", "3 months after" ),
                              labels = c("B3_months","B1_week","A1_week","A3_months"))
gastricbypassL$visit <- as.numeric(gastricbypassL$time) ## convert to numeric
gastricbypassL$baseline <- gastricbypassL$visit<=2
#+END_SRC
#+RESULTS:
rescale the glucagon values
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
gastricbypassL$glucagon <- as.double(scale(gastricbypassL$glucagonAUC))+5
#+END_SRC

#+RESULTS:

and add a group variable:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
gastricbypassL$group <- as.numeric(gastricbypassL$id)%%2
#+END_SRC

#+RESULTS:

\bigskip

_Note:_ the *LMMstar* package is under active development. Newer
package versions may include additional functionalities and fix
previous bugs. The version of the package that is being used is:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
utils::packageVersion("LMMstar")
#+END_SRC

#+RESULTS:
: [1] '0.5.3'

When estimating model coefficients, we will use the internal
optimization routine of the *LMMstar* package (instead of relying on
the =nlme::gls= function, which is the default option):
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
LMMstar.options(optimizer = "FS")
#+END_SRC

#+RESULTS:

\clearpage

* Descriptive statistics
Mean, standard deviation, and other summary statistic can be computed
with respect to a categorical variable (typically time) using the
=summarize= function:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
sss <- summarize(weight+glucagon ~ time, data = gastricbypassL, na.rm = TRUE)
print(sss, digits = 3)
#+END_SRC

#+RESULTS:
:    outcome      time observed missing   mean     sd    min median    max
: 1   weight B3_months       20       0 128.97 20.269 100.90 123.10 173.00
: 2   weight   B1_week       20       0 121.24 18.910  95.70 114.50 162.20
: 3   weight   A1_week       20       0 115.70 18.275  89.90 110.60 155.00
: 4   weight A3_months       20       0 102.36 17.054  78.80  98.50 148.00
: 5 glucagon B3_months       20       0   4.51  0.641   3.61   4.33   6.03
: 6 glucagon   B1_week       19       1   4.39  0.558   3.58   4.23   5.95
: 7 glucagon   A1_week       19       1   6.06  1.044   4.52   5.94   8.27
: 8 glucagon A3_months       20       0   5.06  0.760   3.95   5.03   7.12

Correlation matrices are also ouput when a cluster and ordering
variable have been specified (here respectively =id= and =time=):
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
sss <- summarize(weight ~ time|id, data = gastricbypassL, na.rm = TRUE)
print(sss, digits = 3)
#+END_SRC

#+RESULTS:
#+begin_example
  outcome      time observed missing mean   sd   min median max
1  weight B3_months       20       0  129 20.3 100.9  123.1 173
2  weight   B1_week       20       0  121 18.9  95.7  114.5 162
3  weight   A1_week       20       0  116 18.3  89.9  110.6 155
4  weight A3_months       20       0  102 17.1  78.8   98.5 148

 Pearson's correlation: 
          B3_months B1_week A1_week A3_months
B3_months     1.000   0.990   0.986     0.946
B1_week       0.990   1.000   0.997     0.959
A1_week       0.986   0.997   1.000     0.966
A3_months     0.946   0.959   0.966     1.000
#+end_example

\clearpage

* Linear mixed model
** Covariance patterns
Fit a linear model with *identity* structure:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eId.lmm <- lmm(weight ~ time + glucagon,
               repetition = ~time|id, structure = "ID",
               data = gastricbypassL)
eId.lmm
cat(" covariance structure: \n");getVarCov(eId.lmm)
#+END_SRC

#+RESULTS:
#+begin_example
     Linear regression 

 outcome/cluster/time: weight/id/time 
 data                : 78 observations and distributed in 20 clusters 
 parameters          : 5 mean ((Intercept) timeB1_week timeA1_week timeA3_months glucagon) 
                       1 variance (sigma) 
 log-restr.likelihood: -323.086426918519 
 convergence         : TRUE (1 iterations)
 covariance structure: 
          B3_months  B1_week  A1_week A3_months
B3_months  330.0427   0.0000   0.0000    0.0000
B1_week      0.0000 330.0427   0.0000    0.0000
A1_week      0.0000   0.0000 330.0427    0.0000
A3_months    0.0000   0.0000   0.0000  330.0427
#+end_example

Fit a linear model with *independence* structure:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eInd.lmm <- lmm(weight ~ time + glucagon,
               repetition = ~time|id, structure = "IND",
               data = gastricbypassL)
eInd.lmm
cat(" covariance structure: \n");getVarCov(eInd.lmm)
#+END_SRC

#+RESULTS:
#+begin_example
     Linear regression with heterogeneous residual variance 

 outcome/cluster/time: weight/id/time 
 data                : 78 observations and distributed in 20 clusters 
 parameters          : 5 mean ((Intercept) timeB1_week timeA1_week timeA3_months glucagon) 
                       4 variance (sigma k.B1_week k.A1_week k.A3_months) 
 log-restr.likelihood: -321.457830361849 
 convergence         : TRUE (9 iterations)
 covariance structure: 
          B3_months  B1_week  A1_week A3_months
B3_months  442.6475   0.0000   0.0000    0.0000
B1_week      0.0000 418.9934   0.0000    0.0000
A1_week      0.0000   0.0000 222.8463    0.0000
A3_months    0.0000   0.0000   0.0000  237.2049
#+end_example

\clearpage

Fit a linear mixed model with *compound symmetry* structure:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eCS.lmm <- lmm(weight ~ time + glucagon,
               repetition = ~time|id, structure = "CS",
               data = gastricbypassL)
eCS.lmm
cat(" covariance structure: \n");getVarCov(eCS.lmm)
#+END_SRC

#+RESULTS:
#+begin_example
     Linear Mixed Model with a compound symmetry covariance matrix 

 outcome/cluster/time: weight/id/time 
 data                : 78 observations and distributed in 20 clusters 
 parameters          : 5 mean ((Intercept) timeB1_week timeA1_week timeA3_months glucagon) 
                       1 variance (sigma) 
                       1 correlation (rho) 
 log-restr.likelihood: -243.600523870253 
 convergence         : TRUE (10 iterations)
 covariance structure: 
          B3_months  B1_week  A1_week A3_months
B3_months  355.3062 344.6236 344.6236  344.6236
B1_week    344.6236 355.3062 344.6236  344.6236
A1_week    344.6236 344.6236 355.3062  344.6236
A3_months  344.6236 344.6236 344.6236  355.3062
#+end_example

\clearpage

\noindent Fit a linear mixed model with *stratified compound symmetry* covariance matrix:

#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eSCS.lmm <- lmm(weight ~ time*group,
                repetition = group~time|id, structure = "CS",
                data = gastricbypassL)
eSCS.lmm
cat(" covariance structure: \n");getVarCov(eSCS.lmm)
#+END_SRC

#+RESULTS:
#+begin_example
     Linear Mixed Model with a stratified compound symmetry covariance matrix 

 outcome/cluster/time: weight/id/time 
 data                : 80 observations and distributed in 20 clusters 
 parameters          : 8 mean ((Intercept) timeB1_week timeA1_week timeA3_months group timeB1_week:group timeA1_week:group timeA3_months:group) 
                       2 variance (sigma:0 sigma:1) 
                       2 correlation (rho:0 rho:1) 
 log-restr.likelihood: -233.141302306302 
 convergence         : TRUE (7 iterations)
 covariance structure: 
$`0`
          B3_months  B1_week  A1_week A3_months
B3_months  348.0783 334.7404 334.7404  334.7404
B1_week    334.7404 348.0783 334.7404  334.7404
A1_week    334.7404 334.7404 348.0783  334.7404
A3_months  334.7404 334.7404 334.7404  348.0783

$`1`
          B3_months  B1_week  A1_week A3_months
B3_months  345.1388 340.0877 340.0877  340.0877
B1_week    340.0877 345.1388 340.0877  340.0877
A1_week    340.0877 340.0877 345.1388  340.0877
A3_months  340.0877 340.0877 340.0877  345.1388
#+end_example

\clearpage

\noindent Fit a linear mixed model with *block compound symmetry*
covariance matrix[fn::same as nested random effects]:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eBCS.lmm <- lmm(weight ~ time*group,
                repetition = ~time|id, structure = CS(~baseline, heterogeneous = FALSE),
                data = gastricbypassL)
eBCS.lmm
cat(" covariance structure: \n");getVarCov(eBCS.lmm)
#+END_SRC

#+RESULTS:
#+begin_example
     Linear Mixed Model with a block compound symmetry covariance matrix 

 outcome/cluster/time: weight/id/time 
 data                : 80 observations and distributed in 20 clusters 
 parameters          : 8 mean ((Intercept) timeB1_week timeA1_week timeA3_months group timeB1_week:group timeA1_week:group timeA3_months:group) 
                       1 variance (sigma) 
                       2 correlation (rho(TRUE,TRUE) rho(TRUE,FALSE)) 
 log-restr.likelihood: -234.971305082514 
 convergence         : TRUE (7 iterations)
 covariance structure: 
          B3_months  B1_week  A1_week A3_months
B3_months  346.6085 339.4747 336.3836  336.3836
B1_week    339.4747 346.6085 336.3836  336.3836
A1_week    336.3836 336.3836 346.6085  339.4747
A3_months  336.3836 336.3836 339.4747  346.6085
#+end_example

#+BEGIN_SRC R :exports none :results output :session *R* :cache no
library(lme4)
logLik(lmer(weight ~ time*group + (1|id/baseline), data = gastricbypassL))
#+END_SRC

#+RESULTS:
: 'log Lik.' -234.9713 (df=11)

\noindent Fit a linear mixed model with *block unstructured* covariance matrix:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eBUN.lmm <- lmm(weight ~ time*group,
                repetition = ~time|id, structure = CS(~baseline),
                data = gastricbypassL)
eBUN.lmm
cat(" covariance structure: \n");getVarCov(eBUN.lmm)
#+END_SRC

#+RESULTS:
#+begin_example
     Linear Mixed Model with a block unstructured covariance matrix 

 outcome/cluster/time: weight/id/time 
 data                : 80 observations and distributed in 20 clusters 
 parameters          : 8 mean ((Intercept) timeB1_week timeA1_week timeA3_months group timeB1_week:group timeA1_week:group timeA3_months:group) 
                       2 variance (sigma k.TRUE) 
                       3 correlation (rho(TRUE,TRUE) rho(TRUE,FALSE) rho(FALSE,FALSE)) 
 log-restr.likelihood: -231.80588606934 
 convergence         : TRUE (7 iterations)
 covariance structure: 
          B3_months  B1_week  A1_week A3_months
B3_months  377.4267 372.4602 336.3836  336.3836
B1_week    372.4602 377.4267 336.3836  336.3836
A1_week    336.3836 336.3836 315.7904  306.4892
A3_months  336.3836 336.3836 306.4892  315.7904
#+end_example

\clearpage

\noindent Fit a linear mixed model with *unstructured* covariance matrix:

#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eUN.lmm <- lmm(weight ~ time + glucagon,
               repetition = ~time|id, structure = "UN",
               data = gastricbypassL)
eUN.lmm
cat(" covariance structure: \n");getVarCov(eUN.lmm)
#+END_SRC

#+RESULTS:
#+begin_example
     Linear Mixed Model with an unstructured covariance matrix 

 outcome/cluster/time: weight/id/time 
 data                : 78 observations and distributed in 20 clusters 
 parameters          : 5 mean ((Intercept) timeB1_week timeA1_week timeA3_months glucagon) 
                       4 variance (sigma k.B1_week k.A1_week k.A3_months) 
                       6 correlation (rho(B3_months,B1_week) rho(B3_months,A1_week) rho(B3_months,A3_months) rho(B1_week,A1_week) rho(B1_week,A3_months) rho(A1_week,A3_months)) 
 log-restr.likelihood: -216.318937004305 
 convergence         : TRUE (23 iterations)
 covariance structure: 
          B3_months  B1_week  A1_week A3_months
B3_months  411.3114 381.9734 352.6400  318.8573
B1_week    381.9734 362.7326 335.4649  304.6314
A1_week    352.6400 335.4649 311.6921  285.8077
A3_months  318.8573 304.6314 285.8077  280.9323
#+end_example

\clearpage

\noindent Fit a linear mixed model with *stratified unstructured* covariance matrix:

#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eSUN.lmm <- lmm(weight ~ time*group + glucagon,
                repetition = group~time|id, structure = "UN",
                data = gastricbypassL)
eSUN.lmm
cat(" covariance structure: \n");getVarCov(eSUN.lmm)
#+END_SRC

#+RESULTS:
#+begin_example
     Linear Mixed Model with a stratified unstructured covariance matrix 

 outcome/cluster/time: weight/id/time 
 data                : 78 observations and distributed in 20 clusters 
 parameters          : 9 mean ((Intercept) timeB1_week timeA1_week timeA3_months group glucagon timeB1_week:group timeA1_week:group timeA3_months:group) 
                       8 variance (sigma:0 sigma:1 k.B1_week:0 k.A1_week:0 k.A3_months:0 k.B1_week:1 k.A1_week:1 k.A3_months:1) 
                       12 correlation (rho(B3_months,B1_week):0 rho(B3_months,A1_week):0 rho(B3_months,A3_months):0 rho(B1_week,A1_week):0 rho(B1_week,A3_months):0 rho(A1_week,A3_months):0 rho(B3_months,B1_week):1 rho(B3_months,A1_week):1 rho(B3_months,A3_months):1 rho(B1_week,A1_week):1 rho(B1_week,A3_months):1 rho(A1_week,A3_months):1) 
 log-restr.likelihood: -197.171312062212 
 convergence         : TRUE (51 iterations)
 covariance structure: 
$`0`
          B3_months  B1_week  A1_week A3_months
B3_months  417.3374 382.8829 362.5674  301.7430
B1_week    382.8829 364.4515 346.4039  292.7507
A1_week    362.5674 346.4039 331.1789  282.9301
A3_months  301.7430 292.7507 282.9301  253.3324

$`1`
          B3_months  B1_week  A1_week A3_months
B3_months  383.8877 363.6405 336.5771  350.0416
B1_week    363.6405 347.9898 321.5908  331.5182
A1_week    336.5771 321.5908 297.5329  308.1345
A3_months  350.0416 331.5182 308.1345  334.8267
#+end_example

\clearpage

** Model output

The =summary= method can be used to display the main information
relative to the model fit:
#+BEGIN_SRC R :exports code :results none :session *R* :cache no
summary(eUN.lmm)
#+END_SRC

#+BEGIN_SRC R :exports results :results output :session *R* :cache no
summary(eUN.lmm, hide.mean = TRUE)
#+END_SRC
#+RESULTS:
#+begin_example
           Linear Mixed Model 
 
Dataset: gastricbypassL 

  - 20 clusters 
  - 78 observations were analyzed, 2 were excluded because of missing values 
  - between 3 and 4 observations per cluster 

Summary of the outcome and covariates: 

    $ weight  : num  127 165 110 146 113 ...
    $ time    : Factor w/ 4 levels "B3_months","B1_week",..: 1 1 1 1 1 1 1 1 1 1 ...
    $ glucagon: num  4.03 5.24 4.93 4.32 4.38 ...
    reference level: time=B3_months 

Estimation procedure 

  - Restricted Maximum Likelihood (REML) 
  - log-likelihood :-216.3189
  - parameters: mean = 5, variance = 4, correlation = 6
  - convergence: TRUE (23 iterations) 
    largest |score| = 7.034631e-05 for k.A1_week
            |change|= 1.097373342418e-06 for (Intercept)
 
Residual variance-covariance: unstructured 

  - correlation structure: ~time 
              B3_months B1_week A1_week A3_months
    B3_months     1.000   0.989   0.985     0.938
    B1_week       0.989   1.000   0.998     0.954
    A1_week       0.985   0.998   1.000     0.966
    A3_months     0.938   0.954   0.966     1.000

  - variance structure: ~time 
              standard.deviation     ratio
    B3_months           20.28081 1.0000000
    B1_week             19.04554 0.9390916
    A1_week             17.65480 0.8705176
    A3_months           16.76104 0.8264480
#+end_example

\clearpage

#+BEGIN_SRC R :exports results :results output :session *R* :cache no
oo <- capture.output(summary(eUN.lmm, hide.fit = TRUE, hide.data = TRUE, hide.cor = TRUE, hide.var = TRUE, hide.sd = TRUE))
cat(sapply(oo[-(1:2)],paste0,"\n"))
#+END_SRC

#+RESULTS:
#+begin_example
Fixed effects: weight ~ time + glucagon 
 
               estimate    se     df   lower   upper p.value    
 (Intercept)     132.98 4.664 19.758 123.243 142.717 < 0.001 ***
 timeB1_week     -7.882 0.713 19.171  -9.374   -6.39 < 0.001 ***
 timeA1_week    -11.788 1.018 21.644   -13.9  -9.676 < 0.001 ***
 timeA3_months  -26.122 1.656  18.84 -29.591 -22.654 < 0.001 ***
 glucagon        -0.888 0.242 13.708  -1.408  -0.369 0.00257  **
 
 Uncertainty was quantified using model-based standard errors (column se). 
 Degrees of freedom were computed using a Satterthwaite approximation (column df). 
 The columns lower and upper indicate a 95% confidence interval for each coefficient.
#+end_example

_Note:_ the calculation of the degrees of freedom, especially when
using the observed information can be quite slow. Setting the
arguments =df= to =FALSE= and =type.information= to ="expected"= when
calling =lmm= should lead to a more reasonnable computation time.

** Extract estimated coefficients
The value of the estimated coefficients can be output using =coef=:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
coef(eUN.lmm)
#+END_SRC

#+RESULTS:
:   (Intercept)   timeB1_week   timeA1_week timeA3_months      glucagon 
:   132.9801355    -7.8822331   -11.7879545   -26.1223908    -0.8883081

Variance coefficients can be output by specifying the =effects= argument:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
coef(eUN.lmm, effects = "variance")
#+END_SRC

#+RESULTS:
:       sigma   k.B1_week   k.A1_week k.A3_months 
:  20.2808131   0.9390916   0.8705176   0.8264480


It is possible to apply specific transformation on the variance
coefficients, for instance to obtain the residual variance relative to
each outcome:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
coef(eUN.lmm, effects = "variance", transform.k = "sd")
#+END_SRC

#+RESULTS:
: sigma:B3_months   sigma:B1_week   sigma:A1_week sigma:A3_months 
:        20.28081        19.04554        17.65480        16.76104

The marginal means at each timepoint can be obtained using =dummy.coef=:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
dummy.coef(eUN.lmm)
#+END_SRC

#+RESULTS:
:        time estimate       se       df     lower    upper
: 1 B3_months 128.5386 4.536445 18.97584 119.04289 138.0343
: 2   B1_week 120.6564 4.261691 19.04078 111.73783 129.5749
: 3   A1_week 116.7506 3.956964 19.04925 108.47007 125.0312
: 4 A3_months 102.4162 3.747908 19.05531  94.57328 110.2591

** Extract estimated coefficient and associated uncertainty

The uncertainty about the mean coefficients can be obtained using the
=model.tables= method [fn:: it is equivalent to =confint= method
except that by default it also outputs =se= and =p.value=]:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
model.tables(eUN.lmm)
#+END_SRC

#+RESULTS:
:               estimate    se   df  lower   upper  p.value
: (Intercept)    132.980 4.664 19.8 123.24 142.717 0.00e+00
: timeB1_week     -7.882 0.713 19.2  -9.37  -6.390 9.27e-10
: timeA1_week    -11.788 1.018 21.6 -13.90  -9.676 9.55e-11
: timeA3_months  -26.122 1.656 18.8 -29.59 -22.654 2.62e-12
: glucagon        -0.888 0.242 13.7  -1.41  -0.369 2.57e-03

Values for the all correlation parameters can be displayed
too, by specifying @@latex:\texttt{effect="all"}@@:
#+BEGIN_SRC R :exports code :results output :session *R* :cache no
model.tables(eUN.lmm, effect = "all") ## not shown
#+END_SRC

#+RESULTS:
#+begin_example
                         estimate     se   df   lower   upper  p.value
(Intercept)               132.980 4.6642 19.8 123.243 142.717 0.00e+00
timeB1_week                -7.882 0.7132 19.2  -9.374  -6.390 9.27e-10
timeA1_week               -11.788 1.0175 21.6 -13.900  -9.676 9.55e-11
timeA3_months             -26.122 1.6564 18.8 -29.591 -22.654 2.62e-12
glucagon                   -0.888 0.2416 13.7  -1.408  -0.369 2.57e-03
sigma                      20.281 0.1622 17.9  14.423  28.519       NA
k.B1_week                   0.939 0.0342 19.3   0.874   1.009 0.00e+00
k.A1_week                   0.871 0.0408 20.3   0.800   0.948 0.00e+00
k.A3_months                 0.826 0.0797 19.5   0.700   0.976 3.94e-12
rho(B3_months,B1_week)      0.989 0.2298 32.8   0.972   0.996 7.78e-13
rho(B3_months,A1_week)      0.985 0.2307 26.3   0.961   0.994 5.78e-11
rho(B3_months,A3_months)    0.938 0.2303 23.6   0.847   0.976 1.15e-07
rho(B1_week,A1_week)        0.998 0.2355 27.0   0.994   0.999 3.73e-14
rho(B1_week,A3_months)      0.954 0.2302 24.7   0.886   0.982 1.78e-08
rho(A1_week,A3_months)      0.966 0.2296 27.9   0.915   0.987 1.45e-09
Note: estimates and confidence intervals for sigma, k, rho have been back-transformed. 
      standard errors are not back-transformed.
#+end_example

Because these parameters are constrained (e.g. strictly positive),
they uncertainty is by default computed after transformation
(e.g. =log=) and then backtransformed. 

** Extract estimated residual variance-covariance structure

The method =getVarCov= can be used to output the covariance structure of the residuals:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
getVarCov(eUN.lmm)
#+END_SRC

#+RESULTS:
:           B3_months  B1_week  A1_week A3_months
: B3_months  411.3114 381.9734 352.6400  318.8573
: B1_week    381.9734 362.7326 335.4649  304.6314
: A1_week    352.6400 335.4649 311.6921  285.8077
: A3_months  318.8573 304.6314 285.8077  280.9323

It can also be specific to a "known" individual:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
getVarCov(eUN.lmm, individual = 5)
#+END_SRC

#+RESULTS:
:           B3_months  A1_week A3_months
: B3_months  411.3114 352.6400  318.8573
: A1_week    352.6400 311.6921  285.8077
: A3_months  318.8573 285.8077  280.9323

or for a new individual:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
newdata <- data.frame(id = "X", time = c("B3_months","B1_week","A1_week","A3_months"))
getVarCov(eUN.lmm, individual = newdata)
#+END_SRC

#+RESULTS:
:           B3_months  B1_week  A1_week A3_months
: B3_months  411.3114 381.9734 352.6400  318.8573
: B1_week    381.9734 362.7326 335.4649  304.6314
: A1_week    352.6400 335.4649 311.6921  285.8077
: A3_months  318.8573 304.6314 285.8077  280.9323


\clearpage

** Model diagnostic

The method =plot= can be used to display diagnostic plots about:
- misspecification of the mean structure
#+BEGIN_SRC R :file ./figures/diag-scatterplot.pdf :results graphics file :session *R* :cache no
plot(eUN.lmm, type = "scatterplot")
#+END_SRC

#+RESULTS:
[[file:./figures/diag-scatterplot.pdf]]

#+ATTR_LaTeX: :width 0.4\textwidth :placement [!h]
[[./figures/diag-scatterplot.pdf]]

- misspecification of the variance structure
#+BEGIN_SRC R :file ./figures/diag-scatterplot2.pdf :results graphics file :session *R* :cache no
plot(eUN.lmm, type = "scatterplot2")
#+END_SRC

#+RESULTS:
[[file:./figures/diag-scatterplot2.pdf]]

#+ATTR_LaTeX: :width 0.4\textwidth :placement [!h]
[[./figures/diag-scatterplot2.pdf]]

\clearpage

- misspecification of the correlation structure

#+BEGIN_SRC R :exports code :results output :session *R* :cache no
plot(eUN.lmm, type = "correlation", type.residual = "response")
plot(eUN.lmm, type = "correlation", type.residual = "normalized")
#+END_SRC

#+RESULTS:

#+BEGIN_SRC R :exports none :results output raw drawer :session *R* :cache no
library(ggpubr)
gg <- ggarrange(plot(eUN.lmm, type = "correlation", type.residual = "response", plot = FALSE)$plot,
                plot(eUN.lmm, type = "correlation", type.residual = "normalized", plot = FALSE)$plot,
                common.legend = TRUE)
ggsave(gg, filename = "./figures/diag-correlation.pdf", width = 12)
#+END_SRC


#+RESULTS:
:results:
Saving 12 x 6.38 in image
:end:

#+ATTR_LaTeX: :width 0.6\textwidth :placement [!h]
[[./figures/diag-correlation.pdf]]

- residual distribution vs. normal distribution [fn::see cite:oldford2016self for guidance
  about how to read quantile-quantile plots.]:
#+BEGIN_SRC R :file ./figures/diag-qqplot.pdf :results graphics file :session *R* :cache no
plot(eUN.lmm, type = "qqplot", engine.qqplot = "qqtest")
## Note: the qqtest package to be installed to use the argument engine.plot = "qqtest" 
#+END_SRC

#+RESULTS:
[[file:./figures/diag-qqplot.pdf]]

#+ATTR_LaTeX: :width 0.5\textwidth :placement [!h]
[[./figures/diag-qqplot.pdf]]

\clearpage

The method =residuals= returns the residulas in the wide format:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eUN.diagW <- residuals(eUN.lmm, type = "normalized", format = "wide")
colnames(eUN.diagW) <- gsub("normalized.","",colnames(eUN.diagW))
head(eUN.diagW)
#+END_SRC

#+RESULTS:
:   cluster r.B3_months  r.B1_week   r.A1_week r.A3_months
: 1       1  -0.2897365 -0.2027622 -1.16864038   0.3258573
: 2       2   0.8603117 -1.6492164  0.62578801   1.7370660
: 3       3   0.7273066 -0.4155171 -0.68266741  -0.8510316
: 4       4  -1.6403082 -0.5128368  0.06806206   1.1725813
: 5       5   0.4755409         NA -0.18736415  -0.8634200
: 6       6   1.7801675  1.2847703  2.63004812   0.3505542

or in the long format:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eUN.diagL <- residuals(eUN.lmm, type = "normalized", format = "long")
head(eUN.diagL)
#+END_SRC

#+RESULTS:
: [1] -0.2897365  0.8603117  0.7273066 -1.6403082  0.4755409  1.7801675

Various type of residuals can be extract but the normalized one are
recommanded when doing model checking.

** Model fit

The fitted values can be displayed via the =plot= method or using the =emmeans= package:

#+BEGIN_SRC R :file ./figures/fit-autoplot.pdf :results graphics file :session *R* :cache no
library(ggplot2) ## left panel
plot(eUN.lmm, type = "fit", color = "id", ci.alpha = NA, size.text = 20)
#+END_SRC

#+RESULTS:
[[file:./figures/fit-autoplot.pdf]]

#+BEGIN_SRC R :file ./figures/fit-emmip.pdf :results graphics file :session *R* :cache no
library(emmeans) ## right panel
emmip(eUN.lmm, ~time) + theme(text = element_text(size=20))
#+END_SRC

#+RESULTS:
[[file:./figures/fit-emmip.pdf]]

#+latex: \begin{minipage}{0.45\linewidth}
#+ATTR_LaTeX: :width \textwidth :placement [!h]
[[./figures/fit-autoplot.pdf]]
#+latex: \end{minipage}
#+latex: \begin{minipage}{0.45\linewidth}
#+ATTR_LaTeX: :width \textwidth :placement [!h]
[[./figures/fit-emmip.pdf]]
#+latex: \end{minipage}

# ## ggsave(emmip(eUN.lmm, ~time) + theme(text = element_text(size=20)), filename = "figures/fit-emmip.pdf")
# ## ggsave(autoplot(eUN.lmm, color = "id", plot = FALSE)$plot + theme(text = element_text(size=20)), filename = "figures/fit-autoplot.pdf")

In the first case each possible curve is displayed while in the latter
the average curve (over glucagon values). With the =plot= method,
it is possible to display a curve specific to a glucagon value via the
argument =at=:
#+BEGIN_SRC R :exports code :results output :session *R* :cache no
plot(eUN.lmm, type = "fit", at = data.frame(glucagon = 10), color = "glucagon")
#+END_SRC

#+RESULTS:

It is also possible to display the observed values along with the
fitted values by setting the argument =obs.alpha= to a strictly
positive value below or equal to 1. This argument controls the
transparency of the color used to display the observed values:
#+BEGIN_SRC R :exports code :results output :session *R* :cache no
gg <- plot(eUN.lmm, type = "fit", obs.alpha = 0.2, ci = FALSE,plot = FALSE)$plot
gg <- gg + facet_wrap(~id, labeller = label_both)
gg <- gg + theme(axis.text.x=element_text(angle = 90, hjust = 0))
gg
#+END_SRC

#+RESULTS:
: Advarselsbesked:
: Removed 2 rows containing missing values (geom_point).

# ## ggsave(gg + theme(text = element_text(size=20)), filename = "figures/fit-autoplot-indiv.pdf", width = 12)

#+ATTR_LaTeX: :width \textwidth :placement [!h]
[[./figures/fit-autoplot-indiv.pdf]]

\clearpage

Partial residuals can also be displayed via the =plot= method:
#+header: :width 12 :height 5
#+BEGIN_SRC R :file ./figures/fit-pres.pdf :results graphics file :session *R* :cache no
gg1 <- plot(eUN.lmm, type = "partial", var = "glucagon", plot = FALSE)$plot
gg2 <- plot(eUN.lmm, type = "partial", var = c("(Intercept)","glucagon"), plot = FALSE)$plot
ggarrange(gg1,gg2)
#+END_SRC

#+RESULTS:
[[file:./figures/fit-pres.pdf]]

#+ATTR_LaTeX: :width \textwidth :placement [!h]
[[./figures/fit-pres.pdf]]

Their value can be extracted via the =residuals= method, e.g.:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
df.pres <- residuals(eUN.lmm, type = "partial", var = "glucagon", keep.data = TRUE)
m.pres <- gastricbypassL$weight - model.matrix(~time,gastricbypassL) %*% coef(eUN.lmm)[1:4]
range(df.pres$r.partial - m.pres, na.rm = TRUE)
#+END_SRC

#+RESULTS:
: [1] -1.065814e-14  1.154632e-14

#+BEGIN_SRC R :exports none :results output :session *R* :cache no
df.pres2 <- residuals(eUN.lmm, type = "partial-center", var = "glucagon", keep.data = TRUE)
m.pres2 <- m.pres - mean(gastricbypassL$glucagon, na.rm = TRUE) * coef(eUN.lmm)["glucagon"]
range(df.pres2$r.partial - m.pres2, na.rm = TRUE)
#+END_SRC

#+RESULTS:
: [1] -1.421085e-14  1.110223e-14

\clearpage

** Statistical inference (linear)

The =anova= method can be use to test one or several linear
combinations of the model coefficients using Wald tests. By default,
it will simultaneously test all parameters associated to a variable:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
anova(eUN.lmm)
#+END_SRC

#+RESULTS:
: 
: 	|| mean coefficients || 
: 
:  - Multivariate Wald test (global null hypothesis)
:          statistic df.num df.denom    p.value    
: time        86.743      3   19.005 2.8424e-11 ***
: glucagon    13.518      1   13.708  0.0025716  **
: ---
:  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note that here the p-values are not adjust for multiple comparisons
over variables. It is possible to specify a null hypothesis to be
test: e.g. is there a change in average weight just after taking the
treatment:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
anova(eUN.lmm, effects = c("timeA1_week-timeB1_week=0"))
#+END_SRC

#+RESULTS:
: 
: 	|| User-specified linear hypotheses || 
: 
:  - Multivariate Wald test (global null hypothesis)
:  statistic df.num df.denom    p.value    
:     43.141      1   17.875 3.7234e-06 ***
: ---
:  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

One can also simulateneously tests several null hypotheses:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
e.anova <- anova(eUN.lmm, effects = c("timeA1_week-timeB1_week=0",
                                      "timeA3_months-timeB1_week=0"))
summary(e.anova)
#+END_SRC

#+RESULTS:
#+begin_example

	|| User-specified linear hypotheses || 

 - Multivariate Wald test (global null hypothesis)
 statistic df.num df.denom    p.value    
    98.651      2    18.62 1.2338e-10 ***

 - Univariate Wald test (individual null hypotheses) 
                             estimate        se        df     lower    upper p.value    
timeA1_week - timeB1_week    -3.90572   0.59464  17.87453  -5.31903  -2.4924   2e-05 ***
timeA3_months - timeB1_week -18.24016   1.32283  19.02810 -21.38419 -15.0961  <1e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Standard errors: model-based
(CIs/p-values adjusted for multiple comparisons -- max-test adjustment)
CIs/p-values computed using 1e+05 samples.
#+end_example

or return all pairwise comparisons for a given factor using the =mcp=
function of the multcomp package:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
library(multcomp)
summary(anova(eUN.lmm, effects = mcp(time = "Tukey")))
#+END_SRC

#+RESULTS:
#+begin_example
Singular contrast matrix: contrasts "A1_week - B1_week" "A3_months - B1_week" "A3_months - A1_week" have been removed. 


	|| User-specified linear hypotheses || 

 - Multivariate Wald test (global null hypothesis)
 statistic df.num df.denom    p.value    
    86.743      3   19.005 2.8424e-11 ***

 - Univariate Wald test (individual null hypotheses) 
                       estimate        se        df     lower    upper p.value    
B1_week - B3_months    -7.88223   0.71318  19.17147  -9.81006  -5.9544  <1e-05 ***
A1_week - B3_months   -11.78795   1.01751  21.64404 -14.53843  -9.0375  <1e-05 ***
A3_months - B3_months -26.12239   1.65641  18.84049 -30.59989 -21.6449  <1e-05 ***
A1_week - B1_week      -3.90572   0.59464  17.87453  -5.51312  -2.2983   2e-05 ***
A3_months - B1_week   -18.24016   1.32283  19.02810 -21.81595 -14.6644  <1e-05 ***
A3_months - A1_week   -14.33444   1.05650  20.26658 -17.19030 -11.4786  <1e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Standard errors: model-based
(CIs/p-values adjusted for multiple comparisons -- max-test adjustment)
CIs/p-values computed using 1e+05 samples.
#+end_example


When testing transformed variance or correlation parameters,
parentheses (as in =log(k).B1_week=) cause problem for recognizing
parameters:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
try(
  anova(eUN.lmm,
        effects = c("log(k).B1_week=0","log(k).A1_week=0","log(k).A3_months=0"))
)
#+END_SRC

#+RESULTS:
: Error in .anova_Wald(object, effects = effects, robust = robust, rhs = rhs,  : 
:   Possible mispecification of the argument 'effects' as running mulcomp::glht lead to the following error: 
: Error in parse(text = ex[i]) : <text>:1:7: uventet symbol
: 1: log(k).B1_week
:           ^

It is then advised to build a contrast matrix, e.g.:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
name.coef <- rownames(confint(eUN.lmm, effects = "all"))
name.varcoef <- grep("^k",name.coef, value = TRUE)
C <- matrix(0, nrow = 3, ncol = length(name.coef), dimnames = list(name.varcoef, name.coef))
diag(C[name.varcoef,name.varcoef]) <- 1
C
#+END_SRC

#+RESULTS:
#+begin_example
            (Intercept) timeB1_week timeA1_week timeA3_months glucagon sigma k.B1_week k.A1_week
k.B1_week             0           0           0             0        0     0         1         0
k.A1_week             0           0           0             0        0     0         0         1
k.A3_months           0           0           0             0        0     0         0         0
            k.A3_months rho(B3_months,B1_week) rho(B3_months,A1_week) rho(B3_months,A3_months)
k.B1_week             0                      0                      0                        0
k.A1_week             0                      0                      0                        0
k.A3_months           1                      0                      0                        0
            rho(B1_week,A1_week) rho(B1_week,A3_months) rho(A1_week,A3_months)
k.B1_week                      0                      0                      0
k.A1_week                      0                      0                      0
k.A3_months                    0                      0                      0
#+end_example

And then call the =anova= method specifying the null hypothesis via the
contrast matrix:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
anova(eUN.lmm, effects = C)
#+END_SRC

#+RESULTS: 
: 
: 	|| User-specified linear hypotheses || 
: 
:  - Multivariate Wald test (global null hypothesis)
:  statistic df.num df.denom   p.value   
:     6.2032      3   17.995 0.0044171 **
: ---
:  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note that using the approach of cite:pipper2012versatile it is also
possible to adjust for multiple testing across several =lmm=
objects. To do so, one first fit the mixed models, then use the
=anova= method to indicate which hypotheses are being tested, and
combine them using =rbind=.

\clearpage

Here is a (very artificial) example:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
Manova <- rbind(anova(eInd.lmm, effects = "glucagon = 0"),
                anova(eCS.lmm, effects = "glucagon = 0"),
                anova(eUN.lmm, effects = "glucagon = 0"))
summary(Manova) 
#+END_SRC

#+RESULTS:
#+begin_example

	|| User-specified linear hypotheses || 

 - Multivariate Wald test (global null hypothesis)
 statistic df.num df.denom    p.value    
    8.8925      3      Inf 6.8788e-06 ***

 - Univariate Wald test (individual null hypotheses) 
      estimate        se        df     lower   upper p.value   
[1,]  -8.27006   2.57880  34.20071 -14.88175 -1.6584 0.01173 * 
[2,]   0.82179   0.61997  53.80983  -0.76772  2.4113 0.47205   
[3,]  -0.88831   0.24161  13.70759  -1.50776 -0.2689 0.00399 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Standard errors: model-based
(CIs/p-values adjusted for multiple comparisons -- max-test adjustment)
CIs/p-values computed using 1e+05 samples.
#+end_example

\clearpage

** Statistical inference (non-linear)

The =estimate= function can be used to test one or several non-linear
combinations of model coefficients, using a first order delta method
to quantify uncertainty. The combination has to be specified via a
function (argument =f=). To illustrate its use consider an ANCOVA
analysis:
#+BEGIN_EXPORT latex
\[ Y_{i1} = \textcolor{\darkred}{\alpha} + \textcolor{\darkblue}{\beta} Y_{i,0} + \textcolor{\darkgreen}{\gamma} X_{i} + e_{i} \]
#+END_EXPORT

#+BEGIN_SRC R :exports both :results output :session *R* :cache no
gastricbypassW <- reshape(gastricbypassL[,c("id","time","weight","group")],
                          direction = "wide",
                          timevar = "time", idvar = c("id","group"))
e.ANCOVA <- lm(weight.A1_week ~ weight.B1_week + group, data = gastricbypassW)
summary(e.ANCOVA)$coef
#+END_SRC

#+RESULTS:
:                  Estimate Std. Error    t value     Pr(>|t|)
: (Intercept)    -1.4823022 2.31781138 -0.6395267 5.310047e-01
: weight.B1_week  0.9654917 0.01803988 53.5198489 2.156258e-20
: group           0.2521714 0.66499945  0.3792054 7.092302e-01

We can replicate this analysis by first fitting a mixed model:
#+BEGIN_EXPORT latex
\[ Y_{ij} = \alpha_j + \gamma_j X_{i} + \varepsilon_{i,j} \text{ where } \varepsilon_i \sim \Gaus \left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \sigma^2_1 & \rho \sigma_1 \sigma_2 \\ \rho \sigma_1 \sigma_2 & \sigma^2_2 \end{bmatrix} \right) \]
#+END_EXPORT
#+BEGIN_SRC R :exports code :results output :session *R* :cache no
e.lmmANCOVA <- lmm(weight ~ time+time:group, repetition = ~time|id,
                   data = gastricbypassL[gastricbypassL$visit %in% 2:3,])
#+END_SRC

#+RESULTS:
: Factor variable(s) with empty level: "time"
:  The empty level(s) will be remove internally. Consider applying droplevels to avoid this warning.

and then perform a delta-method:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
lava::estimate(e.lmmANCOVA, f = function(p){
  c(Y1 = as.double(p["rho(B1_week,A1_week)"]*p["k.A1_week"]),
    X1 = as.double(p["timeA1_week:group"]-p["rho(B1_week,A1_week)"]*p["k.A1_week"]*p["timeB1_week:group"]))
})
#+END_SRC

#+RESULTS:
:     estimate         se       df      lower    upper   p.value
: Y1 0.9654917 0.01753161 15.96769  0.9283203 1.002663 0.0000000
: X1 0.2521714 0.64626331 15.00349 -1.1252784 1.629621 0.7018731

Indeed:
#+BEGIN_EXPORT latex
\begin{align*}
\Esp[Y_{i2}|Y_{i1},X_{i}] &= \alpha_2 + \gamma_2 X_{i} + \rho \frac{\sigma_2}{\sigma_1}\left(Y_{i1} - \alpha_1 - \gamma_1 X_{i}\right) \\
                         &= \textcolor{\darkred}{\alpha_2 - \rho \frac{\sigma_2}{\sigma_1} \alpha_1}
                         + \textcolor{\darkblue}{\rho \frac{\sigma_2}{\sigma_1}Y_{i1}}
                         + \textcolor{\darkgreen}{\left(\gamma_2 - \rho \frac{\sigma_2}{\sigma_1} \gamma_1\right)  X_{i} }
\end{align*}
#+END_EXPORT

We obtain identical estimate but different standard-errors/degrees of
freedom compared to the univariate linear model approach. The later is
to be prefer as it does not rely on approximation. The former is
nevertheless useful as it can handle missing data in the outcome
variable.

\clearpage

** Baseline adjustment

The =lmm= contains an "experimental" feature to drop non-identifiable
effects from the model. For instance, let us define two (artifical) groups of
patients:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
gastricbypassL$group <- c("1","2")[as.numeric(gastricbypassL$id) %in% 15:20 + 1]
#+END_SRC
#+RESULTS:
We would like to model group differences only after baseline
(i.e. only at 1 week and 3 months after). For this we will define a
treatment variable being the group variable except before baseline where
it is ="none"=:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
gastricbypassL$treat <- baselineAdjustment(gastricbypassL, variable = "group",
                                           repetition = ~time|id, constrain = c("B3_months","B1_week"),
                                           new.level = "none")
table(treat = gastricbypassL$treat, time = gastricbypassL$time, group = gastricbypassL$group)
#+END_SRC

#+RESULTS:
#+begin_example
, , group = 1

      time
treat  B3_months B1_week A1_week A3_months
  none        14      14       0         0
  1            0       0      14        14
  2            0       0       0         0

, , group = 2

      time
treat  B3_months B1_week A1_week A3_months
  none         6       6       0         0
  1            0       0       0         0
  2            0       0       6         6
#+end_example

Here we will be able to estimate a total of 6 means and therefore can
at most identify 6 effects. However the design matrix for the
interaction model:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
colnames(model.matrix(weight ~ treat*time, data = gastricbypassL))
#+END_SRC

#+RESULTS:
:  [1] "(Intercept)"          "treat1"               "treat2"               "timeB1_week"         
:  [5] "timeA1_week"          "timeA3_months"        "treat1:timeB1_week"   "treat2:timeB1_week"  
:  [9] "treat1:timeA1_week"   "treat2:timeA1_week"   "treat1:timeA3_months" "treat2:timeA3_months"

contains 12 parameters (i.e. 6 too many). The =lmm= function will
internally remove the one that cannot be identified and fit a
simplified model:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eC.lmm <- lmm(weight ~ treat*time, data = gastricbypassL,
              repetition = ~time|id, structure = "UN")
#+END_SRC

#+RESULTS:
: Constant values in the design matrix in interactions "treat:time"
:  Coefficients "treat1" "treat2" "timeA1_week" "timeA3_months" "treat1:timeB1_week" "treat2:timeB1_week" have been removed.

with the following coefficients:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
coef(eC.lmm, effects = "mean")
#+END_SRC

#+RESULTS:
:          (Intercept)          timeB1_week   treat1:timeA1_week   treat2:timeA1_week 
:            128.97000             -7.73000            -12.83949            -14.27452 
: treat1:timeA3_months treat2:timeA3_months 
:            -27.07620            -25.50553

One can vizualize the baseline adjustment via the =autoplot= function:
#+BEGIN_SRC R :file ./figures/gg-baseAdj.pdf :results graphics file :session *R* :cache no
autoplot(eC.lmm, color = "group", ci = FALSE, size.text = 20) 
#+END_SRC

#+RESULTS:
[[file:./figures/gg-baseAdj.pdf]]

#+ATTR_LaTeX: :width 0.4\textwidth :placement [!h]
[[./figures/gg-baseAdj.pdf]]

To more easily compare the two groups, one could set the baseline
treatment to the treatment in the control arm by omitting the argument
=new.level=:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
gastricbypassL$treat2 <- baselineAdjustment(gastricbypassL, variable = "group",
                                            repetition = ~time|id, constrain = c("B3_months","B1_week"))
table(treat = gastricbypassL$treat2, time = gastricbypassL$time, group = gastricbypassL$group)
#+END_SRC

#+RESULTS:
#+begin_example
, , group = 1

     time
treat B3_months B1_week A1_week A3_months
    1        14      14      14        14
    2         0       0       0         0

, , group = 2

     time
treat B3_months B1_week A1_week A3_months
    1         6       6       0         0
    2         0       0       6         6
#+end_example

Fitting the model
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eC2.lmm <- suppressWarnings(lmm(weight ~ treat2*time, data = gastricbypassL,
                                repetition = ~time|id, structure = "UN"))
#+END_SRC

#+RESULTS:
: Constant values in the design matrix in interactions "treat2:time"
:  Coefficients "treat22" "treat22:timeB1_week" have been removed.

will directly output group differences (last two coefficients):
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
model.tables(eC2.lmm)
#+END_SRC
#+RESULTS:
:                       estimate    se   df  lower  upper  p.value
: (Intercept)             128.97 4.532 19.0 119.48 138.46 0.00e+00
: timeB1_week              -7.73 0.697 19.0  -9.19  -6.27 1.00e-09
: timeA1_week             -12.84 0.865 20.5 -14.64 -11.04 2.02e-12
: timeA3_months           -27.08 1.724 21.4 -30.66 -23.50 3.20e-13
: treat22:timeA1_week      -1.44 0.621 16.3  -2.75  -0.12 3.43e-02
: treat22:timeA3_months     1.57 2.463 16.3  -3.64   6.78 5.32e-01

It is also possible to get the estimated mean at each timepoint, using
an equivalent mean structure:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eC3.lmm <- suppressWarnings(lmm(weight ~ 0+treat2:time, data = gastricbypassL,
                                repetition = ~time|id, structure = "UN"))
model.tables(eC3.lmm) ## equivalent to dummy.coef(eC2.lmm)
#+END_SRC

#+RESULTS:
: Constant values in the design matrix in interactions "treat2:time"
:  Coefficients "treat22:timeB3_months" "treat22:timeB1_week" have been removed.
:                       estimate   se   df lower upper p.value
: treat21:timeB3_months      129 4.53 19.0 119.5   138       0
: treat21:timeB1_week        121 4.23 19.0 112.4   130       0
: treat21:timeA1_week        116 4.11 19.1 107.5   125       0
: treat22:timeA1_week        115 4.13 19.4 106.1   123       0
: treat21:timeA3_months      102 3.87 20.2  93.8   110       0
: treat22:timeA3_months      103 4.17 25.2  94.9   112       0

or the baseline mean and the change since baseline:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eC4.lmm <- suppressWarnings(lmm(weight ~ treat2:time, data = gastricbypassL,
                                repetition = ~time|id, structure = "UN"))
model.tables(eC4.lmm)
#+END_SRC

#+RESULTS:
: Constant values in the design matrix in interactions "treat2:time"
:  Coefficients "treat22:timeB3_months" "treat22:timeB1_week" have been removed.
:                       estimate    se   df  lower  upper  p.value
: (Intercept)             128.97 4.532 19.0 119.48 138.46 0.00e+00
: treat21:timeB1_week      -7.73 0.697 19.0  -9.19  -6.27 1.00e-09
: treat21:timeA1_week     -12.84 0.865 20.5 -14.64 -11.04 2.02e-12
: treat22:timeA1_week     -14.27 0.950 26.3 -16.23 -12.32 2.02e-14
: treat21:timeA3_months   -27.08 1.724 21.4 -30.66 -23.50 3.20e-13
: treat22:timeA3_months   -25.51 2.323 22.6 -30.32 -20.69 1.60e-10

** Marginal means

The =emmeans= package can be used to output marginal means. Consider
the following model:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
e.group <- lmm(weight ~ time*group, data = gastricbypassL,
               repetition = ~time|id, structure = "UN")
#+END_SRC

#+RESULTS:

We can for instance compute the average value over time /assuming balanced groups/:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
emmeans(e.group, specs=~time)
#+END_SRC

#+RESULTS:
: NOTE: Results may be misleading due to involvement in interactions
:  time      emmean   SE   df lower.CL upper.CL
:  B3_months    130 5.05 18.0    119.3      141
:  B1_week      122 4.69 18.0    112.5      132
:  A1_week      117 4.55 18.0    107.0      126
:  A3_months    104 4.20 18.1     94.9      113
: 
: Results are averaged over the levels of: group 
: Confidence level used: 0.95

This differs from the average value over time over the whole sample:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
df.pred <- cbind(gastricbypassL, predict(e.group, newdata = gastricbypassL))
summarize(formula = estimate~time, data = df.pred)
#+END_SRC

#+RESULTS:
:    outcome      time observed missing    mean       sd      min   median    max
: 1 estimate B3_months       20       0 128.970 2.270212 127.5214 127.5214 132.35
: 2 estimate   B1_week       20       0 121.240 2.726942 119.5000 119.5000 125.30
: 3 estimate   A1_week       20       0 115.700 2.014981 114.4143 114.4143 118.70
: 4 estimate A3_months       20       0 102.365 3.146729 100.3571 100.3571 107.05

as the groups are not balanced:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
table(group = gastricbypassL$group, time = gastricbypassL$time)
#+END_SRC

#+RESULTS:
:      time
: group B3_months B1_week A1_week A3_months
:     1        14      14      14        14
:     2         6       6       6         6

The "emmeans" approach gives equal "weight" to the expected value of
both group 2:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
mu.group1 <-  as.double(coef(e.group)["(Intercept)"])
mu.group2 <-  as.double(coef(e.group)["(Intercept)"] + coef(e.group)["group2"])
p.group1 <- 14/20          ; p.group2 <- 6/20
c(emmeans = (mu.group1+mu.group2)/2, predict = mu.group1 * p.group1 + mu.group2 * p.group2)
#+END_SRC

#+RESULTS:
:  emmeans  predict 
: 129.9357 128.9700

Which one is relevant depends on the application. The =emmeans=
function can also be used to display expected value in each group over
time:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
emmeans.group <- emmeans(e.group, specs = ~group|time)
emmeans.group
#+END_SRC

#+RESULTS:
#+begin_example
time = B3_months:
 group emmean   SE   df lower.CL upper.CL
 1        128 5.53 18.0    115.9      139
 2        132 8.45 18.0    114.6      150

time = B1_week:
 group emmean   SE   df lower.CL upper.CL
 1        120 5.14 18.0    108.7      130
 2        125 7.85 18.0    108.8      142

time = A1_week:
 group emmean   SE   df lower.CL upper.CL
 1        114 4.99 18.0    103.9      125
 2        119 7.62 18.0    102.7      135

time = A3_months:
 group emmean   SE   df lower.CL upper.CL
 1        100 4.60 18.1     90.7      110
 2        107 7.03 18.1     92.3      122

Confidence level used: 0.95
#+end_example

Using the =pair= function displays the differences:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
epairs.group <- pairs(emmeans.group, reverse = TRUE)
epairs.group
#+END_SRC

#+RESULTS:
#+begin_example
time = B3_months:
 contrast estimate    SE   df t.ratio p.value
 2 - 1        4.83 10.10 18.0   0.478  0.6383

time = B1_week:
 contrast estimate    SE   df t.ratio p.value
 2 - 1        5.80  9.38 18.0   0.618  0.5441

time = A1_week:
 contrast estimate    SE   df t.ratio p.value
 2 - 1        4.29  9.11 18.0   0.471  0.6435

time = A3_months:
 contrast estimate    SE   df t.ratio p.value
 2 - 1        6.69  8.40 18.1   0.797  0.4361
#+end_example

One can adjust for multiple comparison via the =adjust= argument and
display confidence intervals setting the argument =infer= to =TRUE=:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
summary(epairs.group, by = NULL, adjust = "mvt", infer = TRUE)
#+END_SRC

#+RESULTS:
:  contrast time      estimate    SE   df lower.CL upper.CL t.ratio p.value
:  2 - 1    B3_months     4.83 10.10 18.0    -18.0     27.7   0.478  0.7496
:  2 - 1    B1_week       5.80  9.38 18.0    -15.4     27.0   0.618  0.6481
:  2 - 1    A1_week       4.29  9.11 18.0    -16.3     24.9   0.471  0.7551
:  2 - 1    A3_months     6.69  8.40 18.1    -12.3     25.7   0.797  0.5285
: 
: Confidence level used: 0.95 
: Conf-level adjustment: mvt method for 4 estimates 
: P value adjustment: mvt method for 4 tests

This should also work when doing baseline adjustment (because of
baseline adjustment no difference is expected at the first two
timepoints):
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
summary(pairs(emmeans(eC2.lmm , specs = ~treat2|time), reverse = TRUE), by = NULL)
#+END_SRC

#+RESULTS:
: Note: adjust = "tukey" was changed to "sidak"
: because "tukey" is only appropriate for one set of pairwise comparisons
:  contrast time      estimate    SE   df t.ratio p.value
:  2 - 1    B3_months     0.00 0.000  NaN     NaN     NaN
:  2 - 1    B1_week       0.00 0.000  NaN     NaN     NaN
:  2 - 1    A1_week      -1.44 0.621 16.2  -2.311  0.1303
:  2 - 1    A3_months     1.57 2.463 16.3   0.638  0.9522
: 
: P value adjustment: sidak method for 4 tests

** Predictions

Two types of predictions can be performed with the =predict= method:
- *static predictions* that are only conditional on the covariates:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
news <- gastricbypassL[gastricbypassL$id==1,]
news$glucagon <- 0
predict(eUN.lmm, newdata = news)
#+END_SRC

#+RESULTS:
:   estimate       se       df     lower    upper
: 1 132.9801 4.664247 19.75815 123.24305 142.7172
: 2 125.0979 4.388294 19.91418 115.94155 134.2543
: 3 121.1922 4.214230 20.55331 112.41660 129.9678
: 4 106.8577 3.942058 20.95499  98.65871 115.0568

\clearpage

which can be computing by creating a design matrix:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
X.12 <- model.matrix(formula(eUN.lmm), news)
X.12
#+END_SRC

#+RESULTS:
#+begin_example
   (Intercept) timeB1_week timeA1_week timeA3_months glucagon
1            1           0           0             0        0
21           1           1           0             0        0
41           1           0           1             0        0
61           1           0           0             1        0
attr(,"assign")
[1] 0 1 1 1 2
attr(,"contrasts")
attr(,"contrasts")$time
[1] "contr.treatment"
#+end_example

and then multiplying it with the regression coefficients:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
X.12 %*% coef(eUN.lmm)
#+END_SRC

#+RESULTS:
:        [,1]
: 1  132.9801
: 21 125.0979
: 41 121.1922
: 61 106.8577

- *dynamic predictions* that are conditional on the covariates and the
  outcome measured at other timepoints. Consider two subjects for who
  we would like to predict the weight 1 week before the intervention
  based on the weight 3 months before the intervention:
  
#+ATTR_LATEX: :options otherkeywords={}, deletekeywords={}
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
newd <- rbind(
  data.frame(id = 1, time = "B3_months", weight = coef(eUN.lmm)["(Intercept)"], glucagon = 0),
  data.frame(id = 1, time = "B1_week", weight = NA, glucagon = 0),
  data.frame(id = 2, time = "B3_months", weight = 100, glucagon = 0),
  data.frame(id = 2, time = "B1_week", weight = NA, glucagon = 0)
)
predict(eUN.lmm, newdata = newd, type = "dynamic", keep.newdata = TRUE)
#+END_SRC

#+RESULTS:
:   id      time   weight glucagon  estimate        se  df     lower    upper
: 1  1 B3_months 132.9801        0        NA        NA  NA        NA       NA
: 2  1   B1_week       NA        0 125.09790 0.6362754 Inf 123.85083 126.3450
: 3  2 B3_months 100.0000        0        NA        NA  NA        NA       NA
: 4  2   B1_week       NA        0  94.47017 7.2279385 Inf  80.30367 108.6367
  
The first subjects has the average weight while the second has a much
  lower weight. The predicted weight for the first subject is then the
  average weight one week before while it is lower for the second
  subject due to the positive correlation over time. The predicted
  value is computed using the formula of the conditional mean for a
  Gaussian vector:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
mu1 <- coef(eUN.lmm)[1]
mu2 <- sum(coef(eUN.lmm)[1:2])
Omega_11 <- getVarCov(eUN.lmm)["B3_months","B3_months"]
Omega_21 <- getVarCov(eUN.lmm)["B1_week","B3_months"]
as.double(mu2 + Omega_21 * (100 - mu1) / Omega_11)
#+END_SRC

#+RESULTS:
: [1] 94.47017


\clearpage

** Missing values and imputation


We now consider the glucagon level as an outcome. The =summarize=
function can be used to describe the amount of missing data at each
repetition:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
sss <- summarize(glucagon ~ time, data = gastricbypassL, na.rm = TRUE)
cbind(sss[,1:4], pc = paste0(100 * sss$missing / (sss$missing + sss$observed), "%"))
#+END_SRC

#+RESULTS:
:    outcome      time observed missing pc
: 1 glucagon B3_months       20       0 0%
: 2 glucagon   B1_week       19       1 5%
: 3 glucagon   A1_week       19       1 5%
: 4 glucagon A3_months       20       0 0%

Further description of the missing data patterns rely on function
outside the LMMstar package, e.g. appropriate call to =tapply= and
=table=:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
vec.pattern <- tapply(as.numeric(is.na(gastricbypassL$glucagon)),
                      INDEX = gastricbypassL$id,
                      FUN = paste, collapse=".")
table(vec.pattern)
#+END_SRC

#+RESULTS:
: vec.pattern
: 0.0.0.0 0.0.1.0 0.1.0.0 
:      18       1       1

Linear mixed model can handle missing value in the outcome variable,
assuming that missigness is random conditional on the covariate and
observed outcome values. The =lmm= function can be used "as usual":
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eUN.lmmNA <- lmm(glucagon ~ time,
                 repetition = ~time|id, structure = "UN",
                 data = gastricbypassL)
summary(eUN.lmmNA, hide.fit = TRUE,
        hide.cor = TRUE, hide.sd = TRUE, hide.mean = TRUE)
#+END_SRC

#+RESULTS:
#+begin_example
           Linear Mixed Model 
 
Dataset: gastricbypassL 

  - 20 clusters 
  - 78 observations were analyzed, 2 were excluded because of missing values 
  - between 3 and 4 observations per cluster 

Summary of the outcome and covariates: 

    $ glucagon: num  4.03 5.24 4.93 4.32 4.38 ...
    $ time    : Factor w/ 4 levels "B3_months","B1_week",..: 1 1 1 1 1 1 1 1 1 1 ...
    reference level: time=B3_months
#+end_example

The visible difference in the summary is when describing the dataset:
we can see that some repetitions (here 2) have been ignored as the
outcome was missing. So for some clusters only 3 values were analyzed
instead of 4. It is possible to extract the most likely value for
these missing observation using the =fitted= function with argument
=impute=TRUE=:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
fitted(eUN.lmmNA, impute = TRUE)
#+END_SRC

#+RESULTS:
: [1] 4.256984 6.497856

When using the argument =keep.newdata=TRUE=, the missing outcome value
has been replaced by its most likely value (which is the same as the
dynamic prediction, describedy previously):
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
eData <- fitted(eUN.lmmNA, impute = TRUE, keep.newdata = TRUE)
eData[eData$id %in% eData[eData$imputed,"id"],]
#+END_SRC

#+RESULTS:
:    id visit      time weight glucagon baseline group treat treat2 imputed estimate
: 5   5     1 B3_months  113.1 4.383738     TRUE     1  none      1   FALSE 4.514352
: 15 15     1 B3_months  115.0 4.098741     TRUE     2  none      1   FALSE 4.514352
: 25  5     2   B1_week  105.6 4.256984     TRUE     1  none      1    TRUE 4.386762
: 35 15     2   B1_week  109.7 4.509697     TRUE     2  none      1   FALSE 4.386762
: 45  5     3   A1_week   99.9 6.430376    FALSE     1     1      1   FALSE 6.078985
: 55 15     3   A1_week  103.5 6.497856    FALSE     2     2      2    TRUE 6.078985
: 65  5     4 A3_months   87.7 5.275118    FALSE     1     1      1   FALSE 5.057642
: 75 15     4 A3_months   94.1 6.259632    FALSE     2     2      2   FALSE 5.057642

Visually:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
ggplot(eData, aes(x=time,y=glucagon, group=id)) + geom_line() + geom_point(aes(color=imputed))
#+END_SRC

#+RESULTS:

#+BEGIN_SRC R :exports none :results output raw drawer :session *R* :cache no
gg <- ggplot(eData, aes(x=time,y=glucagon, group=id)) + geom_line() + geom_point(aes(color=imputed), size=2)
ggsave(gg+ theme(text = element_text(size=20)), filename = "./figures/imputation.pdf", width = 12)
#+END_SRC
#+RESULTS:
:results:
Saving 12 x 6.38 in image
:end:

#+ATTR_LaTeX: :width 1\textwidth :options trim={0 0 0 0} :placement [!h]
[[./figures/imputation.pdf]]

It is possible to sample from the estimated distribution of the
missing value instead of using the most likely value, e.g. accounting
for residual variance and uncertainty related to parameter estimation:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
set.seed(10)
fitted(eUN.lmmNA, impute = TRUE, se = "total")
fitted(eUN.lmmNA, impute = TRUE, se = "total")
fitted(eUN.lmmNA, impute = TRUE, se = "total")
#+END_SRC

#+RESULTS:
: [1] 4.262434 6.305287
: [1] 3.858267 5.871642
: [1] 4.342624 6.905246

\clearpage

* Data generation
Simulate some data in the wide format:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
set.seed(10) ## ensure reproductibility
n.obs <- 100
n.times <- 4
mu <- rep(0,4)
gamma <- matrix(0, nrow = n.times, ncol = 10) ## add interaction
gamma[,6] <- c(0,1,1.5,1.5)
dW <- sampleRem(n.obs, n.times = n.times, mu = mu, gamma = gamma, format = "wide")
head(round(dW,3))
#+END_SRC

#+RESULTS:
:   id X1 X2 X3 X4 X5     X6     X7     X8    X9    X10     Y1     Y2     Y3     Y4
: 1  1  1  0  1  1  0 -0.367  1.534 -1.894 1.729  0.959  1.791  2.429  3.958  2.991
: 2  2  1  0  1  2  0 -0.410  2.065  1.766 0.761 -0.563  2.500  4.272  3.002  2.019
: 3  3  0  0  2  1  0 -1.720 -0.178  2.357 1.966  1.215 -3.208 -5.908 -4.277 -5.154
: 4  4  0  0  0  1  0  0.923 -2.089  0.233 1.307 -0.906 -2.062  0.397  1.757 -1.380
: 5  5  0  0  2  1  0  0.987  5.880  0.385 0.028  0.820  7.963  7.870  7.388  8.609
: 6  6  0  0  1  1  2 -1.075  0.479  2.202 0.900 -0.739  0.109 -1.602 -1.496 -1.841

Simulate some data in the long format:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
set.seed(10) ## ensure reproductibility
dL <- sampleRem(n.obs, n.times = n.times, mu = mu, gamma = gamma, format = "long")
head(dL)
#+END_SRC

#+RESULTS:
:   id visit        Y X1 X2 X3 X4 X5         X6       X7        X8        X9        X10
: 1  1     1 1.791444  1  0  1  1  0 -0.3665251 1.533815 -1.894425 1.7288665  0.9592499
: 2  1     2 2.428570  1  0  1  1  0 -0.3665251 1.533815 -1.894425 1.7288665  0.9592499
: 3  1     3 3.958350  1  0  1  1  0 -0.3665251 1.533815 -1.894425 1.7288665  0.9592499
: 4  1     4 2.991198  1  0  1  1  0 -0.3665251 1.533815 -1.894425 1.7288665  0.9592499
: 5  2     1 2.500179  1  0  1  2  0 -0.4097541 2.065413  1.765841 0.7613348 -0.5630173
: 6  2     2 4.272357  1  0  1  2  0 -0.4097541 2.065413  1.765841 0.7613348 -0.5630173

\clearpage

* Modifying default options
The =LMMstar.options= method enable to get and set the default options
used by the package. For instance, the default option for the information matrix is:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
LMMstar.options("type.information")
#+END_SRC

#+RESULTS:
: $type.information
: [1] "observed"

To change the default option to "expected" (faster to compute but less accurate p-values and confidence intervals in small samples) use:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
LMMstar.options(type.information = "expected")
#+END_SRC

#+RESULTS:

To restore the original default options do:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
LMMstar.options(reinitialise = TRUE)
#+END_SRC

#+RESULTS:

\clearpage

* R session
Details of the R session used to generate this document:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
sessionInfo()
#+END_SRC

#+RESULTS:
#+begin_example
R version 4.1.1 (2021-08-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=Danish_Denmark.1252  LC_CTYPE=Danish_Denmark.1252    LC_MONETARY=Danish_Denmark.1252
[4] LC_NUMERIC=C                    LC_TIME=Danish_Denmark.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] lme4_1.1-27.1   Matrix_1.4-0    LMMstar_0.5.0   nlme_3.1-153    ggpubr_0.4.0    multcomp_1.4-18
 [7] TH.data_1.1-0   MASS_7.3-54     survival_3.2-13 mvtnorm_1.1-3   qqtest_1.2.0    emmeans_1.7.2  
[13] ggplot2_3.3.5  

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.8          lattice_0.20-45     tidyr_1.1.4         listenv_0.8.0      
 [5] zoo_1.8-9           assertthat_0.2.1    digest_0.6.29       utf8_1.2.2         
 [9] parallelly_1.30.0   R6_2.5.1            plyr_1.8.6          backports_1.4.1    
[13] coda_0.19-4         pillar_1.6.5        rlang_0.4.12        minqa_1.2.4        
[17] nloptr_1.2.2.3      car_3.0-12          textshaping_0.3.6   labeling_0.4.2     
[21] splines_4.1.1       stringr_1.4.0       munsell_0.5.0       broom_0.7.11       
[25] compiler_4.1.1      numDeriv_2016.8-1.1 systemfonts_1.0.3   pkgconfig_2.0.3    
[29] mgcv_1.8-38         globals_0.14.0      tidyselect_1.1.1    gridExtra_2.3      
[33] tibble_3.1.6        codetools_0.2-18    fansi_1.0.2         future_1.23.0      
[37] crayon_1.4.2        dplyr_1.0.7         withr_2.4.3         grid_4.1.1         
[41] xtable_1.8-4        gtable_0.3.0        lifecycle_1.0.1     DBI_1.1.2          
[45] magrittr_2.0.1      scales_1.1.1        estimability_1.3    future.apply_1.8.1 
[49] stringi_1.7.6       carData_3.0-5       farver_2.1.0        ggsignif_0.6.3     
[53] reshape2_1.4.4      ragg_1.2.1          ellipsis_0.3.2      generics_0.1.1     
[57] vctrs_0.3.8         cowplot_1.1.1       boot_1.3-28         sandwich_3.0-1     
[61] lava_1.6.10         tools_4.1.1         glue_1.6.1          purrr_0.3.4        
[65] abind_1.4-5         parallel_4.1.1      colorspace_2.0-2    rstatix_0.7.0
#+end_example

\clearpage

* References
:PROPERTIES:
:UNNUMBERED: t
:END:

#+BEGIN_EXPORT latex
\begingroup
\renewcommand{\section}[2]{}
#+END_EXPORT

bibliographystyle:apalike
[[bibliography:bibliography.bib]]

#+BEGIN_EXPORT latex
\endgroup
#+END_EXPORT

\clearpage

#+BEGIN_EXPORT LaTeX
\appendix
\titleformat{\section}
{\normalfont\Large\bfseries}{Appendix~\thesection}{1em}{}

\renewcommand{\thefigure}{\Alph{figure}}
\renewcommand{\thetable}{\Alph{table}}
\renewcommand{\theequation}{\Alph{equation}}

\setcounter{figure}{0}    
\setcounter{table}{0}    
\setcounter{equation}{0}    
#+END_EXPORT

* Likelihood in a linear mixed model
:PROPERTIES:
:CUSTOM_ID: SM:likelihood
:END:

Denote by \(\VY\) a vector of \(m\) outcomes, \(\VX\) a vector of
\(p\) covariates, \(\mu(\Vparam,\VX)\) the modeled mean, and
\(\Omega(\Vparam,\VX)\) the modeled residual variance-covariance. We
consider \(n\) replicates (i.e. \(\VY_1,\ldots,\VY_n)\) and
\(VX_1,\ldots,\VX_n\)) along with a vector of weights
\(\omega=(w_1,\ldots,w_n)\), which are by default all equal to 1.

** Log-likelihood

The restricted log-likelihood in a linear mixed model can then be
written:
 #+BEGIN_EXPORT LaTeX
\begin{align}
\Likelihood(\Vparam|\VY,\VX) =& \textcolor{\darkred}{ \frac{p}{2} \log(2\pi)-\frac{1}{2} \log\left(\left|\sum_{i=1}^n w_i \VX_i \Omega_i^{-1}(\Vparam) \trans{\VX}_i\right|\right)} \notag \\
& + \sum_{i=1}^{n} w_i \left(\textcolor{\darkblue}{-\frac{m}{2} \log(2\pi) - \frac{1}{2} \log\left|\Omega_i(\Vparam)\right| - \frac{1}{2} (\VY_i-\mu(\Vparam,\VX_i)) \Omega_i(\Vparam)^{-1} \trans{(\VY_i-\mu(\Vparam,\VX_i))}} \right)  \label{eq:log-likelihood}
\end{align}
 #+END_EXPORT
 
 This is what the =logLik= method is computing for the REML
 criteria. The red term is specific to the REML criteria and prevents
 from computing individual contributions to the likelihood[fn::The REML is the
 likelihood of the observations divided by the prior on the estimated
 mean parameters \(\VparamHat_{\mu} \sim \Gaus(\mu,\left(\VX
 \Omega^{-1}(\Vparam) \trans{\VX}\right)^{-1})\). This corresponds to
 \(\frac{1}{\sqrt{2\pi}^p \left|\left(\sum_{i=1}^n \VX_i
 \Omega_i^{-1}(\Vparam) \trans{\VX}_i\right)^{-1}\right|}
 \exp\left(-(\VparamHat_{\mu}-\mu)\left(2\sum_{i=1}^n \VX_i
 \Omega_i^{-1}(\Vparam)
 \trans{\VX}_i\right)^{-1})\trans{(\VparamHat_{\mu}-\mu)}\right)\)
 Since \(\mu\) will be estimated to be \(\Vparam_{\mu}\), the
 exponential term equals 1 and thus does not contribute to the
 log-likelihood. One divided by the other term gives \(\sqrt{2\pi}^p
 \left(\left|\sum_{i=1}^n \VX_i \Omega_i^{-1}(\Vparam)
 \trans{\VX}_i\right|\right)^{-1}\). The log of this term equals the red
 term]. The blue term is what =logLik= outputs for the ML criteria
 when setting the argument =indiv= to =TRUE=.

\bigskip

** Score

 Using that \(\partial \log(\det(X))=tr(X^{-1}\partial(X))\), the
score is obtained by derivating once the log-likelihood, i.e., for
\(\theta \in \Vparam\):
#+BEGIN_EXPORT LaTeX
\begin{align*}
   \Score(\theta) =& \dpartial[\Likelihood(\Vparam|\VY,\VX)][\theta]
= \textcolor{\darkred}{ \frac{1}{2} tr \left( \left(\sum_{i=1}^n w_i \VX_i \Omega_i^{-1}(\Vparam) \trans{\VX}_i\right)^{-1} \left(\sum_{i=1}^n w_i \VX_i \Omega_i^{-1}(\Vparam) \dpartial[\Omega_i(\Vparam)][\theta] \Omega_i(\Vparam)^{-1} \trans{\VX}_i\right)  \right) } \\
&+ \sum_{i=1}^n w_i \left( \textcolor{\darkblue}{ -\frac{1}{2} tr\left(\Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta]\right) + \dpartial[\mu(\Vparam,\VX_i)][\theta] \Omega_i(\Vparam)^{-1} \trans{(\VY_i-\mu(\Vparam,\VX_i))} } \right. \\
 & \qquad \qquad \left. \textcolor{\darkblue}{ + \frac{1}{2} (\VY_i-\mu(\Vparam,\VX_i)) \Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta] \Omega_i(\Vparam)^{-1} \trans{(\VY_i-\mu(\Vparam,\VX_i))} } \right).
\end{align*}
#+END_EXPORT

 This is what the =score= method is computing for the REML
 criteria. The red term is specific to the REML criteria and prevents
 from computing the score relative to each cluster. The blue term is
 what =score= outputs for the ML criteria when setting the argument
 =indiv= to =TRUE=.

\bigskip

\clearpage

** Hessian

Derivating a second time the log-likelihood gives the hessian, \(\Hessian(\Vparam)\), with element[fn::if one is relative to the mean and the other to the variance then they are respectively \(\theta\) and \(\theta'\)]:
#+BEGIN_EXPORT LaTeX
\begin{align*}
& \Hessian(\theta,\theta^{\prime}) = \ddpartial[\Likelihood(\Vparam|\VY,\VX)][\theta][\theta^{\prime}] = \dpartial[\Score(\theta)][\theta^{\prime}] \\
=& \textcolor{\darkred}{\frac{1}{2} tr \left( \left(\sum_{i=1}^n w_i \VX_i \Omega_i^{-1}(\Vparam) \trans{\VX}_i\right)^{-1} \left\{ \sum_{i=1}^n w_i \VX_i \Omega_i^{-1}(\Vparam) \left(\ddpartial[\Omega_i(\Vparam)][\theta][\theta^{\prime}] - 2 \dpartial[\Omega_i(\Vparam)][\theta] \Omega_i^{-1}(\Vparam) \dpartial[\Omega_i(\Vparam)][\theta^{\prime}]\right)\Omega_i(\Vparam)^{-1} \trans{\VX}_i \right.  \right.}  \\
& \textcolor{\darkred}{ \left. \left. + \left(\sum_{i=1}^n w_i \VX_i \Omega_i^{-1}(\Vparam) \dpartial[\Omega_i(\Vparam)][\theta] \Omega_i(\Vparam)^{-1} \trans{\VX}_i\right) \left(\sum_{i=1}^n w_i \VX_i\Omega_i^{-1}(\Vparam) \trans{\VX}_i \right)^{-1} \left(\sum_{i=1}^n w_i \VX_i \Omega_i^{-1}(\Vparam) \dpartial[\Omega_i(\Vparam)][\theta^{\prime}] \Omega_i(\Vparam)^{-1} \trans{\VX}_i\right) \right\} \right) } \\
& +\sum_{i=1}^n w_i \left( \textcolor{\darkblue}{ \frac{1}{2} tr\left(\Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta^{\prime}] \Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta] - \Omega_i(\Vparam)^{-1} \ddpartial[\Omega_i(\Vparam)][\theta][\theta^{\prime}] \right) } \right.\\
& \qquad \textcolor{\darkblue}{ -  \dpartial[\mu(\Vparam,\VX_i)][\theta] \Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta^{\prime}] \Omega_i(\Vparam)^{-1} \trans{\Vvarepsilon_i(\Vparam)} - \dpartial[\mu(\Vparam,\VX_i)][\theta] \Omega_i(\Vparam)^{-1} \trans{\dpartial[\mu(\Vparam,\VX_i)][\theta^{\prime}]} } \\
& \qquad \left. \textcolor{\darkblue}{ + \frac{1}{2} \Vvarepsilon_i(\Vparam) \Omega_i(\Vparam)^{-1} \left(\ddpartial[\Omega_i(\Vparam)][\theta][\theta^{\prime}] - \dpartial[\Omega_i(\Vparam)][\theta^{\prime}] \Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta] - \dpartial[\Omega_i(\Vparam)][\theta] \Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta^{\prime}] \right) \Omega_i(\Vparam)^{-1} \trans{\Vvarepsilon_i(\Vparam)} } \right).
\end{align*}
#+END_EXPORT
where \(\Vvarepsilon_i(\Vparam) = \VY_i-\mu(\Vparam,\VX_i)\).

\bigskip

The =information= method will (by default) return the (observed)
information which is the opposite of the hessian. So multiplying the
previous formula by -1 gives what =information= output for the REML
criteria. The red term is specific to the REML criteria and prevents
from computing the information relative to each cluster. The blue term
is what =information= outputs for the ML criteria (up to a factor -1)
when setting the argument =indiv= to =TRUE=.

\bigskip

A possible simplification is to use the expected hessian at the maximum likelihood. Indeed for
any deterministic matrix \(A\):
- \(\Esp[A \trans{(\VY_i-\mu(\Vparam,\VX_i))}|\VX_i] = 0\)
- \(\Esp[(\VY_i-\mu(\Vparam,\VX_i)) A \trans{(\VY_i-\mu(\Vparam,\VX_i))}||\VX_i] = tr(A \Var(\VY_i-\mu(\Vparam,\VX_i)))\)
when \(\Esp[\VY_i-\mu(\Vparam,\VX_i)]=0\). This leads to:
#+BEGIN_EXPORT LaTeX
\begin{align}
 & \Esp[\Hessian(\theta,\theta^{\prime})|\VX] \notag\\ 
 &= \textcolor{\darkred}{ \frac{1}{2} tr \left( \left(\sum_{i=1}^n w_i \VX_i \Omega_i^{-1}(\Vparam) \trans{\VX}_i\right)^{-1}  \left\{ \sum_{i=1}^n w_i \VX_i \Omega_i^{-1}(\Vparam) \left( \ddpartial[\Omega_i(\Vparam)][\theta][\theta^{\prime}] - 2 \dpartial[\Omega_i(\Vparam)][\theta]  \Omega_i^{-1}(\Vparam) \dpartial[\Omega_i(\Vparam)][\theta^{\prime}]\right) \Omega_i(\Vparam)^{-1} \trans{\VX}_i \right.  \right.} \notag \\
 & \textcolor{\darkred}{ \left. \left. +  \left(\sum_{i=1}^n w_i \VX_i \Omega_i^{-1}(\Vparam) \dpartial[\Omega_i(\Vparam)][\theta] \Omega_i(\Vparam)^{-1} \trans{\VX}_i\right) \left(\sum_{i=1}^n w_i \VX_i \Omega_i^{-1}(\Vparam) \trans{\VX}_i \right)^{-1} \left(\sum_{i=1}^n w_i \VX_i \Omega_i^{-1}(\Vparam) \dpartial[\Omega_i(\Vparam)][\theta^{\prime}] \Omega_i(\Vparam)^{-1} \trans{\VX}_i\right) \right\} \right) } \notag\\
 & + \sum_{i=1}^n w_i \left( \textcolor{\darkblue}{
- \frac{1}{2} tr\left(\Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta^{\prime}] \Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta]\right)
 - \dpartial[\mu(\Vparam,\VX_i)][\theta] \Omega_i(\Vparam)^{-1} \trans{\dpartial[\mu(\Vparam,\VX_i)][\theta^{\prime}]}
 } \right) \label{eq:expectedInfo} 
\end{align}
#+END_EXPORT

This is what =information= output when the argument =type.information=
is set to ="expected"= (up to a factor -1).

\clearpage

** Degrees of freedom

Degrees of freedom are computed using a Satterthwaite approximation,
i.e. for an estimate coefficient \(\widehat{\beta}\in\widehat{\Vparam}\) with standard
error \(\sigma_{\widehat{beta}}\), the degree of freedom is:
#+begin_export latex
\begin{align*}
df\left(\sigma_{\widehat{\beta}}\right) = \frac{2 \sigma^4_{\widehat{\beta}}}{\Var[\widehat{\sigma}_{\widehat{\beta}}]}
\end{align*}
#+end_export
Using a first order Taylor expansion we can approximate the variance term as:
#+begin_export latex
\begin{align*}
\Var[\widehat{\sigma}_{\widehat{\beta}}] & \approx \dpartial[\widehat{\sigma}_{\widehat{\beta}}][\Vparam] \Sigma_{\Vparam}  \trans{\dpartial[\widehat{\sigma}_{\widehat{\beta}}][\Vparam]} \\
& \approx c_{\beta} \left(\widehat{\Information}_{\widehat{\Vparam}}\right)^{-1} \dpartial[\widehat{\Information}_{\widehat{\Vparam}}][\Vparam] \left(\widehat{\Information}_{\widehat{\Vparam}}\right)^{-1} \trans{c_{\beta}} \Sigma_{\Vparam} \trans{c_{\beta}} \left(\widehat{\Information}_{\widehat{\Vparam}}\right)^{-1} \trans{\dpartial[\widehat{\Information}_{\widehat{\Vparam}}][\Vparam]} \left(\widehat{\Information}_{\widehat{\Vparam}}\right)^{-1} c_{\beta}
\end{align*}
#+end_export

  where \(\Sigma_{\Vparam}\) is the variance-covariance matrix of all
  model coefficients, \(\Information_{\Vparam}\) the information
  matrix for all model coefficients, \(c_{\beta}\) a matrix used to
  select the element relative to \(\beta\) in the first derivative of
  the information matrix, and \(\dpartial[.][\Vparam]\) denotes the
  vector of derivatives with respect to all model coefficients.

\bigskip

The derivative of the information matrix (i.e. negative hessian) can
then be computed using numerical derivatives or using analytical
formula. To simplify the derivation of the formula we will only derive
them at the maximum likelihood, i.e. when
\(\Esp\left[\dpartial[\Hessian(\theta,\theta^{\prime}|\VX)][\theta^{\prime\prime}]\right]=\frac{\partial
\Esp[\Hessian(\theta,\theta^{\prime}|\VX)]}{\partial
\theta^{\prime\prime}}\) where the expectation is taken over
\(\VX\). We can therefore take the derivative of formula
eqref:eq:expectedInfo. We first note that its derivative with respect
to the mean parameters is 0. So we just need to compute the derivative
with respect to a variance parameter \(\theta^{\prime\prime}\):
#+BEGIN_EXPORT LaTeX
\begin{align*}
 & \frac{\partial\Esp[\Hessian(\theta,\theta^{\prime})|\VX]}{\partial \theta^{\prime\prime}} \notag\\ 
% &= \textcolor{\darkred}{ \frac{1}{2} tr \left( \left(\sum_{i=1}^n \VX_i \Omega_i^{-1}(\Vparam) \trans{\VX}_i\right)^{-1}  \left\{ \sum_{i=1}^n \VX_i \Omega_i^{-1}(\Vparam) \left( \ddpartial[\Omega_i(\Vparam)][\theta][\theta^{\prime}] - 2 \dpartial[\Omega_i(\Vparam)][\theta]  \Omega_i^{-1}(\Vparam) \dpartial[\Omega_i(\Vparam)][\theta^{\prime}]\right) \Omega_i(\Vparam)^{-1} \trans{\VX}_i \right.  \right.} \notag \\
% & \textcolor{\darkred}{ \left. \left. \qquad +  \left(\sum_{i=1}^n \VX_i \Omega_i^{-1}(\Vparam) \dpartial[\Omega_i(\Vparam)][\theta] \Omega_i(\Vparam)^{-1} \trans{\VX}_i\right) \left(\sum_{i=1}^n \VX_i \Omega_i^{-1}(\Vparam) \trans{\VX}_i \right) \left(\sum_{i=1}^n \VX_i \Omega_i^{-1}(\Vparam) \dpartial[\Omega_i(\Vparam)][\theta^{\prime}] \Omega_i(\Vparam)^{-1} \trans{\VX}_i\right) \right\} \right) } \notag\\
 & + \sum_{i=1}^n w_i \left( \textcolor{\darkblue}{
- \frac{1}{2} tr\left(
-2\Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta^{\prime\prime}] \Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta^{\prime}] \Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta] \right. } \right. \\
& \qquad \qquad \textcolor{\darkblue}{\left. + \Omega_i(\Vparam)^{-1} \ddpartial[\Omega_i(\Vparam)][\theta^{\prime}][\theta^{\prime\prime}] \Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta]
+ \Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta^{\prime}] \Omega_i(\Vparam)^{-1} \ddpartial[\Omega_i(\Vparam)][\theta][\theta^{\prime\prime}]
\right)} \\
& \qquad \qquad  \textcolor{\darkblue}{\left. + \dpartial[\mu(\Vparam,\VX_i)][\theta] \Omega_i(\Vparam)^{-1} \dpartial[\Omega_i(\Vparam)][\theta^{\prime\prime}] \Omega_i(\Vparam)^{-1}   \trans{\dpartial[\mu(\Vparam,\VX_i)][\theta^{\prime}]}
 \right)}
\end{align*}
#+END_EXPORT





\clearpage
  
* Likelihood ratio test with the REML criterion
:PROPERTIES:
:CUSTOM_ID: SM:LRT-REML
:END:

The blue term of autoref:eq:log-likelihood in the log-likelihood is
invariant to re-parameterisation while the red term is not. This means
that a re-parametrisation of \(X\) into \(\tilde{X} = B X\) with \(B\)
invertible would not change the likelihood when using ML but would
decrease the log-likelihood by \(\log(|B|)\) when using REML. 
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
LMMstar.options(optimizer = "FS",
                param.optimizer = c(n.iter = 1000, tol.score = 1e-3, tol.param = 1e-5))
#+END_SRC

#+RESULTS:

\bigskip

Let's take an example:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
## data(gastricbypassL, package = "LMMstar")
dfTest <- gastricbypassL
dfTest$glucagon2 <- dfTest$glucagon*2
#+END_SRC

#+RESULTS:

where we multiply one column of the design matrix by 2. As mentionned
previously this does not affect the log-likelihood when using ML:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
logLik(lmm(weight ~ glucagon, data = dfTest, structure = UN(~time|id), method = "ML"))
logLik(lmm(weight ~ glucagon2, data = dfTest, structure = UN(~time|id), method = "ML"))
#+END_SRC

#+RESULTS:
: [1] -245.7909
: [1] -245.7909

but it does when using REML:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
logLik(lmm(weight ~ glucagon, data = dfTest, structure = UN(~time|id), method = "REML"))
logLik(lmm(weight ~ glucagon2, data = dfTest, structure = UN(~time|id), method = "REML"))
log(2)
#+END_SRC

#+RESULTS:
: [1] -245.0382
: [1] -245.7313
: [1] 0.6931472

Therefore, when comparing models with different mean effects there is
a risk that the difference (or part of it) in log-likelihood is due to
a new parametrisation and no only to a difference in model fit. This
would typically be the case when adding an interaction where we can
have a smaller restricted log-likehood when considering a more complex
model:

#+BEGIN_SRC R :exports both :results output :session *R* :cache no
set.seed(1)
dfTest$ff <- rbinom(NROW(dfTest), size = 1, prob = 0.5)
logLik(lmm(weight ~ glucagon, data = dfTest, structure = UN(~time|id), method = "REML"))
logLik(lmm(weight ~ glucagon*ff, data = dfTest, structure = UN(~time|id), method = "REML"))
#+END_SRC

#+RESULTS:
: [1] -245.0382
: [1] -245.3555

This is quite counter-intuitive as more complex model should lead to
better fit and would never happen when using ML:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
logLik(lmm(weight ~ glucagon, data = dfTest, structure = UN(~time|id), method = "ML"))
logLik(lmm(weight ~ glucagon*ff, data = dfTest, structure = UN(~time|id), method = "ML"))
#+END_SRC

#+RESULTS:
: [1] -245.7909
: [1] -245.3593

This is why, unless one knows what he/she is doing, it is not
recommanded to use likelihood ratio test to assess relevance of mean
parameters in mixed models estimated with REML.

* CONFIG                                                           :noexport:
#+LANGUAGE:  en
#+LaTeX_CLASS: org-article
#+LaTeX_CLASS_OPTIONS: [12pt]
#+OPTIONS:   title:t author:t toc:nil todo:nil
#+OPTIONS:   H:3 num:t 
#+OPTIONS:   TeX:t LaTeX:t
** Display of the document
# ## space between lines
#+LATEX_HEADER: \RequirePackage{setspace} % to modify the space between lines - incompatible with footnote in beamer
#+LaTeX_HEADER:\renewcommand{\baselinestretch}{1.1}
# ## margins
#+LaTeX_HEADER: \geometry{a4paper, left=10mm, right=10mm, top=10mm}
# ## personalize the prefix in the name of the sections
#+LaTeX_HEADER: \usepackage{titlesec}
# ## fix bug in titlesec version
# ##  https://tex.stackexchange.com/questions/299969/titlesec-loss-of-section-numbering-with-the-new-update-2016-03-15
#+LaTeX_HEADER: \usepackage{etoolbox}
#+LaTeX_HEADER: 
#+LaTeX_HEADER: \makeatletter
#+LaTeX_HEADER: \patchcmd{\ttlh@hang}{\parindent\z@}{\parindent\z@\leavevmode}{}{}
#+LaTeX_HEADER: \patchcmd{\ttlh@hang}{\noindent}{}{}{}
#+LaTeX_HEADER: \makeatother
** Color
# ## define new colors
#+LATEX_HEADER: \RequirePackage{colortbl} % arrayrulecolor to mix colors
#+LaTeX_HEADER: \definecolor{myorange}{rgb}{1,0.2,0}
#+LaTeX_HEADER: \definecolor{mypurple}{rgb}{0.7,0,8}
#+LaTeX_HEADER: \definecolor{mycyan}{rgb}{0,0.6,0.6}
#+LaTeX_HEADER: \newcommand{\lightblue}{blue!50!white}
#+LaTeX_HEADER: \newcommand{\darkblue}{blue!80!black}
#+LaTeX_HEADER: \newcommand{\darkgreen}{green!50!black}
#+LaTeX_HEADER: \newcommand{\darkred}{red!50!black}
#+LaTeX_HEADER: \definecolor{gray}{gray}{0.5}
# ## change the color of the links
#+LaTeX_HEADER: \hypersetup{
#+LaTeX_HEADER:  citecolor=[rgb]{0,0.5,0},
#+LaTeX_HEADER:  urlcolor=[rgb]{0,0,0.5},
#+LaTeX_HEADER:  linkcolor=[rgb]{0,0,0.5},
#+LaTeX_HEADER: }
** Font
# https://tex.stackexchange.com/questions/25249/how-do-i-use-a-particular-font-for-a-small-section-of-text-in-my-document
#+LaTeX_HEADER: \newenvironment{note}{\small \color{gray}\fontfamily{lmtt}\selectfont}{\par}
#+LaTeX_HEADER: \newenvironment{activity}{\color{orange}\fontfamily{qzc}\selectfont}{\par}
** Symbols
# ## valid and cross symbols
#+LaTeX_HEADER: \RequirePackage{pifont}
#+LaTeX_HEADER: \RequirePackage{relsize}
#+LaTeX_HEADER: \newcommand{\Cross}{{\raisebox{-0.5ex}%
#+LaTeX_HEADER:		{\relsize{1.5}\ding{56}}}\hspace{1pt} }
#+LaTeX_HEADER: \newcommand{\Valid}{{\raisebox{-0.5ex}%
#+LaTeX_HEADER:		{\relsize{1.5}\ding{52}}}\hspace{1pt} }
#+LaTeX_HEADER: \newcommand{\CrossR}{ \textcolor{red}{\Cross} }
#+LaTeX_HEADER: \newcommand{\ValidV}{ \textcolor{green}{\Valid} }
# ## warning symbol
#+LaTeX_HEADER: \usepackage{stackengine}
#+LaTeX_HEADER: \usepackage{scalerel}
#+LaTeX_HEADER: \newcommand\Warning[1][3ex]{%
#+LaTeX_HEADER:   \renewcommand\stacktype{L}%
#+LaTeX_HEADER:   \scaleto{\stackon[1.3pt]{\color{red}$\triangle$}{\tiny\bfseries !}}{#1}%
#+LaTeX_HEADER:   \xspace
#+LaTeX_HEADER: }
# # R Software
#+LATEX_HEADER: \newcommand\Rlogo{\textbf{\textsf{R}}\xspace} % 
** Code
:PROPERTIES:
:ID: 2ec77c4b-f83d-4612-9a89-a96ba1b7bf70
:END:
# Documentation at https://org-babel.readthedocs.io/en/latest/header-args/#results
# :tangle (yes/no/filename) extract source code with org-babel-tangle-file, see http://orgmode.org/manual/Extracting-source-code.html 
# :cache (yes/no)
# :eval (yes/no/never)
# :results (value/output/silent/graphics/raw/latex)
# :export (code/results/none/both)
#+PROPERTY: header-args :session *R* :tangle yes :cache no ## extra argument need to be on the same line as :session *R*
# Code display:
#+LATEX_HEADER: \RequirePackage{fancyvrb}
#+LATEX_HEADER: \DefineVerbatimEnvironment{verbatim}{Verbatim}{fontsize=\small,formatcom = {\color[rgb]{0.5,0,0}}}
# ## change font size input (global change)
# ## doc: https://ctan.math.illinois.edu/macros/latex/contrib/listings/listings.pdf
# #+LATEX_HEADER: \newskip kipamount    kipamount =6pt plus 0pt minus 6pt
# #+LATEX_HEADER: \lstdefinestyle{code-tiny}{basicstyle=\ttfamily\tiny, aboveskip =  kipamount, belowskip =  kipamount}
# #+LATEX_HEADER: \lstset{style=code-tiny}
# ## change font size input (local change, put just before BEGIN_SRC)
# ## #+ATTR_LATEX: :options basicstyle=\ttfamily\scriptsize
# ## change font size output (global change)
# ## \RecustomVerbatimEnvironment{verbatim}{Verbatim}{fontsize=\tiny,formatcom = {\color[rgb]{0.5,0,0}}}
** Lists
#+LATEX_HEADER: \RequirePackage{enumitem} % better than enumerate
** Image and graphs
#+LATEX_HEADER: \RequirePackage{epstopdf} % to be able to convert .eps to .pdf image files
#+LATEX_HEADER: \RequirePackage{capt-of} % 
#+LATEX_HEADER: \RequirePackage{caption} % newlines in graphics
#+LaTeX_HEADER: \RequirePackage{tikz-cd} % graph
# ## https://tools.ietf.org/doc/texlive-doc/latex/tikz-cd/tikz-cd-doc.pdf
** Table
#+LATEX_HEADER: \RequirePackage{booktabs} % for nice lines in table (e.g. toprule, bottomrule, midrule, cmidrule)
** Inline latex
# @@latex:any arbitrary LaTeX code@@
** Algorithm
#+LATEX_HEADER: \RequirePackage{amsmath}
#+LATEX_HEADER: \RequirePackage{algorithm}
#+LATEX_HEADER: \RequirePackage[noend]{algpseudocode}
** Math
#+LATEX_HEADER: \RequirePackage{dsfont}
#+LATEX_HEADER: \RequirePackage{amsmath,stmaryrd,graphicx}
#+LATEX_HEADER: \RequirePackage{prodint} % product integral symbol (\PRODI)
# ## lemma
# #+LaTeX_HEADER: \RequirePackage{amsthm}
# #+LaTeX_HEADER: \newtheorem{theorem}{Theorem}
# #+LaTeX_HEADER: \newtheorem{lemma}[theorem]{Lemma}
*** Template for shortcut
#+LATEX_HEADER: \usepackage{ifthen}
#+LATEX_HEADER: \usepackage{xifthen}
#+LATEX_HEADER: \usepackage{xargs}
#+LATEX_HEADER: \usepackage{xspace}
#+LATEX_HEADER: \newcommand\defOperator[7]{%
#+LATEX_HEADER:	\ifthenelse{\isempty{#2}}{
#+LATEX_HEADER:		\ifthenelse{\isempty{#1}}{#7{#3}#4}{#7{#3}#4 \left#5 #1 \right#6}
#+LATEX_HEADER:	}{
#+LATEX_HEADER:	\ifthenelse{\isempty{#1}}{#7{#3}#4_{#2}}{#7{#3}#4_{#1}\left#5 #2 \right#6}
#+LATEX_HEADER: }
#+LATEX_HEADER: }
#+LATEX_HEADER: \newcommand\defUOperator[5]{%
#+LATEX_HEADER: \ifthenelse{\isempty{#1}}{
#+LATEX_HEADER:		#5\left#3 #2 \right#4
#+LATEX_HEADER: }{
#+LATEX_HEADER:	\ifthenelse{\isempty{#2}}{\underset{#1}{\operatornamewithlimits{#5}}}{
#+LATEX_HEADER:		\underset{#1}{\operatornamewithlimits{#5}}\left#3 #2 \right#4}
#+LATEX_HEADER: }
#+LATEX_HEADER: }
#+LATEX_HEADER: \newcommand{\defBoldVar}[2]{	
#+LATEX_HEADER:	\ifthenelse{\equal{#2}{T}}{\boldsymbol{#1}}{\mathbf{#1}}
#+LATEX_HEADER: }
**** Probability
#+LATEX_HEADER: \newcommandx\Esp[2][1=,2=]{\defOperator{#1}{#2}{E}{}{\lbrack}{\rbrack}{\mathbb}}
#+LATEX_HEADER: \newcommandx\Prob[2][1=,2=]{\defOperator{#1}{#2}{P}{}{\lbrack}{\rbrack}{\mathbb}}
#+LATEX_HEADER: \newcommandx\Qrob[2][1=,2=]{\defOperator{#1}{#2}{Q}{}{\lbrack}{\rbrack}{\mathbb}}
#+LATEX_HEADER: \newcommandx\Var[2][1=,2=]{\defOperator{#1}{#2}{V}{ar}{\lbrack}{\rbrack}{\mathbb}}
#+LATEX_HEADER: \newcommandx\Cov[2][1=,2=]{\defOperator{#1}{#2}{C}{ov}{\lbrack}{\rbrack}{\mathbb}}
#+LATEX_HEADER: \newcommandx\Binom[2][1=,2=]{\defOperator{#1}{#2}{B}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\Gaus[2][1=,2=]{\defOperator{#1}{#2}{N}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\Wishart[2][1=,2=]{\defOperator{#1}{#2}{W}{ishart}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\Likelihood[2][1=,2=]{\defOperator{#1}{#2}{L}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\logLikelihood[2][1=,2=]{\defOperator{#1}{#2}{\ell}{}{(}{)}{}}
#+LATEX_HEADER: \newcommandx\Information[2][1=,2=]{\defOperator{#1}{#2}{I}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\Hessian[2][1=,2=]{\defOperator{#1}{#2}{H}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\Score[2][1=,2=]{\defOperator{#1}{#2}{S}{}{(}{)}{\mathcal}}
**** Operators
#+LATEX_HEADER: \newcommandx\Vois[2][1=,2=]{\defOperator{#1}{#2}{V}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\IF[2][1=,2=]{\defOperator{#1}{#2}{IF}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\Ind[1][1=]{\defOperator{}{#1}{1}{}{(}{)}{\mathds}}
#+LATEX_HEADER: \newcommandx\Max[2][1=,2=]{\defUOperator{#1}{#2}{(}{)}{min}}
#+LATEX_HEADER: \newcommandx\Min[2][1=,2=]{\defUOperator{#1}{#2}{(}{)}{max}}
#+LATEX_HEADER: \newcommandx\argMax[2][1=,2=]{\defUOperator{#1}{#2}{(}{)}{argmax}}
#+LATEX_HEADER: \newcommandx\argMin[2][1=,2=]{\defUOperator{#1}{#2}{(}{)}{argmin}}
#+LATEX_HEADER: \newcommandx\cvD[2][1=D,2=n \rightarrow \infty]{\xrightarrow[#2]{#1}}
#+LATEX_HEADER: \newcommandx\Hypothesis[2][1=,2=]{
#+LATEX_HEADER:         \ifthenelse{\isempty{#1}}{
#+LATEX_HEADER:         \mathcal{H}
#+LATEX_HEADER:         }{
#+LATEX_HEADER: 	\ifthenelse{\isempty{#2}}{
#+LATEX_HEADER: 		\mathcal{H}_{#1}
#+LATEX_HEADER: 	}{
#+LATEX_HEADER: 	\mathcal{H}^{(#2)}_{#1}
#+LATEX_HEADER:         }
#+LATEX_HEADER:         }
#+LATEX_HEADER: }
#+LATEX_HEADER: \newcommandx\dpartial[4][1=,2=,3=,4=\partial]{
#+LATEX_HEADER: 	\ifthenelse{\isempty{#3}}{
#+LATEX_HEADER: 		\frac{#4 #1}{#4 #2}
#+LATEX_HEADER: 	}{
#+LATEX_HEADER: 	\left.\frac{#4 #1}{#4 #2}\right\rvert_{#3}
#+LATEX_HEADER: }
#+LATEX_HEADER: }
#+LATEX_HEADER: \newcommandx\dTpartial[3][1=,2=,3=]{\dpartial[#1][#2][#3][d]}
#+LATEX_HEADER: \newcommandx\ddpartial[3][1=,2=,3=]{
#+LATEX_HEADER: 	\ifthenelse{\isempty{#3}}{
#+LATEX_HEADER: 		\frac{\partial^{2} #1}{\partial #2^2}
#+LATEX_HEADER: 	}{
#+LATEX_HEADER: 	\frac{\partial^2 #1}{\partial #2\partial #3}
#+LATEX_HEADER: }
#+LATEX_HEADER: } 
**** General math
#+LATEX_HEADER: \newcommand\Real{\mathbb{R}}
#+LATEX_HEADER: \newcommand\Rational{\mathbb{Q}}
#+LATEX_HEADER: \newcommand\Natural{\mathbb{N}}
#+LATEX_HEADER: \newcommand\trans[1]{{#1}^\intercal}%\newcommand\trans[1]{{\vphantom{#1}}^\top{#1}}
#+LATEX_HEADER: \newcommand{\independent}{\mathrel{\text{\scalebox{1.5}{$\perp\mkern-10mu\perp$}}}}
#+LaTeX_HEADER: \newcommand\half{\frac{1}{2}}
#+LaTeX_HEADER: \newcommand\normMax[1]{\left|\left|#1\right|\right|_{max}}
#+LaTeX_HEADER: \newcommand\normTwo[1]{\left|\left|#1\right|\right|_{2}}
#+LATEX_HEADER: \newcommand\Veta{\boldsymbol{\eta}}

** Notations

#+LaTeX_HEADER:\newcommand{\Model}{\mathcal{M}}
#+LaTeX_HEADER:\newcommand{\ModelHat}{\widehat{\mathcal{M}}}

#+LaTeX_HEADER:\newcommand{\param}{\Theta}
#+LaTeX_HEADER:\newcommand{\paramHat}{\widehat{\param}}
#+LaTeX_HEADER:\newcommand{\paramCon}{\widetilde{\param}}

#+LaTeX_HEADER:\newcommand{\Vparam}{\boldsymbol{\param}}
#+LaTeX_HEADER:\newcommand{\VparamT}{\Vparam_0}
#+LaTeX_HEADER:\newcommand{\VparamHat}{\boldsymbol{\paramHat}}
#+LaTeX_HEADER:\newcommand{\VparamCon}{\boldsymbol{\paramCon}}

#+LaTeX_HEADER:\newcommand{\X}{X}
#+LaTeX_HEADER:\newcommand{\x}{x}
#+LaTeX_HEADER:\newcommand{\VX}{\boldsymbol{X}}
#+LaTeX_HEADER:\newcommand{\Vx}{\boldsymbol{x}}

#+LaTeX_HEADER:\newcommand{\Y}{Y}
#+LaTeX_HEADER:\newcommand{\y}{y}
#+LaTeX_HEADER:\newcommand{\VY}{\boldsymbol{Y}}
#+LaTeX_HEADER:\newcommand{\Vy}{\boldsymbol{y}}
#+LaTeX_HEADER:\newcommand{\Vvarepsilon}{\boldsymbol{\varepsilon}}


