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

Pool estimates from one or several linear mixed models. confint is the prefered method as this is meant for internal use.

Usage

## S3 method for class 'Wald_lmm'
aggregate(
  x,
  method,
  rhs = NULL,
  columns = NULL,
  qt = NULL,
  level = 0.95,
  df = NULL,
  tol = 1e-10,
  ...
)

Arguments

Details

Use a first order delta method to evaluate the standard error. When method="p.rejection" the p-value and distribution under the null is evaluated by simulating data using a Gaussian or t-copula, which can be time consumming. The number of sample is controlled by the argument n.sampleCopula in LMMstar.options.

See Also

confint.Wald_lmm, confint.rbindWald_lmm, confint.mlmm

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,
  type.information = NULL,
  robust = NULL,
  df = NULL,
  univariate = TRUE,
  multivariate = TRUE,
  p = NULL,
  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):

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

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

• model (optional): linear mixed model used to generate the Wald tests.

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 = NULL,
  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", ...)

Arguments

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,
  type = NULL,
  at,
  limits = c(-1, 1.00001),
  low = "blue",
  mid = "white",
  high = "red",
  midpoint = 0,
  labeller = "label_value",
  size.text = 16,
  obs.size = NULL,
  digits = 2,
  ...
)

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

## Visualize partial correlation
data(gastricbypassW, package = "LMMstar")
e.pCor <- partialCor(c(weight4,glucagonAUC4)~weight1+glucagonAUC1,
                     data = gastricbypassW)
plot(e.pCor)

## with repeated measurements
data(gastricbypassL, package = "LMMstar")
e.mpCor <- partialCor(c(weight,glucagonAUC)~time, repetition = ~visit|id,
                     data = gastricbypassL)
plot(e.mpCor)
}

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

dtS2 <- summarize(glucagonAUC + weight ~ time|id, data = gastricbypassL, na.rm = TRUE,
                  columns = add("correlation"))
plot(dtS2, variable = "glucagonAUC", type = "correlation", size.text = 1)
plot(dtS2, variable = "weight", 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",
  ordering = NULL,
  backtransform = NULL,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  transform.names = TRUE,
  simplify = TRUE,
  ...
)

Arguments

Value

A numeric vector

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,
  backtransform = NULL,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  transform.names = TRUE,
  simplify = TRUE,
  ...
)

Arguments

Value

A numeric vector

Description

Extract estimated value of linear contrasts involved in Wald tests.

Usage

## S3 method for class 'Wald_lmm'
coef(
  object,
  effects = "Wald",
  method = "none",
  backtransform = NULL,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = 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,
  level = 0.95,
  df = NULL,
  method = NULL,
  columns = NULL,
  ordering = NULL,
  backtransform = NULL,
  ...
)

Arguments

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

Description

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

Usage

## S3 method for class 'resample'
confint(
  object,
  parm = NULL,
  level = 0.95,
  null = NULL,
  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 of linear contrasts involved in Wald tests, along with corresponding p-values. 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 be used to pool estimates.

Usage

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

Arguments

Details

Multiple testing methods:

• "none": no adjustment

• "bonferroni": bonferroni adjustment

• "single-step": max-test adjustment performed by multcomp::glht() finding 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.

• "single-step2": max-test adjustment performed using Monte Carlo integration. 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.

• "holm", "hochberg", "hommel", "BH", "BY", "fdr": other adjustments performed by stats::p.adjust(). No confidence interval is computed.

• "free", "Westfall", "Shaffer": other adjustments performed by multcomp::glht(). No confidence interval is computed.

When degrees-of-freedom differs between individual hypotheses, "single-step2" is recommended over "single-step".

Pooling methods:

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

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. In presence of a categorical covariate varying within cluster, the correlation structure is structured in blocks. For instance with 2 categories:

• $A$: pairs of observations both at level 0.

• $B$: pairs of observations both at level 1.

• $C$: pairs of observations, one with level 0 and the other with level 1.

Usage

CS(formula, var.cluster, twin = TRUE, cross = "CS", var.time, add.time)

Arguments

Details

A typical formula would be ~1. It will model:

• $\sigma^2_{j}=\sigma^2$: a variance constant over repetitions.

• $\rho_{j,j'}=\rho$: a correlation constant over repetitions.

With 4 repetitions ($(j,j') \in \{1,2,3,4\}^2$ and $j\neq j'$) this leads to the following residuals variance-covariance matrix:

$\begin{bmatrix} \sigma^2 & \rho\sigma^2 & \rho\sigma^2 & \rho\sigma^2 \\ \rho\sigma^2 & \sigma^2 & \rho\sigma^2 & \rho\sigma^2 \\ \rho\sigma^2 & \rho\sigma^2 & \sigma^2 & \rho\sigma^2 \\ \rho\sigma^2 & \rho\sigma^2 & \rho\sigma^2 & \sigma^2 \end{bmatrix}$

With positive correlation, this is reparametrisation of a random intercept model.

Covariate constant within-cluster: it is possibly to stratify the covariance pattern on a categorical variable (e.g. sex) using a formula such as sex~1. This will lead to:

$\begin{bmatrix} \sigma_{\text{female}}^2 & \rho_{\text{female}}\sigma_{\text{female}}^2 & \rho_{\text{female}}\sigma_{\text{female}}^2 & \rho_{\text{female}}\sigma_{\text{female}}^2 \\ \rho_{\text{female}}\sigma_{\text{female}}^2 & \sigma_{\text{female}}^2 & \rho_{\text{female}}\sigma_{\text{female}}^2 & \rho_{\text{female}}\sigma_{\text{female}}^2 \\ \rho_{\text{female}}\sigma_{\text{female}}^2 & \rho_{\text{female}}\sigma_{\text{female}}^2 & \sigma_{\text{female}}^2 & \rho_{\text{female}}\sigma_{\text{female}}^2 \\ \rho_{\text{female}}\sigma_{\text{female}}^2 & \rho_{\text{female}}\sigma_{\text{female}}^2 & \rho_{\text{female}}\sigma_{\text{female}}^2 & \sigma_{\text{female}}^2 \end{bmatrix} \text{and} \begin{bmatrix} \sigma_{\text{male}}^2 & \rho_{\text{male}} \sigma^2_{\text{male}} & \rho_{\text{male}} \sigma^2_{\text{male}} & \rho_{\text{male}} \sigma^2_{\text{male}} \\ \rho_{\text{male}} \sigma^2_{\text{male}} & \sigma_{\text{male}}^2 & \rho_{\text{male}} \sigma^2_{\text{male}} & \rho_{\text{male}} \sigma^2_{\text{male}} \\ \rho_{\text{male}} \sigma^2_{\text{male}} & \rho_{\text{male}} \sigma^2_{\text{male}} & \sigma_{\text{male}}^2 & \rho_{\text{male}} \sigma^2_{\text{male}} \\ \rho_{\text{male}} \sigma^2_{\text{male}} & \rho_{\text{male}} \sigma^2_{\text{male}} & \rho_{\text{male}} \sigma^2_{\text{male}} & \sigma_{\text{male}}^2 \end{bmatrix}$

By default, the argument twin=TRUE meaning that ~sex will only make the correlation sex dependent:

$\begin{bmatrix} \sigma^2 & \rho_{\text{female}}\sigma^2 & \rho_{\text{female}}\sigma^2 & \rho_{\text{female}}\sigma^2 \\ \rho_{\text{female}}\sigma^2 & \sigma^2 & \rho_{\text{female}}\sigma^2 & \rho_{\text{female}}\sigma^2 \\ \rho_{\text{female}}\sigma^2 & \rho_{\text{female}}\sigma^2 & \sigma^2 & \rho_{\text{female}}\sigma^2 \\ \rho_{\text{female}}\sigma^2 & \rho_{\text{female}}\sigma^2 & \rho_{\text{female}}\sigma^2 & \sigma^2 \end{bmatrix} \text{and} \begin{bmatrix} \sigma^2 & \rho_{\text{male}} \sigma^2 & \rho_{\text{male}} \sigma^2 & \rho_{\text{male}} \sigma^2 \\ \rho_{\text{male}} \sigma^2 & \sigma^2 & \rho_{\text{male}} \sigma^2 & \rho_{\text{male}} \sigma^2 \\ \rho_{\text{male}} \sigma^2 & \rho_{\text{male}} \sigma^2 & \sigma^2 & \rho_{\text{male}} \sigma^2 \\ \rho_{\text{male}} \sigma^2 & \rho_{\text{male}} \sigma^2 & \rho_{\text{male}} \sigma^2 & \sigma^2 \end{bmatrix}$

To obtain a sex-dependent variance one can either set the argument twin=FALSE or specify a separate formula for the variance and correlation using a list: list(~sex,~sex)

$\begin{bmatrix} \sigma^2 & \rho_{\text{female}}\sigma^2 & \rho_{\text{female}}\sigma^2 & \rho_{\text{female}}\sigma^2 \\ \rho_{\text{female}}\sigma^2 & \sigma^2 & \rho_{\text{female}}\sigma^2 & \rho_{\text{female}}\sigma^2 \\ \rho_{\text{female}}\sigma^2 & \rho_{\text{female}}\sigma^2 & \sigma^2 & \rho_{\text{female}}\sigma^2 \\ \rho_{\text{female}}\sigma^2 & \rho_{\text{female}}\sigma^2 & \rho_{\text{female}}\sigma^2 & \sigma^2 \end{bmatrix} \text{and} \begin{bmatrix} \sigma^2 k_{\text{male}}^2 & \rho_{\text{male}} \sigma^2 k^2_{\text{male}} & \rho_{\text{male}} \sigma^2 k^2_{\text{male}} & \rho_{\text{male}} \sigma^2 k^2_{\text{male}} \\ \rho_{\text{male}} \sigma^2 k^2_{\text{male}} & \sigma^2 k_{\text{male}}^2 & \rho_{\text{male}} \sigma^2 k_{\text{male}}^2 & \rho_{\text{male}} \sigma^2 k^2_{\text{male}} \\ \rho_{\text{male}} \sigma^2 k^2_{\text{male}} & \rho_{\text{male}} \sigma^2 k^2_{\text{male}} & \sigma^2 k^2_{\text{male}} & \rho_{\text{male}} \sigma^2 k^2_{\text{male}} \\ \rho_{\text{male}} \sigma^2 k^2_{\text{male}} & \rho_{\text{male}} \sigma^2 k^2_{\text{male}} & \rho_{\text{male}} \sigma^2 k^2_{\text{male}} & \sigma^2 k^2_{\text{male}} \end{bmatrix}$

This is a reparametrization of the stratified CS structure.

Covariate varying within-cluster: it is not possible to stratify on covariates whose value is not constant across repetition within clusters. But one can make the variance and correlation dependent on such variable. By default (twin=TRUE) two correlation parameters are used: one when the variable values are same between the two repetitions ($\rho_W$) and another when the variable values differ between the two repetitions ($\rho_B$), e.g. if baseline refers to the firt two repetitions then using the formula ~baseline will model:

$\begin{bmatrix} \sigma^2 & \rho_W\sigma^2 & \rho_B\sigma^2 & \rho_B\sigma^2 \\ \rho_W\sigma^2 & \sigma^2 & \rho_B\sigma^2 & \rho_B\sigma^2 \\ \rho_B\sigma^2 & \rho_B\sigma^2 & \sigma^2 & \rho_W\sigma^2 \\ \rho_B\sigma^2 & \rho_B\sigma^2 & \rho_W\sigma^2 & \sigma^2 \end{bmatrix}$

With positive correlation, this is reparametrisation of a nested random effect model.

Setting twin=FALSE enables to obtain a correlation parameter specific to each variable value and each difference in variable value:

$\begin{bmatrix} \sigma^2 & \rho_{\text{baseline}}\sigma^2 & \rho_{\text{baseline,follow-up}}\sigma^2 k_{\text{follow-up}} & \rho_{\text{baseline,follow-up}}\sigma^2 k_{\text{follow-up}} \\ \rho_{\text{baseline}}\sigma^2 & \sigma^2 & \rho_{\text{baseline,follow-up}}\sigma^2 k_{\text{follow-up}} & \rho_{\text{baseline,follow-up}}\sigma^2 k_{\text{follow-up}} \\ \rho_{\text{baseline,follow-up}}\sigma^2 k_{\text{follow-up}} & \rho_{\text{baseline,follow-up}}\sigma^2 k_{\text{follow-up}} & \sigma^2 k^2_{\text{follow-up}} & \rho_{\text{follow-up}}\sigma^2 k^2_{\text{follow-up}} \\ \rho_{\text{baseline,follow-up}}\sigma^2 k_{\text{follow-up}} & \rho_{\text{baseline,follow-up}}\sigma^2 k_{\text{follow-up}} & \rho_{\text{follow-up}}\sigma^2 k^2_{\text{follow-up}} & \sigma^2 k^2_{\text{follow-up}} \end{bmatrix}$

Setting cross="ID" only considers correlation between observations with identical covariate values. Consider clusters of two subjects and modeling: within-subject correlation and within-time correlation for subjects from the same cluster but assuming independence between observations measured at different timepoints on different subjects:

$\begin{bmatrix} \sigma^2 & \rho_{\text{subject}}\sigma^2 & \rho_{\text{subject}}\sigma^2 & \rho_{\text{subject}}\sigma^2 & \rho_{\text{time}}\sigma^2 & 0 & 0 & 0 \\ \rho_{\text{subject}}\sigma^2 & \sigma^2 & \rho_{\text{subject}}\sigma^2 & \rho_{\text{subject}}\sigma^2 & 0 & \rho_{\text{time}}\sigma^2 & 0 & 0 \\ \rho_{\text{subject}}\sigma^2 & \rho_{\text{subject}}\sigma^2 & \sigma^2 & \rho_{\text{subject}}\sigma^2 & 0 & 0 & \rho_{\text{time}}\sigma^2 & 0 \\ \rho_{\text{subject}}\sigma^2 & \rho_{\text{subject}}\sigma^2 & \rho_{\text{subject}}\sigma^2 & \sigma^2 & 0 & 0 & 0 & \rho_{\text{time}}\sigma^2 \\ \rho_{\text{time}}\sigma^2 & 0 & 0 & 0 & \sigma^2 & \rho_{\text{subject}}\sigma^2 & \rho_{\text{subject}}\sigma^2 & \rho_{\text{subject}}\sigma^2\\ 0 & \rho_{\text{time}}\sigma^2 & 0 & 0 & \rho_{\text{subject}}\sigma^2 & \sigma^2 & \rho_{\text{subject}}\sigma^2 & \rho_{\text{subject}}\sigma^2 \\ 0 & 0 & \rho_{\text{time}}\sigma^2 & 0 & \rho_{\text{subject}}\sigma^2 & \rho_{\text{subject}}\sigma^2 & \sigma^2 & \rho_{\text{subject}}\sigma^2\\ 0 & 0 & 0 & \rho_{\text{time}}\sigma^2 & \rho_{\text{subject}}\sigma^2 & \rho_{\text{subject}}\sigma^2 & \rho_{\text{subject}}\sigma^2 & \sigma^2 \\ \end{bmatrix}$

Value

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

Examples

## Not run: 
data(gastricbypassL, package = "LMMstar")
gastricbypassL$sex <- factor(as.numeric(gastricbypassL$id) %% 2,
                             labels = c("female", "male"))
gastricbypassL$baseline <- factor(gastricbypassL$time < 0,
                             labels = c("baseline", "follow-up"))
gastricbypassL$family <- paste0("F",(as.numeric(gastricbypassL$id)-1) %/% 2)

#### default: constant variance and correlation ####
eCS.default <- lmm(glucagonAUC ~ visit, repetition = ~visit|id,
                 structure = CS, data = gastricbypassL)
sigma(eCS.default)
paramCS.default <- model.tables(eCS.default, effects = "all")
paramCS.default
paramCS.default["sigma","estimate"]^2 * paramCS.default["rho(id)","estimate"]

#### Covariate constant within-cluster #### 
## stratified: strata specific variance and correlation
eCS.strata <- lmm(glucagonAUC ~ visit, repetition = ~visit|id,
                 structure = CS(sex~1), data = gastricbypassL)
sigma(eCS.strata)
model.tables(eCS.strata, effects = "all")

## covariate: constant variance and sex-dependent correlation
eCS.sexCor <- lmm(glucagonAUC ~ visit, repetition = ~visit|id,
                  structure = CS(~sex), data = gastricbypassL)
sigma(eCS.sexCor)

## covariate: sex-dependent variance and correlation
eCS.sex <- lmm(glucagonAUC ~ visit, repetition = ~visit|id,
               structure = CS(~sex, twin = TRUE), data = gastricbypassL)
sigma(eCS.sex)

eCS.sex2 <- lmm(glucagonAUC ~ visit, repetition = ~visit|id,
               structure = CS(list(~sex, ~sex)), data = gastricbypassL)
sigma(eCS.sex2)

#### Covariate varying within-cluster ####

## twin: within (rho id/basline) and between (rho id) correlation
eBCS.default <- lmm(glucagonAUC ~ visit, repetition = ~visit|id,
                 structure = CS(~baseline), data = gastricbypassL)
sigma(eBCS.default)
model.tables(eBCS.default, effects = "all")

## heterogeneous: within (rho id/basline) and between (rho id) correlation
eBCS.hetero <- lmm(glucagonAUC ~ visit, repetition = ~visit|id,
                    structure = CS(~baseline, twin = FALSE), data = gastricbypassL)
sigma(eBCS.hetero)
model.tables(eBCS.hetero, effects = "all")

#### Cross compound symmetry ####
eCS.cross <- lmm(glucagonAUC ~ visit, repetition = ~visit|id,
                 structure = CS(~baseline, cross = "ID"),
                 data = gastricbypassL)
sigma(eCS.cross)
model.tables(eCS.cross, effects = "all")

## artificially create a two level structure family-id
## index family member: first, second
gastricbypassL <- gastricbypassL[order(gastricbypassL$id),]
gastricbypassL$member <- repetition(~id|family, data = gastricbypassL,
                                    type = "consecutive", label.rep = c("I","II"))
gastricbypassL$visit.member <- interaction(gastricbypassL[,c("visit","member")])

## constant correlation within subject, constant correlation within time,
## no correlation within family for different subjects at different times
eCS.family <- lmm(glucagonAUC ~ visit, repetition = ~visit.member|family,
                  structure = CS(~id+visit, cross = "ID"),
                  data = gastricbypassL)
sigma(eCS.family)
model.tables(eCS.family, effects = "all")

eHCS.family <- lmm(glucagonAUC ~ visit, repetition = ~visit.member|family,
                  structure = CS(~member+visit, twin = TRUE, cross = "ID"),
                  data = gastricbypassL) ## lack of convergence
sigma(eHCS.family)
model.tables(eHCS.family, effects = "all")

## End(Not run)

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