Type:	Package
Title:	Repeated Measurement Models for Discrete Times
Version:	1.1.0
Date:	2024-05-12
Description:	Companion R package for the course "Statistical analysis of correlated and repeated measurements for health science researchers" taught by the section of Biostatistics of the University of Copenhagen. It implements linear mixed models where the model for the variance-covariance of the residuals is specified via patterns (compound symmetry, toeplitz, unstructured, ...). Statistical inference for mean, variance, and correlation parameters is performed based on the observed information and a Satterthwaite approximation of the degrees of freedom. Normalized residuals are provided to assess model misspecification. Statistical inference can be performed for arbitrary linear or non-linear combination(s) of model coefficients. Predictions can be computed conditional to covariates only or also to outcome values.
License:	GPL-3
Encoding:	UTF-8
URL:	https://github.com/bozenne/LMMstar
BugReports:	https://github.com/bozenne/LMMstar/issues
Depends:	R (≥ 3.5.0)
Imports:	copula, doParallel, foreach, ggplot2, grid, lava, Matrix, multcomp, nlme, numDeriv, parallel, rlang,
Suggests:	asht, data.table, ggh4x, ggpubr, lattice, mvtnorm, lme4, lmerTest, mice, nlmeU, optimx, pbapply, psych, Publish, qqtest, R.rsp, reshape2, rmcorr, scales, testthat
VignetteBuilder:	R.rsp
RoxygenNote:	7.3.1
Collate:	'0-onload.R' 'BuyseTest_discreteRoot.R' 'LMMstar-package.R' 'LMMstar.options.R' 'add.R' 'anova.R' 'autoplot.R' 'backtransform.R' 'baselineAdjustment.R' 'coef.R' 'confint.R' 'constrain.R' 'df.R' 'discreteRoot.R' 'doc-data.R' 'effects.R' 'estimate.R' 'findPatterns.R' 'fitted.R' 'formula.R' 'iid.R' 'information.R' 'levels.R' 'lmm.R' 'lmmCC.R' 'logLik.R' 'manifest.R' 'mlmm.R' 'model.frame.R' 'model.matrix.R' 'model.tables.R' 'moments.R' 'mt.test.R' 'multcomp.R' 'namePatterns.R' 'nobs.R' 'operator.R' 'partialCor.R' 'plot.R' 'precompute.R' 'predict.R' 'print.R' 'profile.R' 'proportion.R' 'ranef.R' 'rbind.R' 'reformat.R' 'remove.R' 'reparametrize.R' 'resample.R' 'residuals.R' 'restaureNA.R' 'sampleRem.R' 'scatterplot.R' 'score.R' 'sigma.R' 'structure-calc_Omega.R' 'structure-calc_d2Omega.R' 'structure-calc_dOmega.R' 'structure-initialization.R' 'structure-skeleton.R' 'structure-skeletonK.R' 'structure-skeletonRho.R' 'structure-skeletonSigma.R' 'structure-update.R' 'structure.R' 'summarize.R' 'summarizeNA.R' 'summary.R' 'terms.R' 'utils.R' 'vcov.R' 'weights.R'
NeedsCompilation:	no
Packaged:	2024-05-12 21:11:58 UTC; bozenne
Author:	Brice Ozenne [aut, cre], Julie Forman [aut]
Maintainer:	Brice Ozenne <brice.mh.ozenne@gmail.com>
Repository:	CRAN
Date/Publication:	2024-05-12 21:43:11 UTC

LMMstar package: repeated measurement models for discrete times

Description

Companion R package for the course "Statistical analysis of correlated and repeated measurements for health science researchers" taught by the section of Biostatistics of the University of Copenhagen. It implements linear mixed models where the model for the variance-covariance of the residuals is specified via patterns (compound symmetry, toeplitz, unstructured, ...). Statistical inference for mean, variance, and correlation parameters is performed based on the observed information and a Satterthwaite approximation of the degrees of freedom. Normalized residuals are provided to assess model misspecification. Statistical inference can be performed for arbitrary linear or non-linear combination(s) of model coefficients. Predictions can be computed conditional to covariates only or also to outcome values.

Details

Notations: the linear mixed model estimated by lmm is denoted:

\mathbf{Y}_{i} = \mathbf{X}_{i}\beta+\boldsymbol{\varepsilon}_i

where

• \mathbf{Y}=(Y_1,\ldots,Y_m): vector of outcomes.

• \mathbf{X}=(X_1,\ldots,X_p): design matrix (extractor: model.matrix.lmm).

• \boldsymbol{\varepsilon}: vector of residuals with 0-mean and variance \Omega_i (extractor: residuals.lmm).

• \beta: estimated mean coefficients relative to X (extractor: coef.lmm).

• \Omega: the modeled variance-covariance of the residuals with diagonal elements \sigma^2_{j} (extractor: sigma.lmm).

• i indexes the cluster (level where replicates are assumed independent).

• j indexes the repetitions, e.g. the variance of \varepsilon_{ij} is \sigma^2_{ij}.

Covariance patterns: \Omega can be parametrized as:

• ID: identity (no correlation, constant variance).

• IND: independent (no correlation, time-specific variance).

• CS: compound symmetry (constant correlation and variance). Can also be used to specify a nested random effect structure or a block specific correlation and variance.

• RE: random effects.

• TOEPLITZ: toeplitz (lag-specific correlation, time-specific variance).

• UN: unstructured (time-specific correlation, time-specific variance).

It possible to stratify each structure with respect to a categorical variable.

Optimizer: the default optimizer, "FS", implements a fisher scoring algorithm descent with back-tracking in case of decreasing or undefined log-likelihood. It does not constrain \Omega to be positive definite which may cause problem in small sample or complex models. It is possible to use other optimizer inferfaced by optimx::optimx.

Keywords: documented methods/functions are classified according to the following keywords

• models: function fitting a statistical model based on a dataset (e.g. lmm, lmmCC, mlmm, mt.test, partialCor)

• htest: methods performing statistical inference based on an existing model (e.g. anova.lmm, estimate.lmm, effects.lmm, profile.lmm, proportion, rbind.Wald_lmm, resample.lmm)

• methods: extractors (e.g. coef.lmm, confint.lmm, df.residual.lmm, fitted.lmm, iid.lmm, information.lmm, levels.lmm, logLik.lmm, manifest.lmm, model.frame.lmm, model.matrix.lmm, model.tables.lmm, nobs.lmm, predict.lmm, ranef.lmm, residuals.lmm, score.lmm, sigma.lmm, summary.lmm, vcov.lmm, weights.Wald_lmm)

• utilities: function used to facilitate the user interface (e.g. add, baselineAdjustment, LMMstar.options, remove, scatterplot, summarize, summarizeNA)

• datasets: dataset stored in the package (e.g. abetaW)

• hplot: graphical display (e.g. autoplot.lmm or plot.lmm)

• datagen: function for generating data sets (e.g. sampleRem)

• multivariate: covariance patterns (e.g. ID, IND, CS, RE, TOEPLITZ, UN, CUSTOM)

Author(s)

Maintainer: Brice Ozenne brice.mh.ozenne@gmail.com (ORCID)

Authors:

• Julie Forman jufo@sund.ku.dk (ORCID)

See Also

Useful links:

• https://github.com/bozenne/LMMstar

• Report bugs at https://github.com/bozenne/LMMstar/issues

Compound Symmetry Structure

Description

Variance-covariance structure where the variance and correlation of the residuals is constant within covariate levels. Can be stratified on a categorical variable. The default has no covariate and therefore the variance and correlation are constant within cluster.

Usage

CS(formula, var.cluster, var.time, type = NULL, group.type = NULL, add.time)

Arguments

Details

A typical formula would be ~1, indicating a variance constant over time and the same correlation between all pairs of times.

Value

An object of class CS that can be passed to the argument structure of the lmm function.

Examples

## no covariates
CS(~1, var.cluster = "id", var.time = "time")
CS(gender~1, var.cluster = "id", var.time = "time")

## covariates
CS(~time, var.cluster = "id", var.time = "time")
CS(gender~time, var.cluster = "id", var.time = "time")
CS(list(~time,~1), var.cluster = "id", var.time = "time")
CS(list(gender~time,gender~1), var.cluster = "id", var.time = "time")

Custom Structure

Description

Variance-covariance structure specified by the user.

Usage

CUSTOM(
  formula,
  var.cluster,
  var.time,
  FCT.sigma,
  dFCT.sigma = NULL,
  d2FCT.sigma = NULL,
  init.sigma,
  FCT.rho,
  dFCT.rho = NULL,
  d2FCT.rho = NULL,
  init.rho,
  add.time
)

Arguments

Value

An object of class CUSTOM that can be passed to the argument structure of the lmm function.

Examples

## Compound symmetry structure
CUSTOM(~1,
       FCT.sigma = function(p,n.time,X){rep(p,n.time)},
       init.sigma = c("sigma"=1),
       dFCT.sigma = function(p,n.time,X){list(sigma = rep(1,n.time))},  
       d2FCT.sigma = function(p,n.time,X){list(sigma = rep(0,n.time))},  
       FCT.rho = function(p,n.time,X){
           matrix(p,n.time,n.time)+diag(1-p,n.time,n.time)
       },
       init.rho = c("rho"=0.5),
       dFCT.rho = function(p,n.time,X){
            list(rho = matrix(1,n.time,n.time)-diag(1,n.time,n.time))
       },
       d2FCT.rho = function(p,n.time,X){list(rho = matrix(0,n.time,n.time))}
)

## 2 block structure
rho.2block <- function(p,n.time,X){
   rho <- matrix(0, nrow = n.time, ncol = n.time)
   rho[1,2] <- rho[2,1] <- rho[4,5] <- rho[5,4] <- p["rho1"]
   rho[1,3] <- rho[3,1] <- rho[4,6] <- rho[6,4] <- p["rho2"]
   rho[2,3] <- rho[3,2] <- rho[5,6] <- rho[6,5] <- p["rho3"]
   rho[4:6,1:3] <- rho[1:3,4:6] <- p["rho4"]
   return(rho)
}
drho.2block <- function(p,n.time,X){
   drho <- list(rho1 = matrix(0, nrow = n.time, ncol = n.time),
                rho2 = matrix(0, nrow = n.time, ncol = n.time),
                rho3 = matrix(0, nrow = n.time, ncol = n.time),
                rho4 = matrix(0, nrow = n.time, ncol = n.time))   
   drho$rho1[1,2] <- drho$rho1[2,1] <- drho$rho1[4,5] <- drho$rho1[5,4] <- 1
   drho$rho2[1,3] <- drho$rho2[3,1] <- drho$rho2[4,6] <- drho$rho2[6,4] <- 1
   drho$rho3[2,3] <- drho$rho3[3,2] <- drho$rho3[5,6] <- drho$rho3[6,5] <- 1
   drho$rho4[4:6,1:3] <- drho$rho4[1:3,4:6] <- 1
   return(drho)
}
d2rho.2block <- function(p,n.time,X){
   d2rho <- list(rho1 = matrix(0, nrow = n.time, ncol = n.time),
                 rho2 = matrix(0, nrow = n.time, ncol = n.time),
                 rho3 = matrix(0, nrow = n.time, ncol = n.time),
                 rho4 = matrix(0, nrow = n.time, ncol = n.time))   
   return(d2rho)
}

CUSTOM(~variable,
       FCT.sigma = function(p,n.time,X){rep(p,n.time)},
       dFCT.sigma = function(p,n.time,X){list(sigma=rep(1,n.time))},
       d2FCT.sigma = function(p,n.time,X){list(sigma=rep(0,n.time))},
       init.sigma = c("sigma"=1),
       FCT.rho = rho.2block,
       dFCT.rho = drho.2block,
       d2FCT.rho = d2rho.2block,
       init.rho = c("rho1"=0.25,"rho2"=0.25,"rho3"=0.25,"rho4"=0.25))

identity Structure

Description

Variance-covariance structure where the residuals are independent and identically distributed. Can be stratified on a categorical variable.

Usage

ID(formula, var.cluster, var.time, add.time)

Arguments

Details

A typical formula would be ~1.

Value

An object of class IND that can be passed to the argument structure of the lmm function.

Examples

ID(NULL, var.cluster = "id", var.time = "time")
ID(~1, var.cluster = "id", var.time = "time")
ID(~gender, var.cluster = "id", var.time = "time")
ID(gender~1, var.cluster = "id", var.time = "time")

Independence Structure

Description

Variance-covariance structure where the residuals are independent but may have different variance. Can be stratified on a categorical variable.

Usage

IND(formula, var.cluster, var.time, add.time)

Arguments

Details

A typical formula would be either ~1 indicating constant variance or ~time indicating a time dependent variance.

Value

An object of class IND that can be passed to the argument structure of the lmm function.

Examples

IND(NULL, var.cluster = "id", var.time = "time", add.time = TRUE)
IND(~1, var.cluster = "id", var.time = "time")
IND(gender~1, var.cluster = "id", var.time = "time")
IND(gender~time, var.cluster = "id", var.time = "time")
IND(~gender+time, var.cluster = "id", var.time = "time")

Global options for LMMstar package

Description

Update or select global options for the LMMstar package.

Usage

LMMstar.options(..., reinitialise = FALSE)

Arguments

Details

The options are:

• backtransform.confint [logical]: should variance/covariance/correlation estimates be back-transformed when they are transformed on the log or atanh scale. Used by confint.

• columns.anova [character vector]: columns to ouput when using anova with argument ci=TRUE.

• columns.confint [character vector]: columns to ouput when using confint.

• columns.summary [character vector]: columns to ouput when displaying the model coefficients using summary.

• df [logical]: should approximate degrees of freedom be computed for Wald and F-tests. Used by lmm, anova, predict, and confint.

• drop.X [logical]: should columns causing non-identifiability of the model coefficients be dropped from the design matrix. Used by lmm.

• effects [character]: parameters relative to which estimates, score, information should be output.

• min.df [integer]: minimum possible degree of freedom. Used by confint.

• method.fit [character]: objective function when fitting the Linear Mixed Model (REML or ML). Used by lmm.

• method.numDeriv [character]: type used to approximate the third derivative of the log-likelihood (when computing the degrees of freedom). Can be "simple" or "Richardson". See numDeriv::jacobian for more details. Used by lmm.

• n.sampleCopula [integer]: number of samples used to compute confidence intervals and p-values adjusted for multiple comparisons via "single-step2". Used by confint.Wald_lmm.

• optimizer [character]: method used to estimate the model parameters. Either "FS", an home-made fisher scoring algorithm, or a method from optimx:optimx like "BFGS"or Nelder-Mead.

• param.optimizer [numeric vector]: default option for the FS optimization routine: maximum number of gradient descent iterations (n.iter), maximum acceptable score value (tol.score), maximum acceptable change in parameter value (tol.param).

• precompute.moments [logical]: Should the cross terms between the residuals and design matrix be pre-computed. Useful when the number of subject is substantially larger than the number of mean paramters.

• sep [character vector]: character used to combined two strings of characters in various functions (lp: .vcov.model.matrix, k.cov/k.strata: .skeletonK, pattern: .findUpatterns, rho.name/rho.strata: .skeletonRho, reformat: .reformat ).

• trace [logical]: Should the progress of the execution of the lmm function be displayed?

• tranform.sigma, tranform.k, tranform.rho: transformation used to compute the confidence intervals/p-values for the variance and correlation parameters. See the detail section of the coef function for more information. Used by lmm, anova and confint.

• type.information [character]: Should the expected or observed information ("expected" or "observed") be used to perform statistical inference? Used by lmm, anova and confint.

Value

A list containing the default options.

Random Effect Structure

Description

Variance-covariance structure parametrized via random effects. Can be stratified on a categorical variable.

Usage

RE(formula, var.cluster, var.time, ranef = NULL, add.time)

Arguments

Details

A typical formula would be ~1, indicating a variance constant over time and the same correlation between all pairs of times.

Value

An object of class CS that can be passed to the argument structure of the lmm function.

Examples

RE(~1, var.cluster = "id", var.time = "time")
RE(~gender, var.cluster = "id", var.time = "time")
RE(gender~(1|id), var.time = "time")

Toeplitz Structure

Description

Variance-covariance structure where the correlation depends on time elapsed between two repetitions. Can be stratified on a categorical variable.

Usage

TOEPLITZ(formula, var.cluster, var.time, type = "LAG", add.time)

Arguments

Details

formula: there can only be at most one covariate for the correlation structure. A typical formula would be ~1, indicating a variance constant over time and a correlation specific to each gap time.

type: for a binary covariate the correlation matrix can be decomposed into four blocs: A, B, B, C. A correspond the correlation within level 0 of the covariate, C within level 1, and B between level 0 and 1. Different correlation structures can be specified:

• "UN": unstructured matrix except for the diagonal elements of C which are constrained to be equal.

• "LAG": Toeplitz structure within A, B, and C, i.e. correlation specific to each time lag and covariate level.

• "CS": block-specific value except for C which has a different value for its diagonal elements.

Value

An object of class TOEPLITZ that can be passed to the argument structure of the lmm function.

Examples

## no covariate
TOEPLITZ(~time, var.cluster = "id", var.time = "time")
TOEPLITZ(gender~time, var.cluster = "id", var.time = "time")
TOEPLITZ(list(~time,~time), var.cluster = "id", var.time = "time")
TOEPLITZ(list(gender~time,gender~time), var.cluster = "id", var.time = "time")

## with covariates
TOEPLITZ(~side, var.cluster = "id", type = "UN",
         var.time = "time", add.time = TRUE)
TOEPLITZ(~side, var.cluster = "id", type = "LAG",
         var.time = "time", add.time = TRUE)
TOEPLITZ(~side, var.cluster = "id", type = "CS",
         var.time = "time", add.time = TRUE)
TOEPLITZ(gender~side, var.cluster = "id", type = "CS",
         var.time = "time", add.time = TRUE)

Unstructured Structure

Description

Variance-covariance structure where the residuals have time-specific variance and correlation. Can be stratified on a categorical variable.

Usage

UN(formula, var.cluster, var.time, add.time)

Arguments

Details

A typical formula would be ~1, indicating a time-specific variance parameter and a correlation parameter specific to each pair of times.

Value

An object of class UN that can be passed to the argument structure of the lmm function.

Examples

UN(NULL, var.cluster = "id", var.time = "time", add.time = TRUE)
UN(~gender, var.cluster = "id", var.time = "time", add.time = TRUE)
UN(gender ~ 1, var.cluster = "id", var.time = "time", add.time = TRUE)
UN(list(~gender,~1), var.cluster = "id", var.time = "time", add.time = TRUE)
UN(list(gender~age,gender~1), var.cluster = "id", var.time = "time", add.time = TRUE)

Data From The Bland Altman Study (Long Format)

Description

Extract data from a longitudinal case control study including 87 patients newly diagnosed with bipolar disorder and 44 age and sex matched healthy controls. Contains demographic data and lifestyle factors at baseline, as well as measures of psychosocial functioning at baseline and 1 year follow-up. This dataset is in the long format (i.e. one line per measurement).

• id: study participant.

• sex: male (M) or female (F).

• age: age in years.

• group: bipolar disorder (BD) or healthy control (HC).

• episode: whether the patient experience an affective episode during follow-up.

• visit: index of time at which pss, fast, and qol measurements where performed.

• year: time at which pss, fast, and qol measurements where performed.

• pss: perceived stress score.

• fast: functioning assessment short test.

• qol: WHO quality of life score.

• educationyears: years of education including basic school.

• alcohol: daily alcohol consumption.

• missingreason: reason of drop out or missed visit.

Usage

data(abetaL)

References

Pech, Josefine, et al. The impact of a new affective episode on psychosocial functioning, quality of life and perceived stress in newly diagnosed patients with bipolar disorder: A prospective one-year case-control study.Journal of Affective Disorders 277 (2020): 486-494.

Data From The abeta Study (Wide Format)

Description

Extract data from a longitudinal case control study including 87 patients newly diagnosed with bipolar disorder and 44 age and sex matched healthy controls. Contains demographic data and lifestyle factors at baseline, as well as measures of psychosocial functioning at baseline and 1 year follow-up. This dataset is in the wide format (i.e. one line per participant).

• id: study participant.

• sex: male (M) or female (F).

• age: age in years.

• group: bipolar disorder (BD) or healthy control (HC).

• episode: whether the patient experience an affective episode during follow-up.

• fast0,fast1: functioning assessment short test at baseline and follow-up.

• qol0,qol1: WHO quality of life score at baseline and follow-up.

• pss0,pss1: perceived stress score at baseline and follow-up.

• educationyears: years of education including basic school.

• alcohol: daily alcohol consumption.

• missingreason: reason of drop out or missed visit.

Usage

data(abetaW)

References

Pech, Josefine, et al. "The impact of a new affective episode on psychosocial functioning, quality of life and perceived stress in newly diagnosed patients with bipolar disorder: A prospective one-year case-control study."Journal of Affective Disorders 277 (2020): 486-494.

Add Columns to Output

Description

Auxiliary function that can be used when specifying the argument columns (e.g. calling confint.lmm) to add columns.

Usage

add(...)

Arguments

Value

A character vector

Examples

set.seed(10)
dW <- sampleRem(25, n.times = 1, format = "long")
e.lmm <- lmm(Y~X1, data = dW)

confint(e.lmm, columns = add("statistic"))

Multivariate Tests For Linear Mixed Model

Description

Simultaneous tests of linear combinations of the model paramaters using Wald tests or Likelihood Ratio Test (LRT).

Usage

## S3 method for class 'lmm'
anova(
  object,
  effects = NULL,
  robust = FALSE,
  multivariate = TRUE,
  rhs = NULL,
  df = !is.null(object$df),
  ci = TRUE,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  transform.names = TRUE,
  ...
)

Arguments

Details

By default adjustment of confidence intervals and p-values for multiple comparisons is based on the distribution of the maximum-statistic. This is refered to as a single-step Dunnett multiple testing procedures in table II of Dmitrienko et al. (2013). It is performed using the multcomp package with the option test = adjusted("single-step") with equal degrees of freedom or by simulation using a Student's t copula with unequal degrees of freedom (more in the note of the details section of confint.Wald_lmm).

Value

A data.frame (LRT) or a list of containing the following elements (Wald):

• multivariate: data.frame containing the multivariate Wald test. The column df.num refers to the degrees of freedom for the numerator (i.e. number of hypotheses) wherease the column df.denum refers to the degrees of freedom for the denominator (i.e. Satterthwaite approximation).

• univariate: data.frame containing each univariate Wald test.

• glht: used internally to call functions from the multcomp package.

• object: list containing key information about the linear mixed model.

• vcov: variance-covariance matrix associated to each parameter of interest (i.e. hypothesis).

• iid: matrix containing the influence function relative to each parameter of interest (i.e. hypothesis).

• args: list containing argument values from the function call.

References

Dmitrienko, A. and D'Agostino, R., Sr (2013), Traditional multiplicity adjustment methods in clinical trials. Statist. Med., 32: 5172-5218. https://doi.org/10.1002/sim.5990.

See Also

summary.Wald_lmm or confint.Wald_lmm for a summary of the results.
autoplot.Wald_lmm for a graphical display of the results.
rbind.Wald_lmm for combining result across models and adjust for multiple comparisons.

Examples

#### simulate data in the long format ####
set.seed(10)
dL <- sampleRem(100, n.times = 3, format = "long")

#### fit Linear Mixed Model ####
eUN.lmm <- lmm(Y ~ visit + X1 + X2 + X5,
               repetition = ~visit|id, structure = "UN", data = dL)

#### Multivariate Wald test ####
## F-tests
anova(eUN.lmm)
anova(eUN.lmm, effects = "all")
anova(eUN.lmm, robust = TRUE, df = FALSE)
summary(anova(eUN.lmm), method = "bonferroni")

## user defined F-test
summary(anova(eUN.lmm, effects = c("X1=0","X2+X5=10")))

## chi2-tests
anova(eUN.lmm, df = FALSE)

## with standard contrast
if(require(multcomp)){
amod <- lmm(breaks ~ tension, data = warpbreaks)
e.amod <- anova(amod, effect = mcp(tension = "Tukey"))
summary(e.amod)
}

#### Likelihood ratio test ####
eUN0.lmm <- lmm(Y ~ X1 + X2, repetition = ~visit|id, structure = "UN", data = dL)
anova(eUN.lmm, eUN0.lmm) 

eCS.lmm <- lmm(Y ~ X1 + X2 + X5, repetition = ~visit|id, structure = "CS", data = dL)
anova(eUN.lmm, eCS.lmm)

Graphical Display For Linear Hypothesis Test

Description

Graphical Display For Linear Hypothesis Test

Usage

## S3 method for class 'Wald_lmm'
autoplot(object, type = "forest", size.text = 16, add.args = NULL, ...)

## S3 method for class 'Wald_lmm'
plot(x, ...)

Arguments

Details

Argument add.args: parameters specific to the forest plot:

• color: [logical] should the estimates be colored by global null hypothesis, e.g. when testing the effect of a 3 factor covariate, the two corresponding coefficient will have the same color. Alternatively a vector of positive integers giving the color with which each estimator should be displayed.

• color: [logical] should the estimates be represented by a different shape per global null hypothesis, e.g. when testing the effect of a 3 factor covariate, the two corresponding coefficient will have the same type of point. Alternatively a vector of positive integers describing the shape to be used for each estimator.

• ci: [logical] should confidence intervals be displayed?

• size.estimate: [numeric, >0] size of the dot used to display the estimates.

• size.ci: [numeric, >0] thickness of the line used to display the confidence intervals.

• width.ci: [numeric, >0] width of the line used to display the confidence intervals.

• size.null: [numeric, >0] thickness of the line used to display the null hypothesis.

Parameters specific to the heatmap plot:

• limits: [numeric vector of length 2] minimum and maximum value of the colorscale relative to the correlation.

• low, mid, high: [character] color for the the colorscale relative to the correlation.

• midpoint: [numeric] correlation value associated with the color defined by argument mid

Value

A list with two elements

• data: data used to create the graphical display.

• plot: ggplot object.

Functions

• plot(Wald_lmm): Graphical Display For Linear Hypothesis Test

Examples

## From the multcomp package
if(require(datasets) && require(ggplot2)){

## only tests with 1 df
ff <- Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality
e.lmm <- lmm(ff, data = swiss)
e.aovlmm <- anova(e.lmm)

autoplot(e.aovlmm, type = "forest")
autoplot(e.aovlmm, type = "heat") ## 3 color gradient
autoplot(e.aovlmm, type = "heat", add.args = list(mid = NULL)) ## 2 color gradient

## test with more than 1 df
e.lmm2 <- lmm(breaks ~ tension + wool, data = warpbreaks)
e.aovlmm2 <- anova(e.lmm2)
autoplot(e.aovlmm2)
autoplot(e.aovlmm2, add.args = list(color = FALSE, shape = FALSE))
}

Graphical Display For Linear Mixed Models

Description

Display fitted values or residual plot for the mean, variance, and correlation structure. Can also display quantile-quantile plot relative to the normal distribution.

Usage

## S3 method for class 'lmm'
autoplot(
  object,
  type = "fit",
  type.residual = NULL,
  obs.alpha = 0,
  obs.size = NULL,
  facet = NULL,
  facet_nrow = NULL,
  facet_ncol = NULL,
  scales = "fixed",
  labeller = "label_value",
  at = NULL,
  time.var = NULL,
  color = NULL,
  position = NULL,
  ci = TRUE,
  ci.alpha = NULL,
  ylim = NULL,
  mean.size = c(3, 1),
  size.text = 16,
  position.errorbar = "identity",
  ...
)

## S3 method for class 'lmm'
plot(x, ...)

Arguments

Value

A list with two elements

• data: data used to create the graphical display.

• plot: ggplot object.

Functions

• plot(lmm): Graphical Display For Linear Mixed Models

See Also

plot.lmm for other graphical display (residual plots, partial residual plots).

Examples

if(require(ggplot2)){

#### simulate data in the long format ####
set.seed(10)
dL <- sampleRem(100, n.times = 3, format = "long")
dL$X1 <- as.factor(dL$X1)

#### fit Linear Mixed Model ####
eCS.lmm <- lmm(Y ~ visit + X1 + X6,
               repetition = ~visit|id, structure = "CS", data = dL, df = FALSE)

#### model fit ####
plot(eCS.lmm, type = "fit", facet =~X1)
## customize display
gg <- autoplot(eCS.lmm, type = "fit", facet =~X1)$plot
gg + coord_cartesian(ylim = c(0,6))
## restrict to specific covariate value
plot(eCS.lmm, type = "fit", at = data.frame(X6=1), color = "X1")

#### qqplot ####
plot(eCS.lmm, type = "qqplot")
plot(eCS.lmm, type = "qqplot", engine.qqplot = "qqtest")

#### residual correlation ####
plot(eCS.lmm, type = "correlation")

#### residual trend ####
plot(eCS.lmm, type = "scatterplot")

#### residual heteroschedasticity ####
plot(eCS.lmm, type = "scatterplot2")

#### partial residuals ####
plot(eCS.lmm, type = "partial", type.residual = "visit") 
plot(eCS.lmm, type = "partial", type.residual = c("(Intercept)","X1","visit"))
plot(eCS.lmm, type = "partial", type.residual = c("(Intercept)","X1","visit"),
facet = ~X1)
}

Graphical Display For Partial Correlation

Description

Extract and display the correlation modeled via the linear mixed model.

Usage

## S3 method for class 'partialCor'
autoplot(
  object,
  size.text = 16,
  limits = c(-1, 1.00001),
  low = "blue",
  mid = "white",
  high = "red",
  midpoint = 0,
  ...
)

## S3 method for class 'partialCor'
plot(x, ...)

Arguments

Value

A list with two elements

• data: data used to create the graphical display.

• plot: ggplot object.

Functions

• plot(partialCor): Graphical Display For Partial Correlation

Examples

if(require(ggplot2)){
data(gastricbypassL, package = "LMMstar")

e.pCor <- partialCor(c(weight,glucagonAUC)~time, repetition = ~visit|id,
                     data = gastricbypassL)
plot(e.pCor)
}

Graphical Display of Profile Likelihood

Description

Graphical representation of the profile likelihood from a linear mixed model

Usage

## S3 method for class 'profile_lmm'
autoplot(
  object,
  type = "logLik",
  quadratic = TRUE,
  ci = FALSE,
  size = c(3, 2, 1, 1),
  linetype = c("dashed", "dashed", "dashed"),
  shape = 19,
  scales = "free",
  nrow = NULL,
  ncol = NULL,
  ...
)

## S3 method for class 'profile_lmm'
plot(x, ...)

Arguments

Value

A list with three elements

• data.fit: data containing the quadratice approximation of the log-likelihood

• data.ci: data containing the confidence intervals.

• plot: ggplot object.

Functions

• plot(profile_lmm): Display Contour of the log-Likelihood

Graphical Display of the Residuals

Description

Graphical representation of the residuals from a linear mixed model. Require a long format (except for the correlation where both format are accepted) and having exported the dataset along with the residual (argument keep.data when calling residuals.lmm).

Usage

## S3 method for class 'residuals_lmm'
autoplot(
  object,
  type = NULL,
  type.residual = NULL,
  time.var = NULL,
  facet = NULL,
  facet_nrow = NULL,
  facet_ncol = NULL,
  engine.qqplot = "ggplot2",
  add.smooth = TRUE,
  digits.cor = 2,
  size.text = 16,
  color = NULL,
  obs.size = NULL,
  mean.size = c(3, 1),
  ci.alpha = 0.25,
  position = NULL,
  scales = "fixed",
  labeller = "label_value",
  ...
)

## S3 method for class 'residuals_lmm'
plot(x, ...)

Arguments

Value

A list with two elements

• data: data used to generate the plot.

• plot: ggplot object.

Functions

• plot(residuals_lmm): Graphical Display of the Residuals

Graphical Display of the Descriptive Statistics

Description

Graphical representation of the descriptive statistics.

Usage

## S3 method for class 'summarize'
autoplot(
  object,
  type = "mean",
  variable = NULL,
  size.text = 16,
  linewidth = 1.25,
  size = 3,
  ...
)

## S3 method for class 'summarize'
plot(x, ...)

Arguments

Value

A list with two elements

• data: data used to generate the plot.

• plot: ggplot object.

Functions

• plot(summarize): Graphical Display of Missing Data Pattern

Examples

data(gastricbypassL, package = "LMMstar")
dtS <- summarize(weight ~ time, data = gastricbypassL)
plot(dtS)
dtS <- summarize(glucagonAUC + weight ~ time|id, data = gastricbypassL, na.rm = TRUE)
plot(dtS, variable = "glucagonAUC")
plot(dtS, variable = "glucagonAUC", type = "correlation", size.text = 1)

Graphical Display of Missing Data Pattern

Description

Graphical representation of the possible missing data patterns in the dataset.

Usage

## S3 method for class 'summarizeNA'
autoplot(
  object,
  variable = NULL,
  size.text = 16,
  add.missing = " missing",
  order.pattern = NULL,
  ...
)

## S3 method for class 'summarizeNA'
plot(x, ...)

Arguments

Value

A list with two elements

• data: data used to create the graphical display.

• plot: ggplot object.

Functions

• plot(summarizeNA): Graphical Display of Missing Data Pattern

Perform Baseline Adjustment

Description

Create a new variable based on a time variable and a group variable where groups are constrained to be equal at specific timepoints.

Usage

baselineAdjustment(
  object,
  variable,
  repetition,
  constrain,
  new.level = NULL,
  collapse.time = NULL
)

Arguments

Value

A vector of length the number of rows of the dataset.

Examples

data(ncgsL, package = "LMMstar")
ncgsL$group <- relevel(ncgsL$group, "placebo")

## baseline adjustment 1
ncgsL$treat <- baselineAdjustment(ncgsL, variable = "group",
               repetition= ~ visit|id, constrain = 1)
table(treat = ncgsL$treat, visit = ncgsL$visit, group = ncgsL$group)

ncgsL$treattime <- baselineAdjustment(ncgsL, variable = "group",
                   repetition= ~ visit|id, constrain = 1, collapse.time = ".")
table(treattime = ncgsL$treattime, visit = ncgsL$visit, group = ncgsL$group)

## baseline adjustment 2
ncgsL$treat2 <- baselineAdjustment(ncgsL, variable = "group",
                 new.level = "baseline",
                 repetition= ~ visit|id, constrain = 1)
table(treat = ncgsL$treat2, visit = ncgsL$visit, group = ncgsL$group)

ncgsL$treattime2 <- baselineAdjustment(ncgsL, variable = "group",
                   new.level = "baseline",
                   repetition= ~ visit|id, constrain = 1, collapse.time = ".")
table(treattime = ncgsL$treattime2, visit = ncgsL$visit, group = ncgsL$group)

Data From The Bland Altman Study (Long Format)

Description

Data From The Bland Altman Study where two methods to measure the peak expiratory flow rate (PEFR) where compared. This dataset is in the long format (i.e. one line per measurement).

• id: patient identifier.

• replicate: index of the measurement (first or second).

• method: device used to make the measurement (Wright peak flow meter or mini Wright peak flow meter).

• pefr: measurement (peak expiratory flow rate).

Usage

data(blandAltmanL)

References

Bland & Altman, Statistical methods for assessing agreement between two methods of clinical measurement, Lancet, 1986; i: 307-310.

Data From The Bland Altman Study (Wide Format)

Description

Data From The Bland Altman Study where two methods to measure the peak expiratory flow rate (PEFR) where compared. This dataset is in the wide format (i.e. one line per patient).

• id: patient identifier.

• wright1: first measurement made with a Wright peak flow meter.

• wright2: second measurement made with a Wright peak flow meter.

• mini1: first measurement made with a mini Wright peak flow meter.

• mini2: second measurement made with a mini Wright peak flow meter.

Usage

data(blandAltmanW)

References

Bland & Altman, Statistical methods for assessing agreement between two methods of clinical measurement, Lancet, 1986; i: 307-310.

Data From The Blood Pressure Study (Long Format)

Description

Data from a cross-over trial comparing the impact of three formulations of a drug on the blood pressure. The study was conducted on 12 male volunteers randomly divided into tree groups and receiving each of the three formulations with a wash-out period of one week.

• id: patient identifier.

• sequence: sequence of treatment .

• treatment: formulation of the treatment A (50 mg tablet) B (100 mg tablet) C (sustained-release formulation capsule)

• period: time period (in weeks).

• duration: duration of the drug (in hours).

Usage

data(bloodpressureL)

References

TO ADD

Data From The Calcium Supplements Study (Long Format)

Description

Data from a randomized study including 112 girls at age 11 investigate the effect of a calcium supplement (n=55) vs. placebo (n=57) on bone mineral density over a 2 year follow-up. The clinical question is: does a calcium supplement help to increase bone gain in adolescent women? This dataset is in the long format (i.e. one line per measurement).

• girl: patient identifier.

• grp: treatment group: calcium supplement (coded C) or placebo (coded P).

• visit: visit index.

• bmd: bone mineral density (mg/cm3).

• time.obs: visit time (in years).

• time.num: scheduled visit time (numeric variable, in years).

• time.fac: scheduled visit time (factor variable).

Usage

data(calciumL)

References

TO ADD

Data From The Calcium Supplements Study (Wide Format)

Description

Data from a randomized study including 112 girls at age 11 investigate the effect of a calcium supplement (n=55) vs. placebo (n=57) on bone mineral density over a 2 year follow-up. The clinical question is: does a calcium supplement help to increase bone gain in adolescent women? This dataset is in the wide format (i.e. one line per patient).

• girl: patient identifier

• grp: treatment group: calcium supplement (coded C) or placebo (coded P).

• obstime1: time after the start of the study at which the first visit took place (in years).

• obstime2: time after the start of the study at which the second visit took place (in years).

• obstime3: time after the start of the study at which the third visit took place (in years).

• obstime4: time after the start of the study at which the fourth visit took place (in years).

• obstime5: time after the start of the study at which the fifth visit took place (in years).

• bmd1: bone mineral density measured at the first visit (in mg/cm3).

• bmd2: bone mineral density measured at the second visit (in mg/cm3).

• bmd3: bone mineral density measured at the third visit (in mg/cm3).

• bmd4: bone mineral density measured at the fourth visit (in mg/cm3).

• bmd5: bone mineral density measured at the fifth visit (in mg/cm3).

Usage

data(calciumW)

References

Vonesh and Chinchilli 1997. Linear and Nonlinear models for the analysis of repeated measurement (Table 5.4.1 on page 228). New York: Marcel Dekker.

CKD long

Description

TODO

• id: patient identifier.

• allocation:

• sex:

• age:

• visit:

• time:

• pwv:

• aix:

• dropout:

Usage

data(ckdL)

References

TO ADD

CKD wide

Description

TODO

• id: patient identifier.

• allocation:

• sex:

• age:

• pwv0:

• pwv12:

• pwv24:

• aix0:

• aix12:

• aix24:

• dropout:

Usage

data(ckdW)

References

TO ADD

Extract Coefficients From a Linear Mixed Model

Description

Extract coefficients from a linear mixed model.

Usage

## S3 method for class 'lmm'
coef(
  object,
  effects = NULL,
  p = NULL,
  transform.sigma = "none",
  transform.k = "none",
  transform.rho = "none",
  transform.names = TRUE,
  ...
)

Arguments

Details

transform.sigma:

• "none" ouput residual standard error.

• "log" ouput log-transformed residual standard error.

• "square" ouput residual variance.

• "logsquare" ouput log-transformed residual variance.

transform.k:

• "none" ouput ratio between the residual standard error of the current level and the reference level.

• "log" ouput log-transformed ratio between the residual standard errors.

• "square" ouput ratio between the residual variances.

• "logsquare" ouput log-transformed ratio between the residual variances.

• "sd" ouput residual standard error of the current level.

• "logsd" ouput residual log-transformed standard error of the current level.

• "var" ouput residual variance of the current level.

• "logvar" ouput residual log-transformed variance of the current level.

transform.rho:

• "none" ouput correlation coefficient.

• "atanh" ouput correlation coefficient after tangent hyperbolic transformation.

• "cov" ouput covariance coefficient.

When using a (pure) compound symmetry covariance structure (structure = "CS"), estimated random effects can be extracted by setting argument effects to "ranef".

Value

A vector with the value of the model coefficients.

See Also

confint.lmm or model.tables.lmm for a data.frame containing estimates with their uncertainty.

Examples

## simulate data in the long format
set.seed(10)
dL <- sampleRem(100, n.times = 3, format = "long")

## fit linear mixed model
eUN.lmm <- lmm(Y ~ X1 + X2 + X5, repetition = ~visit|id, structure = "UN", data = dL, df = FALSE)

## output coefficients
coef(eUN.lmm)
coef(eUN.lmm, effects = "mean")
coef(eUN.lmm, transform.sigma = "none", transform.k = "none", transform.rho = "none")

Extract Coefficients From a Linear Mixed Model

Description

Extract coefficients from a linear mixed model.

Usage

## S3 method for class 'mlmm'
coef(object, effects = "contrast", ordering = "parameter", ...)

Arguments

Confidence Intervals for Multivariate Wald Tests

Description

Compute confidence intervals for linear hypothesis tests, possibly with adjustment for multiple comparisons.

Usage

## S3 method for class 'Wald_lmm'
confint(
  object,
  parm,
  level = 0.95,
  method = NULL,
  columns = NULL,
  backtransform = NULL,
  ...
)

Arguments

Details

Adjustment for multiple comparisons: available methods are:

• "none", "bonferroni", "single-step2"

• "holm", "hochberg", "hommel", "BH", "BY", "fdr": adjustment performed by [stats::p.adjust()], no confidence interval is computed.

• "single-step", "free", "Westfall", "Shaffer": adjustment performed by [multcomp::glht()], for all but the first method no confidence interval is computed.

Note: method "single-step" adjusts for multiple comparisons using equicoordinate quantiles of the multivariate Student's t-distribution over all tests, instead of the univariate quantiles. It assumes equal degrees of freedom in the marginal and is described in section 7.1 of Dmitrienko et al. (2013) under the name single-step Dunnett procedure. The name "single-step" is borrowed from the multcomp package. In the book Bretz et al. (2010) written by the authors of the package, the procedure is refered to as max-t tests which is the terminology adopted in the LMMstar package.
When degrees of freedom differs between individual hypotheses, method "single-step2" is recommended. It simulates data using copula whose marginal distributions are Student's t-distribution (with possibly different degrees of freedom) and elliptical copula with parameters the estimated correlation between the test statistics (via the copula package). It then computes the frequency at which the simulated maximum exceed the observed maximum and appropriate quantile of simulated maximum for the confidence interval.

Pooling estimates: available methods are:

• "average": average estimates

• "pool.fixse": weighted average of the estimates, with weights being the inverse of the squared standard error. The uncertainty about the weights is neglected.

• "pool.se": weighted average of the estimates, with weights being the inverse of the squared standard error. The uncertainty about the weights is computed under independence of the variance parameters between models.

• "pool.gls": weighted average of the estimates, with weights being based on the variance-covariance matrix of the estimates. When this matrix is singular, its spectral decomposition is truncated when the correspodning eigenvalues are below 10^{-12}. The uncertainty about the weights is neglected.

• "pool.gls1": similar to "pool.gls" but ensure that the weights are at most 1 in absolute value by shrinking toward the average.

• "pool.rubin": average of the estimates and compute the uncertainty according to Rubin's rule (Barnard et al. 1999).

References

Barnard and Rubin, Small-sample degrees of freedom with multiple imputation. Biometrika (1999), 86(4):948-955.
Dmitrienko, A. and D'Agostino, R., Sr (2013), Traditional multiplicity adjustment methods in clinical trials. Statist. Med., 32: 5172-5218. Frank Bretz, Torsten Hothorn and Peter Westfall (2010), Multiple Comparisons Using R, CRC Press, Boca Raton.

Statistical Inference for Linear Mixed Model

Description

Compute confidence intervals (CIs) and p-values for the coefficients of a linear mixed model.

Usage

## S3 method for class 'lmm'
confint(
  object,
  parm = NULL,
  level = 0.95,
  effects = NULL,
  robust = FALSE,
  null = NULL,
  columns = NULL,
  df = NULL,
  type.information = NULL,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  transform.names = TRUE,
  backtransform = NULL,
  ...
)

Arguments

Value

A data.frame containing some of the following coefficient (in rows):

• column estimate: the estimate.

• column se: the standard error.

• column statistic: the test statistic.

• column df: the degree of freedom.

• column lower: the lower bound of the confidence interval.

• column upper: the upper bound of the confidence interval.

• column null: the null hypothesis.

• column p.value: the p-value relative to the null hypothesis.

See Also

the function anova to perform inference about linear combinations of coefficients and adjust for multiple comparisons.

coef.lmm for a simpler output (e.g. only estimates).
model.tables.lmm for a more detailed output (e.g. with p-value).

Examples

#### simulate data in the long format ####
set.seed(10)
dL <- sampleRem(100, n.times = 3, format = "long")

#### fit Linear Mixed Model ####
eUN.lmm <- lmm(Y ~ X1 + X2 + X5, repetition = ~visit|id, structure = "UN", data = dL)

#### Confidence intervals ####
## based on a Student's t-distribution with transformation
confint(eUN.lmm, effects = "all")
## based on a Student's t-distribution without transformation
confint(eUN.lmm, effects = "all",
        transform.sigma = "none", transform.k = "none", transform.rho = "none")
## based on a Student's t-distribution transformation but not backtransformed
confint(eUN.lmm, effects = "all", backtransform = FALSE)
## based on a Normal distribution with transformation
confint(eUN.lmm, df = FALSE)

Confidence Intervals for Multiple Linear Mixed Model.

Description

Compute confidence intervals for several linear mixed models.

Usage

## S3 method for class 'mlmm'
confint(
  object,
  parm = NULL,
  level = 0.95,
  method = NULL,
  ordering = "parameter",
  ...
)

Arguments

Details

Statistical inference following pooling is performed according to Rubin's rule whose validity requires the congeniality condition of Meng (1994).

References

Meng X. L.(1994). Multiple-imputation inferences with uncongenial sources of input. Statist. Sci.9, 538–58.

Residuals Degrees of Freedom

Description

Residuals degrees of freedom. Computed as the sum of squared normalized residuals

Usage

## S3 method for class 'lmm'
df.residual(object, ...)

Arguments

Value

A numeric value

Effects Derived For Linear Mixed Model

Description

Estimate average counterfactual outcome or contrast of outcome from linear mixed models.

Usage

## S3 method for class 'lmm'
effects(
  object,
  variable,
  newdata = NULL,
  type = c("identity", "none"),
  conditional = NULL,
  rhs = NULL,
  repetition = NULL,
  multivariate = FALSE,
  prefix.time = NULL,
  prefix.var = TRUE,
  sep.var = ",",
  ...
)

Arguments

Examples

#### simulate data in the long format ####
set.seed(10)
dL <- sampleRem(100, n.times = 3, format = "long")

#### Linear Mixed Model ####
eUN.lmm <- lmm(Y ~ visit + X1 + X2 + X5,
               repetition = ~visit|id, structure = "UN", data = dL)

## outcome
effects(eUN.lmm, variable = "X1")
effects(eUN.lmm, type = "difference", variable = "X1")
effects(eUN.lmm, type = "difference", variable = "X1", repetition = "3")

## change
effects(eUN.lmm, type = "change", variable = "X1")
effects(eUN.lmm, type = "change", variable = "X1", conditional = NULL)
effects(eUN.lmm, type = c("change","difference"), variable = "X1")

## auc
effects(eUN.lmm, type = "auc", variable = "X1")
effects(eUN.lmm, type = c("auc","difference"), variable = "X1")

#### fit Linear Mixed Model with interaction ####
dL$X1.factor <- as.factor(dL$X1)
dL$X2.factor <- as.factor(dL$X2)
eUN.lmmI <- lmm(Y ~ visit * X1.factor + X2.factor + X5,
               repetition = ~visit|id, structure = "UN", data = dL)

## average counterfactual conditional to a categorical covariate
effects(eUN.lmmI, variable = "X1.factor",
        conditional = c("X2.factor"), repetition = "3")
effects(eUN.lmmI, type = "change", variable = "X1.factor",
        conditional = c("X2.factor"), repetition = "3")
effects(eUN.lmmI, type = "auc", variable = "X1.factor",
        conditional = c("X2.factor"), repetition = "3")

## average difference in counterfactual conditional to a categorical covariate
effects(eUN.lmmI, type = "difference", variable = "X1.factor",
        conditional = c("X2.factor"), repetition = "3")
effects(eUN.lmmI, type = c("change","difference"), variable = "X1.factor",
        conditional = c("X2.factor"), repetition = "3")
effects(eUN.lmmI, type = c("auc","difference"), variable = "X1.factor",
        conditional = c("X2.factor"), repetition = "3")

## average difference in counterfactual conditional to a covariate
effects(eUN.lmmI, type = "difference", variable = "X1.factor",
        conditional = list(X5=0:2), repetition = "3")
effects(eUN.lmmI, type = c("difference","change"), variable = "X1.factor",
        conditional = list(X5=0:2))

Extract the Score Function for Multcomp

Description

Extract the Score Function for Multcomp. For internal use.

Usage

## S3 method for class 'lmm'
estfun(x, ...)

Arguments

Value

A matrix containing the score function for each model parameter (columns) relative to each cluster (rows).

Examples

## simulate data in the long format
set.seed(10)
dL <- sampleRem(100, n.times = 3, format = "long")

## fit Linear Mixed Model
eUN.lmm <- lmm(Y ~ X1 + X2 + X5, repetition = ~visit|id, structure = "UN", data = dL, df = FALSE)

## test multiple linear hypotheses
if(require(multcomp)){
LMMstar.options(effects = c("mean"))
e.glht <- multcomp::glht(eUN.lmm)
e.glht$linfct
}

Delta Method for Mixed Models

Description

Perform a first order delta method

Usage

## S3 method for class 'lmm'
estimate(
  x,
  f,
  df = !is.null(x$df),
  robust = FALSE,
  type.information = NULL,
  level = 0.95,
  method.numDeriv = NULL,
  average = FALSE,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  ...
)

Arguments

Examples

if(require(lava) && require(nlme)){

#### Random effect ####
set.seed(10)
dL <- sampleRem(1e2, n.times = 3, format = "long")
e.lmm1 <- lmm(Y ~ X1+X2+X3 + (1|id), repetition = ~visit|id, data = dL)
nlme::ranef(e.lmm1, se = TRUE)
e.ranef <- estimate(e.lmm1, f  = function(p){nlme::ranef(e.lmm1, p = p)})
e.ranef

if(require(ggplot2)){
df.gg <- cbind(index = 1:NROW(e.ranef), e.ranef)
gg.ranef <- ggplot(df.gg, aes(x = index, y=estimate, ymin=lower, ymax = upper))
gg.ranef + geom_point() + geom_errorbar() + ylab("estimated random effect") + xlab("id")
}

#### ANCOVA via mixed model ####
set.seed(10)
d <- sampleRem(1e2, n.time = 2)
e.ANCOVA1 <- lm(Y2~Y1+X1, data = d)

if(require(reshape2)){
   dL2 <- melt(d, id.vars = c("id","X1"),  measure.vars = c("Y1","Y2"),
               value.name = "Y", variable.name = "time")
   dL2$time <- factor(dL2$time, levels = c("Y1","Y2"), labels = c("1","2"))

   ## estimated treatment effect (no baseline constraint)
   e.lmm <- lmm(Y ~ time + time:X1, data = dL2, repetition = ~time|id)

   e.delta <- estimate(e.lmm, function(p){
       c(Y1 = p["rho(1,2)"]*p["k.2"],
         X1 = p["time2:X1"]-p["k.2"]*p["rho(1,2)"]*p["time1:X1"])
   }) ## same estimate and similar standard errors. 
   e.delta ## Degrees of freedom are a bit off though
   cbind(summary(e.ANCOVA1)$coef, df = df.residual(e.ANCOVA1))

   ## estimated treatment effect (baseline constraint)
   dL2$time2 <- as.numeric(dL2$time=="2")
   e.lmmC <- lmm(Y ~ time2 + time2:X1, data = dL2, repetition = ~time|id)
   e.deltaC <- estimate(e.lmmC, function(p){
       c(Y1 = p["rho(1,2)"]*p["k.2"],
         X1 = p["time2:X1"])
   })
   e.deltaC ## Degrees of freedom are a bit more accurate
}

}

Predicted Mean Value For Linear Mixed Model.

Description

Evaluate the expected mean conditional to covariates or the expected outcome values when missing conditional to observed outcome and covariates. Similar to predict.lmm where the values to condition on are, by default, taking from the dataset used to fit the Linear Mixed Model.

Usage

## S3 method for class 'lmm'
fitted(
  object,
  newdata = NULL,
  type = "mean",
  se = NULL,
  df = NULL,
  keep.data = NULL,
  format = "long",
  seed = NULL,
  simplify = TRUE,
  ...
)

Arguments

Value

When format="wide", a data.frame with as many rows as clusters. When format="long", a data.frame with as many rows as observations (keep.data==TRUE) or a vector of length the number of observations (keep.data==TRUE).

Examples

#### single arm trial ####
data(gastricbypassL, package = "LMMstar")
gastricbypassL <- gastricbypassL[order(gastricbypassL$id,gastricbypassL$visit),]
gastricbypassL$weight0 <- unlist(tapply(gastricbypassL$weight,gastricbypassL$id,
function(x){rep(x[1],length(x))}))

eUN.lmm <- lmm(glucagonAUC ~ visit + weight0, repetition = ~visit|id,
               data = gastricbypassL, df = FALSE)

## fitted mean (conditional on covariates only)
fitted(eUN.lmm)
fitted(eUN.lmm, newdata = data.frame(visit = "3", weight0 = 0))
fitted(eUN.lmm, newdata = data.frame(visit = "3", weight0 = 0),
       keep.data = TRUE)

## fitted outcome value (conditional on covariates and covariates)
fitted(eUN.lmm, type = "outcome")
gastricbypassL.O <- fitted(eUN.lmm, type = "outcome", keep.data = TRUE)

if(require(ggplot2)){
gg.outcome <- ggplot(gastricbypassL.O,
                     aes(x=time, y = glucagonAUC, color = impute, group = id))
gg.outcome <- gg.outcome + geom_point() + geom_line()## + facet_wrap(~id)
gg.outcome
}

tapply(gastricbypassL.O$glucagonAUC, gastricbypassL.O$time, mean)
effects(eUN.lmm, variable = NULL)

## fitted change value (conditional on covariates and covariates)
gastricbypassL.C <- fitted(eUN.lmm, type = "change", keep.data = TRUE)

if(require(ggplot2)){
gg.change <- ggplot(gastricbypassL.C,
                    aes(x=time, y = glucagonAUC, color = impute, group = id))
gg.change <- gg.change + geom_point() + geom_line()
gg.change
}

tapply(gastricbypassL.C$glucagonAUC, gastricbypassL.O$time, mean)
effects(eUN.lmm, type = "change", variable = NULL)

## fitted auc (conditional on covariates and covariates)
gastricbypassL.AUC <- fitted(eUN.lmm, type = "auc", keep.data = TRUE)

if(require(ggplot2)){
gg.auc <- ggplot(gastricbypassL.AUC,
                    aes(x = "auc", y = glucagonAUC, color = impute))
gg.auc <- gg.auc + geom_point()
gg.auc
}

mean(gastricbypassL.AUC$glucagonAUC)
effects(eUN.lmm, type = "auc", variable = NULL)

#### two arm trial ####
## Not run: 
if(require(nlmeU) & require(reshape2)){
data(armd.wide, package = "nlmeU")
armd.long <- melt(armd.wide,
                  measure.vars = paste0("visual",c(0,4,12,24,52)),
                  id.var = c("subject","lesion","treat.f","miss.pat"),
                  variable.name = "week",
                  value.name = "visual")

armd.long$week <- factor(armd.long$week, 
                         level = paste0("visual",c(0,4,12,24,52)),
                         labels = c(0,4,12,24,52))

eUN2.lmm <- lmm(visual ~ treat.f*week + lesion,
               repetition = ~week|subject, structure = "UN",
               data = armd.long)

## fitted outcome value (conditional on covariates and covariates)
armd.O <- fitted(eUN2.lmm, type = "outcome", keep.data = TRUE)

gg2.outcome <- ggplot(armd.O,
                     aes(x=week, y = visual, color = impute, group = subject))
gg2.outcome <- gg2.outcome + geom_point() + geom_line() + facet_wrap(~treat.f)
gg2.outcome

aggregate(visual ~ week + treat.f, FUN = mean, data = armd.O)
effects(eUN2.lmm, variable = "treat.f") ## mismatch due to adjustment on lesion

## fitted change value (conditional on covariates and covariates)
armd.C <- fitted(eUN2.lmm, type = "change", keep.data = TRUE)

gg.change <- ggplot(armd.C,
                    aes(x=week, y = visual, color = impute, group = subject))
gg.change <- gg.change + geom_point() + geom_line() + facet_wrap(~treat.f)
gg.change

coef(eUN2.lmm)
effects(eUN2.lmm, type = "change", variable = "treat.f")
effects(eUN2.lmm, type = c("change","difference"), variable = "treat.f")

## fitted auc (conditional on covariates and covariates)
armd.AUC <- fitted(eUN2.lmm, type = "auc", keep.data = TRUE)

gg.auc <- ggplot(armd.AUC, aes(x = treat.f, y = visual, color = impute))
gg.auc <- gg.auc + geom_point()
gg.auc

aggregate(visual ~ treat.f, data = armd.AUC, FUN = "mean")
effects(eUN2.lmm, type = "auc", variable = "treat.f") ## adjusted for lesion
effects(eUN2.lmm, type = c("auc","difference"), variable = "treat.f")
}
## End(Not run)

Data From The Gastric Bypass Study (Long Format)

Description

Data from the gastric bypass study where the bodyweight and serum glucagon (a gut hormone) were measured in 20 obese subjects prior and after gastric bypass surgery. This dataset is in the long format (i.e. one line per measurement).

• id: patient identifier.

• visit: the visit index (factor).

• time: the week at which the visit took place (numeric).

• weight: bodyweight (in kg) measured during the visit.

• glucagonAUC: glucagon measured during the visit.

Usage

data(gastricbypassL)

References

The effect of Roux-en-Y gastric bypass surgery on the gut mucosal gene expression profile and circulating gut hormones. https://easddistribute.m-anage.com/from.storage?image=4iBH9mRQm1kfeEHULC2CxovdlyCtA1EHeVDdoffnZrAUGG9SHTO-U4ItnLU078eVkF1ZUZgYTy7THlTW3KSgFA2

Data From The Gastric Bypass Study (Wide Format)

Description

Data from the gastric bypass study where the bodyweight and serum glucagon (a gut hormone) were measured in 20 obese subjects prior and after gastric bypass surgery. This dataset is in the wide format (i.e. one line per patient).

• id: patient identifier.

• weight1: bodyweight (in kg) 3 months before surgery.

• weight2: bodyweight (in kg) 1 week before surgery.

• weight3: bodyweight (in kg) 1 week after surgery.

• weight4: bodyweight (in kg) 3 months after surgery.

• glucagonAUC1: glucagon (in pmol/l x hours) 3 months before surgery.

• glucagonAUC2: glucagon (in pmol/l x hours) 1 week before surgery.

• glucagonAUC3: glucagon (in pmol/l x hours) 1 week after surgery.

• glucagonAUC4: glucagon (in pmol/l x hours) 3 months after surgery.

Usage

data(gastricbypassW)

References

The effect of Roux-en-Y gastric bypass surgery on the gut mucosal gene expression profile and circulating gut hormones. https://easddistribute.m-anage.com/from.storage?image=4iBH9mRQm1kfeEHULC2CxovdlyCtA1EHeVDdoffnZrAUGG9SHTO-U4ItnLU078eVkF1ZUZgYTy7THlTW3KSgFA2

Extract the Influence Function From a Linear Mixed Model

Description

Extract the influence function from a linear mixed model.

Usage

## S3 method for class 'lmm'
iid(
  x,
  effects = "mean",
  robust = TRUE,
  type.information = NULL,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  transform.names = TRUE,
  ...
)

Arguments

Extract The Information From a Linear Mixed Model

Description

Extract or compute the (expected) second derivative of the log-likelihood of a linear mixed model.

Usage

## S3 method for class 'lmm'
information(
  x,
  effects = NULL,
  newdata = NULL,
  p = NULL,
  indiv = FALSE,
  type.information = NULL,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  transform.names = TRUE,
  ...
)

Arguments

Details

For details about the arguments transform.sigma, transform.k, transform.rho, see the documentation of the coef.lmm function.

Value

When argument indiv is FALSE, a matrix with the value of the infroamtion relative to each pair of coefficient (in rows and columns) and each cluster (in rows). When argument indiv is TRUE, a 3-dimensional array with the value of the information relative to each pair of coefficient (dimension 2 and 3) and each cluster (dimension 1).

Contrasts and Reference Level

Description

Contrasts and reference level used when modeling the mean in a linear mixed modek.

Usage

## S3 method for class 'lmm'
levels(x)

Arguments

Value

a list with two elements

• all: contrast matrix for each categorical or factor variable

• reference: reference level: one value for each categorical variable

Fit Linear Mixed Model

Description

Fit a linear mixed model defined by a mean and a covariance structure.

Usage

lmm(
  formula,
  repetition,
  structure,
  data,
  weights = NULL,
  scale.Omega = NULL,
  method.fit = NULL,
  df = NULL,
  type.information = NULL,
  trace = NULL,
  control = NULL
)

Arguments

Details

Computation time the lmm has not been developped to be a fast function as, by default, it uses REML estimation with the observed information matrix and uses a Satterthwaite approximation to compute degrees of freedom (this require to compute the third derivative of the log-likelihood which is done by numerical differentiation). The computation time can be substantially reduced by using ML estimation with the expected information matrix and no calculation of degrees of freedom: arguments method.fit="ML", type.information="expected", df=FALSE. This will, however, lead to less accurate p-values and confidence intervals in small samples.

By default, the estimation of the model parameters will be made using a Newton Raphson algorithm. This algorithm does not ensure that the residual covariance matrix is positive definite and therefore may sometimes fail. See argument optimizer in LMMstar.options.

Argument control: when using the optimizer "FS", the following elements can be used

• init: starting values for the model parameters.

• n.iter: maximum number of interations of the optimization algorithm.

• tol.score: score value below which convergence has been reached.

• tol.param: difference in estimated parameters from two successive iterations below which convergence has been reached.

• trace: display progress of the optimization procedure.

Argument repetition: when numeric, it will be converted into a factor variable, possibly adding a leading 0 to preserve the ordering. This transformation may cause inconsistency when combining results between different lmm object. This is why the grouping variable should preferably be of type character or factor.

Value

an object of class lmm containing the estimated parameter values, the residuals, and relevant derivatives of the likelihood.

See Also

summary.lmm for a summary of the model fit.
model.tables.lmm for a data.frame containing estimates with their uncertainty.
plot.lmm for a graphical display of the model fit or diagnostic plots.
levels.lmm to display the reference level.
anova.lmm for testing linear combinations of coefficients (F-test, multiple Wald tests)
effects.lmm for evaluating average marginal or counterfactual effects
sigma.lmm for extracting estimated residual variance-covariance matrices.
residuals.lmm for extracting residuals or creating residual plots (e.g. qqplots).
predict.lmm for evaluating mean and variance of the outcome conditional on covariates or other outcome values.

Examples

#### 1- simulate data in the long format ####
set.seed(10)
dL <- sampleRem(100, n.times = 3, format = "long")
dL$X1 <- as.factor(dL$X1)
dL$X2 <- as.factor(dL$X2)

#### 2- fit Linear Mixed Model ####
eCS.lmm <- lmm(Y ~ X1 * X2 + X5, repetition = ~visit|id, structure = "CS", data = dL)

logLik(eCS.lmm) ## -670.9439
summary(eCS.lmm)


#### 3- estimates ####
## reference level
levels(eCS.lmm)$reference
## mean parameters
coef(eCS.lmm)
model.tables(eCS.lmm)
confint(eCS.lmm)

## all parameters
coef(eCS.lmm, effects = "all")
model.tables(eCS.lmm, effects = "all")
confint(eCS.lmm, effects = "all")

## variance-covariance structure
sigma(eCS.lmm)

#### 4- diagnostic plots ####
quantile(residuals(eCS.lmm))
quantile(residuals(eCS.lmm, type = "normalized"))

## Not run: 
if(require(ggplot2)){
  ## investigate misspecification of the mean structure
  plot(eCS.lmm, type = "scatterplot")
  ## investigate misspecification of the variance structure
  plot(eCS.lmm, type = "scatterplot2")
  ## investigate misspecification of the correlation structure
  plot(eCS.lmm, type = "correlation")
  ## investigate misspecification of the residual distribution
  plot(eCS.lmm, type = "qqplot")
}

## End(Not run)

#### 5- statistical inference ####
anova(eCS.lmm) ## effect of each variable
anova(eCS.lmm, effects = "X11-X21=0") ## test specific coefficient
## test several hypothese with adjustment for multiple comparisons
summary(anova(eCS.lmm, effects = c("X11=0","X21=0")))

#### 6- prediction ####
## conditional on covariates
newd <- dL[1:3,]
predict(eCS.lmm, newdata = newd, keep.data = TRUE)
## conditional on covariates and outcome
newd <- dL[1:3,]
newd$Y[3] <- NA
predict(eCS.lmm, newdata = newd, type = "dynamic", keep.data = TRUE)

#### EXTRA ####
if(require(mvtnorm)){
## model for the average over m replicates
## (only works with independent replicates)
Sigma1 <- diag(1,1,1); Sigma5 <- diag(1,5,5)
n <- 100
dfW <- rbind(data.frame(id = 1:n, rep = 5, Y = rowMeans(rmvnorm(n, sigma = Sigma5))),
             data.frame(id = (n+1):(2*n), rep = 1, Y = rmvnorm(n, sigma = Sigma1)))

e.lmmW <- lmm(Y~1, data = dfW, scale.Omega=~rep, control = list(optimizer = "FS"))
e.lmm0 <- lmm(Y~1, data = dfW, control = list(optimizer = "FS"))
model.tables(e.lmmW, effects = "all")
model.tables(e.lmm0, effects = "all")
## TRUE standard error is 1

}

Fit Linear Mixed Model on Complete Case data

Description

Fit a linear mixed model on the complete case data. Mostly useful as a sanity check, to match the results of a univariate analysis on the change.

Usage

lmmCC(object, ...)

## S3 method for class 'formula'
lmmCC(
  object,
  repetition,
  data,
  lm.change = FALSE,
  df = NULL,
  trace = TRUE,
  control = NULL,
  ...
)

## S3 method for class 'lm'
lmmCC(
  object,
  repetition,
  data,
  name.time = "time",
  df = NULL,
  trace = TRUE,
  control = NULL,
  ...
)

Arguments

Value

A lmmCC object, which inherits from lmm.

Examples

#### 1- simulate data in the wide format ####
set.seed(10)
dW <- sampleRem(100, n.times = 3, format = "wide")
dW$Y3[1:10] <- NA
dW$change2 <- dW$Y2 - dW$Y1
dW$change3 <- dW$Y3 - dW$Y1

e.lm2 <- lm(change2 ~ X1 + X2, data = dW)
summary(e.lm2)$coef
e.lm3 <- lm(change3 ~ X1 + X2, data = dW)
summary(e.lm3)$coef

#### 2- complete case LMM from LM ####
e.lmmCC2 <- lmmCC(e.lm2, repetition = change2~Y2-Y1)
model.tables(e.lmmCC2)
e.lmmCC3 <- lmmCC(e.lm3, repetition = change3~Y3-Y1)
model.tables(e.lmmCC3)

#### 3- complete case LMM ####
dL <- reshape(dW[,c("id","X1","X2","Y1","Y2","Y3")],
              direction = "long", 
              varying = c("Y1","Y2","Y3"), sep = "", idvar = "id")
dL$time <- as.character(dL$time)

e.lmm2 <- lmmCC(Y ~ time + time*X1 + time*X2, repetition = ~time|id,
                data = dL[dL$time %in% 1:2,])
model.tables(e.lmm2)
e.lmm3.bis <- lmmCC(Y ~ time + time*X1 + time*X2, repetition = ~time|id,
                data = dL[dL$time %in% c(1,3),], lm.change = TRUE)
model.tables(e.lmm3.bis)
e.lmm3.bis$lm.change

Extract The Log-Likelihood From a Linear Mixed Model

Description

Extract or compute the log-likelihood of a linear mixed model.

Usage

## S3 method for class 'lmm'
logLik(object, newdata = NULL, p = NULL, indiv = FALSE, ...)

Arguments

Details

indiv: only relevant when using maximum likelihood. Must be FALSE when using restricted maximum likelihood.

Value

A numeric value (total logLikelihood) or a vector of numeric values, one for each cluster (cluster specific logLikelihood).

Variables Involved in a Linear Mixed Model

Description

Extract the variables used by the linear mixed model.

Usage

## S3 method for class 'lmm'
manifest(x, effects = "all", original = TRUE, simplify = TRUE, ...)

Arguments

Value

A list of character vectors or a character vector.

Fit Multiple Linear Mixed Model

Description

Fit several linear mixed models, extract relevant coefficients, and combine them into a single table.

Usage

mlmm(
  ...,
  data,
  by,
  contrast.rbind = NULL,
  effects = NULL,
  robust = FALSE,
  df = TRUE,
  ci = TRUE,
  name.short = c(TRUE, TRUE),
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  transform.names = TRUE,
  trace = TRUE
)

Arguments

Details

Grouping variable in argument repetition: when numeric, it will be converted into a factor variable, possibly adding a leading 0 to preserve the ordering. This transformation may cause inconsistency when combining results between different lmm object. This is why the grouping variable should preferably be of type character or factor.

See Also

confint.mlmm for a data.frame containing estimates with their uncertainty.
summary.mlmm for a summary of the model and estimates.
autoplot.Wald_lmm for a graphical display.

Examples

#### univariate regression ####
if(require(lava) && require(multcomp)){

set.seed(10)
d1 <- cbind(sim(lvm(Y~0.5*X1), 25), group = "A")
d2 <- cbind(sim(lvm(Y~0.1*X1), 100), group = "B")
d3 <- cbind(sim(lvm(Y~0.01*X1), 1000), group = "C")
d1$id <- 1:NROW(d1)
d2$id <- 1:NROW(d2)
d3$id <- 1:NROW(d3)

d <- rbind(d1,d2,d3)

e.mlmm <- mlmm(Y~X1, data = d, by = "group", effects = "X1=0")
summary(e.mlmm)
summary(e.mlmm, method = "single-step")
summary(e.mlmm, method = "single-step2")

## re-work contrast
summary(anova(e.mlmm, effects = mcp(X1 = "Dunnett")), method = "none")
## summary(mlmm(Y~X1, data = d, by = "group", effects = mcp(X1="Dunnett")))
}

#### multivariate regression ####
set.seed(10)
dL <- sampleRem(250, n.times = 3, format = "long")

e.mlmm <- mlmm(Y~X1+X2+X6, repetition = ~visit|id, data = dL,
               by = "X4", structure = "CS")
summary(e.mlmm)

e.mlmmX1 <- mlmm(Y~X1+X2+X6, repetition = ~visit|id, data = dL,
               by = "X4", effects = "X1=0", structure = "CS")
summary(e.mlmmX1)
summary(e.mlmmX1, method = "single-step")

Extracting the Model Frame from a Linear Mixed Model

Description

Variables needed to fit the Linear Mixed Model.

Usage

## S3 method for class 'lmm'
model.frame(
  formula,
  newdata = NULL,
  type = NULL,
  add.index = FALSE,
  na.rm = TRUE,
  ...
)

Arguments

Details

Column "XXindexXX" refers to the row of the original dataset (i.e. passed to argument data when calling lmm). When adding rows relative to missing repetitions, since there is no row in the original dataset, a negative sign is used.

Examples

data("armd.wide", package = "nlmeU")
e.lmH <- lmm(visual52 ~ lesion, structure = IND(~treat.f), data = armd.wide)
model.frame(e.lmH)
model.frame(e.lmH, add.index = TRUE)
model.frame(e.lmH, type = "unique")

data("gastricbypassL", package = "LMMstar")
dfL.NNA <- na.omit(gastricbypassL)
e.lmm <- lmm(glucagonAUC ~ time, repetition = ~visit|id, data = dfL.NNA, df = FALSE)
model.frame(e.lmm, type = "unique")
model.frame(e.lmm, type = c("unique","correlation"))
model.frame(e.lmm, type = "add.NA", add.index = TRUE)

Design Matrix for Linear Mixed Model

Description

Extract or construct design matrices for Linear Mixed Model.

Usage

## S3 method for class 'lmm'
model.matrix(
  object,
  newdata = NULL,
  effects = "mean",
  simplify = TRUE,
  drop.X = NULL,
  na.rm = TRUE,
  ...
)

Arguments

Value

When simplify is FALSE, a list with the followin elements:

• mean: design matrix for the mean model

• Y: vector of outcome values

• vcov: list of elements for the variance and correlation models.

• index.cluster: list containing, for each cluster, the location of the corresponding observations in the processed dataset.

• index.clusterTime: list containing, for each cluster, the repetition index corresponding observations.

• index.clusterStrata: list containing, for each cluster, the strata index corresponding observations.

• param: data.frame describing the modle parameters.

• drop.X: logical value indicating whether columns in the design matrix should be dropped if it has not full rank.

• precompute.XX, precompute.XY: moments of X and Y

When simplify is TRUE, this list will be simplified into a list with three elements:

• mean: design matrix for the mean model

• variance: design matrix for the variance model

• correlation: design matrix for the correlation model

or a single design matrixx.

Statistical Inference for Linear Mixed Model

Description

Export estimates, standard errors, degrees of freedom, confidence intervals (CIs) and p-values for the mean coefficients of a linear mixed model.

Usage

## S3 method for class 'lmm'
model.tables(x, columns, ...)

Arguments

Details

This function simply calls confint with a specific value for the argument column.

Multiple Student's t-Test

Description

Perform multiple Student's t-Test via heteroschedastic linear regression and combine the results, possibly adjusted for multiplicity.

Usage

mt.test(formula, data, method = NULL, level = 0.95, trace = TRUE)

Arguments

Details

In presence of missing values, performs a outcome specific complete case analysis.

Value

A data.frame with the estimates, confidence intervals, and p-values relative to each outcome. Depending on the argument method confidence intervals and p-values may be adjusted for multiple comparisons. The data.frame has an attribute mlmm containing the underlying regression models.

Examples

data(calciumW, package = "LMMstar")

t.test(bmd1 ~ grp, data = calciumW)

mt.test(bmd1+bmd2+bmd3+bmd4+bmd5 ~ grp, data = calciumW)
mt.test(bmd1+bmd2+bmd3+bmd4+bmd5 ~ grp|girl, data = calciumW)
mt.test(bmd1+bmd2+bmd3+bmd4+bmd5 ~ grp|girl, data = calciumW, method = "none")

Data From National Cooperative Gallstone Study (Long Format)

Description

Data from the National Cooperative Gallstone Study (NCGS), a randomized study where the level of serum cholesterol was measured at baseline and after intake of high-dose chenondiol (750mg/day) or placebo. This dataset is in the long format (i.e. one line per measurement).

• group: treatment group (highdose or placebo).

• id: patient identifier.

• visit: visit index.

• cholest: cholesterol measurement.

• time: time after the start of the study at which the measurement has been done (in month). Treatment is given at 0+.

Usage

data(ncgsL)

References

Grundy SM, Lan SP, Lachin J. The effects of chenodiol on biliary lipids and their association with gallstone dissolution in the National Cooperative Gallstone Study (NCGS). J Clin Invest. 1984 Apr;73(4):1156-66. doi: 10.1172/JCI111301.

Data From National Cooperative Gallstone Study (Wide Format)

Description

Data from the National Cooperative Gallstone Study (NCGS), a randomized study where the level of serum cholesterol was measured at baseline and after intake of high-dose chenondiol (750mg/day) or placebo. This dataset is in the wide format (i.e. one line per patient).

• group: treatment group (highdose or placebo).

• id: patient identifier.

• cholest1: cholesterol measurement at baseline (before treatment).

• cholest2: cholesterol measurement at 6 months (after treatment).

• cholest3: cholesterol measurement at 12 months (after treatment).

• cholest4: cholesterol measurement at 20 months (after treatment).

• cholest5: cholesterol measurement at 24 months (after treatment).

Usage

data(ncgsW)

References

Grundy SM, Lan SP, Lachin J. The effects of chenodiol on biliary lipids and their association with gallstone dissolution in the National Cooperative Gallstone Study (NCGS). J Clin Invest. 1984 Apr;73(4):1156-66. doi: 10.1172/JCI111301.

Number of Observations from a Linear Mixed Model

Description

Extract the number of observations from a Linear Mixed Model

Usage

## S3 method for class 'lmm'
nobs(object, ...)

Arguments

Value

A vector with 4 elements:

• obs: the number of repetitions with full data

• cluster: the number of clusters with a least one repetition with full data

• missing.obs: the number of repetitions with missing data

• missing.cluster: the number of cluster with only missing data

Partial Correlation

Description

Estimate the partial correlation based on equation 19 of Lloyd et al 2008 (partialCor.lmm) or explicitely modeling the correlation via a linear mixed model (partialCor.list, partialCor.formula). The first option is numerically more efficient and exact with a single observation per cluster. With multiple repetitions, what is being estimated with the first option may not be clear and the second option is therefore preferrable.

Usage

partialCor(object, ...)

## S3 method for class 'list'
partialCor(
  object,
  data,
  repetition = NULL,
  structure = NULL,
  by = NULL,
  effects = NULL,
  rhs = NULL,
  method = "none",
  df = NULL,
  transform.rho = NULL,
  name.short = c(TRUE, FALSE),
  ...
)

## S3 method for class 'formula'
partialCor(object, repetition, ...)

## S3 method for class 'lmm'
partialCor(object, level = 0.95, R2 = FALSE, se = TRUE, df = TRUE, ...)

Arguments

Details

Fit a mixed model to estimate the partial correlation with the following variance-covariance pattern:

• no repetition: unstructure or compound symmetry structure for M observations, M being the number of variables on the left hand side (i.e. outcomes).

• repetition: structure for M*T observations where M being the number of variables (typically 2) and T the number of repetitions. Can be

  • "UN": unstructured (except the off-diagonal containing the correlation parameter which is constant).

  • "PEARSON": same as unstructured except it only uses a single variance parameter per variable, i.e. it assumes constant variance over repetitions.

  • "HLAG": toeplitz by block with variable and repetition specific variance.

  • "LAG": toeplitz by block, i.e. correlation depending on the gap between repetitions and specific to each variable. It assumes constant variance over repetitions.

  • "HCS": heteroschedastic compound symmetry by block, i.e. variable specific correlation constant over repetitions. A specific parameter is used for the off-diagonal crossing the variables at the same repetition (which is the marginal correlation parameter).

  • "CS": compound symmetry by block. It assumes constant variance and correlation over repetitions.

Value

A data.frame with the estimate partial correlation (rho), standard error, degree of freedom, confidence interval, and p-value (test of no correlation). When structure="CS" or structure="HCS" is used with repeated measurements, a second correlation coefficient (r) is output where the between subject variance has been removed (similar to Bland et al. 1995).

References

Bland J M, Altman D G. Statistics notes: Calculating correlation coefficients with repeated observations: Part 1—correlation within subjects BMJ 1995; 310 :446 doi:10.1136/bmj.310.6977.446 Edwards, L.J., Muller, K.E., Wolfinger, R.D., Qaqish, B.F. and Schabenberger, O. (2008), An R2 statistic for fixed effects in the linear mixed model. Statist. Med., 27: 6137-6157. https://doi.org/10.1002/sim.3429

Examples

#### no repetition ####

## example from ppcor::pcor 
y.data <- data.frame(
  hl=c(7,15,19,15,21,22,57,15,20,18),
  disp=c(0.000,0.964,0.000,0.000,0.921,0.000,0.000,1.006,0.000,1.011),
  deg=c(9,2,3,4,1,3,1,3,6,1),
  BC=c(1.78e-02,1.05e-06,1.37e-05,7.18e-03,0.00e+00,0.00e+00,0.00e+00
,  4.48e-03,2.10e-06,0.00e+00)
)
## ppcor::pcor(y.data)

## partial correlation based on a formula
partialCor(c(hl,disp)~BC+deg, data = y.data)
partialCor(hl + disp~BC+deg, data = y.data)
## partial correlation based on a list
partialCor(list(hl~BC+deg,disp~BC+deg), data = y.data)
## via an existing model
e.lm <- lmm(hl~disp+BC+deg, data = y.data)
partialCor(e.lm)

## using a different set of covariates for outcome
partialCor(list(hl~BC+deg, disp~BC), data = y.data)

## statified correlation (using another dataset)
data(gastricbypassW, package = "LMMstar")
gastricbypassW$weight.bin <- gastricbypassW$weight1>=120
partialCor(glucagonAUC1+glucagonAUC2~1, data = gastricbypassW, by = "weight.bin")

## compared correlation between groups
partialCor(glucagonAUC1+glucagonAUC2~1, data = gastricbypassW, by = "weight.bin",
           effects = "Dunnett") 

#### with repetitions ####
## Not run: 
data(gastricbypassL, package = "LMMstar")
## via a mixed model
eUN.lmm <- lmm(weight ~ glucagonAUC+time, repetition =~time|id,
               data = gastricbypassL, structure = "UN")
partialCor(eUN.lmm)

## mean: variable and timepoint specific mean parameter (8)
## variance: variable and timepoint specific variance parameter (8)
## correlation: correlation parameter specific for each variable and time lag (10)
e.cor <- partialCor(weight+glucagonAUC~time, repetition =~time|id,
                    data = gastricbypassL, structure = "LAG")
e.cor
coef(attr(e.cor,"lmm"), effects = "correlation")
if(require(ggplot2)){
autoplot(e.cor)
}

## same except for the mean structure: variable specific mean parameter (2)
e.cor2 <- partialCor(weight+glucagonAUC~time, repetition =~time|id,
                    data = gastricbypassL, structure = "LAG")

## mean: variable and timepoint specific mean parameter (8)
## variance: variable specific variance parameter (2)
## correlation: correlation parameter specific for each variable and some time lag (4)
e.cor3 <- partialCor(weight+glucagonAUC~time, repetition =~time|id,
                     data = gastricbypassL, structure = "CS")
e.cor3
coef(attr(e.cor3,"lmm"), effects = "correlation")
if(require(ggplot2)){
autoplot(e.cor3)
}


## End(Not run)

Data From The Potassium Intake Study (Long Format with intermediate measurements)

Description

Data from the potassium intake study, a randomized placebo-controlled crossover study where the effect of potassium supplement (90 mmol/day) on the renin-angiostensin-aldosteron system (RAAS) was assessed. This dataset is in the long format (i.e. one line per measurement) and contains measurement over 6 timepoints for each time period.

• id: patient identifier.

• sequence: treatment group to which the patient has been randomized.

• period: time period.

• treatment: treatment during the time period.

• time: time within each period.

• aldo: ??

Usage

data(potassiumRepeatedL)

References

Dreier et al. Effect of increased potassium intake on the reninangiotensinaldosterone system and subcutaneous resistance arteries: a randomized crossover study, Nephrol Dial Transplant (2020) 110. doi: 10.1093/ndt/gfaa114

Data From The Potassium Intake Study (Long Format)

Description

Data from the potassium intake study, a randomized placebo-controlled crossover study where the effect of potassium supplement (90 mmol/day) on the renin-angiostensin-aldosteron system (RAAS) was assessed. This dataset is in the long format (i.e. one line per measurement).

• id: patient identifier.

• sequence: treatment group to which the patient has been randomized.

• period: time period.

• treatment: treatment during the time period.

• auc: area under the curve of ?? during the time period.

• bsauc: ??

• aldo: ??

Usage

data(potassiumSingleL)

References

Dreier et al. Effect of increased potassium intake on the reninangiotensinaldosterone system and subcutaneous resistance arteries: a randomized crossover study, Nephrol Dial Transplant (2020) 110. doi: 10.1093/ndt/gfaa114

Data From The Potassium Intake Study (Wide Format)

Description

Data from the potassium intake study, a randomized placebo-controlled crossover study where the effect of potassium supplement (90 mmol/day) on the renin-angiostensin-aldosteron system (RAAS) was assessed. This dataset is in the wide format (i.e. one line per patient).

• id: patient identifier.

• sequence: treatment group to which the patient has been randomized.

• treatment1: treatment during the first time period.

• treatment2: treatment during the second time period.

• auc1: area under the curve of ?? during the first time period.

• auc2: area under the curve of ?? during the second time period.

• bsauc1: ??

• aldo1: ??

• aldo2: ??

Usage

data(potassiumSingleW)

References

Dreier et al. Effect of increased potassium intake on the reninangiotensinaldosterone system and subcutaneous resistance arteries: a randomized crossover study, Nephrol Dial Transplant (2020) 110. doi: 10.1093/ndt/gfaa114

Predicted Mean Value With Uncertainty For Linear Mixed Model

Description

Predicted mean value conditional on covariates or on covariates and other outcome values.

Usage

## S3 method for class 'lmm'
predict(
  object,
  newdata,
  type = "static",
  p = NULL,
  se = NULL,
  robust = FALSE,
  df = NULL,
  level = 0.95,
  keep.data = NULL,
  format = "long",
  export.vcov = FALSE,
  simplify = TRUE,
  ...
)

Arguments

Details

Static prediction are made using the linear predictor X\beta while dynamic prediction uses the conditional normal distribution of the missing outcome given the observed outcomes. So if outcome 1 is observed but not 2, prediction for outcome 2 is obtain by X_2\beta + \sigma_{21}\sigma^{-1}_{22}(Y_1-X_1\beta). In that case, the uncertainty is computed as the sum of the conditional variance \sigma_{22}-\sigma_{21}\sigma^{-1}_{22}\sigma_{12} plus the uncertainty about the estimated conditional mean (obtained via delta method using numerical derivatives).

The model terms are computing similarly to stats::predict.lm, by centering the design matrix around the mean value of the covariates used to fit the model. Then the centered design matrix is multiplied by the mean coefficients and columns assigned to the same variable (e.g. three level factor variable) are summed together.

Value

When format="long", a data.frame containing the following columns:

• estimate: predicted mean.

• se: uncertainty about the predicted mean.

• df: degree of freedom

• lower: lower bound of the confidence interval of the predicted mean

• upper: upper bound of the confidence interval of the predicted mean

When format="wide", a matrix containing the predict means (one line per cluster, one column per timepoint).

Examples

## simulate data in the long format
set.seed(10)
dL <- sampleRem(100, n.times = 3, format = "long")

## fit Linear Mixed Model
eUN.lmm <- lmm(Y ~ visit + X1 + X2 + X5,
               repetition = ~visit|id, structure = "UN", data = dL)

## prediction
newd <- data.frame(X1 = 1, X2 = 2, X5 = 3, visit = factor(1:3, levels = 1:3))
predict(eUN.lmm, newdata = newd)
predict(eUN.lmm, newdata = newd, keep.data = TRUE)
predict(eUN.lmm, newdata = newd, keep.data = TRUE, se = c(TRUE,TRUE))

## dynamic prediction
newd.d1 <- cbind(newd, Y = c(NA,NA,NA))
predict(eUN.lmm, newdata = newd.d1, keep.data = TRUE, type = "dynamic")
newd.d2 <- cbind(newd, Y = c(6.61,NA,NA))
predict(eUN.lmm, newdata = newd.d2, keep.data = TRUE, type = "dynamic")
newd.d3 <- cbind(newd, Y = c(1,NA,NA))
predict(eUN.lmm, newdata = newd.d3, keep.data = TRUE, type = "dynamic")
newd.d4 <- cbind(newd, Y = c(1,1,NA))
predict(eUN.lmm, newdata = newd.d4, keep.data = TRUE, type = "dynamic")

Evaluate Contour of the Log-Likelihood

Description

Display the (restricted) log-likelihood around Maximum Likelihood Estimate (MLE) under specific constrains.

Usage

## S3 method for class 'lmm'
profile(
  fitted,
  effects = NULL,
  profile.likelihood = FALSE,
  maxpts = NULL,
  conf.level = 0.95,
  trace = FALSE,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  transform.names = TRUE,
  ...
)

Arguments

Details

Each parameter defined by the argument effets is treated separately:

• the confidence interval of a parameter is discretized with maxpt points,

• this parameter is set to a discretization value.

• the other parameters are either set to the (unconstrained) MLE (profile.likelihood=FALSE) or to constrained MLE (profile.likelihood=TRUE). The latter case is much more computer intensive as it implies re-running the estimation procedure.

• the (restricted) log-likelihood is evaluated.

Value

A data.frame object containing the log-likelihood for various parameter values.

Examples

data(gastricbypassW, package = "LMMstar")
e.lmm <- lmm(weight2 ~ weight1 + glucagonAUC1,
             data = gastricbypassW, control = list(optimizer = "FS"))

## profile logLiklihood
## Not run: 
e.pro <- profile(e.lmm, effects = "all", maxpts = 10, profile.likelihood = TRUE)
head(e.pro)
plot(e.pro)

## End(Not run)

## along a single parameter axis
e.sliceNone <- profile(e.lmm, effects = "all", maxpts = 10, transform.sigma = "none")
plot(e.sliceNone)
e.sliceLog <- profile(e.lmm, effects = "all", maxpts = 10, transform.sigma = "log")
plot(e.sliceLog)

Proportion of Significant Findings

Description

Evaluate the proportion of test above the statistical significance level

Usage

proportion(object, n.sample, trace, ...)

Arguments

Value

a data.frame with the estimated proportion (estimate column), standard error and confidence interval (when boostrap is used).

Estimate Random Effect From a Linear Mixed Model

Description

Recover the random effects from the variance-covariance parameter of a linear mixed model.

Usage

## S3 method for class 'lmm'
ranef(
  object,
  effects = "mean",
  scale = "absolute",
  se = FALSE,
  df = NULL,
  transform = (effects %in% c("std", "variance")),
  p = NULL,
  newdata = NULL,
  format = "long",
  simplify = TRUE,
  ...
)

Arguments

Details

Consider the following mixed model:

Y = X\beta + \epsilon = X\beta + Z\eta + \xi

where the variance of \epsilon is denoted \Omega, the variance of \eta is denoted \Omega_{\eta}, and the variance of \xi is \sigma^2 I with I is the identity matrix.
The random effets are estimating according to:

E[Y|\eta] = \Omega_{\eta} Z^{t} \Omega^{-1} (Y-X\beta)

Value

A data.frame or a list depending on the argument format.

Examples

if(require(nlme)){
data(gastricbypassL, package = "LMMstar")

## random intercept
e.RI <- lmm(weight ~ time + (1|id), data = gastricbypassL)
ranef(e.RI, effects = "mean")
ranef(e.RI, effects = "mean", se = TRUE)

ranef(e.RI, effects = "variance")
ranef(e.RI, effects = "variance", format = "wide")

}

Linear Hypothesis Testing Across Linear Mixed Models

Description

Linear hypothesis testing accross linear mixed model.

Usage

## S3 method for class 'Wald_lmm'
rbind(model, ..., effects = NULL, rhs = NULL, name = NULL, sep = ": ")

Arguments

Details

WARNING: in presence of measurements from the same cluster across several models, the influence function is used to estimate the covariance between the model parameters. This is why the (robust) standard errors may not match the (model-based) standard error from the linear mixed Setting the argument robust to FALSE when calling anova.lmm will rescale the (robust) standard errors to mimic the original model-based standard errors.

Examples

## simulate data
set.seed(10)
dL <- sampleRem(1e2, n.times = 3, format = "long")

## estimate mixed models
e.lmm1 <- lmm(Y ~ X1+X2+X3, repetition = ~visit|id, data = dL,
              structure = "CS", df = FALSE)
e.lmm2 <- lmm(Y ~ X1+X8+X9, repetition = ~visit|id, data = dL,
              structure = "CS", df = FALSE)

model.tables(e.lmm1) ## model-based standard errors
model.tables(e.lmm1, robust = TRUE) ## robust standard errors

## select null hypotheses & combine (robust standard errors)
AAA <- anova(e.lmm1, ci = TRUE, effect = c("X1|X2,X3"="X1=0","X2|X1,X3"="X2=0"))
BBB <- anova(e.lmm2, ci = TRUE, effect = c("X1|X8,X9"="X1=0"))

ZZZ <- rbind(AAA,BBB)

## select null hypotheses & combine (model-based like standard errors)
AA <- anova(e.lmm1, ci = TRUE, effect = c("X1|X2,X3"="X1=0","X2|X1,X3"="X2=0"),
             robust = FALSE)
BB <- anova(e.lmm2, ci = TRUE, effect = c("X1|X8,X9"="X1=0"),
             robust = FALSE)
ZZ <- rbind(AA,BB)

Remove Columns from Output

Description

Auxiliary function that can be used when specifying the argument columns (e.g. calling confint.lmm) to remove columns.

Usage

remove(...)

Arguments

Value

A character vector

Examples

set.seed(10)
dW <- sampleRem(25, n.times = 1, format = "long")
e.lmm <- lmm(Y~X1, data = dW)

confint(e.lmm, columns = remove("estimate"))

Inference via Resampling for Linear Mixed Model

Description

Non-parametric bootstrap or permutation test for Linear Mixed Models.

Usage

resample(object, type, ...)

## S3 method for class 'lmm'
resample(
  object,
  type,
  effects,
  n.sample = 1000,
  studentized = TRUE,
  level = 0.95,
  correction = TRUE,
  trace = TRUE,
  seed = NULL,
  cpus = 1,
  export.cpus = NULL,
  ...
)

Arguments

Details

All approach are carried at the cluster level:

• Bootstrap: sampling with replacement clusters. If a cluster is picked twice then different cluster id is used for each pick.

• Permutation: permuting covariate values between clusters (this only lead to the null hypothesis when the covariate values are constant within clusters) or assigning new outcome values by, under the null, permuting the normalized residuals, rescaling to residuals, and adding back the permuted fixed effect (any mean effect under H1 would be 0 because of the permutation if the variance-covariance structure is correct). The later procedure originates from Oliver et al (2012).

The studentized bootstrap keep the original estimate and standard error but uses the samples to evaluates the quantiles of the distribution used to form the confidence intervals. The studentized permutation test approximate the distribution of the test statistic under the null (instead of the distribution of the estimate under the null.).

P-values for the bootstrap are computed by centering the bootstrap distribution of the estimate or test statistic around 0 and evaluating the frequency at which it takes values more extreme than the observed estimate or test statistics.

References

Oliver E. Lee and Thomas M. Braun (2012), Permutation Tests for Random Effects in Linear Mixed Models. Biometrics, Journal 68(2).

Examples

## Not run: 

#### univariate linear regression ####
data(gastricbypassW, package = "LMMstar")
## rescale to ease optimization
gastricbypassW$weight1 <- scale(gastricbypassW$weight1)
gastricbypassW$weight2 <- scale(gastricbypassW$weight2)
gastricbypassW$glucagonAUC1 <- scale(gastricbypassW$glucagonAUC1)

e.lm <- lmm(weight2~weight1+glucagonAUC1, data = gastricbypassW)
model.tables(e.lm)

## non-parametric bootstrap
resample(e.lm, type = "boot", effects = c("weight1","glucagonAUC1"), seed = 10)
## permutation test
resample(e.lm, type = "perm-var", effects = "weight1", seed = 10) 
resample(e.lm, type = "perm-var", effects = "glucagonAUC1", seed = 10)
## using multiple cores
resample(e.lm, type = "boot", effects = c("weight1","glucagonAUC1"), cpus = 4)

#### random intercept model ####
data(gastricbypassL, package = "LMMstar")
gastricbypassL$weight <- scale(gastricbypassL$weight)
gastricbypassL$glucagonAUC <- scale(gastricbypassL$glucagonAUC)
gastricbypassL$gender <- as.numeric(gastricbypassL$id) %% 2
gastricbypassLR <- na.omit(gastricbypassL)

eCS.lmm <- lmm(weight~glucagonAUC+gender, data = gastricbypassLR,
               repetition = ~visit|id, structure = "CS")
model.tables(eCS.lmm)

## non-parametric bootstrap
resample(eCS.lmm, type = "boot", effects = c("glucagonAUC","gender"), seed = 10, trace = FALSE)
## permutation test
resample(eCS.lmm, type = "perm-var", effects = "gender", seed = 10)
resample(eCS.lmm, type = "perm-res", effects = "glucagonAUC", seed = 10) 

## End(Not run)

Extract The Residuals From a Linear Mixed Model

Description

Extract or compute the residuals of a linear mixed model.

Usage

## S3 method for class 'lmm'
residuals(
  object,
  type = "response",
  variable = NULL,
  at = NULL,
  newdata = NULL,
  p = NULL,
  format = "long",
  keep.data = FALSE,
  fitted.ci = FALSE,
  simplify = TRUE,
  var,
  ...
)

## S3 method for class 'clmm'
residuals(object, ...)

## S3 method for class 'mlmm'
residuals(object, simplify = TRUE, ...)

Arguments

Details

The argument type defines how the residuals are computed:

• "response": raw residual, i.e. observed outcome minus fitted value \varepsilon_{ij} = Y_{ij} - X_{ij} \hat{\beta}.

• "pearson": each raw residual is divided by its modeled standard deviation \varepsilon_{ij} = \frac{Y_{ij} - X_{ij} \hat{\beta}}{\sqrt{\hat{\omega}_{ij}}}.

• "studentized": same as "pearson" but excluding the contribution of the cluster in the modeled standard deviation \varepsilon_{ij} = \frac{Y_{ij} - X_{ij} \hat{\beta}}{\sqrt{\hat{\omega}_{ij}-\hat{q}_{ij}}}.

• "normalized": raw residuals are multiplied, within clusters, by the inverse of the (upper) Cholesky factor of the modeled residual variance covariance matrix \varepsilon_{ij} = ( Y_{i} - X_{i} \hat{\beta} )\hat{C}^{-1}.

• "normalized2": raw residuals are multiplied, within clusters, by the inverse of the modeled residual variance covariance matrix \varepsilon_{ij} = ( Y_{i} - X_{i} \hat{\beta} )\hat{Omega}^{-1}.

• "scaled": scaled residuals (see PROC MIXED in SAS). Numerically identical to "normalized" but computed by sequentially scaling and centering the residuals, to make them conditionally independent of previous residuals from the same cluster at previous repetitions.

• "partial": partial residuals (\gamma E + \hat{\varepsilon}). A reference level can be also be specified via the attribute "reference" to change the absolute level of the partial residuals. "partial-center": partial residuals with centered continuous covariates (\gamma E + \hat{\varepsilon} where E has been centered, i.e., has 0-mean)

where

• X=(E,W) the design matrix. For partial residuals, it is split according to the variable(s) in argument variable (E) and the rest (W).

• Y the outcome

• \hat{\beta}=(\hat{\gamma},\hat{\delta}) the estimated mean coefficients relative to X=(E,W)

• \hat{\Omega} the modeled variance-covariance of the residuals and \hat{\omega} its diagonal elements

• \hat{C} the upper Cholesky factor of \hat{\Omega}, i.e. upper triangular matrix satisfying \hat{C}^{t} \hat{C} = \hat{\Omega}

• \hat{Q}_i= X_i (X^{t}\hat{\Omega}X)^{-1}X_i^{t} a cluster specific correction factor, approximating the contribution of cluster i to \hat{\Omega}. Its diagonal elements are denoted \hat{q}_i.

• \hat{D}_i the upper Cholesky factor of \hat{\Omega}-\hat{Q}_i

Setting argument fitted.ci to TRUE, simplify to FALSE, format to "long" returns an attribute "grad" containing the first order partial derivatives of the residuals with respect to the model parameters.

Value

lmm: a vector or a data.frame when format="long" (one line per observation, one column per type of residual), a matrix when format="wide" (one line per cluster, one column per timepoint).

Examples

#### simulate data in the long format ####
set.seed(10)
dL <- sampleRem(100, n.times = 3, format = "long")

#### Linear Model ####
e.lm <- lmm(Y ~ visit + X1 + X2 + X6, data = dL)

## partial residuals
pRes <- residuals(e.lm, type = "partial", variable = "X6")
range(residuals(e.lm) + dL$X6 * coef(e.lm)["X6"] - pRes)
e.reslm <- residuals(e.lm, type = "partial", variable = "X6", keep.data = TRUE, simplify = FALSE)
plot(e.reslm)

## partial residuals with specific reference
residuals(e.lm, type = "partial", variable = "X1",
          at = data.frame(visit=factor(2,1:3),X2=0,X6=3))

## partial residuals with centered covariates
residuals(e.lm, type = "partial-center", variable = "X1")

#### Linear Mixed Model ####
eUN.lmm <- lmm(Y ~ visit + X1 + X2 + X5 + X6,
               repetition = ~visit|id, structure = "UN", data = dL)

## residuals
e.resL <- residuals(eUN.lmm, type = "normalized",
                    keep.data = TRUE, simplify = FALSE)
plot(e.resL, type = "qqplot")
plot(e.resL, type = "scatterplot", labeller = ggplot2::label_both)
e.resW <- residuals(eUN.lmm, format = "wide", type = "normalized",
                    simplify = FALSE)
plot(e.resW, type = "correlation")

## residuals and predicted values
residuals(eUN.lmm, type = "all")
residuals(eUN.lmm, type = "all", keep.data = TRUE)

## partial residuals
residuals(eUN.lmm, type = "partial", variable = c("(Intercept)","visit","X6"))

Restaure NA

Description

Restaure NA in the user output that have been excluded when fitting the LMM.

Usage

restaureNA(object, index.na, level, cluster)

Arguments

Sample Longuitudinal Data

Description

Sample longuitudinal data with covariates

Usage

sampleRem(
  n,
  n.times,
  mu = 1:n.times,
  sigma = rep(1, n.times),
  lambda = rep(1, n.times),
  beta = c(2, 1, 0, 0, 0, 1, 1, 0, 0, 0),
  gamma = matrix(0, nrow = n.times, ncol = 10),
  format = "wide",
  latent = FALSE
)

Arguments

Details

The generative model is a latent variable model where each outcome (Y_j) load on the latent variable (\eta) with a coefficient lambda:

Y_j = \mu_j + \lambda_j*\eta + \sigma_j\epsilon_j

The latent variable is related to the covariates (X_1,\ldots,X_(10)):

\eta = \alpha + \beta_1 X_1 + ... + \beta_{10} X_{10} + \xi

\epsilon_j and \xi are independent random variables with standard normal distribution.

Value

a data.frame

Examples

#### generate data in the wide format ####
set.seed(10)
dW <- sampleRem(100, n.times = 3, format = "wide")
head(dW)

#### generate data in the long format ####
set.seed(10)
dL <- sampleRem(100, n.times = 3, format = "long")
head(dL)

Scatterplot for Continuous Variables

Description

Produce a matrix of plot for continuous variables: scatterplots, histograms, correlation and missing values. Inspired from the ggpairs function of the R package GGally.

Usage

scatterplot(
  data,
  formula,
  columns,
  format = NULL,
  group = NULL,
  transform = NULL,
  facet = "grid",
  alpha.point = 1,
  type.diag = "boxplot",
  bins = NULL,
  position.bar = "identity",
  linewidth.density = NULL,
  alpha.area = NULL,
  method.cor = "pearson",
  name.cor = "r",
  size.cor = NULL,
  digits = c(3, 2),
  display.NA = NULL,
  color = NULL,
  xlim = NULL,
  ylim = NULL,
  size.axis = NULL,
  size.legend = NULL,
  size.facet = NULL
)

Arguments

Details

In the long format, the outcome variable contains the numerical values to be displayed. The time variable will be used to spit outcome and display each split separately or jointly with one other split. The identifier links the outcome values across time.

Value

a list of ggplot objects (facet="grid") or a ggplot object (facet="grid2")

Examples

data(gastricbypassL, package = "LMMstar")
gastricbypassL$group <- as.numeric(gastricbypassL$id) %% 3
data(gastricbypassW, package = "LMMstar")

## single group (wide or long format)
scatterplot(gastricbypassL, formula = weight~time|id)
scatterplot(gastricbypassW, columns = paste0("weight",1:4))

## Not run: 
## use histogram instead of boxplot
scatterplot(gastricbypassL, formula = weight~time|id, type.diag = "hist")
scatterplot(gastricbypassL, formula = weight~time|id, type.diag = "hist", bins = 15)

## same scale
scatterplot(gastricbypassL, formula = weight~time|id,
            xlim = "common", ylim = "common")

## transform outcome
scatterplot(gastricbypassL, formula = weight~time|id, transform = "log")

## handling missing values
scatterplot(gastricbypassL, formula = glucagonAUC~time|id)

## coloring per group
scatterplot(gastricbypassL, formula = weight~time|id, group = "group")

## only display NAs
scatterplot(gastricbypassL, formula = glucagonAUC~time|id,
            display.NA = "only", group = "group")
scatterplot(gastricbypassL, formula = glucagonAUC~time|id,
            display.NA = "only", group = "group", size.legend = c(15,2))

## End(Not run)

Simulated Data with 3-level struture (Long Format)

Description

Simulated data a nested structure: Student/Class/School and one outcome.

• school:

• class:

• student:

• outcome:

Usage

data(schoolL)

Extract The Score From a Linear Mixed Model

Description

Extract or compute the first derivative of the log-likelihood of a linear mixed model.

Usage

## S3 method for class 'lmm'
score(
  x,
  effects = "mean",
  newdata = NULL,
  p = NULL,
  indiv = FALSE,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  transform.names = TRUE,
  ...
)

Arguments

Details

For details about the arguments transform.sigma, transform.k, transform.rho, see the documentation of the coef.lmm function.

Value

When argument indiv is FALSE, a vector with the value of the score relative to each coefficient. When argument indiv is TRUE, a matrix with the value of the score relative to each coefficient (in columns) and each cluster (in rows).

Extract The Residuals Variance-Covariance Matrix From a Linear Mixed Model

Description

Extract the unique set of residuals variance-covariance matrices or the one relative to specific clusters.

Usage

## S3 method for class 'lmm'
sigma(
  object,
  cluster = NULL,
  p = NULL,
  chol = FALSE,
  inverse = FALSE,
  simplify = TRUE,
  ...
)

Arguments

Value

A list where each element contains a residual variance-covariance matrix. Can also be directly a matrix when argument is simplify=TRUE and there is a single residual variance-covariance matrix.

Examples

## simulate data in the long format
set.seed(10)
dL <- sampleRem(100, n.times = 3, format = "long")
dL$id.fac <- paste0("id",dL$id)

## fit Linear Mixed Model
eUN.lmm <- lmm(Y ~ X1 + X2 + X5, repetition = ~visit|id.fac,
               structure = "UN", data = dL, df = FALSE)

## extract residuals variance covariance matrix
sigma(eUN.lmm) ## unique patterns
sigma(eUN.lmm, cluster = c("id1","id5")) ## existing clusters
sigma(eUN.lmm, cluster = dL[1:7,,drop=FALSE]) ## new clusters

Compute summary statistics

Description

Compute summary statistics for multiple variables and/or multiple groups and save them in a data frame.

Usage

summarize(
  formula,
  data,
  na.action = stats::na.pass,
  na.rm = FALSE,
  level = 0.95,
  columns = c("observed", "missing", "pc.missing", "mean", "sd", "min", "q1", "median",
    "q3", "max", "correlation"),
  FUN = NULL,
  skip.reference = TRUE,
  digits = NULL,
  filter = NULL,
  ...
)

Arguments

Details

This function is essentially an interface to the stats::aggregate function.
WARNING: it has the same name as a function from the dplyr package. If you have loaded dplyr already, you should use ::: to call summarize i.e. use LMMstar:::summarize.

Confidence intervals (CI) and prediction intervals (PI) for the mean are computed via stats::t.test. Confidence intervals (CI) for the median are computed via asht::medianTest.

Correlation can be assessed when a grouping and ordering variable are given in the formula interface , e.g. Y ~ time|id.

Value

A data frame containing summary statistics (in columns) for each outcome and value of the grouping variables (rows). It has an attribute "correlation" when it was possible to compute the correlation matrix for each outcome with respect to the grouping variable.

Examples

#### simulate data (wide format) ####
set.seed(10)
d <- sampleRem(1e2, n.times = 3)
d$treat <-  sample(LETTERS[1:3], NROW(d), replace=TRUE, prob=c(0.3, 0.3, 0.4) )

## add a missing value
d2 <- d
d2[1,"Y2"] <- NA

#### summarize (wide format) ####
summarize(Y1 ~ 1, data = d)
summarize(Y1 ~ 1, data = d, FUN = quantile, p = c(0.25,0.75))
summarize(Y1+Y2 ~ X1, data = d)
summarize(treat ~ 1, data = d)
summarize(treat ~ 1, skip.reference = FALSE, data = d)

summarize(Y1 ~ X1, data = d2)
summarize(Y1+Y2 ~ X1, data = d2, na.rm = TRUE)
summarize(. ~ treat, data = d2, na.rm = TRUE, filter = "Y")

#### summarize (long format) ####
dL <- reshape(d, idvar = "id", direction = "long",
             v.names = "Y", varying = c("Y1","Y2","Y3"))
summarize(Y ~ time + X1, data = dL)

## compute correlations (single time variable)
e.S <- summarize(Y ~ time + X1 | id, data = dL, na.rm = TRUE)
e.S
attr(e.S, "correlation")

## compute correlations (composite time variable)
dL$time2 <- dL$time == 2
dL$time3 <- dL$time == 3
e.S <- summarize(Y ~ time2 + time3 + X1 | id, data = dL, na.rm = TRUE)
e.S
attr(e.S, "correlation")

Summarize missing data patterns

Description

Summarize missing data patterns.

Usage

summarizeNA(
  data,
  repetition = NULL,
  sep = "",
  newnames = c("variable", "frequency", "missing.pattern", "n.missing"),
  keep.data = TRUE
)

Arguments

Value

a data frame

See Also

autoplot.summarizeNA for a graphical display.

Examples

#### display missing data pattern (wide format) ####
data(gastricbypassW, package = "LMMstar")
e.SNA <- summarizeNA(gastricbypassW) 
plot(e.SNA)
summarizeNA(gastricbypassW, keep.data = FALSE)

#### display missing data pattern (long format) ####
## example 1
data(gastricbypassL, package = "LMMstar")
e.SNAL <- summarizeNA(gastricbypassL, repetition = ~time|id)
plot(e.SNAL, variable = "glucagonAUC")

## example 2
data(calciumL, package = "LMMstar")
mp <- summarizeNA(calciumL, repetition = ~visit|girl)
plot(mp, variable = "bmd")
summarizeNA(calciumL[,c("visit","girl","bmd")], repetition = ~visit|girl)

## example 3
data(vasscoresW, package = "LMMstar")
summarizeNA(vasscoresW)

Summary of Testing for a Linear Mixed Models

Description

Estimates, p-values, and confidence intevals for linear hypothesis testing, possibly adjusted for multiple comparisons.

Usage

## S3 method for class 'Wald_lmm'
summary(
  object,
  print = TRUE,
  seed = NULL,
  columns = NULL,
  legend = TRUE,
  digits = 3,
  digits.df = 1,
  digits.p.value = 3,
  sep = ": ",
  ...
)

Arguments

Details

By default adjustment for multiple comparisons via a single step max-test adjustment, either using the multcomp package (equal degrees of freedom, method="single-step") or the copula package (unequal degrees of freedom, method="single-step2"). See the argument method of confint.Wald_lmm for other adjustments for multiple comparisons.
When multiple multivariate Wald tests are performed, adjustment for multiple comparisons for the univariate Wald tests is performed within each multivariate Wald test. The number of tests ajusted for equal the first degree of freedom of the multivariate Wald statistic.

Adding the value "type" in argument "columns" ensures that the type of parameter that is being test (mean, variance, correlation) is output.

Value

NULL

Summary Output for a Linear Mixed Model

Description

Summary output for a linear mixed model fitted with lmm.

Usage

## S3 method for class 'lmm'
summary(
  object,
  level = 0.95,
  robust = FALSE,
  print = TRUE,
  columns = NULL,
  digits = 3,
  digits.df = 1,
  digits.p.value = 3,
  hide.data = FALSE,
  hide.fit = FALSE,
  hide.cor = NULL,
  type.cor = NULL,
  hide.var = NULL,
  hide.sd = NULL,
  hide.re = NULL,
  hide.mean = FALSE,
  ...
)

Arguments

Value

A list containing elements displayed in the summary:

• correlation: the correlation structure.

• variance: the variance structure.

• sd: the variance structure expressed in term of standard deviations.

• mean: the mean structure.

Summary of Multiple Linear Mixed Models

Description

Estimates, p-values, and confidence intevals for multiple linear mixed models.

Usage

## S3 method for class 'mlmm'
summary(
  object,
  digits = 3,
  method = NULL,
  print = NULL,
  hide.data = FALSE,
  hide.fit = FALSE,
  ...
)

Arguments

Summary for partial correlation

Description

Display estimated partial correlation and associated p-values and confidence intevals.

Usage

## S3 method for class 'partialCor'
summary(object, digits = 3, detail = TRUE, ...)

Arguments

Data From The SWABS Study (Long Format)

Description

Data from the swabs study, where the pneumococcus was studied in 18 families with different space available for the household. This dataset is in the long format (i.e. one line per measurement).

• crowding: space available in the household.

• family: family serial number

• name: type of family member.

• swabs: number of times the swab measurement was positive.

Usage

data(swabsL)

References

TODO

Data From The SWABS Study (Wide Format)

Description

Data from the swabs study, where the pneumococcus was studied in 18 families with different space available for the household. This dataset is in the wide format (i.e. one line per patient).

• crowding: space available in the household.

• family: family serial number

• mother: number of times the swab measurement was positive for the mother.

• father: number of times the swab measurement was positive for the father.

• child1: number of times the swab measurement was positive for the first child.

• child2: number of times the swab measurement was positive for the second child.

• child3: number of times the swab measurement was positive for the third child.

Usage

data(swabsW)

References

Grundy SM, Lan SP, Lachin J. The effects of chenodiol on biliary lipids and their association with gallstone dissolution in the National Cooperative Gallstone Study (SWABS). J Clin Invest. 1984 Apr;73(4):1156-66. doi: 10.1172/JCI111301.

Model Terms For Linear Mixed Models

Description

Model terms for linear mixed models. Used by multcomp::glht.

Usage

## S3 method for class 'lmm'
terms(x, ...)

Arguments

Value

An object of class terms giving a symbolic representation of the mean structure.

Data From The VAS Study (Long Format)

Description

Data from the VAS Study, a randomized controlled clinial trial assessing the healing effect of topical zink sulfate on epidermal wound. The study includes 30 heatlhy volunteers with induced wounds on each buttock which where subsequently treated with a different treatment for each wound. Then the VAS-score (pain sensation on a 0-100mm visual analogue scale) was assessed after each treatment application and summarized by area under the curve. This dataset is in the long format (i.e. one line per measurement).

• id: patient identifier.

• group: treatment group to which the patient has been randomized.

• treat.num:

• vas: VAS-score relative to the wound.

• treatment: Treatment used on the wound. A: active treatment (zink shower gel), B: placebo treatment (shower gel without zink), C: control treatment (demineralized water).

Usage

data(vasscoresL)

References

TODO

Data From The VAS Study (Wide Format)

Description

Data from the VAS Study, a randomized controlled clinial trial assessing the healing effect of topical zink sulfate on epidermal wound. The study includes 30 heatlhy volunteers with induced wounds on each buttock which where subsequently treated with a different treatment for each wound. Then the VAS-score (pain sensation on a 0-100mm visual analogue scale) was assessed after each treatment application and summarized by area under the curve. This dataset is in the wide format (i.e. one line per patient).

• id: patient identifier.

• group: treatment group to which the patient has been randomized.

• vasA: VAS-score when using a zink shower gel.

• vasB: VAS-score when using a placebo treatment (shower gel without zink).

• vasC: VAS-score when using a control treatment with demineralized water.

Usage

data(vasscoresW)

References

TODO

Extract The Variance-Covariance Matrix From a Linear Mixed Model

Description

Extract the variance-covariance matrix of the model coefficients of a linear mixed model.

Usage

## S3 method for class 'lmm'
vcov(
  object,
  effects = "mean",
  robust = FALSE,
  df = FALSE,
  strata = NULL,
  newdata = NULL,
  p = NULL,
  type.information = NULL,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  transform.names = TRUE,
  ...
)

Arguments

Details

For details about the arguments transform.sigma, transform.k, transform.rho, see the documentation of the coef.lmm function.

Value

A matrix with an attribute "df" when argument df is set to TRUE.

Data From The Vitamin Study (Long Format)

Description

Data from the vitamin Study, a randomized study where the growth of guinea pigs was monitored before and after intake of vitamin E/placebo. The weight of each guinea pig was recorded at the end of week 1, 3, 4, 5, 6, and 7. Vitamin E/placebo is given at the beginning of week 5. This dataset is in the long format (i.e. one line per measurement).

• group: treatment group: vitamin or placebo.

• animal: identifier

• weight1: weight (in g) of the pig at the end of week 1 (before treatment).

• weight3: weight (in g) of the pig at the end of week 3 (before treatment).

• weight4: weight (in g) of the pig at the end of week 4 (before treatment).

• weight5: weight (in g) of the pig at the end of week 5 (after treatment).

• weight6: weight (in g) of the pig at the end of week 6 (after treatment).

• weight7: weight (in g) of the pig at the end of week 7 (after treatment).

Usage

data(vitaminL)

References

Crowder and Hand (1990, p. 27) Analysis of Repeated Measures.

Data From The Vitamin Study (Wide Format)

Description

Data from the vitamin Study, a randomized study where the growth of guinea pigs was monitored before and after intake of vitamin E/placebo. The weight of each guinea pig was recorded at the end of week 1, 3, 4, 5, 6, and 7. Vitamin E/placebo is given at the beginning of week 5. This dataset is in the wide format (i.e. one line per patient).

• group: treatment group: vitamin or placebo.

• animal: identifier

• weight1: weight (in g) of the pig at the end of week 1 (before treatment).

• weight3: weight (in g) of the pig at the end of week 3 (before treatment).

• weight4: weight (in g) of the pig at the end of week 4 (before treatment).

• weight5: weight (in g) of the pig at the end of week 5 (after treatment).

• weight6: weight (in g) of the pig at the end of week 6 (after treatment).

• weight7: weight (in g) of the pig at the end of week 7 (after treatment).

Usage

data(vitaminW)

References

TODO

Extract Weights Used to Pool Estimates

Description

Extract weights used to pool estimates.

Usage

## S3 method for class 'Wald_lmm'
weights(object, method, ...)

Arguments

Value

a numeric vector whose elements sum to 1.

Examples

set.seed(10)
dL <- sampleRem(250, n.times = 3, format = "long")

e.mlmm <- mlmm(Y~X1+X2+X6, repetition = ~visit|id, data = dL,
               by = "X4", effects = "X1=0", structure = "CS")
weights(e.mlmm, method = "average")
weights(e.mlmm, method = "pool.fixse")
weights(e.mlmm, method = "pool.se")
weights(e.mlmm, method = "pool.gls")