Constrained inference for linear mixed models.

Description

Constrained inference on linear fixed and mixed models using residual bootstrap. Covariates and random effects are permitted but not required.

Appropriate credit should be given when publishing results obtained using CLME, or when developing other programs/packages based off of this one. Use citation(package="CLME") for Bibtex information.

The work was produced in part with funding from the Intramural Research Program of the NIH, National Institute of Environmental Health Sciences (Z01 ES101744).

Details

This package was introduced in Jelsema and Peddada (2016). The primary function is clme. The other functions in this package may be run separately, but in general are designed for use by clme.

The method which is implemented is the constrained linear mixed effects model described in Farnan, Ivanova, and Peddada (2014). See that paper for more details regarding the method. Here we give a brief overview of the assumed model:

Y = X_{1}\theta_{1} + X_{2}\theta_{2} + U\xi + \epsilon

where

• X_1 is a N \times p_1 design matrix.

• \theta_1 are the coefficients (often treatment effects).

• X_2 is a N \times p_2 matrix of fixed covariates.

• \theta_1 are the coefficients for the covariates.

• U is a N \times c matrix of random effects.

• \xi is a zero-mean random vector with covariance T (see below).

• \epsilon is a zero-mean random vector with covariance \Sigma (see below).

Neither covariates (X_2) nor random effects (U) are required by the model or CLME. The covariance matrix of \xi is given by:

T = diag\left( \tau^{2}_{1}I_{c_{1}}, \tau^{2}_{2}I_{c_{2}} , \dots , \tau^{2}_{q}I_{c_{q}} \right)

The first c_{1} random effects will share a common variance, \tau^{2}_{1}, the next c_{2} random effects will share a common variance, and so on. Note that c = \sum_{i=1}^{q} c_i. Homogeneity of variances in the random effects can be induced by letting q=1 (hence c_{1}=c=ncol(U)).

Similarly, the covariance matrix of \epsilon is given by:

\Sigma = diag\left( \sigma^{2}_{1}I_{n_{1}}, \sigma^{2}_{2}I_{n_{2}} , \dots , \sigma^{2}_{q}I_{n_{k}} \right)

Again, the first n_{1} observations will share a common variance, \sigma^{2}_{1}, the next n_{2} will share a common variance, and so on. Note that N = \sum_{i=1}^{k} n_i. Homogeneity of variances in the residuals can be induced by letting k=1.

The order constraints are defined by the matrix A. This is an r \times p matrix where r is the number of constraints, and p=p_{1}+p_{2} is the dimension of \theta = ( \theta_{1}' , \theta_{2}')'. Formally the hypothesis being tested is:

H_{a}: A\theta > 0

For several default orders (simple, umbrella, simple tree) the A matrix can be automatically generated. Alternatively, the user may define a custom A matrix to test other patterns among the elements of \theta. See create.constraints and clme for more details.

For computational reasons, the implementation is not identical to the model expressed. Particularly, the fixed-effects matrix (or matrices) and the random effects matrix are assumed to be columns in a data frame, not passed as matrices. The A matrix is not r\ times p, but r\ times 2, where each row gives the indices of the constrained coefficients. See create.constraints for further explanation.

The creation of this package CLME, this manual, and the vignette were all supported by the Intramural Research Program of the United States' National Institutes of Health (Z01 ES101744).

Author(s)

Maintainer: Casey M. Jelsema jelsema.casey@gmail.com

Authors:

• Shyamal D. Peddada

References

Jelsema, C. M. and Peddada, S. D. (2016). CLME: An R Package for Linear Mixed Effects Models under Inequality Constraints. Journal of Statistical Software, 75(1), 1-32. doi:10.18637/jss.v075.i01

Farnan, L., Ivanova, A., and Peddada, S. D. (2014). Linear Mixed Efects Models under Inequality Constraints with Applications. PLOS ONE, 9(1). e84778. doi: 10.1371/journal.pone.0084778

See Also

Useful links:

• Report bugs at https://github.com/jelsema/CLME/issues

Akaike information criterion

Description

Calculates the Akaike and Bayesian information criterion for objects of class clme.

Calculates the Akaike and Bayesian information criterion for objects of class clme.

Usage

## S3 method for class 'clme'
AIC(object, ..., k = 2)

## S3 method for class 'summary.clme'
AIC(object, ..., k = 2)

Arguments

Details

The log-likelihood is assumed to be the Normal distribution. The model uses residual bootstrap methodology, and Normality is neither required nor assumed. Therefore the log-likelihood and these information criterion may not be useful measures for comparing models. For k=2, the function computes the AIC. To obtain BIC, set k = log( n/(2*pi) ); which the method BIC.clme does.

Value

Returns the information criterion (numeric).

See Also

CLME-package clme

CLME-package clme

Examples

data( rat.blood )

cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)

AIC( clme.out )
AIC( clme.out, k=log( nobs(clme.out)/(2*pi) ) )

Bayesian information criterion

Description

Calculates the Bayesian information criterion for objects of class clme.

Calculates the Akaike and Bayesian information criterion for objects of class clme.

Usage

## S3 method for class 'clme'
BIC(object, ..., k = log(nobs(object)/(2 * pi)))

## S3 method for class 'summary.clme'
BIC(object, ..., k = log(nobs(object)/(2 * pi)))

Arguments

Details

The log-likelihood is assumed to be the Normal distribution. The model uses residual bootstrap methodology, and Normality is neither required nor assumed. Therefore the log-likelihood and these information criterion may not be useful measures for comparing models. For k=2, the function computes the AIC. To obtain BIC, set k = log( n/(2*pi) ); which the method BIC.clme does.

Value

Returns the Bayesian information criterion (numeric).

See Also

CLME-package clme

CLME-package clme

Examples

data( rat.blood )

cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)

BIC( clme.out )
BIC( clme.out, k=log( nobs(clme.out)/(2*pi) ) )

Variance components

Description

Extracts variance components for objects of class clme.

Usage

VarCorr(x, sigma, rdig)

## S3 method for class 'summary.clme'
VarCorr(x, sigma, rdig)

## S3 method for class 'clme'
VarCorr(x, sigma, rdig)

Arguments

Value

Numeric.

See Also

CLME-package clme

Examples

data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)
                 
VarCorr( clme.out )

Constrained Inference for Linear Mixed Effects Models

Description

Constrained inference for linear fixed or mixed effects models using distribution-free bootstrap methodology

Usage

clme(
  formula,
  data = NULL,
  gfix = NULL,
  constraints = list(),
  tsf = lrt.stat,
  tsf.ind = w.stat.ind,
  mySolver = "LS",
  all_pair = FALSE,
  verbose = c(FALSE, FALSE, FALSE),
  ...
)

Arguments

Details

If any random effects are included, the function computes MINQUE estimates of variance components. After, clme_em is run to obtain the observed values. If nsim>0, a bootstrap test is performed using resid_boot. For the argument levels the first list element should be the column index (in data) of the constrained effect. The second element should be the true order of the levels.

Value

The output of clme is an object of the class clme, which is list with elements:

• theta estimates of \theta coefficients

• theta estimates of \theta_0 coefficients under the null hypothesis

• ssq estimate of residual variance(s), \sigma^{2}_{i}.

• tsq estimate of random effects variance component(s), \tau^{2}_{i}.

• cov.theta the unconstrained covariance matrix of \theta

• ts.glb test statistic for the global hypothesis.

• ts.ind test statistics for each of the constraints.

• mySolver the solver used for isotonization.

• constraints list containing the constraints (A) and the contrast for the global test (B).

• dframe data frame containing the variables in the model.

• residuals matrix containing residuals. For mixed models three types of residuals are given.

• random.effects estimates of random effects.

• gfix group sample sizes for residual variances.

• gran group sizes for random effect variance components.

• gfix_group group names for residual variances.

• formula the formula used in the model.

• call the function call.

• order list describing the specified or estimated constraints.

• P1 the number of constrained parameters.

• nsim the number of bootstrap simulations used for inference.

Note

The argument constraints is a list containing the order restrictions. The elements are order, node, decreasing, A, and B, though not all are necessary. The function can calculate the last two for default orders (simple, umbrella, or simple tree). For default orders, constraints should be a list containing any subset of order, node, and descending. See Figure 1 from Jelsema \& Peddada (2016); the pictured node of the simple tree orders (middle column) is 1, and the node for the umbrella orders (right column) is 3. These may be vectors (e.g. order=('simple','umbrella') ). If any of these three are missing, the function will test for all possible values of the missing element(s), excluding simple tree.

For non-default orders, the elements A and B should be provided. A is an r \times2 matrix (where r is the number of linear constraints, 0 < r. Each row should contain two indices, the first element is the index of the lesser coefficient, the second element is the index of the greater coefficient. So a row of (1,2) corresponds to the constraint \theta_1 \leq \theta_2, and a row (4,3) corresponds to the constraint \theta_4 \leq \theta_3, etc. Element B should hold similar contrasts, specifically those needed for calculating the Williams' type test statistic (B is only needed if tsf=w.stat) The argument tsf is a function to calculate the desired test statistic. The default function calculates likelihood ratio type test statistic. A Williams type test statistic, which is the maximum of the test statistic over the constraints in constraints\$B, is also available, and custom functions may be defined. See w.stat for details. By default, homogeneity of variances is assumed for residuals (e.g., gfix does not define groups) and for each random effect. Some values can be passed to clme that are not used in this function. For instance, seed and nsim can each be passed as an argument here, and summary.clme will use these values.

References

Jelsema, C. M. and Peddada, S. D. (2016). CLME: An R Package for Linear Mixed Effects Models under Inequality Constraints. Journal of Statistical Software, 75(1), 1-32. doi:10.18637/jss.v075.i01

Examples

data( rat.blood )
cons <- list(order="simple", decreasing=FALSE, node=1 )

clme.out <- clme(mcv ~ time + temp + sex + (1|id), data=rat.blood , 
                 constraints=cons, seed=42, nsim=10 )

Constrained EM algorithm for linear fixed or mixed effects models.

Description

clme_em_fixed performs a constrained EM algorithm for linear fixed effects models.

clme_em_mixed performs a constrained EM algorithm for linear mixed effects models.

clme_em is the general function, it will call the others. These Expectation-maximization (EM) algorithms estimate model parameters and compute a test statistic.

Usage

clme_em_fixed(
  Y,
  X1,
  X2 = NULL,
  U = NULL,
  Nks = dim(X1)[1],
  Qs = dim(U)[2],
  constraints,
  mq.phi = NULL,
  tsf = lrt.stat,
  tsf.ind = w.stat.ind,
  mySolver = "LS",
  em.iter = 500,
  em.eps = 1e-04,
  all_pair = FALSE,
  dvar = NULL,
  verbose = FALSE,
  ...
)

clme_em_mixed(
  Y,
  X1,
  X2 = NULL,
  U = NULL,
  Nks = dim(X1)[1],
  Qs = dim(U)[2],
  constraints,
  mq.phi = NULL,
  tsf = lrt.stat,
  tsf.ind = w.stat.ind,
  mySolver = "LS",
  em.iter = 500,
  em.eps = 1e-04,
  all_pair = FALSE,
  dvar = NULL,
  verbose = FALSE,
  ...
)

clme_em(
  Y,
  X1,
  X2 = NULL,
  U = NULL,
  Nks = nrow(X1),
  Qs = ncol(U),
  constraints,
  mq.phi = NULL,
  tsf = lrt.stat,
  tsf.ind = w.stat.ind,
  mySolver = "LS",
  em.iter = 500,
  em.eps = 1e-04,
  all_pair = FALSE,
  dvar = NULL,
  verbose = FALSE,
  ...
)

Arguments

Details

Argument constraints is a list including at least the elements A, B, and Anull. This argument can be generated by function create.constraints.

Value

The function returns a list with the elements:

• theta coefficient estimates.

• theta.null vector of coefficient estimates under the null hypothesis.

• ssq estimate of residual variance term(s).

• tsq estimate of variance components for any random effects.

• cov.theta covariance matrix of the unconstrained coefficients.

• ts.glb test statistic for the global hypothesis.

• ts.ind test statistics for each of the constraints.

• mySolver the solver used for isotonization.

Note

There are few error catches in these functions. If only the EM estimates are desired, users are recommended to run clme setting nsim=0.

By default, homogeneous variances are assumed for the residuals and (if included) random effects. Heterogeneity can be induced using the arguments Nks and Qs, which refer to the vectors (n_{1}, n_{2}, \ldots, n_{k}) and (c_{1}, c_{2}, \ldots, c_{q}) , respectively. See CLME-package for further explanation the model and these values.

See w.stat and lrt.stat for more details on using custom test statistics.

See Also

CLME-package clme create.constraints lrt.stat w.stat

Examples

data( rat.blood )

model_mats <- model_terms_clme( mcv ~ time + temp + sex + (1|id), data = rat.blood )

Y  <- model_mats$Y
X1 <- model_mats$X1
X2 <- model_mats$X2
U  <- model_mats$U

cons <- list(order = "simple", decreasing = FALSE, node = 1 )

clme.out <- clme_em(Y = Y, X1 = X1, X2 = X2, U = U, constraints = cons)

Computes various types of residuals

Description

Computes several types of residuals for objects of class clme.

Usage

clme_resids(formula, data, gfix = NULL)

Arguments

Details

For fixed-effects models Y = X\beta + \epsilon, residuals are given as \hat{e} = Y - X\hat{\beta}. For mixed-effects models Y = X\beta + + U\xi + \epsilon, three types of residuals are available. PA = Y - X\hat{\beta}\ SS = U\hat{\xi}\ FM = Y - X\hat{\beta} - U\hat{\xi}

Value

List containing the elements PA, SS, FM, cov.theta, xi, ssq, tsq. PA, SS, FM are defined above (for fixed-effects models, the residuals are only PA). Then cov.theta is the unconstrained covariance matrix of the fixed-effects coefficients, xi is the vector of random effect estimates, and ssq and tsq are unconstrained estimates of the variance components.

Note

There are few error catches in these functions. If only the EM estimates are desired, users are recommended to run clme setting nsim=0.

By default, homogeneous variances are assumed for the residuals and (if included) random effects. Heterogeneity can be induced using the arguments Nks and Qs, which refer to the vectors (n_{1}, n_{2}, \ldots, n_{k}) and (c_{1}, c_{2}, \ldots, c_{q}) , respectively. See CLME-package for further explanation the model and these values.

See w.stat and lrt.stat for more details on using custom test statistics.

See Also

CLME-package clme

Examples

## Not run: 
data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )

clme.out <- clme_resids(mcv ~ time + temp + sex + (1|id), data = rat.blood )

## End(Not run)

Individual confidence intervals

Description

Calculates confidence intervals for fixed effects parameter estimates in objects of class clme.

Calculates confidence intervals for fixed effects parameter estimates in objects of class clme.

Usage

## S3 method for class 'clme'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'summary.clme'
confint(object, parm, level = 0.95, ...)

Arguments

Details

Confidence intervals are computed using Standard Normal critical values. Standard errors are taken from the covariance matrix of the unconstrained parameter estimates.

Value

Returns a matrix with two columns named lcl and ucl (lower and upper confidence limit).

See Also

CLME-package clme

CLME-package clme

Examples

data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)

confint( clme.out )

Generate common order constraints

Description

Automatically generates the constraints in the format used by clme. Allowed orders are simple, simple tree, and umbrella orders.

Usage

create.constraints(P1, constraints)

Arguments

Details

The elements of constraints are:

• order: string. Currently “simple”, “simple.tree” and “umbrella” are supported.

• node: numeric, the node of the coefficients (unnecessary for simple orders).

• decreasing: logical. For simple orders, is the trend decreasing? For umbrella and simple tree, does the nodal parameter have the greatest value (e.g., the peak, instead of the valley)?

See clme for more information and a depiction of these three elements.

Value

The function returns a list containing the elements of input argument constraints as well as

• A matrix of dimension r \times 2 containing the order constraints, where r is the number of linear constraints.

• B matrix containing the contrasts necessary for computation of the Williams' type test statistic (may be identical to A).

• Anull matrix similar to A which defines all possible constraints. Used to obtain parameter estimates under the null hypothesis.

• order the input argument for constraints\$order.

• node the input argument for constraints\$node.

• decreasing the input argument for constraints\$decreasing

See w.stat for more information on B

Note

The function clme also utilizes the argument constraints. For clme, this argument may either be identical to the argument of this function, or may be the output of create.constraints (that is, a list containing appropriate matrices A, Anull, and if necessary, B).

An example the the A matrix might be:

This matrix defines what CLME describes as a decreasing umbrella order. The first row defines the constraint that \theta_1 \leq \theta_2, the second row defined the constraint \theta_2 \leq \theta_3, the third row defines \theta_4 \leq \theta_3, and so on. The values are indexes, and the left column is the index of the parameter constrained to be smaller.

See Also

clme, w.stat

Examples

## Not run: 
  # For simple order, the node does not matter
  create.constraints( P1 = 5, constraints = list( order='simple' , 
                                                  decreasing=FALSE ))
  
  # Compare constraints against decreasing=TRUE
  create.constraints( P1 = 5, constraints=list( order='simple' , 
                                                decreasing=TRUE ))
  
  # Umbrella order
  create.constraints( P1 = 5, constraints=list( order='umbrella' , node=3
                                                , decreasing=FALSE ))

## End(Not run)

Fibroid Growth Study

Description

This data set contains a subset of the data from the Fibroid Growth Study.

Usage

data(fibroid)

Format

A data frame containing 240 observations on 9 variables.

Details

The response variable lfgr was calculated as the change in log fibroid volume, divided by the length of time between measurements. The growth rates were averaged to produce a single value for each fibroid, which was scaled to represent a 6-month percent change in volume.

References

Peddada, Laughlin, Miner, Guyon, Haneke, Vahdat, Semelka, Kowalik, Armao, Davis, and Baird(2008). Growth of Uterine Leiomyomata Among Premenopausal Black and White Women. Proceedings of the National Academy of Sciences of the United States of America, 105(50), 19887-19892. URL http://www.pnas.org/content/105/50/19887.full.pdf.

Extract fixed effects

Description

Extracts the fixed effects estimates from objects of class clme.

Usage

## S3 method for class 'clme'
fixef(object, ...)

## S3 method for class 'summary.clme'
fixef(object, ...)

## S3 method for class 'clme'
fixef(object, ...)

fixed.effects(object, ...)

## S3 method for class 'summary.clme'
fixed.effects(object, ...)

## S3 method for class 'clme'
fixed.effects(object, ...)

## S3 method for class 'clme'
coefficients(object, ...)

## S3 method for class 'clme'
coef(object, ...)

## S3 method for class 'summary.clme'
coefficients(object, ...)

## S3 method for class 'summary.clme'
coef(object, ...)

Arguments

Value

Returns a numeric vector.

See Also

CLME-package clme

Examples

data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)

fixef( clme.out )

Extract formula

Description

Extracts the formula from objects of class clme.

Usage

## S3 method for class 'clme'
formula(x, ...)

Arguments

Details

The package CLME parametrizes the model with no intercept term. If an intercept was included, it will be removed automatically.

Value

Returns a formula object

See Also

CLME-package clme

Examples

data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)

formula( clme.out )

Constructor method for objects S3 class clme

Description

Test if an object is of class clme or coerce an object to be such.

Usage

is.clme(x)

as.clme(x, ...)

Arguments

Value

Returns an object of the class clme.

See Also

CLME-package clme

Examples

data( rat.blood )

cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)

is.clme( clme.out )
as.clme( clme.out )

Log-likelihood

Description

Computes the log-likelihood of the fitted model for objects of class clme.

Usage

## S3 method for class 'clme'
logLik(object, ...)

## S3 method for class 'summary.clme'
logLik(object, ...)

Arguments

Details

The log-likelihood is computed using the Normal distribution. The model uses residual bootstrap methodology, and Normality is neither required nor assumed. Therefore the log-likelihood may not be a useful measure in the context of CLME.

Value

Numeric.

See Also

CLME-package clme

logLik.clme

Examples

data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)

logLik( clme.out )

Likelihood ratio type statistic (global)

Description

Calculates the likeihood ratio type test statistic (under Normality assumption) for a constrained linear mixed effects model. This is the default test statistic for CLME.

Usage

lrt.stat(theta, theta.null, cov.theta, ...)

Arguments

Value

Output is a numeric value.

Note

This is an internal function, unlikely to be useful outside of CLME-package. To define custom functions, the arguments available are:

theta, theta.null, cov.theta, B, A, Y, X1, X2, U, tsq, ssq, Nks, and Qs.

Of the additional arguments, B and A are identical to those produced by create.constraints. The rest, Y, X1, X2, U, tsq, , ssq, Nks, and Qs, are equivalent to arguments to clme_em.

Custom functions must produce numeric output. Output may have length greater than 1, which corresponds to testing multiple global hypotheses.

See Also

clme_em, w.stat

Examples

data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )

clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)

# Individually compute lrt statistic
lrt.stat(clme.out$theta, clme.out$theta.null, clme.out$cov.theta )

MINQUE Algorithm

Description

Algorithm to obtain MINQUE estimates of variance components of a linear mixed effects model.

Usage

minque(
  Y,
  X1,
  X2 = NULL,
  U = NULL,
  Nks = dim(X1)[1],
  Qs = dim(U)[2],
  mq.eps = 1e-04,
  mq.iter = 500,
  verbose = FALSE,
  ...
)

Arguments

Details

By default, the model assumes homogeneity of variances for both the residuals and the random effects (if included). See the Details in clme_em for more information on how to use the arguments Nks and Qs to permit heterogeneous variances.

Value

The function returns a vector of the form (\tau^{2}_{1}, \tau^{2}_{2}, \ldots, \tau^{2}_{q}, \sigma^{2}_{1},\sigma^{2}_{2},\ldots, \sigma^{2}_{k})'. If there are no random effects, then the output is just (\sigma^{2}_{1},\sigma^{2}_{2},\ldots, \sigma^{2}_{k})'.

Note

This function is called by several other function in CLME to obtain estimates of the random effect variances. If there are no random effects, they will not call minque.

Examples

data( rat.blood )

model_mats <- model_terms_clme( mcv ~ time + temp + sex + (1|id) , 
                                data = rat.blood )
Y  <- model_mats$Y
X1 <- model_mats$X1
X2 <- model_mats$X2
U  <- model_mats$U

# No covariates or random effects
minque(Y = Y, X1 = X1 )

# Include covariates and random effects
minque(Y = Y, X1 = X1, X2 = X2, U = U )

Extract the model design matrix.

Description

Extracts the fixed-effects design matrix from objects of class clme.

Usage

## S3 method for class 'clme'
model.matrix(object, type = "fixef", ...)

## S3 method for class 'summary.clme'
model.matrix(object, ...)

Arguments

Value

Returns a matrix.

See Also

CLME-package clme

model.matrix.clme

Examples

## Not run: 
data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)

model.matrix( clme.out )

## End(Not run)

Create model matrices for clme

Description

Parses formulas to creates model matrices for clme.

Usage

model_terms_clme(formula, data)

Arguments

Value

A list with the elements:

Note

The first term on the right-hand side of the formula should be the fixed effect with constrained coefficients. Random effects are represented with a vertical bar, so for example the random effect U would be included by Y ~ X1 + (1|U).

The intercept is removed automatically. This is done to ensure that parameter estimates are of the means of interest, rather than being expressed as a mean with offsets.

See Also

CLME-package clme

Examples

data( rat.blood )
model_terms_clme( mcv ~ time + temp + sex + (1|id) , data = rat.blood )

Number of observations

Description

Obtains the number of observations used to fit an model for objects of class clme.

Usage

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

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

Arguments

Value

Numeric.

See Also

CLME-package clme

nobs.clme

Examples

data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)

nobs( clme.out )

S3 method to plot objects of class clme

Description

Generates a basic plot of estimated coefficients which are subject to constraints (\theta_1 ). Lines indicate individual constraints (not global tests) and significance.

Usage

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

Arguments

Note

While it is possible to plot the output of a clme fit, this will only plot the fitted means. To indicate significance, plotting must be performed on the summary of a clme fit. This method will change the class so that plot.summary.clme will be called properly.

See Also

CLME-package clme plot.summary.clme

Examples

## Not run: 
  set.seed( 42 )
  data( rat.blood )
  cons <- list(order = "simple", decreasing = FALSE, node = 1 )
  clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                   constraints = cons, seed = 42, nsim = 10)
  plot( clme.out )

## End(Not run)

S3 method to plot objects of class clme

Description

Generates a basic plot of estimated coefficients which are subject to constraints (\theta_1 ). Lines indicate individual constraints (not global tests) and significance.

Usage

## S3 method for class 'summary.clme'
plot(
  x,
  alpha = 0.05,
  legendx = "below",
  inset = 0.01,
  ci = FALSE,
  ylim = NULL,
  cex = 1.75,
  pch = 21,
  bg = "white",
  xlab = expression(paste("Component of ", theta[1])),
  ylab = expression(paste("Estimated Value of ", theta[1])),
  tree = NULL,
  ...
)

Arguments

Details

All of the individual contrasts in the constraints\$A matrix are tested and plotted. The global test is not represented (unless it happens to coincide with an individual contrast). Only the elements of \theta which appear in any constraints (e.g. the elements of \theta_{1}) are plotted. Coefficients for the covariates are not plotted. Solid lines denote no significant difference, while dashed lines denote statistical significance. Significance is determined by the individual p-value being less than or equal to the supplied \alpha threshold. By default a legend denoting the meaning of solid and dashed lines will be placed below the graph. Argument legendx may be set to a legend keyword (e.g. legend=''bottomright'') to place it inside the graph at the specified location. Setting legendx to FALSE or to a non-supported keyword suppresses the legend. Confidence intervals for the coefficients may be plotted. They are individual confidence intervals, and are computed using the covariance matrix of the unconstrained estimates of \theta_{1}. These confidence intervals have higher coverage probability than the nominal value, and as such may appear to be in conflict with the significance tests. Alternate forms of confidence intervals may be provided in future updates.#'

See Also

CLME-package clme

Examples

## Not run: 
  set.seed( 42 )
  data( rat.blood )
  cons <- list(order = "simple", decreasing = FALSE, node = 1 )
  clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                   constraints = cons, seed = 42, nsim = 10)
  clme.out2 <- summary( clme.out )
  plot( clme.out2 )

## End(Not run)

Printout of fitted object.

Description

Prints basic information on a fitted object of class clme.

Usage

## S3 method for class 'clme'
print(x, ...)

Arguments

Value

Text printed to console.

See Also

CLME-package clme

Examples

## Not run: 
data( rat.blood )
set.seed( 42 )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 10)

print( clme.out )

## End(Not run)

S3 method to print a summary for objects of class clme

Description

Summarizes the output of objects of class clme, such as those produced by clme. Prints a tabulated display of global and individual tests, as well as parameter estimates.

Usage

## S3 method for class 'summary.clme'
print(x, alpha = 0.05, digits = 4, ...)

Arguments

Value

NULL, just prints results to the console.

Note

The individual tests are performed on the specified order. If no specific order was specified, then the individual tests are performed on the estimated order.

See Also

CLME-package clme

Examples

## Not run: 
  set.seed( 42 )
  data( rat.blood )
  cons <- list(order = "simple", decreasing = FALSE, node = 1 )
  clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                   constraints = cons, seed = 42, nsim = 10)
  
  summary( clme.out )

## End(Not run)

Printout for variance components

Description

Prints variance components of an objects of clme.

Usage

## S3 method for class 'varcorr_clme'
print(object, rdig = 5, ...)

Arguments

Value

Text printed to console.

See Also

CLME-package clme

Examples

## Not run: 
data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)
                 
print.varcorr_clme( clme.out )

## End(Not run)

Extract random effects

Description

Extract random effects

Extract random effects

Usage

random.effects(object, ...)

## S3 method for class 'summary.clme'
random.effects(object, ...)

Arguments

Extract random effects

Description

Extracts the random effects estimates from objects of class clme.

Usage

## S3 method for class 'clme'
ranef(object, ...)

## S3 method for class 'summary.clme'
ranef(object, ...)

## S3 method for class 'clme'
ranef(object, ...)

## S3 method for class 'clme'
random.effects(object, ...)

Arguments

Value

Returns a numeric vector.

See Also

CLME-package clme

Examples

data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)

ranef( clme.out )

Experiment on mice

Description

This data set contains the data from an experiment on 24 Sprague-Dawley rats from Cora et al (2012).

Usage

data(rat.blood)

Format

A data frame containing 241 observations on 13 variables.

Details

The response variable lfgr was calculated as the change in log fibroid volume, divided by the length of time between measurements. The growth rates were averaged to produce a single value for each fibroid, which was scaled to represent a 6-month percent change in volume.

References

Cora M, King D, Betz L, Wilson R, and Travlos G (2012). Artifactual changes in Sprauge-Dawley rat hematologic parameters after storage of samples at 3 C and 21 C. Journal of the American Association for Laboratory Animal Science, 51(5), 616-621. URL http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3447451/.

Obtain Residual Bootstrap

Description

Generates bootstrap samples of the data vector.

Usage

resid_boot(
  formula,
  data,
  gfix = NULL,
  eps = NULL,
  xi = NULL,
  null.resids = TRUE,
  theta = NULL,
  ssq = NULL,
  tsq = NULL,
  cov.theta = NULL,
  seed = NULL,
  nsim = 1000,
  mySolver = "LS",
  ...
)

Arguments

Details

If any of the parameters theta, ssq, tsq, eps, or xi are provided, the function will use those values in generating the bootstrap samples. They will be estimated if not submitted. Ifnull.resids=TRUE, then theta will be projected onto the space of the null hypothesis ( H_{0}: \theta_1 = \theta_2 = ... = \theta_{p_1}) regardless of whether it is provided or estimated. To generate bootstraps with a specific theta, set null.residuals=FALSE.

Value

Output is N \ times nsim matrix, where each column is a bootstrap sample of the response data Y.

Note

This function is primarily designed to be called by clme.

By default, homogeneous variances are assumed for the residuals and (if included) random effects. Heterogeneity can be induced using the arguments Nks and Qs, which refer to the vectors (n_{1}, n_{2}, \ldots, n_{k}) and (c_{1}, c_{2}, \ldots, c_{q}) , respectively. See clme_em for further explanation of these values.

See Also

clme

Examples

data( rat.blood )
boot_sample <- resid_boot(mcv ~ time + temp + sex + (1|id), nsim = 10, 
                          data = rat.blood, null.resids = TRUE  )

Various types of residuals

Description

Computes several types of residuals for objects of class clme.

Usage

## S3 method for class 'clme'
residuals(object, type = "FM", ...)

## S3 method for class 'summary.clme'
residuals(object, type = "FM", ...)

Arguments

Details

For fixed-effects models Y = X\beta + \epsilon, residuals are given as

\hat{e} = Y - X\hat{\beta}

. For mixed-effects models Y = X\beta + + U\xi + \epsilon, three types of residuals are available. PA = Y - X\hat{\beta}\ SS = U\hat{\xi}\ FM = Y - X\hat{\beta} - U\hat{\xi}

Value

Returns a numeric matrix.

See Also

CLME-package clme

Examples

## Not run: 
data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)
                 
residuals( clme.out, type='PA' )

## End(Not run)

Shiny GUI for CLME

Description

Opens a graphical user interface to run CLME, built from the shiny package.

The UI for the shiny app in CLME

The server for the shiny app in CLME

Usage

shiny_clme()

shinyUI_clme

shinyServer_clme(input, output)

Arguments

Format

An object of class shiny.tag.list (inherits from list) of length 3.

Details

Currently the GUI does not allow specification of custom orders for the alternative hypothesis. Future versions may enable this capability. The data should be a CSV or table-delimited file with the first row being a header. Variables are identified using their column letter or number (e.g., 1 or A). Separate multiple variables with a comma (e.g., 1,2,4 or A,B,D), or select a range of variables with a dash (e.g., 1-4 or A-D). Set to 'None' (default) to indicate no covariates or random effects. If group levels for the constrained effect are character, they may not be read in the proper order. An extra column may contain the ordered group levels (it may therefore have different length than the rest of the dataset).

Note

This function is primarily designed to call clme.

Examples

## Not run:  shiny_clme() 

Residual variance components

Description

Extract residual variance components for objects of class clme.

Usage

## S3 method for class 'clme'
sigma(object, ...)

Arguments

Value

Numeric.

See Also

CLME-package clme

Examples

data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)
         
sigma( clme.out )

Residual variance components

Description

Extract residual variance components for objects of class clme.

Usage

## S3 method for class 'summary.clme'
sigma(object, ...)

Arguments

Value

Numeric.

See Also

CLME-package clme

Examples

data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)
                 
sigma( clme.out )

Produce summary values for objects of class clme

Description

Summarizes the output of objects of class clme, such as those produced by clme.

Usage

## S3 method for class 'clme'
summary(object, nsim = 1000, seed = NULL, verbose = c(FALSE, FALSE), ...)

Arguments

Value

The output of summary.clme is an object of the class summary.clme. This is a list containing the input object (of class clme), along with elements:

See Also

CLME-package clme

Examples

## Not run: 
  set.seed( 42 )
  data( rat.blood )
  cons <- list(order = "simple", decreasing = FALSE, node = 1 )
  clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                   constraints = cons, seed = 42, nsim = 10)
  
  summary( clme.out )

## End(Not run)

Variance-covariance matrix

Description

Extracts variance-covariance matrix for objects of class clme.

Usage

## S3 method for class 'clme'
vcov(object, ...)

## S3 method for class 'summary.clme'
vcov(object, ...)

Arguments

Value

Numeric matrix.

See Also

CLME-package clme

Examples

data( rat.blood )
cons <- list(order = "simple", decreasing = FALSE, node = 1 )
clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
                 constraints = cons, seed = 42, nsim = 0)
                 
vcov( clme.out )

Williams' Type Test Statistic.

Description

Calculates a Williams' type test statistic for a constrained linear mixed effects model.

Usage

w.stat(theta, cov.theta, B, A, ...)

w.stat.ind(theta, cov.theta, B, A, ...)

Arguments

Details

See create.constraints for an example of A. Argument B is similar, but defines the global contrast for a Williams' type test statistic. This is the largest hypothesized difference in the constrained coefficients. So for an increasing simple order, the test statistic is the difference between the two extreme coefficients, \theta_1 and \theta_{p_1}, divided by the standard error (unconstrained). For an umbrella order order, two contrasts are considered, \theta_1 to \theta_{s}, and \theta_{p_1} to \theta_{s}, each divided by the appropriate unconstrained standard error. A general way to express this statistic is:

W = max \theta_{B[i,2]} - \theta_{B[i,1]} / sqrt( VAR( \theta_{B[i,2]} - \theta_{B[i,1]} ) )

where the numerator is the difference in the constrained estimates, and the standard error in the denominator is based on the covariance matrix of the unconstrained estimates.

The function w.stat.ind does the same, but uses the A matrix which defines all of the individual constraints, and returns a test statistic for each constraints instead of taking the maximum.

Value

Output is a numeric value.

Note

See lrt.stat for information on creating custom test statistics.

Examples

theta  <- exp(1:4/4)
th.cov <- diag(4)
X1     <- matrix( 0 , nrow=1 , ncol=4 )
const  <- create.constraints( P1=4 , constraints=list(order='simple' ,
                                                    decreasing=FALSE) )

w.stat( theta , th.cov , const$B , const$A )

w.stat.ind( theta , th.cov , const$B , const$A )