Package 'LMMstar'

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.

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.

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"))

Description

Test linear combinations of parameters from a linear mixed model using Wald test or Likelihood Ratio Test (LRT).

Usage

## S3 method for class 'lmm'
anova(
  object,
  effects = NULL,
  rhs = NULL,
  robust = NULL,
  df = NULL,
  univariate = TRUE,
  multivariate = TRUE,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  simplify = NULL,
  ...
)

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.

• 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)

Description

Evaluate the Area Under the Curve (AUC) based on discrete observations y = f(x), using interpolation of order 0 (step), 1 (trapezoidal), or 3 (natural cubic splines). The AUC is the integral of a function over an interval [from,to].

Usage

approxAUC(
  x,
  y,
  from,
  to,
  method = "trapezoid",
  subdivisions = 100,
  name = "AUC",
  na.rm = FALSE
)

Arguments

Details

This function is a simplified version of the AUC function from the DescTools package.

Value

a numeric value.

Examples

## same examples as DescTools::AUC
approxAUC(x=c(1,3), y=c(1,1), from = 1, to = 3)

approxAUC(x=1:3, y=c(1,2,4), from = 1, to = 3)
approxAUC(x=1:3, y=c(1,2,4), from = 1, to = 3, method = "step")

x <- seq(0, pi, length.out=200)
approxAUC(x=x, y=sin(x), from = 0, to = pi)
approxAUC(x=x, y=sin(x), from = 0, to = pi, method = "spline")

Description

Graphical representation for correlation matrix

Usage

## S3 method for class 'correlate'
autoplot(
  object,
  index,
  size.text = 16,
  name.legend = "Correlation",
  title = NULL,
  scale = ggplot2::scale_fill_gradient2,
  type.cor = "both",
  low = "blue",
  high = "red",
  mid = "white",
  midpoint = mean(limits),
  limits = c(-1, 1),
  ...
)

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

Arguments

Functions

• plot(correlate): Graphical Display For Correlation Matrix

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 a Linear Mixed Model

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)
}

Description

Display estimated linear contrasts applied on parameters from subgroup-specific linear mixed models.

Usage

## S3 method for class 'mlmm'
autoplot(object, type = "forest", ...)

## S3 method for class 'mlmm'
plot(x, facet_nrow = NULL, facet_ncol = NULL, labeller = "label_value", ...)

Functions

• plot(mlmm): Graphical Display For Multiple Linear Mixed Models

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)
}

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

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,
  scales = "fixed",
  labeller = "label_value",
  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,
  ylim = NULL,
  position = NULL,
  ...
)

## 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

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)

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,
  decreasing.order = FALSE,
  title = NULL,
  labeller = "label_value",
  ...
)

## 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

Description

Display estimated linear contrast applied on parameters from a linear mixed model.

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 Wald Tests Applied to a Linear Mixed Model

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))
}

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)

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.

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.

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

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

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.

Description

TODO

• id: patient identifier.

• allocation:

• sex:

• age:

• visit:

• time:

• pwv:

• aix:

• dropout:

Usage

data(ckdL)

References

TO ADD

Description

TODO

• id: patient identifier.

• allocation:

• sex:

• age:

• pwv0:

• pwv12:

• pwv24:

• aix0:

• aix12:

• aix24:

• dropout:

Usage

data(ckdW)

References

TO ADD

Description

Extract estimated parameters from a linear mixed model.

Usage

## S3 method for class 'lmm'
coef(
  object,
  effects = NULL,
  p = NULL,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = "none",
  transform.names = TRUE,
  simplify = 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.

Value

A vector with the value of the model coefficients.
When using simplify=FALSE the character strings in attribute "type" refer to:

• "mu": mean parameters.

• "sigma": standard deviation parameters,

• "k": ratio between standard deviation,

• "rho": correlation parameter

The character strings in attribute "sigma" refer, for each parameter, to a possible corresponding standard deviation parameter. Those in attribute "k.x" and "k.y" refer to the ratio parameter. NA indicates no corresponding standard deviation or ratio parameter.

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 = "variance", transform.k = "sd")
coef(eUN.lmm, effects = "all", simplify = FALSE)

Description

Combine estimated parameters or linear contrasts applied on parameters from group-specific linear mixed models.

Usage

## S3 method for class 'mlmm'
coef(
  object,
  effects = "Wald",
  method = "none",
  p = NULL,
  ordering = "model",
  backtransform = object$args$backtransform,
  transform.sigma = "none",
  transform.k = "none",
  transform.rho = "none",
  transform.names = TRUE,
  simplify = TRUE,
  ...
)

Arguments

Description

Combine estimated values across linear contrasts applied on parameters from different linear mixed models.

Usage

## S3 method for class 'rbindWald_lmm'
coef(
  object,
  effects = "Wald",
  method = "none",
  ordering = NULL,
  transform.names = TRUE,
  backtransform = NULL,
  simplify = TRUE,
  ...
)

Arguments

Details

Argument effects: when evaluating the proportion of rejected hypotheses (effects="p.rejection") a "single-step" method will be used by default to evaluate the critical quantile. This can be changed by adding adjustment method, e.g. effects=c("bonferronin","p.rejection", in the argument.

Description

Extract estimated value of linear contrasts applied on parameters from a linear mixed model.

Usage

## S3 method for class 'Wald_lmm'
coef(
  object,
  effects = "Wald",
  backtransform = NULL,
  transform.names = TRUE,
  simplify = TRUE,
  ...
)

Arguments

Description

Compute confidence intervals (CIs) relative to parameters from 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 degrees-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)

Description

Compute pointwise or simultaneous confidence intervals relative to parameter or linear contrasts of parameters from group-specific linear mixed models. Can also output p-values (corresponding to pointwise confidence intervals) or adjusted p-values (corresponding to simultaneous confidence intervals).

Usage

## S3 method for class 'mlmm'
confint(
  object,
  parm = NULL,
  level = 0.95,
  method = NULL,
  df = NULL,
  columns = NULL,
  backtransform = 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). 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

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

Description

Combine pointwise or simultaneous confidence intervals relative to linear contrasts of parameters from different linear mixed models. Can also output p-values (corresponding to pointwise confidence intervals) or adjusted p-values (corresponding to simultaneous confidence intervals).

Usage

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

Arguments

Details

Argument method: the following pooling are available:

• "average": average estimates

• "pool.se": weighted average of the estimates, with weights being the inverse of the squared standard error.

• "pool.gls": weighted average of the estimates, with weights being based on the variance-covariance matrix of the estimates. When this matrix is singular, the Moore–Penrose inverse is used which correspond to truncate the spectral decomposition for eigenvalues below $10^{-12}$.

• "pool.gls1": similar to "pool.gls" with weights shrinked toward the average whenever they exceed 1 in absolute value.

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

• "p.rejection": proportion of null hypotheses where there is evidence for an effect. By default the critical quantile (defining the level of evidence required) is evaluated using a "single-step" method but this can be changed by adding adjustment method in the argument method, e.g. effects=c("bonferronin","p.rejection").

Description

Compute confidence intervals for linear mixed model using resampling (permutation or bootstrap).

Usage

## S3 method for class 'resample'
confint(
  object,
  parm = NULL,
  null = NULL,
  level = 0.95,
  method = NULL,
  columns = NULL,
  correction = TRUE,
  ...
)

Arguments

Details

Argument correction: if we denote by n.sample the number of permutations that have been performed, having the correction avoids p-values of 0 by ensuring they are at least as.numeric(correction)/(n.sample+as.numeric(correction)), e.g. at least 0.000999001 for a thousand permutations.

Argument correction (bootstrap): percentile confidence intervals are computed based on the quantile of the resampling distribution.
Gaussian confidence intervals and p-values are computed by assuming a normal distribution and estimating its variance based on the bootstrap distribution.
Studentized confidence intervals are computed using quantiles based on the boostrap distribution after it has been centered and rescaled by the point estimate and its standard error Percentile and Studentized p-values are computed by finding the coverage level such that one of the bound of the confidence interval equals the null.

Argument correction (bootstrap): Percentile p-values are computed based on the proportion of times the estimate was more extreme than the permutation estimates.
Gaussian p-values are computed by assuming a normal distribution and estimating its variance based on the permutation distribution.
Studentized p-values are computed based on the proportion of times the test statistic was more extreme than the permutation test statistic.
No confidence intervals are provided.

Description

Compute pointwise or simultaneous confidence intervals relative to linear contrasts of parameters from a linear mixed model. Pointwise confidence intervals have nominal coverage w.r.t. a single contrast whereas simultaneous confidence intervals have nominal coverage w.r.t. to all contrasts. Can also output p-values (corresponding to pointwise confidence intervals) or adjusted p-values (corresponding to simultaneous confidence intervals).

Usage

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

Arguments

Details

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.

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.

Description

Compute correlation matrix for multiple variables and/or multiple groups.

Usage

correlate(
  formula,
  data,
  repetition,
  use = "everything",
  method = "pearson",
  collapse.value = ".",
  collapse.var = ".",
  na.rm
)

Arguments

See Also

as.array to convert the output to an array or as.matrix to convert the output to a matrix (when a single outcome and no grouping variable).
summarize for other summary statistics (e.g. mean, standard deviation, ...).

Examples

#### simulate data (wide format) ####
data(gastricbypassL, package = "LMMstar")

## compute correlations (composite time variable)
e.S <- correlate(weight ~ 1, data = gastricbypassL, repetition = ~time|id)
e.S
as.matrix(e.S)

e.S21 <- correlate(glucagonAUC + weight ~ 1, data = gastricbypassL,
                  repetition = ~time|id, use = "pairwise")
e.S21
as.array(e.S21)
as.matrix(e.S21, index = "weight")

gastricbypassL$sex <- as.numeric(gastricbypassL$id) %% 2
e.S12 <- correlate(weight ~ sex, data = gastricbypassL,
                   repetition = ~time|id, use = "pairwise")
e.S12
as.array(e.S12)

e.S22 <- correlate(glucagonAUC + weight ~ sex, data = gastricbypassL,
                   repetition = ~time|id, use = "pairwise")
e.S22
as.array(e.S22)

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")

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))

Description

Estimate the residual degrees-of-freedom from a linear mixed model.

Usage

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

Arguments

Details

The residual degrees-of-freedom is estimated using the sum of squared normalized residuals.

Value

A numeric value

Description

Combine the residuals degrees-of-freedom from group-specific linear mixed models.

Usage

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

Arguments

Details

The residual degrees-of-freedom is estimated separately for each model using the sum of squared normalized residuals.

Value

A numeric vector with one element for each model.

Description

This expands the mean coefficients of the linear mixed model into one coefficient per level of the original variable, i.e., including the reference level(s) where the fitted coefficients are 0.

Usage

## S3 method for class 'lmm'
dummy.coef(object, use.na = FALSE, ...)

Arguments

Examples

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

## fit mixed model with interaction
eUN.lmm <- lmm(Y ~ visit*group, repetition =~visit|id, data = dL, df = FALSE)
dummy.coef(eUN.lmm)

## fit mixed model with baseline constraint
dL$drug <- dL$group
dL[dL$visit==1,"drug"] <- "d0"
eUN.clmm <- lmm(Y ~ visit + visit:drug, repetition =~visit|id, data = dL, df = FALSE)
dummy.coef(eUN.clmm)

dL$drug <- factor(dL$group, levels = c("none","d0","d1"))
dL[dL$visit==1,"drug"] <- "none"
eUN.clmm.none <- lmm(Y ~ visit:drug, repetition =~visit|id, data = dL, df = FALSE)
dummy.coef(eUN.clmm.none)

Description

Estimate average counterfactual outcome or contrast of outcome based on a linear mixed model.

Usage

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

Arguments

Details

The uncertainty is quantified assuming the contrast matrix to be a-priori known. Said otherwise the standard error does not account for the uncertainty about the covariate distribution.

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
e.YbyX1 <- effects(eUN.lmm, variable = "X1")
e.YbyX1
summary(e.YbyX1)
model.tables(e.YbyX1)
coef(e.YbyX1, type = "contrast")
effects(eUN.lmm, effects = "difference", variable = "X1")
effects(eUN.lmm, effects = "difference", variable = "X1", repetition = "3")

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

## auc
effects(eUN.lmm, type = "auc", variable = "X1")
effects(eUN.lmm, type = "auc", effects = "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 = "X2.factor", repetition = "3")
effects(eUN.lmmI, type = "change", variable = "X1.factor",
        conditional = "X2.factor", repetition = "3")
effects(eUN.lmmI, type = "auc", variable = "X1.factor", conditional = "X2.factor")

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

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

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
}

Description

Estimate standard errors, confidence intervals, and p-values for a smooth transformation of parameters from a linear mixed model.

Usage

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

Arguments

Details

Based a first order delta method to evaluate the variance of the estimate. The derivative of the transformation is evaluated using numerical differentiation (numDeriv::jacobian).

Argument robust: the Satterhwaite approximation for the degrees-of-freedom of robust standard errors is often unreliable. This is why the default is to use the degrees-of-freedom of the modeled based standard error instead.

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(object, p){nlme::ranef(object, 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)

dL2 <- reshape(d, direction = "long", idvar = c("id","X1"), 
               timevar = "time", times = c("1","2"), varying = c("Y1","Y2"), 
               v.names = "Y")

## 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
}

Description

Estimate standard errors, confidence intervals, and p-values for a smooth transformation of parameters from group-specific linear mixed models.

Usage

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

Arguments

Description

Evaluate the expected outcome conditional on covariates or the expected missing outcome value conditional on observed outcomes and covariates based on a 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

Details

Essentially a wrapper function of predict.lmm where the fitted value are stored in the outcome column of the dataset instead of in a separate column.

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)

Description

Evaluate the expected mean conditional to covariates or the expected missing outcome values conditional to observed outcome values and covariates. based on group-specific linear mixed models.

Usage

## S3 method for class 'mlmm'
fitted(object, newdata = NULL, simplify = TRUE, ...)

Arguments

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 (in pmol/l x hours).

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

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

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")

Description

Extract the influence function for an ML or REML estimator of parameters from a linear mixed model.

Usage

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

Arguments

Details

The influence function equals the individual score rescaled by the (inverse) information. With the expected information and for a lmm fitted using ML, the information is block diagonal so the influence function for the mean and variance parameters can be computed separately. Otherwise the information and individual score relative to all model parameters should be considered. The later is probablematic when using REML as the REML term is the ratio of two term linear in the individual contributions which is not itself linear in the individual contributions. As an add-hoc solution, the denominator is treated as fixed so the ratio is decomposed w.r.t. its numerator.

Value

A matrix with one row per observation and one column per parameter.

• df=TRUE: with an attribute "df" containing a numeric vector with one element per parameter.
  

• effects includes "gradient": with an attribute "gradient" containing a 3 dimensional array with dimension the number of parameters.

Examples

data(gastricbypassL)

#### Case WITHOUT cross terms ####
e.lmmREML <- lmm(weight ~ visit, method.fit = "REML", df = FALSE,
                 repetition = ~ visit|id, data = gastricbypassL)
e.lmmML <- lmm(weight ~ visit, method.fit = "ML", df = FALSE,
              repetition = ~ visit|id, data = gastricbypassL)

name.mu <- names(coef(e.lmmREML, effects = "mean"))
name.sigmakrho <- names(coef(e.lmmREML, effects = c("variance","correlation")))
info.REML <- information(e.lmmREML, effects = "all", transform.names = FALSE)
info.ML <- information(e.lmmML, effects = "all", transform.names = FALSE)
info.REML2ML <- information(e.lmmML, p = coef(e.lmmREML, effects = "all"),
                            effects = "all", transform.names = FALSE)

range(info.REML[name.mu,name.sigmakrho])
range(info.ML[name.mu,name.sigmakrho])
range(info.REML[name.mu,]-info.REML2ML[name.mu,])
range(iid(e.lmmREML, REML2ML = TRUE) - iid(e.lmmREML, REML2ML = FALSE))
## neglectable differences

#### Case WITH cross terms ####
e2.lmmREML <- lmm(glucagonAUC ~ visit + weight, method.fit = "REML", df = FALSE,
                 repetition = ~ visit|id, data = gastricbypassL)
e2.lmmML <- lmm(glucagonAUC ~ visit + weight, method.fit = "ML", df = FALSE,
              repetition = ~ visit|id, data = gastricbypassL)

name2.mu <- names(coef(e2.lmmREML, effects = "mean"))
name2.sigmakrho <- names(coef(e2.lmmREML, effects = c("variance","correlation")))
info2.REML <- information(e2.lmmREML, effects = "all", transform.names = FALSE)
info2.ML <- information(e2.lmmML, effects = "all", transform.names = FALSE)
info2.REML2ML <- information(e2.lmmML, p = coef(e2.lmmREML, effects = "all"),
                            effects = "all", transform.names = FALSE)

range(info2.REML[name.mu,]-info2.REML2ML[name.mu,])
## neglectable difference
range(info2.REML[name.mu,name.sigmakrho])
range(info2.ML[name.mu,name.sigmakrho])
range(iid(e2.lmmREML, REML2ML = TRUE) - iid(e2.lmmREML, REML2ML = FALSE))
## non-neglectable differences
diag(crossprod(iid(e2.lmmREML, REML2ML = TRUE)))/diag(vcov(e2.lmmREML))
diag(crossprod(iid(e2.lmmREML, REML2ML = FALSE)))/diag(vcov(e2.lmmREML))

Description

Extract the influence function for ML or REML estimators of parameters from group-specific linear mixed models.

Usage

## S3 method for class 'mlmm'
iid(x, effects = "contrast", p = NULL, ordering = "by", simplify = TRUE, ...)

Description

Extract the influence function of linear contrasts applied to an ML or REML estimator of parameters from a linear mixed model.

Usage

## S3 method for class 'Wald_lmm'
iid(x, effects = "Wald", transform.names = TRUE, ...)

Arguments

Value

A matrix with one row per observation and one column per parameter (effects="Wald" or effects="all") or a logical value (effects="test").

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")

Description

Extract or compute minus the 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).